Question: Python or shell scripting code ro count individual mapped reads of amplicons
0
gravatar for sp
3.1 years ago by
sp20
sp20 wrote:

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 alignment next-gen • 920 views
ADD COMMENTlink modified 3.0 years ago by Zaag720 • written 3.1 years ago by sp20
2
gravatar for Zaag
3.0 years ago by
Zaag720
Amsterdam
Zaag720 wrote:

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
ADD COMMENTlink modified 3.0 years ago • written 3.0 years ago by Zaag720
1

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

ADD REPLYlink written 3.0 years ago by Daniel3.7k

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.

ADD REPLYlink written 3.0 years ago by sp20
0
gravatar for mastal511
3.1 years ago by
mastal5112.0k
mastal5112.0k wrote:

Try samtools idxstats.

ADD COMMENTlink written 3.1 years ago by mastal5112.0k

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

ADD REPLYlink written 3.1 years ago by sp20
0
gravatar for WouterDeCoster
3.1 years ago by
Belgium
WouterDeCoster38k wrote:

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

ADD COMMENTlink written 3.1 years ago by WouterDeCoster38k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1748 users visited in the last hour