samtools mpileup keeps making an empty file

Question

[James Reeve]

I'm trying to call SNPs from a pooled DNA sample using the PoPoolation2 pipeline. I'm stuck on a step that requires me to make a mpileup file from my BAM files. According to multiple tutorials this should work;

samtools mpileup -B input1.bam input2.bam > out.mpileup

However, when I try it makes an empty file.

I've tried adding a reference, re-sorting my BAM files and I checked for errors using Picardtool's ValidateSamFile. Even weirder, I can generate a file using GATK's Pileup command. Any ideas why samtools mpileup doesn't work for me?

[Pierre Lindenbaum]

version of samtools ? any error message ? number of reads per bam ? average mapping quality of the reads ? etc...

[James Reeve]

Samtools v1.9. No error messages. Average 470 million reads per bam. MAPQ = 48.56.

[tbg]

i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be? I asked the same question here but wondering whether you have other ideas...

Answer 1

[gab]

Try to pass -f <reference fasta> as a parameter before the bam files

[Pierre Lindenbaum]

mpileup should produce an output without reference.

[James Reeve]

I tired that, to no avail.

Answer 2

[Pierre Lindenbaum]

show me two alignments from a random BAM please.

[James Reeve]

Alignment 1:

E00434:123:HM2NVCCXY:8:1101:8237:62452  1185    HiC_chromosome_1        1       60      34S117M =       47      197     GATCGTTGCTATATTTGTTCATAGTTTGTTTATTTGTTTTTCCACTGCTTTGCCAAAATAAACGGTAAAAAGTGTTTCCAAAACTTTTTCCTGGATATAACGCACACGATTTTCCAAGCACAAAAACATGGCGTATAACGGCAGGAGAGAC AAFFFAJJJ<JFJJJJJJJFJ7JJJJFJJJJFJ-FJJJJJAFAJ-FF<JJJJJJA7FJJJJJJJJJAJJJAJJ7FJAJ<AFAJJJFFJJAFFJJJFJAJ-AFJJJJ7A---AJJFAAF<JFFJAFFFAA-<A7AJ-----AA)77))-7-7 MD:Z:108C8      PG:Z:MarkDuplicates     RG:Z:Index_12   NM:i:1  AS:i:112        XS:i:72

Alignment 2:

E00434:123:HM2NVCCXY:8:2224:23906:72860 77  *   0   0   *   TTTAGTCTTTTTTGTCTTAATAGGCTCCCCTCGATTTAATTCACTCCCCCCAATTAATGAAAACAACTCAGCAGTTATCAACCCCATCAAACGACTTGCATGAGGGAGTATTCTCGCTGGTCTCTTAATTACTTCAGATCGGGAGAGCACC AAF<A<FFJJJJJFFFFJ-FFJ<JJFFJJJ<<--<JJJJJJ<-A7<FJFJJ-A<J<F--FJJJ-FJ--7<-7<7-7FJ<<J-AJJ-7<FJA7-A7-AAFFFAFJJJ<A--7<AJF<JJ7A--<A-<-<FFA7AA-<FJAFAF)F-AAFFJ) PG:Z:MarkDuplicates RG:Z:Index_12   AS:i:0  XS:i:0

[Pierre Lindenbaum]

alignment 2 is not mapped, alignment 1 has sam flag 1185 = read is PCR or optical duplicate . Is there any non-exotic sam record in your sam file?

[James Reeve]

Here's one that is more normal

E00434:176:HV2G5CCXY:3:1109:27123:67955 97      chrI    1       60      50S79M  =       1       79      CTTCTCCTCTTCCTCTTCTTCCTCTTCCTCCTCCTCTTCTTCCTCTTCCTCTTCTTCTCCTCTTCCTCTTCCTTCTCCCTCTTCCTCCCTGTCGGCGTGTCATCAGATCTGACCAGTGTGTGTGTGTGT       AAFFFJJJJJJJJJJFJJJJJJJJJFJJJJJFJJJJJFJFJFJJJJJJJFFAJJJJFJJJJJFFJJJJJJJJJAF<FJJAJAAJFJFJJJJFJFJJ<7F<JFFJJJJJJJF<JJJFFJJ<7A7<J<FJF       SA:Z:chrXII,8525877,+,17S40M72S,0,2;    MD:Z:79 NM:i:0  AS:i:79 XS:i:27 RG:Z:IDT_64     PG:Z:MarkDuplicates

[Pierre Lindenbaum]

this read is paired but "not properly paired" , pileup will ignore this

[James Reeve]

I checked all my reads. Most reads are flagged as improperly paired. I toggled the option to not skip anomalous read pairs -A and it works. Thanks Pierre for walking me through this.