Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
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.
Attention:
We will be out of the office for a Broad Institute event from Dec 10th to Dec 11th 2019. We will be back to monitor the GATK forum on Dec 12th 2019. In the meantime we encourage you to help out other community members with their queries.
Thank you for your patience!

Issues with Base Quality Score Recalibration

raffsteinraffstein SwitzerlandMember

Hi, I am trying to do a base quality score recalibration of my BAM files. I am working with individual based RADseq data where the reads were cleaned and aligned against a reference genome using bowtie2. The reference genome comes from a relatively closely related taxon (~25 Mio years of divergence), however the species in this taxonomic family underwent a recent genome duplication, hence I am running bowtie2 with the -k option (which retains multiple alignments rather than picking randomly one if multiple matches of equal quality occur) and only retain the uniquely aligned reads for the downstream analyses. The downside to this is, that the MAPQ values become artificial.

After generating a recalibration data file, I run GATK PrintReads (version 3.1-1-g07a4bf8) as follow:

java -Xmx16G -jar GenomeAnalysisTK.jar -T PrintReads -R Ref_genome.fasta -I Individual.bam -BQSR Library_recaldata.grp -o Individual.rc.bam

Which yields the following error after running through 80% of the file:
# ERROR MESSAGE: Exception when processing alignment for BAM index HWI-ST1206:32:C3VGTACXX:5:1104:14976:28487 84b aligned read.

Here the read in question (the one in the middle):

HWI-ST1206:32:C3VGTACXX:5:2313:18281:8513   16  chrUn   536921741   0   84M *   0   0   CCCAGAGTGCAAAGAACCTTGGCAGGACACTAGACAACACCCTGTTGTTTTCTGCAAACATCAAAGCAGTGACTCTGTCCTGCA    [email protected]JJJJJJJJJJJJJHHHHHFF    AS:i:-29    XS:i:-35    XN:i:0  XM:i:5  XO:i:0  XG:i:0  NM:i:5  MD:Z:11T18G1A14G28C7    YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:2314:6453:46729   16  chrUn   536921741   0   84M *   0   0   CCCAGAGTGCAAAGAACCTTGGCAGGACACTAGACAACACCCTGTTGTTTTCTGCAAACATCAAAGCAGTGACTCTGTCCTGCA DDDDEEEEEEFFFEDFHHHHHJJIJJIJIIIIIGJJIFFJJJIJJJJJJJJJJJJJJJJJJJJJJJJJJJIHJJJJJHHHHHFF   AS:i:-29    XS:i:-35    XN:i:0  XM:i:5  XO:i:0  XG:i:0  NM:i:5  MD:Z:11T18G1A14G28C7    YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:1104:14976:28487  16  chrUn   536961523   255 84M *   0   0   ATTGGATTGCAGTGCAGGTACACTTTGTATGGATCGTTGTGTGTGTTTACTGTGACTTCTCTTAAGAAGACGCGTGTCCCTGCA</del>  DDDDDDDDDDDDDDDDDDDDDCEEEEEFFFFFHFFJJJJJJJJJIIHD8JJJJIHIJJJJJJIJJJJJJIJJIJJJJHHHHHFF    AS:i:-9 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:20G27T35   YT:Z:UU RG:Z:70130.satrFamilyX**
HWI-ST1206:32:C3VGTACXX:5:1104:10760:63491  16  chrUn   536961523   255 84M *   0   0   ATTGGATTGCAGTGCAGGTACACTTTGTATGGATCGTTGTGTGTGTTTTCTGTGACTTCTCTTAAGAAGACGCGTGTCCCTGCA    [email protected]JIJJJIJJJJJJJHHHHHFF    AS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:20G63  YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:1104:16496:63907  16  chrUn   536961523   255 84M *   0   0   ATTGGATTGCAGTGCAGGTACACTTTGTATGGATCGTTGTGTGTGTTTTCTGTGACTTCTCTTAAGAAGACGCGTGTCCCTGCA    DDDDDDDDDDDDDDDDDDCCCCDDEEECEEFFFFFHJJJIJJJJJJJJJJJJJIHJJJJJJJJJJIJJJJJJJJJIJHHHHHFF    AS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:20G63  YT:Z:UU RG:Z:70130.satrFamilyX

I subsequently verified my BAM files with Picard, which yields the following error
# "bin field of BAM record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned"
for 84554 reads that are almost consecutive in the same genomic region (chrUn - which is composed of sequences for which no chromosome anchoring information exist). However, the reads flagged by Picard occur about 2 million base pairs before the region in question for GATK and PrintReads seemed to have no problem with them.

For the sake of completeness, here the beginning and end of the section that was flagged by Picard.
It starts with the second of the following reads:

HWI-ST1206:32:C3VGTACXX:5:2315:1142:72568   256 chrUn   997157147   2   84M *   0   0   TGCAGGACCACAGTCTTTTCACCAGCTCCTCAAAATGGTTGAAGAAGTAGTATTTTGGGGGCTCGTATAAACTCGTATGCCGTC    FFHHHHHJJIJJJDFHHGIH[email protected]    AS:i:-50    XS:i:-50    XN:i:0  XM:i:10 XO:i:0  XG:i:0  NM:i:10 MD:Z:2T41T8C7A7G0G1A2G3A2G1 YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:1102:4043:34769   16  chrUn   997184926   255 84M *   0   0   TCGCCTACAAACAAAACCCCACCCATGCACACCTGATTGGTGCACGTAAAACTGGGAAAGAAGTTACCTCCCTGGAAGCCTGCA    [email protected]GCJJIJJJJJIJJHHHHHFF    AS:i:-12    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:41T15G26   YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:1105:2829:22242   16  chrUn   997184926   255 84M *   0   0   TCGCCTACAAACAAAACCCCACCCATGCACACCTGATTGGTGCACGTAAAACTGGGAAAGAAGTTACCTCCCTGGAAGCCTGCA    [email protected]?;DDEE?HHIIGIIHHGJIIJJJJJJJJIIJJJJJJJJJJJJJJJIJIIHE?JIGJJJJJJJJHHHHHFF    AS:i:-12    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:41T15G26   YT:Z:UU RG:Z:70130.satrFamilyX

And stops at the following read (note the second read here is ok):

HWI-ST1206:32:C3VGTACXX:5:2314:10564:44125  16  chrUn   1110931683  255 84M *   0   0   GCTAGCTGACTACTGTAGCCCCAGCCCTTATTAGCCTGCACTGTTATAAGGTGTTATGATTCACACACATTGAAATGGCCTGCA    [email protected]IJIJJJJJJJJJJHHHGHFF    AS:i:-10    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:37A39A6    YT:Z:UU RG:Z:70130.satrFamilyX
HWI-ST1206:32:C3VGTACXX:5:1101:4103:100272  256 chrUn   1110931763  0   84M *   0   0   TGCAGGACAATGCCAGTGAATGGAGCCTAATTACCCTCTGTAATGAAATGTCACATTTTTGTTTCTGGGTCCTGAATAGATGTG    FFHHHHHJJJJJJJJJGIIJJJJIJJJJJJIJJJJJJJJJJJJJJJJJJIJJJJIIJJJJIJJJJJJIHHHFHFFFFFFEEEEE    AS:i:-35    XS:i:-35    XN:i:0  XM:i:6  XO:i:0  XG:i:0  NM:i:6  MD:Z:13A34A5G4C1C11A10  YT:Z:UU RG:Z:70130.satrFamilyX

I personally cannot find something specific that would explain why PrintReads does not work. The pipeline that I use works perfectly when I use an assembly based on my RADseq data, yet not with the reference genome. Thus, I would like to ask, if you could perhaps help me.
Thanks!

Tagged:

Answers

Sign In or Register to comment.