FBN1 Cohort Analysis: Correlation and Synergy
In this tutorial, we analyze a cohort of phenopackets associated with the gene FBN1. The example data is retrieved from the phenopacket store.
Using this dataset, we demonstrate a complete workflow with phenosign, including:
dataset construction
correlation analysis
synergy analysis of phenotypic features
We recommend running this tutorial in a Jupyter notebook for interactive exploration and visualization, although it can also be executed as a standard Python script.
1. Load phenopackets
We start by loading phenopackets for the FBN1 cohort:
[1]:
from ppktstore.registry import configure_phenopacket_registry
registry = configure_phenopacket_registry()
with registry.open_phenopacket_store() as ps:
phenopackets = [
p for p in ps.iter_cohort_phenopackets("FBN1")
]
print(f"Loaded {len(phenopackets)} phenopackets")
Loaded 143 phenopackets
The function load_phenopackets_by_cohort provides convenient access to publicly available datasets from the phenopacket store.
Note：You can also use your own phenopacket data. As long as your data follows the GA4GH Phenopacket schema, it can be directly used in the workflow.
2. Build the dataset
Next, we convert phenopackets into a structured dataset suitable for downstream statistical analysis.
[2]:
from phenosign import PhenotypeDatasetBuilder
dataset = PhenotypeDatasetBuilder(phenopackets).build()
Individuals Processed: 100%|██████████| 143/143 [01:27<00:00, 1.64 individuals/s]
The parameters specify the transcript of interest, the variant class to include. Detailed explanations are provided in the Usage section.
The resulting dataset is now ready for correlation and synergy analysis.
3. Correlation analysis
Next, we compute pairwise correlations between HPO features across individuals in the cohort.
[ ]:
from phenosign import HPOCorrelationAnalyzer
analyzer = HPOCorrelationAnalyzer(dataset=dataset)
results = analyzer.compute_correlation_matrix()
results.results_table.head()
Calculating pairwise correlation: 100%|██████████| 728/728 [00:01<00:00, 394.09it/s]
| HPO_A | HPO_A_label | HPO_B | HPO_B_label | correlation | p_value | adj_p_value | n(A:E/B:E) | n(A:E/B:O) | n(A:O/B:E) | n(A:O/B:O) | n_individuals | n_pmids | pmids | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 467 | HP:0004322 | Short stature | HP:0004942 | Aortic aneurysm | -0.827131 | 2.049015e-14 | 9.773799e-12 | 7 | 34 | 33 | 0 | 74 | 5 | 10756346;11175294;12203992;20375004;21683322 |
| 50 | HP:0000098 | Tall stature | HP:0004942 | Aortic aneurysm | 0.796637 | 1.940118e-13 | 4.627183e-11 | 45 | 7 | 1 | 26 | 79 | 13 | 10756346;11175294;12203992;20375004;20979188;2... |
| 48 | HP:0000098 | Tall stature | HP:0004322 | Short stature | -0.685681 | 5.424742e-12 | 8.625340e-10 | 16 | 33 | 37 | 0 | 86 | 5 | 10756346;11175294;12203992;20375004;21683322 |
| 44 | HP:0000098 | Tall stature | HP:0002616 | Aortic root aneurysm | 0.646903 | 1.766457e-11 | 2.106500e-09 | 70 | 7 | 10 | 26 | 113 | 18 | 10756346;11175294;12203992;20375004;20979188;2... |
| 439 | HP:0002616 | Aortic root aneurysm | HP:0004322 | Short stature | -0.634615 | 2.474657e-10 | 2.360823e-08 | 19 | 33 | 33 | 0 | 85 | 5 | 10756346;11175294;12203992;20375004;21683322 |
This step identifies pairs of HPO terms that tend to co-occur or show mutually exclusive patterns.
In the result table:
HPO_AandHPO_Bare the two HPO terms being comparedcorrelationindicates the strength and direction of associationp_value_correctedprovides the adjusted significance level
Interpretation:
positive correlation → the two phenotypes tend to appear together
negative correlation → the phenotypes tend to occur in different individuals
values near zero → little or no association
Detailed descriptions of parameters and output columns are provided in the Usage section.
4. Visualize correlation results
To better interpret these relationships, we visualize them as a heatmap.
[8]:
import os
from IPython.display import HTML
STATIC_DIR = "_static"
os.makedirs(STATIC_DIR, exist_ok=True)
def plotly_html_link(fig, filename, link_text="View interactive heatmap"):
"""
Generate Plotly HTML file and return an RTD-friendly link.
"""
filepath = os.path.join(STATIC_DIR, filename)
fig.write_html(filepath, include_plotlyjs="cdn", full_html=True)
return HTML(f'<p>To view the interactive plot, click the link below:</p>'
f'<a href="_static/{filename}" target="_blank">{link_text}</a>')
[5]:
fig1 = results.plot_correlation_heatmap_with_significance(
title_name="Cohort FBN1",
)
[6]:
plotly_html_link(fig1, "correlation_heatmap.html", link_text="Click here to view interactive correlation heatmap")
[6]:
To view the interactive plot, click the link below:
Click here to view interactive correlation heatmapThis visualization highlights the strongest and most statistically significant associations between phenotypic features.
Strong positive correlations indicate phenotypes that frequently co-occur
Strong negative correlations indicate phenotypes that tend to occur in different individuals
Note：Hover over the heatmap to see detailed information for each phenotype pair.The thresholds control which interactions are displayed. Lower thresholds include more pairs, while higher thresholds focus on the strongest signals.
5. Inspect available targets
Before running synergy analysis, we inspect the available target variables in the dataset.
[7]:
diseases_df, sex_df, genes_df, variant_effects_df = dataset.describe_conditions()
diseases_df
[7]:
| label | n_individuals | |
|---|---|---|
| disease_id | ||
| OMIM:154700 | Marfan syndrome | 50 (35.0%) |
| OMIM:129600 | Ectopia lentis, familial | 44 (30.8%) |
| OMIM:614185 | Geleophysic dysplasia 2 | 19 (13.3%) |
| OMIM:102370 | Acromicric dysplasia | 13 (9.1%) |
| OMIM:616914 | Marfan lipodystrophy syndrome | 9 (6.3%) |
| OMIM:184900 | Stiff skin syndrome | 8 (5.6%) |
[8]:
sex_df
[8]:
| n_individuals | |
|---|---|
| sex | |
| female | 57 (39.9%) |
| male | 54 (37.8%) |
| unknown | 32 (22.4%) |
[9]:
genes_df
[9]:
| n_individuals | |
|---|---|
| gene_symbol | |
| FBN1 | 143 (100.0%) |
[10]:
variant_effects_df
[10]:
| transcript_id | NM_000138.5 | NM_001406716.1 | NM_001406717.1 |
|---|---|---|---|
| variant_effect | |||
| MISSENSE_VARIANT | 119 (78.8%) | 119 (78.8%) | 8 (88.9%) |
| SPLICE_REGION_VARIANT | 6 (4.0%) | 6 (4.0%) | 1 (11.1%) |
| STOP_GAINED | 4 (2.6%) | 4 (2.6%) | 0 (0.0%) |
| SPLICE_DONOR_VARIANT | 6 (4.0%) | 6 (4.0%) | 0 (0.0%) |
| INFRAME_DELETION | 1 (0.7%) | 1 (0.7%) | 0 (0.0%) |
| INFRAME_INSERTION | 1 (0.7%) | 1 (0.7%) | 0 (0.0%) |
| FRAMESHIFT_VARIANT | 7 (4.6%) | 7 (4.6%) | 0 (0.0%) |
| SPLICE_DONOR_5TH_BASE_VARIANT | 1 (0.7%) | 1 (0.7%) | 0 (0.0%) |
| INTRON_VARIANT | 2 (1.3%) | 2 (1.3%) | 0 (0.0%) |
| SYNONYMOUS_VARIANT | 1 (0.7%) | 1 (0.7%) | 0 (0.0%) |
| SPLICE_DONOR_REGION_VARIANT | 1 (0.7%) | 1 (0.7%) | 0 (0.0%) |
| SPLICE_ACCEPTOR_VARIANT | 2 (1.3%) | 2 (1.3%) | 0 (0.0%) |
This summary shows which target variables (e.g., disease labels or variant conditions) can be used for downstream analysis.
These targets define the conditions under which phenotype–phenotype relationships are evaluated in the synergy analysis.
In other words, synergy analysis asks whether the association between two phenotypes changes across different conditions (e.g., variant classes or disease groups).
6. Synergy analysis
While correlation captures overall associations between phenotypes, it does not account for how these relationships may differ across conditions.
Synergy analysis addresses this by evaluating whether the association between two phenotypes changes depending on a target variable (e.g., variant class or disease group).
We begin by initializing the synergy analyzer:
[3]:
from phenosign import SynergyAnalyzer
synergy_analyzer = SynergyAnalyzer(dataset=dataset)
Next, we compute the synergy matrix for the selected target:
[4]:
from gpsea.model import VariantEffect
from phenosign import has_variant_effect
condition = dataset.get_variant_condition(has_variant_effect(transcript_id = "NM_000138.5",
variant_effect=VariantEffect.MISSENSE_VARIANT)
)
[5]:
results_variant = synergy_analyzer.compute_synergy_matrix(
condition=condition
)
results_variant.results_table.head()
Calculating pairwise synergy: 100%|██████████| 606/606 [03:05<00:00, 3.27it/s]
[5]:
| HPO_A | HPO_A_label | HPO_B | HPO_B_label | synergy | p_value | adj_p_value | n(A:E/B:E)_y0 | n(A:E/B:O)_y0 | n(A:O/B:E)_y0 | n(A:O/B:O)_y0 | N_y0 | n(A:E/B:E)_y1 | n(A:E/B:O)_y1 | n(A:O/B:E)_y1 | n(A:O/B:O)_y1 | N_y1 | n_individuals | n_pmids | pmids | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 39 | HP:0000098 | Tall stature | HP:0002650 | Scoliosis | 0.148600 | 0.0002 | 0.016157 | 2 | 4 | 7 | 2 | 15 | 36 | 1 | 12 | 16 | 65 | 80 | 11 | 10756346;11175294;12203992;20375004;22219643;2... |
| 33 | HP:0000098 | Tall stature | HP:0001382 | Joint hypermobility | 0.106631 | 0.0002 | 0.016157 | 2 | 4 | 6 | 2 | 14 | 58 | 2 | 12 | 17 | 89 | 103 | 12 | 10756346;11175294;12203992;20979188;21594992;2... |
| 31 | HP:0000098 | Tall stature | HP:0001187 | Hyperextensibility of the finger joints | 0.041146 | 0.0002 | 0.016157 | 2 | 3 | 6 | 0 | 11 | 58 | 0 | 12 | 0 | 70 | 81 | 12 | 10756346;11175294;12203992;20979188;21594992;2... |
| 277 | HP:0001187 | Hyperextensibility of the finger joints | HP:0001519 | Disproportionate tall stature | 0.028099 | 0.0002 | 0.016157 | 4 | 5 | 3 | 0 | 12 | 67 | 4 | 0 | 0 | 71 | 83 | 11 | 10756346;11175294;12203992;20979188;21594992;2... |
| 291 | HP:0001382 | Joint hypermobility | HP:0001519 | Disproportionate tall stature | 0.078793 | 0.0002 | 0.016157 | 4 | 5 | 5 | 1 | 15 | 67 | 4 | 9 | 10 | 90 | 105 | 11 | 10756346;11175294;12203992;20979188;21594992;2... |
[9]:
from phenosign import has_disease
condition_disease = dataset.get_condition(has_disease("OMIM:154700"), name="disease:Marfan syndrome")
[10]:
results_disease = synergy_analyzer.compute_synergy_matrix(
condition=condition_disease
)
results_disease.results_table.head()
Calculating pairwise synergy: 100%|██████████| 606/606 [01:44<00:00, 5.78it/s]
[10]:
| HPO_A | HPO_A_label | HPO_B | HPO_B_label | synergy | p_value | adj_p_value | n(A:E/B:E)_y0 | n(A:E/B:O)_y0 | n(A:O/B:E)_y0 | n(A:O/B:O)_y0 | N_y0 | n(A:E/B:E)_y1 | n(A:E/B:O)_y1 | n(A:O/B:E)_y1 | n(A:O/B:O)_y1 | N_y1 | n_individuals | n_pmids | pmids | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 180 | HP:0001519 | Disproportionate tall stature | HP:0004959 | Descending thoracic aorta aneurysm | 0.073524 | 0.0002 | 0.00103 | 45 | 0 | 0 | 0 | 45 | 0 | 1 | 5 | 0 | 6 | 51 | 8 | 10756346;20375004;21594992;21683322;24039054;2... |
| 181 | HP:0001519 | Disproportionate tall stature | HP:0008848 | Moderately short stature | -0.045049 | 0.0002 | 0.00103 | 27 | 9 | 0 | 0 | 36 | 29 | 0 | 21 | 0 | 50 | 86 | 5 | 10756346;11175294;12203992;20375004;21683322 |
| 229 | HP:0003502 | Mild short stature | HP:0003510 | Severe short stature | 0.016259 | 0.0002 | 0.00103 | 16 | 18 | 2 | 0 | 36 | 50 | 0 | 0 | 0 | 50 | 86 | 5 | 10756346;11175294;12203992;20375004;21683322 |
| 163 | HP:0001382 | Joint hypermobility | HP:0003510 | Severe short stature | -0.133008 | 0.0002 | 0.00103 | 11 | 18 | 0 | 0 | 29 | 26 | 0 | 22 | 0 | 48 | 77 | 4 | 10756346;11175294;12203992;21683322 |
| 164 | HP:0001382 | Joint hypermobility | HP:0004322 | Short stature | -0.251557 | 0.0002 | 0.00103 | 0 | 32 | 0 | 0 | 32 | 26 | 0 | 22 | 0 | 48 | 80 | 4 | 10756346;11175294;12203992;21683322 |
The resulting table reports pairwise synergy scores between HPO terms with respect to the selected target.
HPO_AandHPO_Bare the two phenotypes being evaluatedsynergymeasures how much additional information the pair provides about the target compared to individual featuresp_value_correctedindicates statistical significance after multiple testing correction
Interpretation:
positive synergy → the two phenotypes jointly provide additional information about the target
near zero → the phenotypes contribute largely independently
negative synergy → the phenotypes are redundant with respect to the target
Detailed descriptions of all output columns are provided in the Usage section.
7. Visualize synergy results
We can visualize synergy results as a heatmap.
[11]:
fig2 = results_variant.plot_synergy_heatmap(
)
[12]:
plotly_html_link(fig2, "synergy_heatmap.html", link_text="Click here to view interactive synergy heatmap")
[12]:
To view the interactive plot, click the link below:
Click here to view interactive synergy heatmap[13]:
fig3 = results_disease.plot_synergy_heatmap(
)
[14]:
plotly_html_link(fig3, "synergy_heatmap_disease.html", link_text="Click here to view interactive synergy heatmap for Marfan syndrome")
[14]:
To view the interactive plot, click the link below:
Click here to view interactive synergy heatmap for Marfan syndromeThis visualization highlights phenotype pairs that show strong and statistically significant synergy with respect to the selected target.
Note：Hover over the heatmap to see detailed information for each phenotype pair.The thresholds control which interactions are displayed. Lower thresholds include more pairs, while higher thresholds focus on the strongest signals.
Summary
In this tutorial, we:
Loaded phenopacket data
Constructed a structured dataset
Explored phenotype–phenotype relationships using correlation
Identified condition-specific interactions using synergy analysis
Together, these steps provide a workflow for uncovering both global and condition-dependent relationships between phenotypic features.
For additional usage patterns and parameter options, see the Usage section.