Supplementary Figures 5 and 7 correspond to Section 2. Supplementary Figure 6 is first mentioned in the second results section of the paper, but its content primarily relates to Section 3.
These figures are separate for page size reason, to keep a reasonable loading time.
The formatting of the figures may differ slightly from those in the paper, but they display the same data points.
All code cells are folded by default. To view any cell, click “Code” to expand it, or use the code options near the main title above to unfold all at once.
Some code may be repeated, as the original Python notebook was designed for figures to be generated semi-independently.
cross_val_analysis = pd.merge( concat_results_10fold, chrY_df, left_index=True, right_on="filename", suffixes=("", "_DROP") )cross_val_analysis.drop( columns=[c for c in cross_val_analysis.columns if c.endswith("_DROP")], inplace=True)
Define function zscore_per_assay to compute and graph the metric for each assay instead of globally.
Supp. Fig. 5B: Distribution of average z-score signal of epigenomes (dots) over chrY per sex (female in red, male in blue) for each assay individually (showing only the fold change track type for the ChIP datasets, and the two types of WGBS and RNA-seq were merged). Dashed lines represent means, solid lines the medians, boxes the quartiles, and whiskers the farthest points within 1.5× the interquartile range.
C - Female/Male chrY signal z-score cluster separation
Define function merged_assays_separation_distance that computes and graphs the showing separation distance between male/female zscore clusters.
Supp. Fig. 5C: Effect of a prediction score threshold on the aggregated mean (blue) and median (green) sex z-score male/female cluster distances, as well as and corresponding file subset size (red) of ChIP-related assays from panel B.
Images extracted from Epilogos viewer, using specified coordinates (XIST and FIRRE positions), and:
View mode: Paired
Dataset: IHEC
Pairwise: Male VS Female 100 samples
Saliency Metric: S1
Supp. Fig. 5E: Epilogos pairwise comparisons of male (top) vs female (bottom) showing portions of important regions for the Sex classifier, including the XIST (left) and FIRRE (right) genes.
See Annex A for a more detailled Epilogos color legend.
F - Genome browser for biospecimen important regions
Supp. Fig. 5F: Genome browser representation of the important regions shown in Figure 2I.
G - RNA-Seq read handling: Unique vs UniqueMultiple vs Unstranded
The RNA-Seq classifiers were trained on the Unique representation (multi-mapped reads excluded). This panel compares three read-handling variants of the same RNA samples on every non-assay task: Unique, UniqueMultiple (multi-mapped reads included) and Unstranded (summed) (the two Unique stranded tracks summed into one signal per EpiRR). Each fold’s model is applied to its held-out RNA validation set for the three variants, so the comparison is on samples the model did not train on and whose true label is in the classifier’s class space.
The predictions were produced by src/python/epiclass/utils/notebooks/predictions/make_rna_read_handling_predictions.py, which pairs the three variants per sample (resolving md5sums to EpiRRs and restricting to each fold’s validation set) and writes a long-form predictions table. Every prediction carries the Comet experiment key of the classifier that produced it. The figure below reads that table and scores Accuracy and macro-F1 per task, fold and variant.
Code
rna_variations_dir = base_data_dir /"rna_variations"PathChecker.check_directory(rna_variations_dir)# The file starts with a '#'-commented provenance header (metadata source of# true_class, model identity); comment="#" skips it.rna_preds = pd.read_csv( rna_variations_dir /"rna_read_handling_predictions.csv", comment="#")def score_read_handling(group: pd.DataFrame) -> pd.Series:"""Accuracy and macro-F1 of one (task, fold, variant) group vs the true label.""" labels =sorted(group["true_class"].unique())return pd.Series( {"Accuracy": accuracy_score(group["true_class"], group["predicted_class"]),"F1_macro": f1_score( group["true_class"], group["predicted_class"], labels=labels, average="macro", zero_division=0, ), } )rna_metrics = ( rna_preds.groupby(["model", "split", "mapping"]) .apply(score_read_handling) .reset_index())
Code
# Canonical task order + display names; assay is excluded (mRNA/total-RNA confusion# is not relevant to this comparison).TASK_ORDER = ["cancer", "life_stage", "sex", "cell_type", "biomaterial"]TASK_LABELS = {"cancer": "Cancer status","life_stage": "Life stage","sex": "Sex","cell_type": "Biospecimen","biomaterial": "Biomaterial type",}MAPPINGS = ["Unique", "Unstranded (summed)", "UniqueMultiple"]# Okabe-Ito colourblind-safe palette.MAPPING_COLORS = {"Unique": "#0072B2","UniqueMultiple": "#D55E00","Unstranded (summed)": "#009E73",}metrics = rna_metrics[rna_metrics["model"] !="assay"]present =set(metrics["model"])task_order_labels = [TASK_LABELS[t] for t in TASK_ORDER if t in present]fig = make_subplots( rows=1, cols=2, subplot_titles=["Accuracy", "F1_macro"], horizontal_spacing=0.08)for col, metric inenumerate(["Accuracy", "F1_macro"], start=1):for mapping in MAPPINGS: sub = metrics[metrics["mapping"] == mapping] fig.add_trace( go.Box( x=sub["model"].map(TASK_LABELS), y=sub[metric], name=mapping, fillcolor=MAPPING_COLORS[mapping], line=dict(color="black", width=1.2), marker=dict(size=3, color="white", line_width=1), boxmean=True, boxpoints="all", pointpos=0,# Show the fold each point comes from on hover (as in Fig 2A,B). hovertemplate="%{text}", text=[f"{split}: {value:.4f}"for split, value inzip(sub["split"], sub[metric]) ], legendgroup=mapping, showlegend=col ==1, ), row=1, col=col, )fig.update_xaxes(categoryorder="array", categoryarray=task_order_labels)# Semi auto-scaled: fix the top, let the bottom auto-scale to the data.fig.update_yaxes(autorange="min", range=[None, 1.02], row=1, col=1)fig.update_yaxes(autorange="min", range=[None, 1.02], row=1, col=2)fig.update_layout( boxmode="group", template="plotly_white", yaxis_title="Value", height=550, width=1400, legend=dict(orientation="h", y=-0.15),)fig.show()
Supp. Fig. 5G: Accuracy and macro-F1 performance metrics using only RNA-seq datasets across three data types and five classifiers, which were trained exclusively on strand-specific signals. Each fold’s trained model was applied to its own validation samples (Unique, blue), after summing forward and reverse tracks (Unstranded, green), and on a different version of the RNA-seq data downloaded from EpiATLAS where multi-mapped reads were retained (UniqueMultiple, orange) used only for this particular analysis. Changes in Unstranded and UniqueMultiple relative to Unique correspond to the impact of the data type on the inference since the same trained model was applied to all three data types. Each point is the score from one of the 10 cross-validation folds, computed across all signal files in that fold’s validation set. Solid lines are medians, dashed lines means, boxes the quartiles, and whiskers the farthest points within 1.5× the interquartile range.
For full code, since the processing is more complex, see src/python/epiclass/utils/notebooks/paper/confidence_threshold.ipynb (permalink).
Supplementary Figure 6: Impact of prediction score threshold on performance metrics. Impact of prediction score threshold on accuracy, F1-score and number of files for both EpiATLAS cross-validation performance and inference on datasets from other databases with provided or extracted labels. Performance for Assay, Sex, Cancer, Biomaterial type and Life stage classifiers are shown, for EpiATLAS, ENCODE core/non-core, ChIP-Atlas and Recount3 datasets. The number of classes (C) and the number of files analyzed (N) used to calculate the performances are shown at the bottom for each graph. The 11 classes of the Assay classifiers for EpiATLAS correspond to the six ChIP-Seq histone modifications, their control Input file, and two protocols of both RNA-Seq and WGBS, while for the 9 classes of ENCODE the two protocols were grouped, and indeed only the seven ChIP-related assays were used for ChIP-Atlas and RNA-Seq for Recount3. For the Sex classifier the third class corresponds to ‘mixed’, absent for ENCODE. The Cancer classifier is binary (where non-cancer is a mix of healthy and other diseases). For the Biomaterial classifier the ‘primary cell culture’ class is missing from all public sources (but ‘primary cell’, ‘primary tissue’ and ‘cell line’ are present), while the three classes (perinatal, pediatric, adult) were always used for the Life stage classifier.
Supplementary Figure 7 - Biospecimen classifier - ChromScore for high-SHAP regions
For full code, since the processing is more complex, see src/python/epiclass/utils/notebooks/paper/chromscore.ipynb (permalink).
