Commit c942e006 authored by Jeroen F.J. Laros's avatar Jeroen F.J. Laros Committed by Vermaat
Browse files

Add back translator interface

parent a98f48f1
"""
Mutalyzer interface for variant back translation (p. to c.).
"""
from __future__ import unicode_literals
from backtranslate.backtranslate import BackTranslate
from backtranslate.util import protein_letters_3to1
from Bio.Data import CodonTable
from extractor.variant import Allele, DNAVar, ISeq, ISeqList
from . import ncbi, Retriever
from .grammar import Grammar
from .output import Output
def backtranslate(output, description, table_id=1):
"""
Back translation of an amino acid substitution to all possible
causal one-nucleotide substitutions.
:arg object output: Output object for logging.
:arg unicode description: Amino acid substitution in HGVS format.
:arg int table_id: Translation table id.
:returns: List of DNA substitutions in HGVS format.
:rtype: list(unicode)
"""
# TODO: We currently only support NM and NP references, but in principle
# we could also support other references.
# TODO: For NP references where we don't find a link to the corresponding
# NM, we don't check if the specified reference amino acid corresponds
# to the NP reference sequence.
parse_tree = Grammar(output).parse(description)
if not parse_tree:
return []
if parse_tree.RefType != 'p':
output.addMessage(
__file__, 3, 'ENOPROT', 'Reference type is not p. (protein).')
return []
if parse_tree.RawVar.MutationType != 'subst':
output.addMessage(
__file__, 3, 'ENOSUBST', 'Variant is not a substitution.')
return []
if parse_tree.Version:
accession_number = '{}.{}'.format(parse_tree.RefSeqAcc, parse_tree.Version)
else:
accession_number = parse_tree.RefSeqAcc
position = int(parse_tree.RawVar.Main)
reference_amino_acid, amino_acid = parse_tree.RawVar.Args
if len(reference_amino_acid) == 3:
reference_amino_acid = protein_letters_3to1[reference_amino_acid]
if len(amino_acid) == 3:
amino_acid = protein_letters_3to1[amino_acid]
bt = BackTranslate(table_id)
# FIXME: Rancid workaround to silence fatal error raised by `loadrecord`.
output_ = Output('')
retriever = Retriever.GenBankRetriever(output_)
# The genbank retriever does not (yet) support protein references, but we
# cannot reliably distinguish between different reference types from the
# variant description before downloading the reference.
# Therefore, we just try to download the reference. This will succeed for
# transcript references, but fail for protein references.
# As a quick and dirty optimization, we shortcut this for accessions
# starting with 'NP_', of which we know that they are protein references.
# In the future we hope to support protein references directly.
if accession_number.startswith('NP_'):
genbank_record = None
else:
genbank_record = retriever.loadrecord(accession_number)
if not genbank_record:
# Assuming RefSeqAcc is an NP, try to get the corresponding NM.
version = int(parse_tree.Version) if parse_tree.Version else None
try:
transcript = ncbi.protein_to_transcript(parse_tree.RefSeqAcc,
version)
except ncbi.NoLinkError:
pass
else:
if transcript[1] is not None:
accession_number = '{}.{}'.format(*transcript)
else:
output.addMessage(
__file__, 2, 'WNOVERSION',
'Found corresponding nucleotide sequence, but note that '
'the version number is missing.')
accession_number = transcript[0]
genbank_record = retriever.loadrecord(accession_number)
offset = (position - 1) * 3
if genbank_record and genbank_record.molType == 'n':
# Only NM for now.
cds_loc = genbank_record.geneList[0].transcriptList[0].CDS.location
codon = genbank_record \
.seq[cds_loc[0] - 1:cds_loc[1]][offset:offset + 3]
forward_table = CodonTable.unambiguous_dna_by_id[table_id]
found_ref = forward_table.forward_table[unicode(codon)]
if reference_amino_acid != found_ref:
output.addMessage(
__file__, 3, 'EREF',
'{} not found at position {}, found {} instead.'.format(
reference_amino_acid, position, found_ref))
substitutions = bt.with_dna(unicode(codon), amino_acid)
else:
# Assume NP.
output.addMessage(
__file__, 2, 'WNODNA',
'Nucleotide reference sequence could not be found, using '
'protein fallback method.')
accession_number = 'UNKNOWN'
substitutions = bt.without_dna(reference_amino_acid, amino_acid)
if (reference_amino_acid, amino_acid) in bt.improvable():
output.addMessage(
__file__, 2, 'WIMPROVE',
'The back translation for this variant can be improved by '
'supplying a nucleotide reference sequence.')
return ['{}:c.{}'.format(accession_number, v)
for v in subst_to_hgvs(substitutions, offset)]
def subst_to_hgvs(substitutions, offset=0):
"""
Translate a set of substitutions to HGVS.
:arg dict substitutions: Set of single nucleotide substitutions indexed by
position.
:arg int offset: Codon position in the CDS.
:returns: Substitutions in HGVS format.
:rtype: set(Allele)
"""
variants = set()
for position in substitutions:
for substitution in substitutions[position]:
variants.add(Allele([DNAVar(
start=position + offset + 1,
end=position + offset + 1,
sample_start=position + offset + 1,
sample_end=position + offset + 1,
type='subst',
deleted=ISeqList([ISeq(sequence=substitution[0])]),
inserted=ISeqList([ISeq(sequence=substitution[1])])
)]))
return variants
{% extends "base.html" %}
{% set active_page = "back-translator" %}
{% set page_title = "Back Translator" %}
{% block content %}
<p class="alert alert-warning">
Please note that this is an experimental service.
</p>
<p>
Back translation from amino acid substitutions to nucleotide substitutions.
</p>
<p>
The underlying algorithm is also available as a command line utility and
programming library
(<a href="https://github.com/mutalyzer/backtranslate">source</a>,
<a href="https://pypi.python.org/pypi/backtranslate">package</a>).
</p>
<p>
Please supply an amino acid substitution.
</p>
<form class="form" action="{{ url_for('.back_translator') }}" method="get">
<div class="form-group">
<label for="description">Variant description</label>
<input class="form-control form-pre" type="text"
name="description" id="description" value="{{ description }}" placeholder="Variant description using HGVS format">
<p>Examples:
<code class="example-input" data-for="description">NM_003002.3:p.Asp92Tyr</code>,
<code class="example-input" data-for="description">NP_002993.1:p.Asp92Glu</code>,
<code class="example-input" data-for="description">NP_000000.0:p.Asp92Tyr</code>,
<code class="example-input" data-for="description">NP_000000.0:p.Leu92Phe</code>
</p>
</div>
<div class="form-group button-group">
<input type="submit" class="btn btn-primary" value="Back translate">
<a href="https://humgenprojects.lumc.nl/trac/mutalyzer/wiki/BackTranslator" target="new" class="btn btn-default pull-right">Help</a>
</div>
</form>
{% if description %}
<hr>
{% for m in messages %}
{% if m.class == "error" %}
<p class="alert alert-danger" title="{{ m.level }} (origin: {{ m.origin }})">{{ m.description }}</p>
{% elif m.class == "warning" %}
<p class="alert alert-warning" title="{{ m.level }} (origin: {{ m.origin }})">{{ m.description }}</p>
{% elif m.class == "information" %}
<p class="alert alert-info" title="{{ m.level }} (origin: {{ m.origin }})">{{ m.description }}</p>
{% elif m.class == "debug" %}
<p class="alert alert-info" title="{{ m.level }} (origin: {{ m.origin }})">{{ m.description }}</p>
{% endif %}
{% endfor %}
{% if summary == "0 Errors, 0 Warnings." %}
<p class="alert alert-success summary">{{ summary }}</p>
{% else %}
<p>{{summary}}</p>
{% endif %}
{% if not errors %}
<hr>
<h4>Input</h4>
<p><code>{{ description }}</code></p>
<h4>Possible variants</h4>
{% for variant in variants %}
{% if variant.startswith('UNKNOWN:') %}
<p><code>{{ variant }}</code></p>
{% else %}
<p><code><a href="{{ url_for('.name_checker', description=variant) }}">{{ variant }}</a></code></p>
{% endif %}
{% endfor %}
{% endif %}
{% endif %}{# variants #}
{% endblock content %}
......@@ -38,6 +38,7 @@
('website.snp_converter', {}, 'snp-converter', 'SNP Converter', False, False),
('website.name_generator', {}, 'name-generator', 'Name Generator', False, False),
('website.description_extractor', {}, 'description-extractor', 'Description Extractor', False, False),
('website.back_translator', {}, 'back-translator', 'Back Translator', False, False),
('website.reference_loader', {}, 'reference-loader', 'Reference File Loader', False, True),
('website.batch_jobs', {}, 'batch-jobs', 'Batch Jobs', True, False),
......
......@@ -23,8 +23,8 @@ from sqlalchemy.orm.exc import NoResultFound
import extractor
import mutalyzer
from mutalyzer import (announce, File, Retriever, Scheduler, stats,
util, variantchecker)
from mutalyzer import (announce, backtranslator, File, Retriever, Scheduler,
stats, util, variantchecker)
from mutalyzer.config import settings
from mutalyzer.db.models import BATCH_JOB_TYPES
from mutalyzer.db.models import Assembly, BatchJob
......@@ -682,6 +682,33 @@ def reference_loader_submit():
errors=errors)
@website.route('/back-translator')
def back_translator():
"""
Back translator.
"""
output = Output(__file__)
output.addMessage(
__file__, -1, 'INFO',
'Received Back Translate request from {}'.format(request.remote_addr))
stats.increment_counter('back-translator/website')
description = request.args.get('description')
variants = []
if description:
variants = backtranslator.backtranslate(output, description)
errors, warnings, summary = output.Summary()
messages = map(util.message_info, output.getMessages())
output.addMessage(__file__, -1, 'INFO', 'Finished Back Translate request')
return render_template(
'back-translator.html', errors=errors, summary=summary,
description=description or '', messages=messages, variants=variants)
@website.route('/description-extractor')
def description_extractor():
"""
......
-e git+https://github.com/LUMC/magic-python.git#egg=Magic_file_extensions
description-extractor==2.2.1
Flask==0.10.1
Jinja2==2.7.3
MySQL-python==1.2.5
......@@ -8,9 +7,11 @@ SQLAlchemy==0.9.8
Sphinx==1.2.3
Werkzeug==0.9.6
alembic==0.8.2
backtranslate==0.0.5
biopython==1.64
chardet==2.3.0
cssselect==0.9.1
description-extractor==2.2.1
lxml==3.4.0
mock==1.0.1
mockredispy==2.9.0.9
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment