UnmappedReadFilter

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:

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    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?

  • 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!!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    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?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @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

  • 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

    Issue · Github
    by Sheila

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

    @AsJ
    Hi,

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

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    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?

  • AsJAsJ JPMember

    Hi VdAuwera,

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

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

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

Sign In or Register to comment.