Tweet
Hi all,
I'm trying to get the most reliable consensus sequence from raw VCF files generated by samtools and bcftools. I tried "vcfutils.pl" coming with samtools, but I always got warning messages when running the perl script and I'm afraid the output may be not what we want since we are play with haploid data.
1. My strategy is finding the route with highest quality score. For example, there are three possible routes from xxxxx43T to xxxxx47A in the following region:
a. 43 T -> 44 T -> 45 C -> 46 C -> 47 A , the sequence is TTCCA
b. 43 T -> 44-46 TCC deletion -> 47 A, the sequence is TA
c. 43 T -> 44 T -> 45 C -> 46 C deletion -> 47 A, the sequence is TTCA
Obviously, route "a" is the most reliable one according to the overall quality score. I worked on the raw variance calls because even a variance was not that reliable, I would still take it if its quality was better than corresponding reference nucleotide(s).
Question 1. Does this strategy sound reasonable?
2. Sometimes the variance site comes with genotype inferred:
but somethimes not:
Question 2. What's the reason for the difference?
The genotype is somewhat important in my analysis because our data come from haploid sample. Only homogeneous genotypes are expected to get. If a reliable heterogeneous genotype is called, it indicates possible contamination happened, which should be filtered for quality control purpose.
Question 3. In this example, should I consider the quality score of the genotype to calculate the final score of this nucleotide?
The quality score of the variance is 41.1, and the score of the genotype to be 1/1 is 33 when the variance is TRUE. So the final score of this nucleotide to be 'T' should be some number less than 33, right?
3. And, if there is an insertion and it's taken as the most reliable annotation of the sequence (suppose it's true):
then
Question 4. What's the score for each individual nucleotide in the inserted sequence?
I think the overall score should be the Phred transformation of the probability that both the variance (score 30) and the genotype (score 59) are true. But I have no idea how to estimate the quality of each nucleotide when I need to generate FASTQ output of the consensus sequences. Maybe I should assume all the nucleotides are of the same quality in the alternative sequence "AGGGGGG"?
Thanks a lot
I'm trying to get the most reliable consensus sequence from raw VCF files generated by samtools and bcftools. I tried "vcfutils.pl" coming with samtools, but I always got warning messages when running the perl script and I'm afraid the output may be not what we want since we are play with haploid data.
1. My strategy is finding the route with highest quality score. For example, there are three possible routes from xxxxx43T to xxxxx47A in the following region:
Code:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ../sam2bam/GIAMTH101_sorted.bam chr1 147119243 . T . 241 . DP=97;AF1=0;CI95=1.5,0;DP4=42,47,0,0;MQ=20;FQ=-282 PL 0 chr1 147119244 . T . 242 . DP=96;AF1=0;CI95=1.5,0;DP4=42,47,0,0;MQ=20;FQ=-282 PL 0 chr1 147119244 . TCC . 108 . INDEL;DP=96;AF1=0;CI95=1.5,0;DP4=13,16,1,0;MQ=20;FQ=-105;PV4=0.47,1,1,1 PL 0 chr1 147119245 . C . 242 . DP=97;AF1=0;CI95=1.5,0;DP4=41,47,0,0;MQ=20;FQ=-282 PL 0 chr1 147119246 . C . 242 . DP=97;AF1=0;CI95=1.5,0;DP4=41,47,0,0;MQ=20;FQ=-282 PL 0 chr1 147119246 . C . 108 . INDEL;DP=97;AF1=0;CI95=1.5,0;DP4=13,16,1,0;MQ=20;FQ=-105;PV4=0.47,4.7e-07,1,1 PL 0 chr1 147119247 . A . 241 . DP=97;AF1=0;CI95=1.5,0;DP4=43,50,0,0;MQ=20;FQ=-282 PL 0
b. 43 T -> 44-46 TCC deletion -> 47 A, the sequence is TA
c. 43 T -> 44 T -> 45 C -> 46 C deletion -> 47 A, the sequence is TTCA
Obviously, route "a" is the most reliable one according to the overall quality score. I worked on the raw variance calls because even a variance was not that reliable, I would still take it if its quality was better than corresponding reference nucleotide(s).
Question 1. Does this strategy sound reasonable?
2. Sometimes the variance site comes with genotype inferred:
Code:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ../sam2bam/GIAMTH101_sorted.bam chr1 41098739 . A T 41.1 . DP=6;AF1=1;CI95=0.5,1;DP4=0,0,6,0;MQ=20;FQ=-45 GT:PL:GQ 1/1:74,18,0:33
Code:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ../sam2bam/GIAMTH101_sorted.bam chr1 41092204 . G . 239 . DP=157;AF1=0;CI95=1.5,0;DP4=82,70,0,0;MQ=20;FQ=-282 PL 0 chr1 41092204 . G . 38.5 . INDEL;DP=157;AF1=0;CI95=1.5,0;DP4=2,4,0,1;MQ=20;FQ=-35.5;PV4=1,1,1,0.28 PL 0
The genotype is somewhat important in my analysis because our data come from haploid sample. Only homogeneous genotypes are expected to get. If a reliable heterogeneous genotype is called, it indicates possible contamination happened, which should be filtered for quality control purpose.
Question 3. In this example, should I consider the quality score of the genotype to calculate the final score of this nucleotide?
Code:
chr1 41098739 . A T 41.1 . DP=6;AF1=1;CI95=0.5,1;DP4=0,0,6,0;MQ=20;FQ=-45 GT:PL:GQ 1/1:74,18,0:33
3. And, if there is an insertion and it's taken as the most reliable annotation of the sequence (suppose it's true):
Code:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ../sam2bam/GIAMTH101_sorted.bam chr1 147381185 . AGGGGG AGGGGGG 30 . INDEL;DP=134;AF1=0.5;CI95=0.5,0.5;DP4=39,5,7,10;MQ=20;FQ=21.5;PV4=0.00032,9.3e-06,1,1 GT:PL:GQ 1/1:56,0,80:59
Question 4. What's the score for each individual nucleotide in the inserted sequence?
I think the overall score should be the Phred transformation of the probability that both the variance (score 30) and the genotype (score 59) are true. But I have no idea how to estimate the quality of each nucleotide when I need to generate FASTQ output of the consensus sequences. Maybe I should assume all the nucleotides are of the same quality in the alternative sequence "AGGGGGG"?
Thanks a lot
Comment