7
$\begingroup$

I am running samtools mpileup (v1.4) on a bam file with very choppy coverage (ChIP-seq style data). I want to get a first-pass list of positions with SNVs and their frequency as reported by the read counts, but no matter what I do, I keep getting all SNVs filtered out as not passing QC.

What's the magic parameter set for an initial list of SNVs and frequencies?

EDIT: this is a question I posted on "the other" website, but didn't get a reply there.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Would you be able to try another caller just to check? Something like varscan? $\endgroup$
    – nuin
    May 26, 2017 at 14:54
  • $\begingroup$ @nuin I may give it a try, does varscan allow for switching off the filtering? $\endgroup$
    – 719016
    May 26, 2017 at 15:32
  • $\begingroup$ Could you add some examples of commands you tried? $\endgroup$
    – bli
    May 26, 2017 at 16:25
  • 1
    $\begingroup$ @719016 Try this command, and change options like --min-coverage _, _--min-reads2 and --min-avg-qual. This will also help you check if your problem is samtools related or BAM file related $\endgroup$
    – nuin
    May 26, 2017 at 17:40
  • $\begingroup$ am just a bit worried that if the depth is not high and also ChIP-Seq data has biases then it is not high confidence SNPs. The whole point of HC SNP are the read depth as well which is not so much in traditional ChIP-Seq $\endgroup$
    – ivivek_ngs
    May 31, 2017 at 14:55

2 Answers 2

Reset to default
7
$\begingroup$

I used this in the past for ChIP-seq data and it generated SNVs:

samtools mpileup \
--uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \
--fasta-ref ref.fa file.bam \
| bcftools call --consensus-caller \
> out.vcf

This was samtools 1.3 in case that makes a difference.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I got it to work with these parameters. Thanks! $\endgroup$
    – 719016
    May 27, 2017 at 9:08
5
$\begingroup$

Another approach is htsbox. You can get a candidate list with:

htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf

Here, -q sets min mapping quality, -Q sets min base quality, -v outputs variants only -c outputs VCF, -C gives you base counts on both strands and finally -s5 requires at least 5 high-quality bases to call out an allele. It is useful when your data are failing the assumptions made by typical variant callers.

Why not samtools+bcftools or varscan? Transparency and speed. This command line simply counts based on the parameters you use. It does not apply additional operations. And because of this, it is over an order of magnitude faster than samtools mpileup or varscan. It is worth noting that samtools uses BAQ by default, which reduces FPs occasionally. However, BAQ is not quite necessary for longer Illumina reads and it hurts sensitivity at the same time.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I am definitely going to try htsbox, as you say it works by simply counting based on the parameters. $\endgroup$
    – 719016
    May 30, 2017 at 13:24

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.