I have a reference transcriptome of one fungal species and have illumina sequences of a transcriptome from a closely related species. I have mapped these onto my reference transcriptome using BLAT with 50-60% reads mapping.
I have generated SAM/BAM files using samtools. I now want to extract the consensus sequence of this mapping to generate the sequence in order to do phylogentic analyses, however here I run into problems when calling SNPs.
The sequence I want to get will consist of n's at sites were no reads mapped and the read sequence at the positions were the reads mapped. I'm able to create this sequence using samtools mpileup, bcftools, and vcfutils.pl vcf2fq.
I want to retain every single variant, also at sites where only a single read mapped which I know is not the usual case for snp calling. Setting the -A flag in bcftools provides all variant position with a coverage down to two reads, but not the ones with a single read.
Here is an example from the vcf file of two positions with a read depth of one. The variant difference first (pos. 206) and a perfect match afterwards (pos. 207). With a read depth of 1, the variant at position 206 is not called:
contig00004 206 . T . 10.4 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=20;FQ=-30 PL
P 20:1
contig00004 207 . G . 36 . DP=2;AF1=0;AC1=0;DP4=0,2,0,0;MQ=20;FQ=-33 PLP 0:2
with a read depth of 2, the variant (pos. 212) is called and present in output:
contig00004 212 . T C,X 6.02 . DP=2;VDB=0.0347;AF1=1;AC1=2;DP4=0,0,0,2;MQ=20;FQ=-33 GT:PLP:GQ 1/1:36,6,0,36,6,36:2:6
Does anyone know how to tweak the samtools parameters to obtain every possible snp at read depths of a single read (maybe using the -P and -p settings in bcftools view?).
Or am I completely off here and should be using some completely different programs/tools to extract this information from SAM/BAM alignment files?
any help or suggestions will be really appreciated...
I have generated SAM/BAM files using samtools. I now want to extract the consensus sequence of this mapping to generate the sequence in order to do phylogentic analyses, however here I run into problems when calling SNPs.
The sequence I want to get will consist of n's at sites were no reads mapped and the read sequence at the positions were the reads mapped. I'm able to create this sequence using samtools mpileup, bcftools, and vcfutils.pl vcf2fq.
I want to retain every single variant, also at sites where only a single read mapped which I know is not the usual case for snp calling. Setting the -A flag in bcftools provides all variant position with a coverage down to two reads, but not the ones with a single read.
Here is an example from the vcf file of two positions with a read depth of one. The variant difference first (pos. 206) and a perfect match afterwards (pos. 207). With a read depth of 1, the variant at position 206 is not called:
contig00004 206 . T . 10.4 . DP=1;AF1=1;AC1=2;DP4=0,0,0,1;MQ=20;FQ=-30 PL
P 20:1contig00004 207 . G . 36 . DP=2;AF1=0;AC1=0;DP4=0,2,0,0;MQ=20;FQ=-33 PL
P 0:2with a read depth of 2, the variant (pos. 212) is called and present in output:
contig00004 212 . T C,X 6.02 . DP=2;VDB=0.0347;AF1=1;AC1=2;DP4=0,0,0,2;MQ=20;FQ=-33 GT:PL
P:GQ 1/1:36,6,0,36,6,36:2:6Does anyone know how to tweak the samtools parameters to obtain every possible snp at read depths of a single read (maybe using the -P and -p settings in bcftools view?).
Or am I completely off here and should be using some completely different programs/tools to extract this information from SAM/BAM alignment files?
any help or suggestions will be really appreciated...