Hi everyone, I was wondering if you might be able to help me with specifying a minimum read depth when calling consensus sequences?
I'm running samtools-0.1.18 on Ubuntu 10.04. I have produced VCF files for all sites for my bacterial genome sequences (from BWA BAM files), which seem ok (based on using the following commandline: /samtools-0.1.18/samtools mpileup -E -q30 -Q30 -uf REF.fa isolate1.merged.bam | /samtools-0.1.18/bcftools/bcftools view -cg - > isolate1.all_sites.vcf), and am using the following command line in vcftools with the aim of specifying a minimum read depth of 3 when calling the genome consensus sequences:
/samtools-0.1.18/bcfutils/vcfutils.pl vcf2fq -d3 AB10.all_sites_E_q30_Q30.vcf > AB10.all_sites_E_q30_Q30_RD3.fq (which I convert to FASTA using a python script).
Whether or not I specify the option -d3 (I believe 3 is the default anyway) in the commandline, in all cases where the VCF file states that DP=1 or DP=2, rather than not calling the base, it calls a lowercase letter at these sites.
Over the course of my genomes, this amounts to approximately 48 000 bases in each genome which have been called lowercase, which includes approximately 650 variant sites. If these variant sites are of poor quality and include potentially false positive SNPs, this may be adversely affecting my downstream phylogenetic analysis.
My questions are:
a) Does anyone know if this is normal for it to insert a lowercase base rather than calling 'n' at these sites (which is what I was incorrectly assuming it would do!)
b) What are your opinions on whether I should just replace all of these sites with an 'n' to be on the conservative side?
Any help is much appreciated, SeqAnswers has been so great at helping in the past!
Best Wishes,
Laura
Comment