very interesting Pontus, thanks. (And welcome on Biostar)
• 0 views
•
link
written
by
Pierre Lindenbaum
How To Retrieve Gene Variations From Ensembl Using The Perl Api?
Dear all,
A very practical question. I need to retrieve all variations (loose sense) of the MEN1 gene from the human genome at ensembl. Of course, I want to do it programmatically. But, I'm having a hard time using the perl API. I think I didn't grasp the adaptor idea (maybe the nomenclature scheme). Has someone a hint on how to use a gene ID to recover variations?
I've tried some examples from the tutorial. But, something went wrong.
• 11,269 views
•
link
updated
by
Biostar
•
written
by
Jarretinha
•
written
by
Jarretinha
4 answers
Hi,
The numbers you are looking at in the variation table on the gene page in Ensembl are the number of different consequences that affects the transcripts of the gene by variations that fall within it. In other words, one single variation can have multiple consequences on one transcript (i.e. 'intronic' and 'essential splice site') and it can also have different consequences depending on what transcript you're looking at (e.g. 'intronic' in one transcript but 'non-synonymous coding' in another). The table summarizes these numbers.
To get just the unique variations that falls within this gene, the API script provided above will work fine. If you need the variation consequences for each variation on each transcript, you can do something like this:
my $species='human';
my $reg = 'Bio::EnsEMBL::Registry';
$reg->load_registry_from_db(-host => 'ensembldb.ensembl.org',-user => 'anonymous');
my $gene_adaptor =$reg->get_adaptor($species, 'core','Gene');
my $gene = $gene_adaptor->fetch_by_stable_id('ENSG00000133895');
my $tv_adaptor = $reg->get_adaptor($species,'variation','transcriptvariation');
my $tvs = $tv_adaptor->fetch_all_by_Transcripts($gene->get_all_Transcripts());
foreach my $tv (@{$tvs}) {
print join("\t",($tv->variation_feature->variation_name(),$tv->transcript_stable_id,join(",",@{$tv->consequence_type()}))) . "\n";
}
Hope this helps
/Pontus Larsson - Ensembl Variation
PS. Note that in the Ensembl pipeline, we do some consistency checking on the data we import. Prior to release 61, we deleted data for variations that fail these checks from the database but starting from release 61, we just flag them and keep the data. By default, these flagged variations are not displayed in the web browser so they won't be included in the variation consequence table. These flagged variations can be viewed in the location view by turning on the 'Failed variations' track. The default API behaviour is also to not return these variations but this behaviour can be changed (see the variation API tutorial for details). The MySQL query given above will return these flagged variations and you would need to do a left join to the 'failed_variation' table to properly filter them out.
• 27 views
•
link
updated
by
Ram
•
written
by
Pontus Larsson •
•
written
by
Pontus Larsson •
Wonderful! Now, everything makes sense!!! Gonna flag this as the right one for the excelent explanation. But, everyone contributed to made my day!!! YEAH! Variations in hand right now via perl API and SQL!!! Certainly will include a mention to Biostar in the paper acknowledgement!
• 0 views
•
link
written
by
Jarretinha
my $species='human';
my $reg = 'Bio::EnsEMBL::Registry';
$reg->load_registry_from_db(-host => 'ensembldb.ensembl.org',-user => 'anonymous');
my $vfa = $reg->get_adaptor($species, 'variation', 'variationfeature');
my $slice_adaptor = $reg->get_adaptor($species, 'core', 'slice');
my $gene_adaptor =$reg->get_adaptor($species, 'core','Gene');
my $gene = $gene_adaptor->fetch_by_stable_id('ENSG00000133895');
my $slice = $gene->feature_Slice();
my $counter = 0;
foreach my $vf (@{$vfa->fetch_all_by_Slice($slice)}) {
$counter++;
print $vf->seq_region_start(), '-', $vf->seq_region_end(), ' ', $vf->allele_string(), "\n";
}
print "counter: $counterndonen";
• 0 views
•
link
updated
by
Ram
•
written
by
User 6659 •
•
written
by
User 6659 •
That's good too. I've tried smth very similar. Your script returned 474 results. But, I can't understand exactly what's happening. It returned HGMD mutations, but just part of them (6348 in total). Should I recurse over transcripts?
• 0 views
•
link
written
by
Jarretinha
i don't know what information you want to be fair. Do you want just HGMD variants? Do you want the variation_features (i.e. the locations of the variations) or all the transcripts affected by the variants. Do you want SNPs and indels?
• 0 views
•
link
written
by
Andrea_Bio •
i don't know what information you want. Do you want just HGMD variants? Do you want the variation_features (i.e. the locations of the variations) or all the transcripts affected by the variants. Do you want SNPs and indels?
• 0 views
•
link
written
by
User 6659 •
In fact I want all variants of a gene. I need them to model how this gene is evolving in the human population. I don't understand Ensembl numbers. Check this numbers http://www.ensembl.org/Homo_sapiens/Gene/Variation_Gene/Table?g=ENSG00000133895;r=11:64570988-64578766
• 0 views
•
link
updated
by
Ram
•
written
by
Jarretinha
•
written
by
Jarretinha
I need all of them!!! That's the problem. Check these numbers http://www.ensembl.org/Homo_sapiens/Gene/Variation_Gene/Table?g=ENSG00000133895;r=11:64570988-64578766
• 0 views
•
link
updated
by
Ram
•
written
by
Jarretinha
•
written
by
Jarretinha
the api and db return different numbers becase the api only gets dbSNP and hgmd whereas the db gets everything (cosmic too). I don't know where the website gets the figure 6348 from at all. I get 473 vfs in dbSNP and cosmic with 13699 transcript variants with the db, and 474 and 13700 with the api. One has gone missing somewhere though :)
• 0 views
•
link
written
by
User 6659 •
You want somatic variations too?
• 0 views
•
link
written
by
User 6659 •
if you want somatic mutations as well i don't know how to get them with the api. if you just want inherited variants you can use the api code already give but recurse over transcripts (and you'll get 13700). if you do need somati variants then you can use this query in the db
• 0 views
•
link
written
by
User 6659 •
mysql> select count(*) from variation_feature vf \
inner join seq_region sr on sr.seq_region_id = vf.seq_region_id \
inner join variation v on vf.variation_id = v.variation_id \
left outer join failed_variation fv on fv.variation_id = v.variation_id \
inner join transcript_variation tv \
on tv.variation_feature_id = vf.variation_feature_id \
inner join source s on s.source_id = vf.source_id \
where sr.name = '11' and vf.seq_region_start >=64570988 and \
vf.seq_region_end <= 64578766 and fv.variation_id is null;
• 0 views
•
link
updated
by
Ram
•
written
by
User 6659 •
•
written
by
User 6659 •
I'll surely try this. The number differences are really puzzling me. Anyway, MEN1 is a much more complicated gene than I tought.
• 0 views
•
link
written
by
Jarretinha
some bad luck with numbers here actually. I thought that because i got the same number of variants as the previous respondent who filtered by gene, there was only one gene on that chromosome. I should have checked but there is > 1 gene there. The perl and sql given assumes only one gene there. I've modified the perl api code and sql looking for only transcripts in that gene and get figures of about 9000 but not 6000
• 0 views
•
link
written
by
User 6659 •
i deleted some previous comments for clarity. some bad luck with numbers here actually. I thought that because i got the same number of variants with the sql as the previous respondent who filtered by gene, there was only one gene on that chromosome. I should have checked but there is > 1 gene there. The perl given assumes only one gene there. I've modified the perl api code and sql looking for only transcripts in that gene
• 0 views
•
link
written
by
User 6659 •
i deleted some previous comments for clarity. some bad luck with numbers here actually. I thought that because i got the same number of variants with my sql without filtering by gene as the previous respondent (who filtered by gene), there was only one gene on that chromosome. I should have checked but there is > 1 gene there. The perl given assumes only one gene there. I've modified the perl api code and sql looking for only transcripts in that gene
• 0 views
•
link
written
by
User 6659 •
i deleted some previous comments for clarity. some bad luck with numbers here actually. I thought that because i got the same number of variants with my sql without filtering by gene as the previous respondent (whose query seemed to filter by gene), there was only one gene on that chromosome. I should have checked but there is > 1 gene there. The perl given assumes only one gene there. I've modified the perl api code and sql looking for only transcripts in that gene – unknown (google) 0 secs ago
• 0 views
•
link
written
by
User 6659 •
i deleted some previous comments for clarity. some bad luck with numbers here actually. I thought that because i got the same number of variants with my sql without filtering by gene as the previous respondent (whose query seemed to filter by gene), there was only one gene on that chromosome. I should have checked but there is > 1 gene there. The perl given assumes only one gene there. I've modified the perl api code and sql looking for only transcripts in that gene
• 0 views
•
link
written
by
User 6659 •
Just to add to Pontus' answer, if you are also interested in somatic mutations (which are included in the numbers in the table on the website) you can fetch these with the fetch_all_somatic_by_Transcripts() method on the TranscriptVariationAdaptor. The script above could be modified to include these as follows:
my $tvs = $tv_adaptor->fetch_all_by_Transcripts($gene->get_all_Transcripts());
push @$tvs, @{ $tv_adaptor->fetch_all_somatic_by_Transcripts($gene->get_all_Transcripts()) };
And then loop over $tvs as before.
Cheers,
Graham Ritchie - Ensembl Variation
• 0 views
•
link
updated
by
Ram
•
written
by
Graham Ritchie •
•
written
by
Graham Ritchie •
Many thanks for the tip! Be welcome you too to Biostar! Now I've got the gist of the API. Let's see what can I do . . .
• 0 views
•
link
written
by
Jarretinha
I thought it would be interesting to find the mysql query for this request.
I'm not a specialist of the SQL schema for ENSEMBL and I would be interested if someone could improve and/or criticize this SQL query.
> mysql -A -h ensembldb.ensembl.org -u anonymous -P 5306
use homo_sapiens_core_61_37f;
select distinct
XREF.display_label,
VAR.name
from
external_db as EXTDB,
xref as XREF,
gene as GENE,
homo_sapiens_variation_61_37f.variation_feature as VARFEAT,
homo_sapiens_variation_61_37f.variation as VAR
where
EXTDB.db_name="HGNC" and
XREF.external_db_id = EXTDB.external_db_id and
XREF.display_label="MEN1" and
XREF.xref_id = GENE.display_xref_id and
GENE.seq_region_id = VARFEAT.seq_region_id and
VARFEAT.seq_region_start >= GENE.seq_region_start and
VARFEAT.seq_region_end <= GENE.seq_region_end and
VARFEAT.variation_id = VAR.variation_id
order by VARFEAT.seq_region_start
;
+---------------+-------------+
| display_label | name |
+---------------+-------------+
| MEN1 | rs117705251 |
| MEN1 | rs34171410 |
| MEN1 | rs71581758 |
| MEN1 | rs1804850 |
| MEN1 | rs71581759 |
| MEN1 | rs71996491 |
| MEN1 | rs1804848 |
| MEN1 | rs1804849 |
| MEN1 | rs111895237 |
| MEN1 | rs71581749 |
(...)
• 0 views
•
link
updated
by
Ram
•
written
by
Pierre Lindenbaum
•
written
by
Pierre Lindenbaum
Absolutely nice! A SQL approach would be equally suitable! And you're right about ensembl schemas. They're a complexity nightmare. So, let's see how this query behaves...
• 0 views
•
link
written
by
Jarretinha
Hum... I think we hit the same problem with different approaches. When you enter MEN1 in the websearch, it returns 542 variants (mutations+variants in ensembl sense). But, when you browse the human gene there are 9751! Your query returned 621. I'm a bit lost right now. Time to tweak a little ...
• 0 views
•
link
written
by
Jarretinha
That query returns 621 because it includes failed variants. There are 618 without failed variants.
• 0 views
•
link
written
by
User 6659 •
mysql> select * from variation_feature vf \
inner join seq_region sr on sr.seq_region_id = vf.seq_region_id \
inner join variation v on vf.variation_id = v.variation_id left \
outer join failed_variation fv on fv.variation_id = v.variation_id \
where sr.name = '11' and vf.seq_region_start >=64570988 and \
vf.seq_region_end <= 64 578766 and fv.variation_id is null;
• 0 views
•
link
updated
by
Ram
•
written
by
User 6659 •
•
written
by
User 6659 •
if you choose to go with the database schema directly and not use the api you are going to lose all of the 'protection' the api gives you such as sorting the failed variation features from the good ones
• 0 views
•
link
written
by
User 6659 •
but in the web search ~600 (!) SNPs have been duplicated. e.g: I counted "rs653534" 19 times.
• 0 views
•
link
written
by
Pierre Lindenbaum
Log in to answer this question.
More posts like this
-
Ensembl API - Global Ortholog Alignments
written by Nikit •
-
Mapping NCBI gene id with Bioproject id
written by preethibagopi28 •
-
Teaching a course on RESTful microservices; looking for example projects
written by Fukusahachi •
-
How to get Population Genotype Frequency from 1000 genomes Perl API
written by mass2 •
-
News: Upcoming Workshop: Ensembl Perl APIs
written by Denise CS
-
Best method to retrieve genome sequences
written by kay.hellboy •
-
Convert Pmids To Nihmsid And Pmcid Using Eutils
written by Smandape
-
Retrieving All Available Frequency Data For A Snp Using Ensembl Api Tools
written by Krisr
-
Is There Such A Thing As A Ucsc Api?
written by Steve Moss
-
How To Generate The Set Of Sequences At A Given Edit Distance From A Reference?
written by Jarretinha
