To celebrate the release of GATK 4.0, we are giving away free credits for running the GATK4 Best Practices pipelines in FireCloud, our secure online analysis portal. It’s first come first serve, so sign up now to claim your free credits worth $250. Sponsored by Google Cloud. Learn more at https://software.broadinstitute.org/firecloud/documentation/freecredits

PICARD AlignmentSummaryMetrics

Hi
I used the following command to gather summary metrics for my bam file generated via bowtie2 (tophat to be specific):

java -jar /usr/share/picard-tools-1.136/picard.jar CollectAlignmentSummaryMetrics INPUT=Sample_DY10.tophat.bam OUTPUT=tmpmetrics/alignmentmetrics R=/mnt/storage/ref_genome/Homo_sapiens/UCSC_hg19/UCSC/hg19/Sequence/Bowtie2Index/genome.fa

The output file is attached.
The question I have is that the metrics PF_HQ_MEDIAN_MISMATCHES has a very high number (66). When I look at NM tag in the bam file, I see that the median is 1 with max NM = 2
I am wondering how this number is calculated by PICARD.

Any help is appreciated.

Issue · Github
by Sheila

Issue Number
1315
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
sooheelee

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @newbie16
    Hi,

    It says the metric is "The median number of mismatches versus the reference sequence in reads that were aligned to the reference at high quality (i.e. PF_HQ_ALIGNED READS) in this article. However, I am not sure what an appropriate number is for the metric. I will check with the team and get back to you.

    -Sheila

  • shleeshlee CambridgeMember, Broadie, Moderator

    Hi @newbie16,

    We've narrowed down your excessive PF_HQ_MEDIAN_MISMATCHES to three possibilities. Either CIGAR string S bases (softclips) are counted towards mismatches, or CIGAR string N bases (reference-skip bases, e.g. for intronic sequences), or both. Considering these types of bases in your alignment records, does your excessive median mismatches make sense?

    In terms of reads for which this metric is calculated, these I believe have to have MAPQ > 20 (therefore must be aligned) and cannot be supplementary. The tool takes alignment blocks in the record, defined by the CIGAR string, and iterates over each of them to add to the mismatch count by directly comparing the base to the reference. Comparisons are case-insensitive.

  • shleeshlee CambridgeMember, Broadie, Moderator

    @newbie16,

    Someone from the team informs me that the RNA samples have a PF_HQ_MEDIAN_MISMATCHES value typically around 0-2. So what I wrote above may be wrong. Can you post some of your alignment records so we can take a look at the SAM flag values, CIGAR string, etc?

  • Hi
    Thanks for looking into it. I have uploaded a sample bam on google drive with below link. The PF_HQ_MEDIAN_MISMATCHES value for this file was 66.

    https://drive.google.com/open?id=0B2tk2ztVP7NIZHQ1V1pFNGxJeU0

  • Hi @shlee

    The SplitNCigarReads program is emitting following error:

    ERROR MESSAGE: Unsupported CIGAR operator N in read HISEQ:262:C99J2ACXX:8:2206:7109:2528 at chr1:4776766. If you are working with RNA-Seq data, see https://www.broadinstitute.org/gatk/guide/article?id=3891 for guidance. If you choose to disregard those instructions, or for other uses, you have the option of either filtering out all reads with operator N in their CIGAR string (add --filter_reads_with_N_cigar to your command line) or overriding this check (add -U ALLOW_N_CIGAR_READS to your command line). Notice however that the latter is unsupported, so if you use it and encounter any problems, the GATK support team not be able to help you.

    The command I used is this:
    java -jar ~/software/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar -T SplitNCigarReads -R /mnt/storage/ref_genome/mouse_mm10/mouse_mm10/genome.fa -I Sample_14011001.tophat.sub0p005.reorder.bam -o Sample_14011001.tophat.sub0p005.reorder.NoNs.bam

    The aligned read (from sam file) mentioned in the above error is below:

    @DFFFDHDHFBFFG>;EHHGJEGIFHGBGC@FH@GEGBFGHGGGGG@CGBCGGIEGG)=(=@=CG=C>;EEEHBDECCBD?CDCCBBD>A>4:AC<?AA> AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1 RG:Z:Sample_14011001
    HISEQ:262:C99J2ACXX:8:2206:7109:2528 83 chr1 4777618 50 31M4919N69M = 4776766 -5871 TATTAATTTTTGCTTGAAAAGTATCAGCACCCTCTTCAACCAGCTGGACTCCATAATCCCTCTTAAGCGGCTGGATGGTCACACCTCTCCCATTCACAAG @DCCC>;DCA@A;;(CED@CCFCEAHGIIGIFF=HFJIIIHF:IEIFGGGBGIIIHFGDGIGHGHFCFCJIHGBFHFHFFFFFCCC AS:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:100 NM:i:0 XS:A:- NH:i:1 RG:Z:Sample_14011001

    Could you please help

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @newbie16
    Hi,

    In this case, it is okay to use -U ALLOW_N_CIGAR_READS. We added a note in that article the error message points you to :smile:

    -Sheila

  • Thanks @shlee and @Sheila
    Once I got rid of N's the PF_HQ_MEDIAN_MISMATCHES looks ok, i.e. 0

Sign In or Register to comment.