Question: Counting number of reads that start at a given position using bbmap pileup.sh
0
gravatar for ashafik
15 days ago by
ashafik0
ashafik0 wrote:

I'm trying to get the number of reads that start at a given position. I'm using bbmap pileup.sh. This is the script:

./pileup.sh in=sorted.bam basecov=coverage.txt startcov=t

How can I get coordinate information. For example, chrM 1033: 10 (so 10 reads start at that position).

Thank you

alignment • 94 views
ADD COMMENTlink modified 15 days ago by Pierre Lindenbaum133k • written 15 days ago by ashafik0
0
gravatar for GenoMax
15 days ago by
GenoMax95k
United States
GenoMax95k wrote:

Looking at inline help that command should generate the information you are looking for. Is it not doing that? I will have to test and check.

Edit: I checked this and the command should work. You will need to provide more info.

Also take a look at mosdepth (LINK) as a fast alternative.

ADD COMMENTlink modified 15 days ago • written 15 days ago by GenoMax95k

Thanks for the reply. This is the output I get:

#RefName    Pos Coverage
1   0   0
1   1   0
1   2   0
1   3   0
1   4   0
1   5   0
1   6   0
1   7   0
1   8   0
1   9   0
1   10  0
1   11  0
1   12  0
1   13  0
1   14  0
1   15  0
1   16  0
1   17  0
1   18  0
1   19  0
1   20  0
1   21  0
1   22  0
1   23  0
1   24  0
1   25  0
ADD REPLYlink modified 15 days ago by GenoMax95k • written 15 days ago by ashafik0

Since you used

startcov=t          Only track start positions of reads.

should produce counts on reads that start at that position. Is this not working?

ADD REPLYlink written 15 days ago by GenoMax95k

It's not working for me. I just get the output that I showed. The pos column goes to 187588718 and counting. The other two columns don't change.

ADD REPLYlink written 15 days ago by ashafik0

Are you sure sorted.bam has alignments? I tested pileup.sh with a file I had and was able to get counts. Can you show output of samtools view sorted.bam | head -6?

ADD REPLYlink modified 15 days ago • written 15 days ago by GenoMax95k

What do you mean by alignments? Here is the output:

MG01HX02:885:H3YW7CCX2:2:2224:20273:34465   256 1   3000472 1   33M *   0   0   TTTCCACTTGGTTGATTTCAGCTCTGAGTTTGA   AJFFJAJFAFJFJFFJAAAJAFAFFJJJJJFJJ   AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:33 YT:Z:UU NH:i:10
MG01HX02:885:H3YW7CCX2:2:2111:32299:33076   256 1   3014889 1   62M *   0   0   AGTTTGCAAGTCCAATGGGCCTCTATTTGCAGTGATGGCCGACTAGGCCATCTTTTGATACA  JJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:62 YT:Z:UU NH:i:10
MG01HX02:885:H3YW7CCX2:2:2201:20953:1924    272 1   3014930 1   54M *   0   0   ACTAGGCCATCTTTTGATACATATGCAGCTAGAGACAAGAGCTCCGGGGTACTA  F-JFJA7-JFFJFFA77FJ7JJJ<JJJFJFJJFJFAFJF-JJF<FJF-JJJFJJ  AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:54 YT:Z:UU NH:i:10
MG01HX02:885:H3YW7CCX2:2:1118:16701:27433   256 1   3014932 1   2S64M   *   0   0   AATAGGCCATCTTTTGATACATATGCAGCTAGAGACAAGAGCTCCGGGGTACTAGTTAGTTCATAT  JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ  AS:i:-4 ZS:i:-4 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:64 YT:Z:UU NH:i:2
MG01HX02:885:H3YW7CCX2:2:1123:7101:72789    256 1   3016650 1   43M *   0   0   TATCCTTGAGAAGAGTTTTTGCTATCCTCGTTTTTTTGTTATT JJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJ AS:i:0  ZS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:43 YT:Z:UU NH:i:10
MG01HX02:885:H3YW7CCX2:2:1116:23490:48652   0   1   3018672 60  100M    *   0   0   TTTTGTTTTAGGATAAAATGTTCTGTAGATATCTGTCAAGTCCATTTGTTTCATCACTTCTGTTAGTTTCACTGTGTCCCTGTTTAGTTTCTGTTTCCAT    JJFJJJJJJJJFJJJJFJJJJJJFJFJJFJJJJJJJJFJFJJJJJJJJJJJJJJJJFJFJJJJJJJJJJJJFJJJJJJ<JJ<AJJJJJJJJJFJJJJJJJ    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:100    YT:Z:UU NH:i:1
ADD REPLYlink modified 15 days ago by GenoMax95k • written 15 days ago by ashafik0

The alignments look fine. Looking at this you should get something in your coverage file around 3014930 in second column.

If you do grep 3014930 coverage.txt what do you see?

ADD REPLYlink modified 15 days ago • written 15 days ago by GenoMax95k

Yes I see counts now. Thank you for your time.

ADD REPLYlink written 15 days ago by ashafik0
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: 988 users visited in the last hour
_