samtools mpileup with many BAM files

[kz26]

How does samtools mpileup work with hundreds of BAM files? I have a very large, high-coverage dataset that I need to call variants in. I suppose I could just run samtools mpileup ... *.bam, but I feel like that could turn out badly.


Alternatively, I could call variants in each BAM separately and generate many VCF files. Are there any tools out there that will merge VCF files and do the appropriate math with quality scores / read depths / etc.?

[arolfe]

Why do you think it won't work on hundreds of bam files? The files are sorted, so it presumably runs through all the files at once but is only looking at a little data from each file at a given time. I've done work on viral genomes with 20k-30k read depths and had no problems.

mpileup does have some limitations as to how many reads it will use at any given locus. I believe it's 8000 and is hardcoded in; that may bite you more than runtime issues if you really have that much data. If you're running up against this limit, I can send you a patch to samtools but you'll need to compile samtools yourself.

[kz26]

So I found that while samtools does work fine with a large amount of input files, it is much faster to concatenate all the BAM files first, sort the single large file, and run mpileup on that.

[arthurmelo]

I think this command is good ...
samtools mpileup -uD -f <reference>.fasta ./1_sorted.bam ./2_sorted.bam ./3_sorted.bam ./4_sorted.bam | bcftools view -bvcg - > SNPs.raw.bcf
bcftools view SNPs.raw.bcf | vcfutils.pl varFilter -D100 > SNPs.flt.vcf
.
.
.
And to SNP counting:

vcf-stats SNP.flt.vcf > stats.txt
grep 'snp_count' stats.txt > snp_count.txt
awk '{sum+=$1} END {print sum}' snp_count.txt