vcf to consensus call reference instead of N

[spike1985]

Hi all

I have a few genomes sequenced using illumina. I have used samtools and vcfutils to make a consensus for each. All pretty standard stuff. However using vctutils to make the consensus gives me far to many N nucleotides to be happy with.

Is there a way that i can create this consensus but where it currently has provided an N nucleotide it actually inserts the reference nucleotide in its place?

Your help would be much appreciated. Any ideas?

[vivek_]

Have you tried samtools mpileup?

[spike1985]

Sorry thats what ive used in the process. I pretty much used this exact command:

samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

Anything i can change there to get it to call the ref seq rather than N's?

[vivek_]

Normally you get those Ns when you don't specify the reference or specify a different version of the reference file in the command line.

One thing I would check is if the mpileup is producing the reference bases at each position and not Ns without piping the first output to bcftools.

[spike1985]

Thanks for the response i'll test the mpileup tonight.

I'm pretty sure i have the reference sorted. I've used snpEff to analyse the SNPs in the genomes and for each SNP it has had the correct reference sequence.

Could be be a problem of coverage? as in the sequencing reads might not cover this area?

[vivek_]

Could be coverage, but if that's the case then it shouldn't be making a variant call at that site right with a certain degree of confidence right?