This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Python or shell scripting code ro count individual mapped reads of amplicons

I made 100 multiplex amplicon library and sequenced. Aligned it with bwa and converted to sorted-indexed bam file following the instruction of biostars.org/p/41951/ Using "flagstat" option in samtools, I can count the total reads and total mapped reads. To count individual reads, I use "samtools view input.bam chr1:1234-2345 | wc -l" and got the result.

However, my library is composed of 100 amplicons, which means that I need to do this thing 100 times if I do one by one manually. Can anyone kindly help me to share python or shell script to do this job automatically? I have basic knowledge on both tools. Thanks in advance!

sequencing next-gen alignment
• 2,783 views • link
updated by Zaag Zaag • written by sp • sp

3 answers

From the command line:

for i in $(ls *.bam) ; 
do
samtools view $i chr1:1234-2345 | wc -l > $i.individualcount.txt ;
samtools flagstat $i > $i.flagstat.txt
echo "Done with $i"
done ;

should give you for each bam:

  • name.bam.individual count.txt
  • name.bam.flagstat.txt
• 0 views • link
written by Zaag Zaag

FYI, the top line can just be for i in *bam. Less buttons, works the same.

• 0 views • link
written by Daniel Daniel

Thanks Zaag! It does not work for my purpose though. I need iterate coordinates on same bam file. I was digging around and found "bedtools multicov" did the job! I also want to try NGSUtils, but failed to install on my ubuntu.

• 0 views • link
written by sp • sp

Try samtools idxstats.

• 0 views • link
written by mastal511 mastal511

Can you show me a working example? I am stupid biologist.

• 0 views • link
written by sp • sp

You can take a look at the python module pysam pileup to get counts.

• 0 views • link
written by WouterDeCoster WouterDeCoster

Log in to answer this question.

More posts like this