Percent reads mapped in Picard

Question

[phylofun]

What tool in Picard will tell me the number of reads that mapped to my reference sequences? I would like to use this program (not SAMtools) to identify the percentage of mapped reads (i.e. number of reads mapped/number of total reads).

[phylofun]

I'm getting this error when I use that tool:

ERROR: Record 1, Read name M00366:9:000000000-A89UT:1:1103:22880:8023, Alignment start should != 0 because reference name != *.

[finswimmer]

Hello,

which version of picard are you using?

fin swimmer

[phylofun]

Latest release of Picard: version 2.18.27

[finswimmer]

Ok, let's have a look at the reads. Please post the output of:

$ picard ViewSam I=input.bam ALIGNMENT_STATUS=All PF_STATUS=All | grep "M00366:9:000000000-A89UT:1:1103:22880:8023"

(Why is picard making this so complicated? Same thing using samtools: $ samtools view -h input.bam|grep "M00366:9:000000000-A89UT:1:1103:22880:8023")

[phylofun]

I get the same error as before if I use picard ViewSam. However, if I use samtools view, I get:

M00366:9:000000000-A89UT:1:1103:22880:8023  129 uce-6534_wj_3_333991    0   37  251M    uce-7297_wj_3_333991    684 0   TCTACCGTAAATTTGAGACTTGGAAGAGAAACCTTTCCTCATTTTTCTGTGTGTAACGTTTGTAAAATATTGAATAGGAAAACATCAAAGCCAGGGGGAAATATAAGTAACTCACTGCTTACGTTTTAATTCCTCACCTTTCCCCTCTGGGTGAGTTCCCCCTCCTTTTTCCCCCCTGAATTTAAAATGCTGGAAGGAAAAGCAACTGAGGTCCAAGTTTCAGATGCTGAATTCTCTGCTAATTGAAATTA 3>A3AFBAABBFBGGGGGGGGGGHHFHHGHHHHHHHHHHHHGHHHHGGEBGGHHGHGGHHHGHHGHFHHFHHEHFFFGHHHHHHHFFHHFHHHGGCGCGDHGHHHHHGFHHFHHFFHHHFHBGGFBFGHFHHHHGHHHHHHHBGHHGHHGHFGHHHFHGHFEGGGHGHGHHGHHGGGAFFFEHHHHHHHEFG/CGFGHHHHGHHFCCGFFGGHFEGGGFFFFFFFGGGGGGGGGGGF0=FGGGGFE0FFF0 X0:i:1  X1:i:0  MD:Z:0  PG:Z:MarkDuplicates RG:Z:wj_3_333991    XG:i:0  AM:i:37 NM:i:0  SM:i:37 XM:i:1  XO:i:0  XT:A:U
M00366:9:000000000-A89UT:1:1103:22880:8023  65  uce-7297_wj_3_333991    684 37  250M    uce-6534_wj_3_333991    0   0   ATTCAGCGTCTCTGCTCTGCACAAACGGTACACACACGATCCCGAGCATGTGCTACTAAGGAGCATCCGTTCCCTCTGTCAGTTTGAAAATGCATTGATCTCTGTTTTTCTGACTACAAAAAGACATGTCTTGAAGGCACAGAGTGAGAGCTACATATCTCTTCAGCAACCATATATAGAAACATGCCCAGTAATTTCAATTAGCAGAGAATTCAGCATCTGAAACTTGGACCTCAGTTGCTTTTCCTTC  @BBBBFFB@ABFB5GBGGGGGGHHHFEEAGGFGHGGGGGHHGHCEGFGHFHGHGHHHHFGFFHHHHGGGFFFFHGHHFGGHGGHHGHHFFHHHHFHHHHHHGHBGD4FFHF4FGFHHHHHHHGGHHHHHHHE3?GBGFHHHBFHF3G2@@<@GHHHHGFHHHGHG1FGHFHFBDGBGHHEHGHHHHHFHFE1>=GGHHHHGGHF/GDG<GGHHEHHGHHHHFGEHFGFG0CHHHGGFBFBGGGEBFCFFF  X0:i:1  X1:i:0  MD:Z:147A102    PG:Z:MarkDuplicates RG:Z:wj_3_333991    XG:i:0  AM:i:37 NM:i:1  SM:i:37 XM:i:1  XO:i:0  XT:A:U

[phylofun]

I just used samtools flagstat command and was able to get % reads mapped. However, should I expect to get the same or different result using CollectAlignmentSummaryMetrics?

[finswimmer]

Hello again,

you should expect the same result. But your bam lines are somewhat weird. For the first read the POS is 0. This would be valid, if the read is unmapped but the mate could be mapped. But this isn't given by the FLAG values is column 2.

How have you done the alignment?

fin swimmer

[phylofun]

I did the alignment with bwa aln and bwa sampe (a long time ago)

[finswimmer]

Hello again,

to fix the wrong bam file you can do this:

Make sure you have a current version of samtools

$ samtools sort -n input.bam | samtools fixmate - - | samtools sort > fixed.bam

You will receive the warning:

[W::sam_parse1] mapped query cannot have zero coordinate; treated as unmapped

You don't need to care about it, because this is what we expected.

fin swimmer

Answer 1

[finswimmer]

CollectAlignmentSummaryMetrics should give you the information you need.

But why don't you like to use samtools?