<div dir="ltr">Hi Nicolas,<div><br></div><div>You don't need to have the tabix indexed one - this serves a different purpose for the VEP, in that it just finds any overlaps between features in that file and your input variants. Correct me if I'm wrong but I don't think this adds anything to the output that you can't get from the custom cache files.</div><div><br>If you have created cache files from the GTF then this will give you essentially the same location-based annotation and of course the additional sequence-based annotation.</div><div><br></div><div>Will</div></div><div class="gmail_extra"><br><div class="gmail_quote">On 9 October 2014 10:03, Nicolas Thierry-Mieg <span dir="ltr"><<a href="mailto:Nicolas.Thierry-Mieg@imag.fr" target="_blank">Nicolas.Thierry-Mieg@imag.fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Will,<br>
<br>
thanks for your reply.<br>
<br>
So, we need two different GTF files: one with exons and CDS's in transcript order for building the VEP cache, and one in genomic position order (required by tabix and therefore by <a href="http://variant_effect_predictor.pl" target="_blank">variant_effect_predictor.pl</a>, since the latter requires a tabix-indexed bgzipped gff).<br>
<br>
I can confirm that the following works on my small example: transcript.gff is my original file (sorted by increasing genomic positions), and transcript_reorder.gff is similar but sorted in transcript order.<br>
<br>
perl <a href="http://gtf2vep.pl" target="_blank">gtf2vep.pl</a> -i transcript_reorder.gff -f hs.chrom_3.fasta --dir cache/ -d 1 -s homo_sapiens<span class=""><br>
<br>
bgzip -c transcript.gff > transcript.gff.gz<br>
tabix -p gff transcript.gff.gz<br></span>
perl <a href="http://variant_effect_predictor.pl" target="_blank">variant_effect_predictor.pl</a> --custom transcript.gff.gz,MyGene,gff,<u></u>overlap,0 --offline -i snp.vcf --vcf --dir_cache cache/ --cache_version 1 --force --no_stats<br>
<br>
We'll now script sortGFF4VepCache.pl and test the above process on a genome scale. I'll report back to the list with our findings.<br>
<br>
You may want to update the following page, which incorrectly asserts that "The GTF file have its features sorted in chromosomal order" to build a cache:<br>
<a href="http://www.ensembl.org/info/docs/tools/vep/script/vep_cache.html#gtf" target="_blank">http://www.ensembl.org/info/<u></u>docs/tools/vep/script/vep_<u></u>cache.html#gtf</a><br>
<br>
<br>
Thanks again.<br>
<br>
Regards,<br>
Nicolas<span class=""><br>
<br>
<br>
On 10/08/2014 03:04 PM, Will McLaren wrote:<br>
</span><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">
Hi Nicolas,<br>
<br>
Apologies for the delay in replying to this.<br>
<br></span>
The <a href="http://gtf2vep.pl" target="_blank">gtf2vep.pl</a> <<a href="http://gtf2vep.pl" target="_blank">http://gtf2vep.pl</a>> parser at the moment is somewhat<span class=""><br>
limited, and it has very strict expectations on the format of the data<br>
coming in.<br>
<br>
One of these is that the exon and CDS lines appear in _transcript_ order<br>
(i.e. 5' -> 3'), even if this means the first exon listed for a<br>
transcript in a gene has a genomic position larger (more 3' on the<br>
genome) than the last.<br>
<br>
If I change the order of your input to the following (the start_codon<br>
and stop_codon lines are actually ignored by the parser):<br>
<br></span>
3protein_<u></u>codingexon113077591113077946.-<u></u>.gene_id toto; transcript_id<br>
toto_a; exon_number 1;<br>
3protein_<u></u>codingCDS113077591113077746.-<u></u>0gene_id toto; transcript_id<span class=""><br>
toto_a; product_id toto_a; exon_number 1;<br></span>
3protein_<u></u>codingexon113063155113063561.-<u></u>.gene_id toto; transcript_id<br>
toto_a; exon_number 2;<br>
3protein_<u></u>codingCDS113063355113063561.-<u></u>0gene_id toto; transcript_id<span class=""><br>
toto_a; product_id toto_a; exon_number 2;<br>
<br>
it works OK, and I get the expected stop_gained consequence for your<br>
input variant:<br>
<br>
#Uploaded_variation Location Allele Gene Feature<br>
Feature_type Consequence cDNA_position CDS_position<br>
Protein_position Amino_acids Codons Existing_variation Extra<br>
3_113063450_G/A 3:113063450 A toto toto_a Transcript<br>
stop_gained 468 268 90 R/* Cga/Tga - STRAND=-1<br>
<br>
The parser is something we intend to revisit and improve sometime soon,<br>
so thanks for your patience in the meantime!<br>
<br>
Regards<br>
<br>
Will McLaren<br>
Ensembl Variation<br>
<br>
<br>
On 1 October 2014 12:08, Nicolas Thierry-Mieg<br></span><span class="">
<<a href="mailto:Nicolas.Thierry-Mieg@imag.fr" target="_blank">Nicolas.Thierry-Mieg@imag.fr</a> <mailto:<a href="mailto:Nicolas.Thierry-Mieg@imag.fr" target="_blank">Nicolas.Thierry-Mieg@<u></u>imag.fr</a>>> wrote:<br>
<br>
Hello developers,<br>
<br>
we would like to use VEP with custom transcript definitions. To this<br>
end we are building a VEP cache, following instructions found here:<br></span>
<a href="http://www.ensembl.org/info/__docs/tools/vep/script/vep___cache.html" target="_blank">http://www.ensembl.org/info/__<u></u>docs/tools/vep/script/vep___<u></u>cache.html</a><div><div class="h5"><br>
<<a href="http://www.ensembl.org/info/docs/tools/vep/script/vep_cache.html" target="_blank">http://www.ensembl.org/info/<u></u>docs/tools/vep/script/vep_<u></u>cache.html</a>><br>
<br>
For some transcripts things seem to work fine, but for others the<br>
SNP effects reported by VEP are incorrect.<br>
<br>
We have constructed a small example that exhibits the problem. I am<br>
attaching the GFF but also copying it inline, in case attachments<br>
don't work on this ML:<br>
3 protein_coding exon 113063155 113063561 .<br>
- . gene_id toto; transcript_id toto_a; exon_number 2;<br>
3 protein_coding stop_codon 113063352 113063354<br>
. - . gene_id toto; transcript_id toto_a;<br>
product_id toto_a;<br>
3 protein_coding CDS 113063355 113063561 .<br>
- 0 gene_id toto; transcript_id toto_a; product_id<br>
toto_a; exon_number 2;<br>
3 protein_coding exon 113077591 113077946 .<br>
- . gene_id toto; transcript_id toto_a; exon_number 1;<br>
3 protein_coding CDS 113077591 113077746 .<br>
- 0 gene_id toto; transcript_id toto_a; product_id<br>
toto_a; exon_number 1;<br>
3 protein_coding start_codon 113077744 113077746<br>
. - . gene_id toto; transcript_id toto_a;<br>
product_id toto_a;<br>
<br>
This is a (fabricated) two-exon transcript on the reverse strand of<br>
chromosome 3, using hg19 coordinates.<br>
I have checked and re-checked, the two CDS features indeed represent<br>
a CDS (no STOPs in-frame).<br>
For reference, the CDS sequences from exons 1 and 2 (hg19, chrom3,<br>
reverse strand) are respectively:<br></div></div>
atgaacagatcttcttatctgcaggatttt<u></u>__<u></u>gatgtagattcccaaatccgtgctgagatg<u></u>__<u></u>cacagaaaaacagctttcaaaattcaacaa<u></u>__<u></u>gtggaaaaggaattagcttgggaaaaagag<u></u>__<u></u>aaacatgaactcggcctaatgaagctaaag<u></u>__aatcgg<br>
agatttcgagatccactggaaagtgatact<u></u>__<u></u>attgtggttcatgccatactgagtgaccac<u></u>__<u></u>aagatatcctcttacaggctggtgcagccc<u></u>__<u></u>tctaagtactccaaattcaaacgagctagt<u></u>__<u></u>caatcagagagaaaaccaagcaaattggac<u></u>__<u></u>aggtttgaaaaagagggacctggaagaaag<u></u>__gacagccagagagatgcaggtagccta<span class=""><br>
<br>
<br>
Using this GFF, I build the cache with:<br>
bgzip -c transcript.gff > transcript.gff.gz<br>
tabix -p gff transcript.gff.gz<br></span>
perl <a href="http://gtf2vep.pl" target="_blank">gtf2vep.pl</a> <<a href="http://gtf2vep.pl" target="_blank">http://gtf2vep.pl</a>> -i transcript.gff -f<span class=""><br>
hs.chrom_3.fasta --dir cache/ -d 1 -s homo_sapiens<br>
<br>
<br>
I then have a single SNV in a VCF file:<br>
#CHROM POS ID REF ALT QUAL FILTER INFO<br>
3 113063450 . G A . PASS<br>
<br>
This SNV falls within the CDS, it is a nonsense mutation<br>
(stop-gained). However, when I run VEP it sees it as a<br>
"3_prime_UTR_variant".<br>
<br>
<br>
Specifically, I run VEP with:<br>
perl <a href="http://variant_effect_predictor.pl" target="_blank">variant_effect_predictor.pl</a><br></span>
<<a href="http://variant_effect_predictor.pl" target="_blank">http://variant_effect_<u></u>predictor.pl</a>> --custom<br>
transcript.gff.gz,MyGene,gff,_<u></u>_overlap,0 --offline -i snp.vcf --vcf<span class=""><br>
--dir_cache cache/ --cache_version 1<br>
<br>
I obtain:<br>
##VEP=v76 cache=cache//homo_sapiens/1 db=.<br></span>
##INFO=<ID=CSQ,Number=.,Type=_<u></u>_String,Description="__<u></u>Consequence<span class=""><br>
type as predicted by VEP. Format:<br></span>
Allele|Gene|Feature|Feature___<u></u>type|Consequence|cDNA___<u></u>position|CDS_position|Protein_<u></u>__position|Amino_acids|Codons|<u></u>__Existing_variation|DISTANCE|<u></u>__STRAND"><br>
##INFO=<ID=MyGene,Number=.,__<u></u>Type=String,Description="__<u></u>transcript.gff.gz<span class=""><br>
(overlap)"><br>
#CHROM POS ID REF ALT QUAL FILTER INFO<br>
3 113063450 . G A . PASS<br></span>
CSQ=A|toto|toto_a|Transcript|_<u></u>_3_prime_UTR_variant|468||||||<u></u>|__-1;MyGene=CDS_3:113063355-_<u></u>_113063561,exon_3:113063155-__<u></u>113063561<span class=""><br>
<br>
<br>
<br>
This is with the latest release of VEP (76).<br>
<br>
<br>
Please let me know if more info is needed to debug this. It would be<br>
great if VEP could be used with custom annotations!<br>
<br>
<br>
Regards,<br>
Nicolas Thierry-Mieg<br>
<br>
<br>
<br>
--<br></span>
------------------------------<u></u>__----------------------------<u></u>-<span class=""><br>
Nicolas Thierry-Mieg<br>
Laboratoire TIMC-IMAG/BCM, CNRS UMR 5525<br>
Pavillon Taillefer, Faculte de Medecine<br>
38700 La Tronche, France<br></span>
tel: <a href="tel:%28%2B33%29456.520.067" value="+33456520067" target="_blank">(+33)456.520.067</a> <tel:%28%2B33%29456.520.067>, fax:<br>
<a href="tel:%28%2B33%29456.520.055" value="+33456520055" target="_blank">(+33)456.520.055</a> <tel:%28%2B33%29456.520.055><br>
------------------------------<u></u>__----------------------------<u></u>--<br>
<br>
<br>
______________________________<u></u>_________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a> <mailto:<a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a>><span class=""><br>
Posting guidelines and subscribe/unsubscribe info:<br>
<a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/<u></u>mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
<br>
<br>
<br>
<br>
______________________________<u></u>_________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a><br>
Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/<u></u>mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
<br>
</span></blockquote><div class="HOEnZb"><div class="h5">
<br>
-- <br>
------------------------------<u></u>-----------------------------<br>
Nicolas Thierry-Mieg<br>
Laboratoire TIMC-IMAG/BCM, CNRS UMR 5525<br>
Pavillon Taillefer, Faculte de Medecine<br>
38700 La Tronche, France<br>
tel: <a href="tel:%28%2B33%29456.520.067" value="+33456520067" target="_blank">(+33)456.520.067</a>, fax: <a href="tel:%28%2B33%29456.520.055" value="+33456520055" target="_blank">(+33)456.520.055</a><br>
------------------------------<u></u>------------------------------<br>
<br>
______________________________<u></u>_________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a><br>
Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/<u></u>mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
</div></div></blockquote></div><br></div>