Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Reassign all mapping qualities above a certain value

ryanfriedman22ryanfriedman22 Washington University in St. LouisMember

I'm trying to do variant calling on RNA-seq data aligned with Tophat. I know that I need to use SplitNCigarReads, but I'm having some trouble reassigning mapping qualities. I can successfully reassign MQ of 255 as suggested in the best practices, but am still getting errors for MQ values that are above 60. Is there any way that I can reassign all mapping qualities above a certain value? It seems like there's only tools to do this for one specific MQ and not an entire range.

Best Answers

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Can you post the errors you're getting please?

  • ryanfriedman22ryanfriedman22 Washington University in St. LouisMember

    I'm running a bunch of jobs in parallel, so my error messages are all intermingled. However, here's an example of one message; it only occurs for values larger than 60.

    SAM/BAM/CRAM file [email protected]49330a8 appears to be using the wrong encoding for quality scores: we encountered an extremely high quality score of 65. Please see https://www.broadinstitute.org/gatk/guide?id=6470 for more details and options related to this error.

    This is appearing when I run the command:

    java -Xmx8000M -Xms8000M -jar ~ryan.friedman/GenomeAnalysisTK.jar -T PrintReads -R ${genomeDir}/crNeoH99.fasta -I accepted_hits.reordered.bam -o accepted_hits.remap.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -- filter_reads_with_N_cigar

  • ryanfriedman22ryanfriedman22 Washington University in St. LouisMember

    Now I'm getting a different error:

    Bad input: while fixing mis-encoded base qualities we encountered a read that was correctly encoded; we cannot handle such a mixture of reads so unfortunately the BAM must be fixed with some other tool

    using the command:

    java -Xmx8000M -Xms8000M -jar ~ryan.friedman/GenomeAnalysisTK.jar -T PrintReads -R ${genomeDir}/crNeoH99.fasta -I accepted_hits.reordered.bam -o accepted_hits.remap.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 --filter_reads_with_N_cigar --fix_misencoded_quality_scores

    Can I not both reassign MQ and fix misencoded scores? Alternatively, if I'm fixing the misencoded scores, do I even need to use PrintReads in the first place or can I just add the --fix_misencoded_quality_scores to the following command?

    java -Xmx8000M -Xms8000M -jar ~ryan.friedman/GenomeAnalysisTK.jar -T SplitNCigarReads -R ${genomeDir}/crNeoH99.fasta -I accepted_hits.remap.bam -o split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

  • SheilaSheila Broad InstituteMember, Broadie admin

    @ryanfriedman22
    Hi,

    Have a look at this thread for some more help.

    -Sheila

  • ryanfriedman22ryanfriedman22 Washington University in St. LouisMember

    Hi Sheila,

    Are you suggesting I should use the seqtk seq tool to realign the base quals in my original FASTQ files, then rerun Tophat and pass that through GATK to split reads? I have 521 different sets of RNA-seq data I'm using and only some of them are causing base quality errors, so I just want to make sure I'm interpreting that thread correctly before I re-run 521 jobs.

    -Ryan

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Ryan,

    If you only have a few samples producing these errors and all samples were processed the same way, then it's possible that it's not the encoding that is at fault, but a few aberrant records. You should check which is your case. If it's just a few aberrant records, you may be able to use the argument for ignoring misencoded quality scores. You shouldn't have to remap anything.

  • ryanfriedman22ryanfriedman22 Washington University in St. LouisMember

    I'm getting this error in at least 164 of my samples, so I don't think it's that. What else could it be if I don't need to remap anything? After I run Tophat on each sample, I use Picard to add read group information and then Samtools to sort and remove duplicate reads. I'd like to keep this if possible since I'm using Samtools to do the same with DNA-seq data.

    I should mention that these aren't paired-end reads, but I can't imagine that would make a difference in any of this.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Best thing to do at this point is check the encoding and check the distribution of base quality values. Base quals are determined by the sequencer software, and should not have been touched by any of the subsequent tools until GATK.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah yes you should only use that argument on samples that throw the error.

Sign In or Register to comment.