Counting number of reads that start at a given position using bbmap pileup.sh

15 days ago by

GenoMax ♦ 95k

United States

GenoMax ♦ 95k 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 COMMENT • link modified 15 days ago • written 15 days ago by GenoMax ♦ 95k

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 REPLY • link modified 15 days ago by GenoMax ♦ 95k • written 15 days ago by ashafik • 0

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 REPLY • link written 15 days ago by GenoMax ♦ 95k

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 REPLY • link written 15 days ago by ashafik • 0

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 REPLY • link modified 15 days ago • written 15 days ago by GenoMax ♦ 95k

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 REPLY • link modified 15 days ago by GenoMax ♦ 95k • written 15 days ago by ashafik • 0

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 REPLY • link modified 15 days ago • written 15 days ago by GenoMax ♦ 95k

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

ADD REPLY • link written 15 days ago by ashafik • 0

Similar posts • Search »