7
$\begingroup$

I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.

I am heavily filtering reads that are not paired and of a certain orientation. Applying the filtering, I get a few thousands of reads:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

Gives me 89502 reads.

If I then pipe this into samtools mpileup, I get no results:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 0.

I tried different combinations of filtering, and when I do both -f 16 and -f 32 returns empty, but if I do either of those, then it works:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 1056.

Any ideas why? My thinking was that it would work with --excl-flags 0.

EDIT: substituting mpileup for depth does work, and prints out each position and the depth as expected.

EDIT2: adding -q 0 to mpileup gives the same empty result.

Thanks in advance

Improve this question
$\endgroup$
5
  • $\begingroup$ Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference? $\endgroup$
    – Kamil S Jaron
    May 24, 2017 at 16:52
  • $\begingroup$ Are you aware that -f 16 combined with -f 32 (btw, you can just -f 48 and skip the piping) are filtering for discordant reads? What happens if you add -q 0 to mpileup? $\endgroup$
    – Devon Ryan
    May 24, 2017 at 17:26
  • $\begingroup$ Since depth and mpileup share a fair bit of code, I have to wonder if you simply have no variants there. What happens if you specify -a and filter out the positions with no coverage? $\endgroup$
    – Devon Ryan
    May 24, 2017 at 23:02
  • $\begingroup$ BTW, I suspect this will be a case where you'll need to post the BAM file somewhere. There are a few (currently commented out) printf()s in the samtools code that could probably be used to track down what's actually happening here. $\endgroup$
    – Devon Ryan
    May 25, 2017 at 19:01
  • $\begingroup$ Devon already mentioned it in passing but it’s too important to get lost: Samtools flags are bit masks that you need to combine! Don’t call samtools repeatedly; call it once, with the combined bit flag. You can use the BROAD online tool to compute the combined flag if you don’t want to do it manually. $\endgroup$
    – Konrad Rudolph
    May 26, 2017 at 10:29

1 Answer 1

Reset to default
2
$\begingroup$

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:

...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):

...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Bingo! Adding -A, --count-orphans do not discard anomalous read pairs did the trick! $\endgroup$
    – 719016
    May 25, 2017 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.