Entering edit mode
9.4 years ago
wangyang703092
▴
120
Hi,
I have a set of reads mapped against a reference genome and need take consensus sequence.So I used the Samtools-Bcftools-Vcfutils.pl pipeline ,I used the command:
samtools mpileup -uf <reference.fasta> <file.sorted.bam> | \
bcftools view -cg - | \
vcfutils.pl vcf2fq | \
perl fastq2fasta.pl cns.fasta
In the vcf file,there actually exists indels, but in the final consensus sequence, there is no indels.
For example, in the vcf
Supercontig_1.1 2953301 . T TCGTGGCCT 178 . INDEL;IS=2,0.014493;DP=138;VDB=2.191453e-01;AF1=1;AC1=2;DP4=0,0,75,55;MQ=20;FQ=-290 GT:PL:GQ 1/1:219,255,0:99
Someone have any suggestions?