|
Hello,
I am trying to use "samtools mpileup --no-BAQ" to retrieve variants in reads at certain positions and their original BAQ values. But it looks like "--no-BAQ" does not work because I get different numbers from what I see in BAMs and IGV. When I change from "--no-BAQ" to "--redo-BAQ" I get the same numbers. Could you please give me some insight on it? Perhaps I miss something or it is a bug? I use samtools/1.11
I run the following command: samtools mpileup -Q 0 --no-BAQ -r ${CHROM}:${POS}-${POS} --output-QNAME
--output-BP --fasta-ref $GENOME $BAM_PATH/$P_BAM
chr1 91888322 T 43 ,$G..Gg,,gg.Gg,g.g,g,gGG,G,,GGG...G,gGG.,...
?]??????!??5?????????55?5??555???5??55????? 150,148,145,131,128,128,128,118,116,113,110,105,103,91,85,84,83,82,81,80,76,63,61,61,60,53,52,50,49,49,43,39,36,34,25,17,13,13,11,9,5,3,3 ST-E00270:172:HMNWGCCXX:7:2215:20851:26536,ST-E00270:172:HMNWGCCXX:7:1205:5670:45523,ST-E00270:172:HMNWGCCXX:7:1116:12073:30246,ST-E00270:172:HMNWGCCXX:7:2209:16376:63261,ST-E00270:172:HMNWGCCXX:7:1204:11657:51008,ST-E00270:172:HMNWGCCXX:7:2205:8572:55719,ST-E00270:172:HMNWGCCXX:7:1106:21440:4526,ST-E00270:172:HMNWGCCXX:7:2122:28990:52555,ST-E00270:172:HMNWGCCXX:7:1205:5670:45523,ST-E00270:172:HMNWGCCXX:7:1201:25317:47316,ST-E00270:172:HMNWGCCXX:7:2104:24860:4825,ST-E00270:172:HMNWGCCXX:7:2114:29011:2979,ST-E00270:172:HMNWGCCXX:7:2201:5792:10714,ST-E00270:172:HMNWGCCXX:7:1121:13606:18379,ST-E00270:172:HMNWGCCXX:7:1114:31335:40407,ST-E00270:172:HMNWGCCXX:7:2206:19512:41163,ST-E00270:172:HMNWGCCXX:7:2105:21704:35115,ST-E00270:172:HMNWGCCXX:7:2206:3721:33744,ST-E00270:172:HMNWGCCXX:7:2115:22841:2417,ST-E00270:172:HMNWGCCXX:7:2114:29640:64650,ST-E00270:172:HMNWGCCXX:7:1103:20963:18502,ST-E00270:172:HMNWGCCXX:7:1208:20821:43607,ST-E00270:172:HMNWGCCXX:7:1210:23622:7673,ST-E00270:172:HMNWGCCXX:7:1214:9607:67779,ST-E00270:172:HMNWGCCXX:7:2209:29061:20717,ST-E00270:172:HMNWGCCXX:7:1216:4980:42077,ST-E00270:172:HMNWGCCXX:7:2216:21054:9695,ST-E00270:172:HMNWGCCXX:7:1102:5426:40776,ST-E00270:172:HMNWGCCXX:7:2220:21602:2487,ST-E00270:172:HMNWGCCXX:7:1202:4533:14687,ST-E00270:172:HMNWGCCXX:7:2118:14123:14248,ST-E00270:172:HMNWGCCXX:7:2103:25804:13246,ST-E00270:172:HMNWGCCXX:7:2122:5345:22282,ST-E00270:172:HMNWGCCXX:7:2213:6299:57477,ST-E00270:172:HMNWGCCXX:7:1111:19877:47351,ST-E00270:172:HMNWGCCXX:7:1115:23399:24163,ST-E00270:172:HMNWGCCXX:7:1101:19177:20331,ST-E00270:172:HMNWGCCXX:7:2202:26585:6759,ST-E00270:172:HMNWGCCXX:7:1222:22079:69098,ST-E00270:172:HMNWGCCXX:7:1120:16559:61151,ST-E00270:172:HMNWGCCXX:7:2120:23399:25358,ST-E00270:172:HMNWGCCXX:7:2206:17320:26361,ST-E00270:172:HMNWGCCXX:7:1219:20141:12191
For example, G with quality "]" at position 148 in read ST-E00270:172:HMNWGCCXX:7:1205:5670:45523.
When I grep the read from BAM, I get the quality sequence with no "]" at 148 position.
samtools view $BAM_PATH/$P_BAM | grep ST-E00270:172:HMNWGCCXX:7:1205:5670:45523
ST-E00270:172:HMNWGCCXX:7:1205:5670:45523 99 chr1 91888175 60 150M = 91888207 182 ATCCACGTGATATGTGACTTGCTCCTCCTTACCTTCTGCCATGATTGTGAGGTCTCCCCAGCCATGTGGAACTGTAAGTCCAATAAACCTTTCTTTTGTAAATTGCCCAGTCTCGGGTATGTCTTTATCAGCAGCATGAAAACAGACGAA
??????5???????????????????????????????????????????????????????????????????????????????????????????????????????????5??????????????????????????????????? MC:Z:150M AS:i:142 XS:i:98 MD:Z:93C53T2 NM:i:2 RG:Z:0
ST-E00270:172:HMNWGCCXX:7:1205:5670:45523 147 chr1 91888207 60 150M = 91888175 -182 CTTCTGCCATGATTGTGAGGTCTCCCCAGCCATGTGGAACTGTAAGTCCAATAAACCTTTCTTTTGTAAATTGCCCAGTCTCGGGTATGTCTTTATCAGCAGCATGAAAACAGACGAATACAGTATTTATGCCTTGTTAACAAGTTTTAG
??????????5+?????????????????????????????????????????????????????????????????????5????????????????????????????????5??????????????????????????????????? MC:Z:150M AS:i:140 XS:i:98 MD:Z:61C53T34 NM:i:2 RG:Z:0
Thank you!
Best regards, Mariya. |
_______________________________________________ Samtools-help mailing list Samtools-help@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/samtools-help
