SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem generating consensus files EmmaB Bioinformatics 2 11-09-2012 02:30 PM
New consensus from VCF dbrami Bioinformatics 9 09-20-2012 11:35 PM
Consensus from a VCF file szilva Bioinformatics 1 09-20-2012 10:01 PM
Generating vcf files with samtools Brittd Bioinformatics 1 06-06-2012 12:41 AM
Differences in Consensus Sequences bhan Genomic Resequencing 0 08-01-2011 08:36 AM

Reply
 
Thread Tools
Old 07-29-2011, 02:58 PM   #1
sulicon
Member
 
Location: Los Angeles

Join Date: Aug 2010
Posts: 41
Default Questions about generating consensus sequences from the VCF file

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:
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
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:
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
but somethimes not:
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
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?
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
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):
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
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
sulicon is offline   Reply With Quote
Old 07-24-2012, 07:47 AM   #2
nupurgupta
Member
 
Location: New Jersey

Join Date: Aug 2010
Posts: 29
Default

Hi,
Did you finally figure this out, and write a script for it?
Thanks
nupurgupta is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:57 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO