I have a set of illumina paired end reads (rone.fq and rtwo.fq), and I have mapped those reads to a reference (ref.fa).
bwa aln ref.fa rone.fq > rtwo.sai bwa aln ref.fa rtwo.fq > rtwo.sai bwa sampe ref.fa rone.sai rtwo.sai rone.fq rtwo.fq > aln.sam #convert sam to bam samtools view -bS aln.sam > aln.bam #sort bam file samtools view aln.bam aln_sorted
Now I am trying to get SNP and INDEL information using the following script
samtools mpileup -ugf ref.fa aln_sorted.bam >aln.bcf #convert to vcf for viewing bcftools view aln.bcf > aln.vcf
I get the following VCF (this is just one line of the VCF I actually get)
##fileformat=VCFv4.1 ##samtoolsVersion=0.1.17 (r973:277) #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT aln_sorted.bam gi|110227054|gb|AE004091.2| 32 . A C,X 0 . DP=14;I16=9,1,1,0,352,12854,23,529,600,36000,29,841,159,3395,21,441 PL 0,10,189,30,192,202
I am assuming that samtools thinks this is a SNP, because the REF=A and ALT=C,X. My question is: How can I know the coverage depth of this SNP given this VCF? If this VCF does not have information enough to show the SNP coverage depth then how can I creat a VCF that shows that information? Thanks