This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Percent reads mapped in Picard

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).

reads mapping next-generation sequencing picard
• 3,313 views • link
written by phylofun 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 != *.

• 0 views • link
written by phylofun phylofun

Hello,

which version of picard are you using?

fin swimmer

• 0 views • link
written by finswimmer finswimmer

Latest release of Picard: version 2.18.27

• 0 views • link
written by phylofun phylofun

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")

• 0 views • link
written by finswimmer finswimmer

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
• 0 views • link
written by phylofun 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?

• 0 views • link
written by phylofun phylofun

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

• 0 views • link
written by finswimmer finswimmer

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

• 0 views • link
written by phylofun phylofun

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

• 0 views • link
written by finswimmer finswimmer

1 answer

CollectAlignmentSummaryMetrics should give you the information you need.

But why don't you like to use samtools?

• 0 views • link
written by finswimmer finswimmer

Log in to answer this question.

More posts like this