number of aligned reads inconsistency between sam2tsv and samtools mpileup

Huanle commented 3 years ago

Subject of the issue

Hi @lindenb,

Thanks very much for developing sam2tsv and other bioinformatics java tool-kits.

I have a question regarding the inconsistency in terms of the number of mapped reads I observed in one of my sample.

My environment

the version of sam2tsv I am using: 0558422;
version of java: openjdk 11.0.2 2019-01-15; OpenJDK Runtime Environment 18.9 (build 11.0.2+9); OpenJDK 64-Bit Server VM 18.9 (build 11.0.2+9, mixed mode)
samtools: Version: 1.9 (using htslib 1.9)

Steps to reproduce

run sam2tsv

samtools view -h test.filt.bam | java -jar sam2tsv.jar -R ref.fa | awk '($7 ==1 && $NF !~ /S/) || /READ/'
#READ_NAME      FLAG    CHROM   READ_POS        BASE    QUAL    REF_POS REF     OP
09fa5020-9343-47d6-9f75-56b7c8be0ba8    16      YAL005C 146     A       @       1       A       M
415554f1-d76b-4f6d-8ade-c01ef795e926    16      YAL005C 153     A       >       1       A       M
db9f9880-117e-4239-9bb3-ce22ae59de2e    16      YAL005C 183     A       C       1       A       M
aaed3e60-38a1-4e7a-a469-87aa0f735f95    16      YAL005C 18      A       /       1       A       M
45955b97-24c0-45c3-85a2-f5d7661d4ba2    16      YAL005C 170     A       =       1       A       M
45777d56-1a8e-4a3c-a28c-90dd17ef829c    16      YAL005C 188     A       B       1       A       M
f51a1929-b5a4-443e-8623-4bad8f01bbf5    16      YAL005C 77      A       <       1       A       M
3d57b5c4-5877-4061-b14f-cdaf64a32c11    16      YAL005C 18      A       -       1       A       M
dd617aec-b243-4af4-b47b-cc1edf099568    16      YAL005C 191     A       C       1       A       M
b2fc369a-57d7-40ec-8ee5-83e37b113f25    16      YAL005C 206     A       .       1       A       M
aae4689b-fe64-4553-b139-59944034a676    16      YAL005C 140     A       5       1       A       M
4c4676a5-4c2f-471c-a5d0-0939f25143b2    16      YAL005C 171     A       C       1       A       M
c2ed1eb7-6d19-43d6-95a3-a2d4af3751f6    16      YAL005C 16      A       $       1       A       M
dd821987-409f-432b-b43c-6f28ceb886f6    16      YAL005C 83      A       8       1       A       M
51d35b1c-abf7-4b9d-a7dd-2ea2e5530d5e    16      YAL007C 21      A       ;       1       A       M
9fc9a0fd-92ea-4e15-aebc-7d82eee78c84    16      YAL007C 33      A       A       1       A       M
5c9b0e6a-a3ea-4287-8970-6683f051bb96    16      YAL007C 1       A       7       1       A       M
ea638c17-fa28-43e1-9cc3-2cd22a5fc291    16      YAL007C 95      A       1       1       A       M
9e96a945-460b-47a2-8a8b-efed62d8602c    16      YAL007C 33      A       ;       1       A       M
7633e6dd-47eb-4c1e-adb8-83263672ac37    16      YAL010C 137     A       D       1       A       M

So according to this, there are 20 reads mapped to YAL010C :1.

run samtools mpileup

samtools  mpileup -6 -f ref.fa -q 0 -Q 0 test.filt.bam|head -n1
[mpileup] 1 samples in 1 input files
YAL005C 1       A       14      ^],^],^],^],^],^T,^],^],^],^],^N,^],^\,^],      !!!!!!!!!!!!!!

It shows that there are 14 reads mapped to the position.

What I observed from IGV the number of reads mapped to this position is also 14.

Expected behaviour

I expect no inconsistency in terms of mapped reads between the two.

Actual behaviour

I do not know what is the underlying cause of this. Can you help me figure it out? I attached the relevant files data.tar.gz used in the above commands. Thanks a lot.

lindenb commented 3 years ago

$ java -jar dist/sam2tsv.jar ~/tmp/data/test.filt.bam | awk '($4=="YAL005C" && $8==1)'  | cat -n
     1  09fa5020-9343-47d6-9f75-56b7c8be0ba8    16  60  YAL005C 146 A   @   1   N   M
     2  415554f1-d76b-4f6d-8ade-c01ef795e926    16  60  YAL005C 153 A   >   1   N   M
     3  db9f9880-117e-4239-9bb3-ce22ae59de2e    16  60  YAL005C 183 A   C   1   N   M
     4  aaed3e60-38a1-4e7a-a469-87aa0f735f95    16  60  YAL005C 18  A   /   1   N   M
     5  45955b97-24c0-45c3-85a2-f5d7661d4ba2    16  60  YAL005C 170 A   =   1   N   M
     6  45777d56-1a8e-4a3c-a28c-90dd17ef829c    16  51  YAL005C 188 A   B   1   N   M
     7  f51a1929-b5a4-443e-8623-4bad8f01bbf5    16  60  YAL005C 77  A   <   1   N   M
     8  3d57b5c4-5877-4061-b14f-cdaf64a32c11    16  60  YAL005C 18  A   -   1   N   M
     9  dd617aec-b243-4af4-b47b-cc1edf099568    16  60  YAL005C 191 A   C   1   N   M
    10  b2fc369a-57d7-40ec-8ee5-83e37b113f25    16  60  YAL005C 206 A   .   1   N   M
    11  aae4689b-fe64-4553-b139-59944034a676    16  45  YAL005C 140 A   5   1   N   M
    12  4c4676a5-4c2f-471c-a5d0-0939f25143b2    16  60  YAL005C 171 A   C   1   N   M
    13  c2ed1eb7-6d19-43d6-95a3-a2d4af3751f6    16  59  YAL005C 16  A   $   1   N   M
    14  dd821987-409f-432b-b43c-6f28ceb886f6    16  60  YAL005C 83  A   8   1   N   M
    15  0c366064-d758-41fe-990b-8f67a7b346b2    16  60  YAL005C 131 T   %   1   N   S
    16  39773790-1ffc-4461-8073-f882bf5bf6b3    16  60  YAL005C 172 A   ,   1   N   S
    17  a328222d-831f-44d1-a6e0-8cf7820b37e6    16  60  YAL005C 172 T   %   1   N   S
    18  9e1338a5-6aca-4fc4-9213-ea25ab4c38ff    16  60  YAL005C 20  G   -   1   N   S

I see 18 reads at this position. 4 are soft clipped operator==S. 18-4=14

Huanle commented 3 years ago

Hi @lindenb , Thanks a lot for your quick reply. I assumed the alignment file has only one reference. How careless i was!

lindenb / jvarkit