Verify IndelRealigner

Biorunner88Biorunner88 Member
edited February 2016 in Ask the GATK team

Hi

I'm using IndelRealigner, and it seems to take long time, from 5 to 8h, I have a strange sensation that I may not put the right commands but I'm not sure. Previously of what I show below I have Marked duplicates with MarkDuplicates.

/home/Programas/samtools-0.1.19/samtools index /local/ktroule/SNPiR/Panc047/accepted_hits_RG_Reordered_deduplicated.bam

java -Xmx20g -jar /home/Programas/GATK-3.1-1-g07a4bf8/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /local/Referencias/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome.fa -I /local/SNPiR/Panc042/accepted_hits_RG_Reordered_deduplicated.bam --unsafe ALLOW_N_CIGAR_READS --known /local/Referencias/HG19/Mills_and_1000G_gold_standard.indels.hg19.vcf -o /local/SNPiR/Panc042/Panc042_indel.intervals

java -Xmx20g -jar /home/Programas/GATK-3.1-1-g07a4bf8/GenomeAnalysisTK.jar -T IndelRealigner --targetIntervals /local/SNPiR/Panc042/Panc042_indel.intervals  -R /local/Referencias/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome.fa --unsafe ALLOW_N_CIGAR_READS -I /local/SNPiR/Panc042/accepted_hits_RG_Reordered_deduplicated.bam -o /local/SNPiR/Panc042/accepted_hits_RG_Reordered_deduplicated_Realigned.bam

This is not really about an error, but about If the program takes that long and about if the parameters I have used are the rights ones, as I have seen that few ones use the downsampling.

Finally, if I have not read wrongly InderRealigner does not support multithread, at least not directly from the program.

thanks for yor time.

Post edited by Geraldine_VdAuwera on
Tagged:

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @Biorunner88
    Hi,

    Your commands look fine. I wonder if the program is simply stalling on a section where there are a lot of reads. Is the progress meter showing a particular position that seems to take a really long time?

    -Sheila

  • I'm just following a protocol specified for SNPiR program.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    If you mean the protocol outlined in http://lilab.stanford.edu/SNPiR/readme it is terribly out of date and will produce inferior results. Please have a look at our best practices for RNAseq.

  • Yes I meant that protocol. It's giving me painful times. i just wanted to try it among other three which includes GATK's best practices.

    In the meantime. Can I use the ouput of theses commands for any other RNAseq variant calling program besides the stated in GATK's best practices? Such as RVBoost?

    java -Xmx20g -jar /home/Programas/GATK/GenomeAnalysisTK.jar -T SplitNCigarReads -U ALLOW_N_CIGAR_READS -R /local/ktroule/Referencias/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome.fa -I /local/ktroule/SNPiR/Panc042/accepted_hits_RG_Reordered_deduplicated.bam -o /local/ktroule/SNPiR/Panc042/accepted_hits_RG_Reordered_deduplicated_SplitNCigar.bam

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Frankly I think that protocol is a waste of your time now, but it's up to you of course.

    I can't speak to the requirements of other programs but it should be usable in that way, yes.

  • Just another question related to this. I'm following GATK's best practices protocol for RNAseq. And there is one step that confuses me.

    Split'N'Trim + ReassignMappingQuality

    I have done
    java -Xmx20g -jar /home/Programas/GATK/GenomeAnalysisTK.jar -T SplitNCigarReads -U ALLOW_N_CIGAR_READS -R /local/Referencias/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome.fa -I /local/SNPiR/Panc047/accepted_hits_RG_Reordered_deduplicated.bam -o /local/SNPiR/Panc047/accepted_hits_RG_Reordered_deduplicated_SplitNCigar.bam

    I'm not sure if I can go directly to the indel realignment or I must do something related to the ReassignMappingQuality. I dont knwo if this last is part of the SplitNCigarReads or it's another step that I have not been able to find.

    thanks again for your time.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    It's done with a read filter that you add to your command. Which protocol document are you following?

  • Understood, so I should add ReassignOneMappingQuality -RMQF 255 -RMQT 60 to the previous command, right?

    I was looking in this documentation and no reference to MapQ is done. https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_rnaseq_SplitNCigarReads.php

    And I saw that now change in MapQ was done after the SplitNCigarReads, just 50, 3,1 and 0s for MapQ (Tophat2's MapQ )

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    You have to preface that with -rf to work.

    We have an article that details the full RNAseq best practices in the methods section of the documentation guide.

  • @Geraldine_VdAuwera said:
    You have to preface that with -rf to work.

    We have an article that details the full RNAseq best practices in the methods section of the documentation guide.

    Yes indeed, I'm using in the command, just didn't paste it here.

    Thanks for your time.

Sign In or Register to comment.