Indel Realigner - Exception counting mismatches AND Missing read group or read group not defined

Hello,
For the past several months, I have consistently been receiving 2 error messages depending on which version of GATK I am using.

If I use an older version (such as 1.0.5588), I am able to successfully generate the intervals file using RealignerTargetCreator with no problems. When I then run IndelRealigner, I receive messages for individual reads about Exception counting mismatches. I typically remove these reads by hand (even though there appears to be nothing wrong with them) and carry on. When I then move to the UnifiedGenotyper to call SNPs, I use the realigned output file from the IndelRealigner and now get the message about missing read group or read group not defined. However, the input bam file for the IndelRealigner worked just fine. But the output from the IndelRealigner now has a read group issue. I don't understand

If I use a newer version (2.3.9 or 2.4.9), I immediately get the message about missing read group or read group not defined when I first run RealignerTargetCreator. When I look at the headers from the bam file, everything seems correct.

I have double checked my reference file, its index, my bam file and the index bam file. Nothing seems to have an error.

Any suggestions on how to move forward? I have been troubleshooting for months with other colleagues and have gotten nowhere. Thank you for your time and recommendations.

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Lee,

    In future, please do not hesitate to post questions or issues here when you encounter them -- we can probably save you a lot of time!

    First, you should definitely not be using versions 1.x, as they have many bugs that we have fixed since, and we have made many improvements to the tools as well. Also, we cannot provide support for such old versions.

    Second, regarding your read group issue: can you confirm that your reads actually have read group information? The GATK requires read group info to be attached to every single reads, otherwise the tools will not function properly. If your reads do not have read groups attached, you can add that information using Picard's AddOrReplaceReadGroup tool. If you're not sure, please post the header of your bam file and I'll have a look.

  • LEELEE Member

    Hi Geraldine,
    Thanks for the response! I am attaching the header of the bam file. I think it is okay which is why I am confused. But maybe I am missing something. Thanks so much!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    OK, you have RG tags in the header, that's a good sign. But it's possible that some of the reads in the bam file are missing a read group. To find out if that is the case, you can try this command:

    samtools view mybam.bam | head -n 10000 | grep -v RG: | wc -l
    

    This command will tell you the number of reads that are missing a read group.

    Just replace "mybam.bam" by your actual bam file's name. Note that this makes a couple of assumptions: that a read without a read group occurs in the first 10000 reads, and the string "RG:" occurs in no other tags.

  • LEELEE Member

    Hi Geraldine,
    Thanks! The output from samtools was 0 reads. Ideas for a next step? I really appreciate this!

  • LEELEE Member

    I also increased the n to 100000 and 1000000 and still found 0 reads that were missing a read group.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm. Can you please post your original RealignerTargetCreator command here for me, then run it again, this time with -l DEBUG added? This will output a lot of information to the console including genome positions; so when it errors out it should give us an idea of the position at which the bad reads occur (if they do). I'll need to look at the last stack of info that gets output before the error.

  • LEELEE Member

    And n = 10000000, still 0 reads with a missing read group. Just trying to give you some more information to work with. Thanks!

  • LEELEE Member

    java -XMx2g -jar /programs/GenomeAnalysisTK-2.4.9/GenomeAnalysisTK.jar -T RealignerTargetCreator -R clustered_transcriptome_200.fasta -o merged_output.intervals -I merged.bam --minReadsAtLocus 3

  • LEELEE Member

    Here is output from the debugging step.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Did you try attaching a file? I can't see it. It may have been refused by the forum if it was too large. Feel free to just paste the last page or so's worth into a comment.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, there we go: it's actually in the error message, read HWI-ST1085:123:D1AM3ACXX:3:2216:18167:96607 is the culprit. And it looks like it's mapped to either contig28871 or the one after that. Since this is a non-standard reference I'm not sure how far in that would be but this should at least help you narrow it down. You should try to track down where that read comes from. Based on your filenames it looks like you merged some data together from different sources; maybe one of them was not processed properly.

  • LEELEE Member

    Should I remove the read? What if there are multiple reads that are a problem? I have also successfully made it past this point with other samples and received a message about the file as a whole missing BAM headers, without a specific read mentioned.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Put simply, the GATK will refuse to proceed anytime it encounters a read without a read group, so this needs to be fixed. If it is a single read that somehow "escaped" getting tagged, you can safely remove it from the file. But if there are more of these, then it's going to be very painful to identify and remove them one by one. It really would be better to track down where the problem comes from. You may need to examine the original data files. If you have separate unmerged bams you can use the same samtools command from earlier on those and see if it's just one of them that didn't get tagged. If so you can add a tag using Picard tools.

  • LEELEE Member

    Hi Geraldine,
    I double checked the original unmerged files and all of the reads have read groups. Additionally, in the merged file, the read that GATK is indicating as the first problem read (I am assuming there will be more from previous experience), has a read group. So I don't know how to fix a problem of adding a read group when a read group is present. I have been trying to track down the source of this problem for months with no success. I don't know how to fix a problem of an error message indicating no read group for a read when there is a read group present for all reads.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I understand your frustration, but we'll get to the bottom of this. Can you just tell me what is the read group indicated for that read? And is that read group represented in the bam header of the merged file?

    (feel free to just post the full line from the bam for that "trouble" read)

  • LEELEE Member

    Apologizing in advance for this question - it shows what a novice I am. I know how to remove a read manually from a SAM file. Is there a way to remove a specific read from the BAM file?

    From the BAM file:

    @[email protected]>;<=BBD?CC XT:A:U NM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:68 RG:Z:CS1_pooled

    Additional information - this is the very first read in the file.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    OK, this is just weird, because that read group is definitely in your header. Can you please make a bam snippet including that read and upload it to our FTP? Instructions here: http://www.broadinstitute.org/gatk/guide/article?id=1894

    I need to test this locally. If it uses a custom reference I'll also need a copy of that.

Sign In or Register to comment.