samtools mpileup keeps making an empty file
2
0
Entering edit mode
5.8 years ago
James Reeve ▴ 130

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?

samtools pileup • 4.8k views
ADD COMMENT
2
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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...

ADD REPLY
0
Entering edit mode
5.8 years ago
gab ▴ 20

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

ADD COMMENT
1
Entering edit mode

mpileup should produce an output without reference.

ADD REPLY
0
Entering edit mode

I tired that, to no avail.

ADD REPLY
0
Entering edit mode
5.8 years ago

show me two alignments from a random BAM please.

ADD COMMENT
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 2087 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6