Hello!
I have a vcf file that details the alleles of a few hundred individuals' exomes.
I would like to merge this with a human reference genome (exome) in order to see how my individuals' alleles differ from the reference, rather than differ amongst themselves.
Is this possible? I am not sure where to start... I have vcftools ready for analysis.
assuming your file is tab delimited and you have a reference genome indexed with samtools faidx with the correct chromosome names.
grep -v CHROM input.txt | while IFS='' read L; do echo -ne "${L}\t" && echo "${L}" |awk '{printf("%s:%s-%s\n",$1,$2,$2);}' | xargs samtools faidx /path/to/human_hg19.fasta | grep -v '^>' | tr -d '\n'; echo ; done
UPDATE: from a real VCF and using bioalcidaejdk:
java -jar dist/bioalcidaejdk.jar -e 'stream().forEach(V->{println(V.getContig()+"\t"+V.getStart()+"\t"+V.getReference().getDisplayString()+"\t"+V.getAlleles().size()+"\t"+V.getGenotypes().stream().flatMap(G->G.getAlleles().stream()).filter(A->A.isCalled()).count()+"\t"+V.getGenotypes().stream().flatMap(G->G.getAlleles().stream()).filter(A->A.isCalled()).map(A->A.getDisplayString()).collect(Collectors.groupingBy(Function.identity(), Collectors.counting())));});' input.vcf
rotavirus 51 A 2 8 {A=6, G=2}
rotavirus 91 A 2 8 {A=7, T=1}
rotavirus 130 T 2 8 {C=1, T=7}
rotavirus 232 T 2 8 {A=1, T=7}
rotavirus 267 C 2 8 {C=7, G=1}
rotavirus 424 A 2 8 {A=7, G=1}
rotavirus 520 T 2 8 {A=1, T=7}
rotavirus 536 A 2 8 {A=6, T=2}
rotavirus 562 A 2 8 {A=7, G=1}
rotavirus 583 G 2 8 {C=1, G=7}
rotavirus 661 T 2 8 {A=1, T=7}
rotavirus 693 T 2 8 {T=6, G=2}
rotavirus 738 T 2 8 {A=1, T=7}
rotavirus 799 A 2 8 {A=6, C=2}
rotavirus 812 G 2 8 {T=2, G=6}
rotavirus 833 G 2 8 {A=2, G=6}
rotavirus 916 A 2 8 {A=6, T=2}
rotavirus 946 C 2 8 {A=1, C=7}
rotavirus 961 T 2 8 {A=1, T=7}
rotavirus 1044 A 2 8 {A=6, T=2}
rotavirus 1045 C 2 8 {C=6, G=2}
rotavirus 1054 C 2 8 {C=6, G=2}
rotavirus 1064 G 2 8 {A=2, G=6}
Thank you very much for the effort, Pierre!
It turns out I was just confused. The source of my vcf data claimed that the first column in the --freq or --counts output was NOT the reference column, but it actually was. So there was no problem to begin with...
Log in to answer this question.
the REF column should answer the question. Or give us an example of input/output.
I don't think there is a reference... having a reference is optional I think?
(There is mention of a reference 'Homo_sapiens_assembly19' in the header lines though)
Example:
so you're wrong; this is NOT a vcf file.
The example I pasted is from the output of the '--counts' option in vcftools, acting on the vcf file.
The actual vcf file is huge... what part would you like to see?
So you have a vcf, which does contain a REF field?
In addition, variant calling is usually performed by comparing to the reference genome. Sounds like you already got what you want.
Yes it appears I already had what I wanted... the source of my vcf data claimed otherwise, and I have little experience with these things, so was easily confused!