Python or shell scripting code ro count individual mapped reads of amplicons

Question

[sp]

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!

Answer 1

[Zaag]

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

[Daniel]

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

[sp]

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.

Answer 2

[mastal511]

Try samtools idxstats.

[sp]

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

Answer 3

[WouterDeCoster]

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