I'm trying to convert VCF to BED, loosing as little information as possible, but the closer I look the more problems arise.
A simple heterozygote SNP is not a problem (Skipped some information for readability):
VCF
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT IND_1 Ind_2
1 209949474 rs41303065 G A 4448.64 PASS AC=1;... GT:... 0/1:... 0/0:...
is converted to
BED
#Chrom Start Stop Ref Alt ID FORMAT Ind_1 Ind_2
1 209949474 209949474 G A rs41303065 GT:... 0/1:... 0/0:...
A multiallelic call should be split into the same number of rows as there are alternatives?
At least this is the way that the well used convert2annovar.pl does it. We have:
VCF
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT IND_1 Ind_2
1 209969979 rs5780538 A AAC,AACAC 1556.66 PASS AC=2,0;... GT:... 0/1:... 1/2:...
BED
#Chrom Start Stop Ref Alt ID FORMAT Ind_1 Ind_2
1 209969980 209969982 - AC rs41303065 GT:... 0/1:... 1/2:...
1 209969980 209969984 - ACAC rs41303065 GT:... 0/1:... 1/2:...
My questions are:
Is the dbSNP-ID relevant for all alternatives in one position?
How do we preserve the genotype call here?
Individual 1 is genotyped A/AC and individual 2 AAC/ACAC, how do we show this when the call is split?
If we want to annotate the frequency of this variant by using some frequency database, will the different alternatives get different frequencies?
I know that most people skip these positions since it is most likely that we have false positives when multiallelic calls, I would really like to not skip things without closer examination.
We learn from GATK that the larger number of individuals called simultaneously the better quality our call set will have, also the more multiallelic calls we will see.
So when annotating our callset for future analysis it is importans that these calls can be handeled.
I really hope that some of you can contribute with some ideas and solutions.
Best regards,
Måns
A simple heterozygote SNP is not a problem (Skipped some information for readability):
VCF
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT IND_1 Ind_2
1 209949474 rs41303065 G A 4448.64 PASS AC=1;... GT:... 0/1:... 0/0:...
is converted to
BED
#Chrom Start Stop Ref Alt ID FORMAT Ind_1 Ind_2
1 209949474 209949474 G A rs41303065 GT:... 0/1:... 0/0:...
A multiallelic call should be split into the same number of rows as there are alternatives?
At least this is the way that the well used convert2annovar.pl does it. We have:
VCF
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT IND_1 Ind_2
1 209969979 rs5780538 A AAC,AACAC 1556.66 PASS AC=2,0;... GT:... 0/1:... 1/2:...
BED
#Chrom Start Stop Ref Alt ID FORMAT Ind_1 Ind_2
1 209969980 209969982 - AC rs41303065 GT:... 0/1:... 1/2:...
1 209969980 209969984 - ACAC rs41303065 GT:... 0/1:... 1/2:...
My questions are:
Is the dbSNP-ID relevant for all alternatives in one position?
How do we preserve the genotype call here?
Individual 1 is genotyped A/AC and individual 2 AAC/ACAC, how do we show this when the call is split?
If we want to annotate the frequency of this variant by using some frequency database, will the different alternatives get different frequencies?
I know that most people skip these positions since it is most likely that we have false positives when multiallelic calls, I would really like to not skip things without closer examination.
We learn from GATK that the larger number of individuals called simultaneously the better quality our call set will have, also the more multiallelic calls we will see.
So when annotating our callset for future analysis it is importans that these calls can be handeled.
I really hope that some of you can contribute with some ideas and solutions.
Best regards,
Måns