We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

UnmappedReadFilter

AsJAsJ JPMember ✭
in Ask the GATK team

I would like to know how UnmappedReadFilter identify umapped reads from bam files?
one more thing, the read count is much higher than number of reads i have in the original sam file.So, how mutect2 counts the number of reads?

Mutect2 log:
MicroScheduler - 123005 reads were filtered out during the traversal out of approximately 767207040 total reads
but the original sam file has 30,777,253 reads
Any idea?
Thank you

Tagged:
· Share on Twitter Share on Google+

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Unmapped reads are tagged with a flag by the aligner. We just read the flag. See https://software.broadinstitute.org/gatk/blog?id=7019 for an explanation of SAM flags.

    How are you counting the number of reads in your file?

    · Share on Twitter Share on Google+
  • AsJAsJ JPMember ✭

    Thnks for the reply. I think there is something strange in counting the number of reads.
    I counted the reads using wc -l file.sam and samtools flagstat
    I got same number of reads using both wc -l file.sam and samtools flagstat

    ====

    samtools flagstat file.bam
    30777253+ 0 in total (QC-passed reads + QC-failed reads)
    ....
    ....

    =====

    By any means the MuTect2 count should never be higher than wc -l. It can be lower due to filteration of reads,..etc but not higher.
    Do you know why i got this large number of reads ?

    Thank you!!

    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    That is indeed very odd. What happens if you just run on a small region, does it give you correct counts, or are they also inflated?

    · Share on Twitter Share on Google+
  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @AsJ
    Hi,

    Can you please try using CountReads with -drf MalformedReadFilter -drf BadCigarFilter? I wonder if there is some discrepancy between Samtools and GATK.

    Thanks,
    Sheila

    · Share on Twitter Share on Google+
  • AsJAsJ JPMember ✭

    Hi,
    Thank you for the answer.
    I have tried running CountReads and i got the same number of reads as samtools or wc -l file.sam.
    This means that the log message by Mutect2 was not correct.

    ------

    java -jar GenomeAnalysisTK.jar -R reference.fa -T CountReads -I sortedfile.bam

    Output:
    INFO ...... CountReads - CountReads counted 30777253 reads in the traversal

    · Share on Twitter Share on Google+

    Issue · Github
    by Sheila

    Issue Number
    1830
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @AsJ
    Hi,

    Can you confirm this happens with the latest version/nightly build? I will check with the team what they think.

    -Sheila

    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @AsJ, our developer pointed out that the program is probably counting the combined read counts from both the tumor and the normal bam. Can you confirm you are running this with both tumor and normal bams?

    · Share on Twitter Share on Google+
  • AsJAsJ JPMember ✭

    Hi VdAuwera,

    Yes, I am running this with both tumor and normal bams.

    · Share on Twitter Share on Google+
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ok, then that makes sense. We'll see if we can make that clearer in the output.

    · Share on Twitter Share on Google+
Sign In or Register to comment.