I want to get the number of each covering each position in the genome using samtools. I used samtools depth, mpileup and view and they all give different results. I tried with different versions of samtools (0.1.1 and 1.6) and using different options, but the results don't match. Is there something that I am missing?
samtools depth file.bam -aa -q 0 -Q 0 | head
chr 1 1083
chr 2 1349
chr 3 1485
chr 4 1669
chr 5 1789
chr 6 1834
chr 7 2137
chr 8 2318
chr 9 2402
chr 10 2762
samtools mpileup file.com -aa -q 0 -Q 0 -f genome.fas | head | cut -f1,2,4
chr 1 T 1081
chr 2 T 1347
chr 3 A 1483
chr 4 A 1667
chr 5 A 1787
chr 6 A 1832
chr 7 C 2135
chr 8 A 2316
chr 9 G 2410
chr 10 C 2765
samtools view file.bam chr:0-1 -c
2947
•
written
by
george •
1 answer
bam2depth.c skips/ignores some reads:https://github.com/samtools/samtools/blob/3a7d0f02e8693620014cb59aea17028c5ab427e2/bam2depth.c#L65
if ( b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP) ) continue;
if ( (int)b->core.qual < aux->min_mapQ ) continue;
if ( aux->min_len && bam_cigar2qlen(b->core.n_cigar, bam_get_cigar(b)) < aux->min_len ) continue;
not the same filters, you will find things related to PROPER_PAIR : https://github.com/samtools/samtools/blob/master/bam_plcmd.c#L282
else if ((ma->conf->flag&MPLP_NO_ORPHAN) && (b->core.flag&BAM_FPAIRED) && !(b->core.flag&BAM_FPROPER_PAIR)) skip = 1;
see all under: https://github.com/samtools/samtools/blob/master/bam_plcmd.c#L229
Log in to answer this question.
