diff --git a/docs/report/tbx5_all_variants.html b/docs/report/tbx5_all_variants.html new file mode 100644 index 000000000..a96181629 --- /dev/null +++ b/docs/report/tbx5_all_variants.html @@ -0,0 +1,583 @@ + + + + + Cohort + + + + +

GPSEA cohort analysis: All variant alleles

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+

Variant alleles

+ A total of 53 unique alleles were identified in the cohort. +
Variant keyVariant (cDNA)Variant (protein)EffectsCount
12_114385521_114385521_C_Tc.710G>Ap.Arg237Glnmissense22
12_114401830_114401830_C_Tc.238G>Ap.Gly80Argmissense20
12_114385563_114385563_G_Ac.668C>Tp.Thr223Metmissense8
12_114398675_114398675_G_Tc.408C>Ap.Tyr136Terstop gained6
12_114399514_114399514_A_Cc.361T>Gp.Trp121Glymissense, splice region5
12_114403792_114403792_C_CGc.106_107insCp.Ser36ThrfsTer25frameshift5
12_114398682_114398682_C_CGc.400dupp.Arg134ProfsTer49frameshift5
12_114385474_114385474_A_Gc.755+2T>CNonesplice donor4
12_114398656_114398656_C_CGc.426dupp.Ala143ArgfsTer40frameshift4
12_114366360_114366360_C_Tc.787G>Ap.Val263Metmissense4
12_114385522_114385522_G_Ac.709C>Tp.Arg237Trpmissense4
12_114403798_114403798_G_GCc.100dupp.Ala34GlyfsTer27frameshift4
12_114399613_114399613_T_Ac.262A>Tp.Lys88Terstop gained3
12_114401827_114401827_T_Ac.241A>Tp.Arg81Trpmissense, splice region3
12_114366312_114366312_G_Ac.835C>Tp.Arg279Terstop gained3
12_114366366_114366366_T_Ac.781A>Tp.Ser261Cysmissense3
12_114403798_114403799_GC_Gc.100delp.Ala34ProfsTer32frameshift3
12_114385475_114385475_C_Tc.755+1G>ANonesplice donor3
12_114401853_114401853_G_Tc.215C>Ap.Thr72Lysmissense3
12_114385521_114385521_C_Gc.710G>Cp.Arg237Promissense2
12_114398568_114398568_C_Ac.510+5G>TNonesplice donor 5th base, intronic2
12_114366207_114366208_GC_Gc.939delp.Gln315ArgfsTer79frameshift2
12_114366274_114366274_G_Tc.873C>Ap.Tyr291Terstop gained2
12_114366267_114366267_C_Ac.880G>Tp.Glu294Terstop gained2
12_114398578_114398579_CA_Cc.504delp.Phe168LeufsTer6frameshift2
12_114394762_114394763_CA_Cc.641delp.Val214GlyfsTer12frameshift2
12_114398666_114398667_TG_Tc.416delp.Pro139GlnfsTer11frameshift2
12_114403754_114403754_G_Tc.145C>Ap.Gln49Lysmissense, splice region2
12_114385550_114385550_A_AATTATTCTCAGc.680_681insCTGAGAATAATp.Ile227_Glu228insTerinframe insertion, stop retainined2
12_114366348_114366349_CT_Cc.798delp.Val267TrpfsTer127frameshift1
12_114399559_114399559_T_Cc.316A>Gp.Ile106Valmissense1
12_114356064_114356065_TA_Tc.1024delp.Tyr342ThrfsTer52frameshift1
12_114355755_114355756_TG_Tc.1333delp.His445MetfsTer137frameshift1
12_114398632_114398632_G_Ac.451C>Tp.Gln151Terstop gained1
12_114401907_114401907_A_Gc.161T>Cp.Ile54Thrmissense1
12_114355723_114355723_G_Ac.1366C>Tp.Gln456Terstop gained1
12_114398602_114398602_T_Gc.481A>Cp.Thr161Promissense1
12_114398708_114398709_GC_Gc.374delp.Gly125AlafsTer25frameshift1
12_114385553_114385553_C_Ac.678G>Tp.Lys226Asnmissense1
12_114399633_114399633_C_Gc.243-1G>CNonesplice acceptor1
12_114401873_114401874_TA_Tc.194delp.Leu65GlnfsTer10frameshift1
12_114399594_114399594_A_Cc.281T>Gp.Leu94Argmissense1
12_114399622_114399622_G_Tc.253C>Ap.Pro85Thrmissense1
12_114394743_114394746_TGTG_Tc.658_660delp.His220delinframe deletion1
12_114394820_114394820_C_Gc.584G>Cp.Gly195Alamissense1
12_114401846_114401846_C_Gc.222G>Cp.Met74Ilemissense1
12_114366241_114366242_CT_Cc.905delp.Gln302ArgfsTer92frameshift1
12_114355784_114355785_CA_Cc.1304delp.Leu435ArgfsTer147frameshift1
12_114398626_114398627_CG_Cc.456delp.Val153SerfsTer21frameshift1
12_114403859_114403859_G_Tc.40C>Ap.Pro14Thrmissense1
12_114399625_114399629_ACATC_Ac.246_249delp.Met83PhefsTer6frameshift1
12_114394817_114394817_G_Cc.587C>Gp.Ser196Terstop gained1
12_114401921_114401921_C_Gc.148-1G>CNonesplice acceptor1
+ + + + + \ No newline at end of file diff --git a/docs/tutorial.rst b/docs/tutorial.rst index 8917ce46a..e448e8c44 100644 --- a/docs/tutorial.rst +++ b/docs/tutorial.rst @@ -103,7 +103,13 @@ We loaded the patient data into a `cohort` which is ready for the next steps. Explore cohort ^^^^^^^^^^^^^^ -We can now explore the cohort to see how many patients are included. +GPSEA helps with gaining insight into the cohort by providing + + +Show cohort summary +------------------- + +The summary report provides an overview about the HPO terms, variants, diseases, and variant effects that occurr most frequently: >>> from gpsea.view import CohortViewable >>> viewer = CohortViewable(hpo) @@ -121,6 +127,10 @@ We can now explore the cohort to see how many patients are included. from IPython.display import HTML, display display(HTML(report)) + +Plot distribution of variants with respect to the protein sequence +------------------------------------------------------------------ + Now we can show the distribution of variants with respect to the encoded protein. We first obtain `tx_coordinates` (:class:`~gpsea.model.TranscriptCoordinates`) and `protein_meta` (:class:`~gpsea.model.ProteinMetadata`) @@ -154,6 +164,27 @@ and we follow with plotting the diagram of the mutations on the protein: :width: 600px +.. _show-cohort-variants: + +Summarize all variant alleles +----------------------------- + +We can prepare a table of all variant alleles that occurr in the cohort. +Each table row corresponds to a single allele and lists the variant key, +the predicted effect on the transcript (*cDNA*) and protein of interest, +the variant effects, and the number of patients who present +with one or more variant alleles (*Count*): + +>>> from gpsea.view import CohortVariantViewer +>>> viewer = CohortVariantViewer(tx_id=tx_id) +>>> report = viewer.process(cohort=cohort) +>>> with open('docs/report/tbx5_all_variants.html', 'w') as fh: # doctest: +SKIP +... _ = fh.write(report) + +.. raw:: html + :file: report/tbx5_all_variants.html + + Prepare genotype and phenotype predicates ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/gpsea/model/_variant_effects.py b/src/gpsea/model/_variant_effects.py index e15452e08..9567592d4 100644 --- a/src/gpsea/model/_variant_effects.py +++ b/src/gpsea/model/_variant_effects.py @@ -1,7 +1,10 @@ -from enum import Enum +import enum +import typing +import hpotk -class VariantEffect(Enum): + +class VariantEffect(enum.Enum): """ `VariantEffect` represents consequences of a variant on transcript that are supported by GPSEA. @@ -47,7 +50,7 @@ class VariantEffect(Enum): THREE_PRIME_UTR_VARIANT = "SO:0001624" NON_CODING_TRANSCRIPT_EXON_VARIANT = "SO:0001792" INTRON_VARIANT = "SO:0001627" - NMD_TRANSCRIPT_VARIANT = "SO:0001621", + NMD_TRANSCRIPT_VARIANT = "SO:0001621" NON_CODING_TRANSCRIPT_VARIANT = "SO:0001619" UPSTREAM_GENE_VARIANT = "SO:0001631" DOWNSTREAM_GENE_VARIANT = "SO:0001632" @@ -62,6 +65,43 @@ class VariantEffect(Enum): INTERGENIC_VARIANT = "SO:0001628" SEQUENCE_VARIANT = "SO:0001060" + def to_display(self) -> str: + """ + Get a concise name of the variant effect that is suitable for showing to humans. + + Example + ^^^^^^^ + + >>> from gpsea.model import VariantEffect + >>> VariantEffect.MISSENSE_VARIANT.to_display() + 'missense' + >>> VariantEffect.SPLICE_DONOR_5TH_BASE_VARIANT.to_display() + 'splice donor 5th base' + + :returns: a `str` with the name or `'n/a'` if the variant effect was not assigned a concise name. + """ + return effect_to_display.get(self, "n/a") + + @staticmethod + def structural_so_id_to_display(so_term: typing.Union[hpotk.TermId, str]) -> str: + """ + Get a `str` with a concise name for representing a Sequence Ontology (SO) term identifier. + + Example + ^^^^^^^ + + >>> from gpsea.model import VariantEffect + >>> VariantEffect.structural_so_id_to_display('SO:1000029') + 'chromosomal deletion' + + :param so_term: a CURIE `str` or a :class:`~hpotk.TermId` with the query SO term. + :returns: a `str` with the concise name for the SO term or `'n/a'` if a name has not been assigned yet. + """ + if isinstance(so_term, hpotk.TermId): + so_term = so_term.value + + return so_id_to_display.get(so_term, "n/a") + def __init__(self, curie: str): self._curie = curie @@ -75,3 +115,54 @@ def curie(self) -> str: def __str__(self) -> str: return self.name.lower() + + +effect_to_display = { + VariantEffect.TRANSCRIPT_ABLATION: "transcript ablation", + VariantEffect.SPLICE_ACCEPTOR_VARIANT: "splice acceptor", + VariantEffect.SPLICE_DONOR_VARIANT: "splice donor", + VariantEffect.STOP_GAINED: "stop gained", + VariantEffect.FRAMESHIFT_VARIANT: "frameshift", + VariantEffect.STOP_LOST: "stop lost", + VariantEffect.START_LOST: "start lost", + VariantEffect.TRANSCRIPT_AMPLIFICATION: "transcript amplification", + VariantEffect.INFRAME_INSERTION: "inframe insertion", + VariantEffect.INFRAME_DELETION: "inframe deletion", + VariantEffect.MISSENSE_VARIANT: "missense", + VariantEffect.PROTEIN_ALTERING_VARIANT: "protein altering", + VariantEffect.SPLICE_REGION_VARIANT: "splice region", + VariantEffect.SPLICE_DONOR_5TH_BASE_VARIANT: "splice donor 5th base", + VariantEffect.SPLICE_DONOR_REGION_VARIANT: "splice donor", + VariantEffect.SPLICE_POLYPYRIMIDINE_TRACT_VARIANT: "splice polypyrimidine", + VariantEffect.INCOMPLETE_TERMINAL_CODON_VARIANT: "incomplete terminal codon", + VariantEffect.START_RETAINED_VARIANT: "start retained", + VariantEffect.STOP_RETAINED_VARIANT: "stop retainined", + VariantEffect.SYNONYMOUS_VARIANT: "synonymous", + VariantEffect.CODING_SEQUENCE_VARIANT: "coding sequence", + VariantEffect.MATURE_MIRNA_VARIANT: "mature miRNA", + VariantEffect.FIVE_PRIME_UTR_VARIANT: "5UTR", + VariantEffect.THREE_PRIME_UTR_VARIANT: "3UTR", + VariantEffect.NON_CODING_TRANSCRIPT_EXON_VARIANT: "non-coding transcript exon", + VariantEffect.INTRON_VARIANT: "intronic", + VariantEffect.NMD_TRANSCRIPT_VARIANT: "NMD transcript", + VariantEffect.NON_CODING_TRANSCRIPT_VARIANT: "non-coding transcript", + VariantEffect.UPSTREAM_GENE_VARIANT: "upstream of gene", + VariantEffect.DOWNSTREAM_GENE_VARIANT: "downstream of gene", + VariantEffect.TFBS_ABLATION: "TFBS ablation", + VariantEffect.TFBS_AMPLIFICATION: "TFBS amplification", + VariantEffect.TF_BINDING_SITE_VARIANT: "TFBS binding site", + VariantEffect.REGULATORY_REGION_ABLATION: "regulatory region ablation", + VariantEffect.REGULATORY_REGION_AMPLIFICATION: "regulatory region amplification", + VariantEffect.FEATURE_ELONGATION: "feature elongation", + VariantEffect.REGULATORY_REGION_VARIANT: "regulatory region", + VariantEffect.FEATURE_TRUNCATION: "feature truncation", + VariantEffect.INTERGENIC_VARIANT: "intergenic", + VariantEffect.SEQUENCE_VARIANT: "sequence variant", +} + +so_id_to_display = { + "SO:1000029": "chromosomal deletion", + "SO:1000037": "chromosomal duplication", + "SO:1000030": "chromosomal_inversion", + "SO:1000044": "chromosomal_translocation", +} diff --git a/src/gpsea/view/__init__.py b/src/gpsea/view/__init__.py index 7bb12e8e3..d6a6c22ec 100644 --- a/src/gpsea/view/__init__.py +++ b/src/gpsea/view/__init__.py @@ -1,4 +1,5 @@ from ._cohort import CohortViewable +from ._cohort_variant_viewer import CohortVariantViewer from ._disease import DiseaseViewable from ._phenotype_analysis import summarize_hpo_analysis from ._protein_viewer import ProteinViewable @@ -9,6 +10,7 @@ from ._formatter import VariantFormatter __all__ = [ + 'CohortVariantViewer', 'CohortViewable', 'ProteinVisualizer', 'ProteinVisualizable', 'ProteinViewable', 'DiseaseViewable', diff --git a/src/gpsea/view/_cohort_variant_viewer.py b/src/gpsea/view/_cohort_variant_viewer.py new file mode 100644 index 000000000..0933ce02e --- /dev/null +++ b/src/gpsea/view/_cohort_variant_viewer.py @@ -0,0 +1,158 @@ +import typing + +from jinja2 import Environment, PackageLoader +from collections import namedtuple, defaultdict + +from gpsea.model import Cohort, Variant, VariantEffect +from ._formatter import VariantFormatter + + +ToDisplay = namedtuple('ToDisplay', ['hgvs_cdna', 'hgvsp', 'variant_effects']) + +VariantData = namedtuple('VariantData', ['variant_key', 'hgvs_cdna', 'hgvsp', 'variant_effects']) + + +class CohortVariantViewer: + """ + `CohortVariantViewer` creates an HTML report with the cohort variants. + + The report can be either written into an HTML file or displayed in a Jupyter notebook. + + See :ref:show-cohort-variants: for an example usage. + """ + + def __init__( + self, + tx_id: str + ): + """ + Args: + tx_id (str): The transcript identifier (Usually, the MANE RefSeq transcript, that should start with "NM_") + """ + environment = Environment(loader=PackageLoader('gpsea.view', 'templates')) + self._cohort_template = environment.get_template("all_variants.html") + self._var_formatter = VariantFormatter(tx_id) + if not tx_id.startswith("NM"): + print(f"[WARNING] Non-RefSeq transcript id: {tx_id}") + self._transcript_id = tx_id + + def process( + self, + cohort: Cohort, + only_hgvs: bool = True + ) -> str: + """ + Create an HTML that should be shown with ``display(HTML(..))`` of the ipython package. + + Args: + cohort (Cohort): The cohort being analyzed in the current notebook. + only_hgvs (bool): Do not show the transcript ID part of the HGVS annotation, just the annotation. + + Returns: + str: an HTML string with parameterized template for rendering + """ + context = self._prepare_context(cohort, only_hgvs=only_hgvs) + return self._cohort_template.render(context) + + def _prepare_context( + self, + cohort: Cohort, + only_hgvs: bool, + ) -> typing.Mapping[str, typing.Any]: + variant_count_dictionary = defaultdict(int) + variant_effect_count_dictionary = defaultdict(int) + variant_key_to_variant = dict() + for var in cohort.all_variants(): + var_key = var.variant_info.variant_key + vdata = self._get_variant_data(var, only_hgvs) + variant_key_to_variant[var_key] = vdata + variant_count_dictionary[var_key] += 1 + for v_eff in vdata.variant_effects: + variant_effect_count_dictionary[v_eff] += 1 + variant_counts = list() + variant_effect_counts = list() + for var_key, count in variant_count_dictionary.items(): + var_data = variant_key_to_variant[var_key] + variant_counts.append( + { + "variant_key": var_data.variant_key, + "variant": var_data.hgvs_cdna, + "variant_name": var_data.hgvs_cdna, + "protein_name": var_data.hgvsp, + "variant_effects": ", ".join(var_data.variant_effects), + "count": count, + } + ) + for v_effect, count in variant_effect_count_dictionary.items(): + variant_effect_counts.append( + { + "effect": v_effect, + "count": count + } + ) + variant_counts = sorted(variant_counts, key=lambda row: row.get("count"), reverse=True) + variant_effect_counts = sorted(variant_effect_counts, key=lambda row: row.get("count"), reverse=True) + + # The following dictionary is used by the Jinja2 HTML template + return { + "has_transcript": False, + "variant_count_list": variant_counts, + "variant_effect_count_list": variant_effect_counts, + "total_unique_allele_count": len(variant_counts) + } + + def _get_variant_data( + self, + variant: Variant, + only_hgvs: bool + ) -> VariantData: + """ + Get user-friendly strings (e.g., HGVS for our target transcript) to match to the chromosomal strings + Args: + variant (Variant): The variant to be formatted. + only_hgvs (bool): do not show the transcript ID part of the HGVS annotation, just the annotation. + + Returns: + VariantData: a named tuple with variant data formatted for human consumption. + """ + variant_key = variant.variant_info.variant_key + if variant.variant_info.has_sv_info(): + sv_info = variant.variant_info.sv_info + gene_symbol = sv_info.gene_symbol + display = f"SV involving {gene_symbol}" + effect = VariantEffect.structural_so_id_to_display(so_term=sv_info.structural_type) + return VariantData( + variant_key=variant_key, + hgvs_cdna=display, + hgvsp="p.?", + variant_effects=[effect], + ) + else: + variant_key = variant.variant_info.variant_key + display = self._var_formatter.format_as_string(variant) + tx_annotation = variant.get_tx_anno_by_tx_id(self._transcript_id) + if tx_annotation is not None: + hgvsp = tx_annotation.hgvsp + var_effects = [var_eff.to_display() for var_eff in tx_annotation.variant_effects] + else: + hgvsp = None + var_effects = [] + if only_hgvs: + # do not show the transcript id + fields_dna = display.split(":") + if len(fields_dna) > 1: + display_hgvs_cDNA = fields_dna[1] + else: + display_hgvs_cDNA = fields_dna[0] + + fields_ps = hgvsp.split(":") if hgvsp is not None else (None,) + if len(fields_ps) > 1: + hgvsp = fields_ps[1] + else: + hgvsp = fields_ps[0] + return VariantData( + variant_key=variant_key, + hgvs_cdna=display_hgvs_cDNA, + hgvsp=hgvsp, + variant_effects=var_effects, + ) diff --git a/src/gpsea/view/templates/all_variants.html b/src/gpsea/view/templates/all_variants.html new file mode 100644 index 000000000..9640bd1af --- /dev/null +++ b/src/gpsea/view/templates/all_variants.html @@ -0,0 +1,117 @@ + + + + + Cohort + + + + +

GPSEA cohort analysis: All variant alleles

+ + + + + + + + + + + + + + + {% for vdata in variant_count_list %} + + + + + + + + {% endfor %} + +
+

Variant alleles

+ A total of {{ total_unique_allele_count }} unique alleles were identified in the cohort. +
Variant keyVariant (cDNA)Variant (protein)EffectsCount
{{ vdata.variant_key }}{{ vdata.variant_name }}{{ vdata.protein_name }}{{ vdata.variant_effects }}{{ vdata.count }}
+ + + + + diff --git a/tests/conftest.py b/tests/conftest.py index 64ed2e817..907b7d838 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -111,12 +111,19 @@ def suox_cohort( @pytest.fixture(scope='session') -def suox_gt_predicate() -> GenotypePolyPredicate: +def suox_mane_tx_id() -> str: + return 'NM_001032386.2' + + +@pytest.fixture(scope='session') +def suox_gt_predicate( + suox_mane_tx_id: str, +) -> GenotypePolyPredicate: # To bin the patients to a group with >1 MISSENSE variant or 0 MISSENSE variants. - suox_mane_tx_id = 'NM_001032386.2' + return boolean_predicate( variant_predicate=VariantPredicates.variant_effect( - effect=VariantEffect.MISSENSE_VARIANT, + effect=VariantEffect.MISSENSE_VARIANT, tx_id=suox_mane_tx_id ) ) diff --git a/tests/view/test_variant_viewer.py b/tests/view/test_variant_viewer.py new file mode 100644 index 000000000..31fe12bbd --- /dev/null +++ b/tests/view/test_variant_viewer.py @@ -0,0 +1,16 @@ +import pytest + +from gpsea.model import Cohort +from gpsea.view import CohortVariantViewer + + +@pytest.mark.skip("For manual run only") +def test_viewer( + suox_mane_tx_id: str, + suox_cohort: Cohort, +): + viewer = CohortVariantViewer(tx_id=suox_mane_tx_id) + html = viewer.process(suox_cohort) + + with open("all_variants.html", "w") as fh: + fh.write(html)