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

Question

[ashafik]

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

Answer 1

[GenoMax]

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.

[ashafik]

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

[GenoMax]

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?

[ashafik]

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.

[GenoMax]

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?

[ashafik]

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

[GenoMax]

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?

[ashafik]

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