diff --git a/bin/__plot_selectionfeatures.py b/bin/__plot_selectionfeatures.py deleted file mode 100755 index 51afae8d..00000000 --- a/bin/__plot_selectionfeatures.py +++ /dev/null @@ -1,899 +0,0 @@ -#!/usr/bin/env python - - -import pandas as pd -import json -import os, sys -import matplotlib.pyplot as plt -import seaborn as sns -import numpy as np -from matplotlib.ticker import MaxNLocator -from mpl_toolkits.axes_grid1.inset_locator import inset_axes -from read_utils import custom_na_values - -pd.set_option('display.max_columns', None) - -# Suppress warnings -import warnings -pd.options.mode.chained_assignment = None -warnings.filterwarnings("ignore", message="FixedFormatter should only be used together with FixedLocator") - - -import pandas as pd -import os, sys -import matplotlib.pyplot as plt -import seaborn as sns -import numpy as np -from scipy.stats import linregress,norm -import tabix -import matplotlib.cm as cm -import matplotlib.colors as mcolors -from matplotlib.patches import Patch -from matplotlib.lines import Line2D - -import sys, os -import pandas as pd -import numpy as np -import matplotlib.pyplot as plt - - -def plot_single_needle(gene, snvs_maf_obs, - seq_info_df, - sample = '' - ): - - snvs_maf_obs_gene = snvs_maf_obs[(snvs_maf_obs["canonical_SYMBOL"] == gene) - & (snvs_maf_obs["TYPE"] == "SNV")] - - # Subset gene - seq_info_df_gene = seq_info_df[seq_info_df["Gene"] == gene] - gene_len = len(seq_info_df_gene.Seq.values[0]) - pos_gene = pd.DataFrame({"Pos": range(1, gene_len + 1)}) - - # Get per-position SNV mutations count - obs_snv_count_gene = snvs_maf_obs_gene.groupby("canonical_Protein_position").size().reset_index(name='Count').astype(int) - obs_snv_count_gene.columns = ["Pos", "Count"] - obs_snv_count_gene = pos_gene.merge(obs_snv_count_gene, how="left", on="Pos") - - - # Determine the maximum y-value for this gene - max_y_gene = obs_snv_count_gene["Count"].max() - - - # Create a single figure with two subplots - fig, axs = plt.subplots(1, 1, figsize=(10, 3), sharex=True) - - - # Plot the data for observed and randomized mutations in separate subplots - plotting_needle_from_counts(obs_snv_count_gene, axs, max_y=max_y_gene) - - # Set titles for each subplot - axs.set_title(f'{gene} {sample} - {int(obs_snv_count_gene["Count"].sum()):,} observed SNVs') - - # Adjust layout to prevent overlap - plt.tight_layout() - - return fig - - -def plotting_needle_from_counts(data_gene, - ax = None, - max_y=None, - col_pos_track='#003366', - label_pos_track='observed', - col_hv_lines='grey'): - - # Determine max_y if not provided - if max_y is None: - max_y = data_gene["Count"].max() - - # Calculate the precise margin to add to the y-axis - marker_size = 60 - marker_radius = np.sqrt(marker_size / np.pi) - marker_margin = marker_radius * 0.5 # Convert marker radius to a suitable margin - - # Add the margin to max_y - max_y += marker_margin - - - # If ax is not provided, create a new figure and axis - if ax is None: - fig, ax = plt.subplots(figsize=(10, 3)) - else: - fig = None # No need to create a figure if ax is provided - - - ax.vlines(data_gene["Pos"], ymin=0, ymax=data_gene["Count"], color=col_hv_lines, lw=1, zorder=1, alpha=0.5) - ax.scatter(data_gene["Pos"], data_gene["Count"], color='white', zorder=3, lw=1, ec="white") # To cover the overlapping needle top part - ax.scatter(data_gene["Pos"].values, data_gene["Count"].values, color=col_pos_track, zorder=4, - alpha=0.7, lw=0.1, ec="black", s=30, label=label_pos_track) - - # Remove the right and top spines - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - - # Set the y-axis limit - ax.set_ylim(0, max_y) - - # Add labels - ax.set_xlabel('Position') - ax.set_ylabel('Count') - - # If a new figure was created, show it - if fig is not None: - plt.show() - - return fig - - -def manager(mutations_file, o3d_seq_file, sample_name, sample_name_out): - # Load your MAF DataFrame (raw_annotated_maf) - maf = pd.read_csv(mutations_file, sep = "\t", header = 0, na_values = custom_na_values) - o3d_seq_df = pd.read_csv(o3d_seq_file, sep = "\t", header = 0) - - maf_f = maf[maf["canonical_Protein_position"] != '-'].reset_index(drop = True) - del maf - - gene_order = sorted(pd.unique(maf_f["canonical_SYMBOL"])) - - - os.makedirs(f"{sample_name_out}") - - # Loop over each gene to plot - for geneeee in gene_order: - print(geneeee) - fig_needles = plot_single_needle(geneeee, maf_f, - o3d_seq_df, - sample=sample_name - ) - - fig_needles.savefig(f"{sample_name_out}/{geneeee}.needles.pdf", bbox_inches='tight') - plt.close() - - - - -# @click.command() -# @click.option('--sample_name', type=str, help='Name of the sample being processed.') -# @click.option('--mut_file', type=click.Path(exists=True), help='Input mutation file') -# @click.option('--out_maf', type=click.Path(), help='Output MAF file') -# @click.option('--json_filters', type=click.Path(exists=True), help='Input mutation filtering criteria file') -# @click.option('--req_plots', type=click.Path(exists=True), help='Column names to output') -# # @click.option('--plot', is_flag=True, help='Generate plot and save as PDF') - -# def main(sample_name, mut_file, out_maf, json_filters, req_plots): # , plot): -# click.echo(f"Subsetting MAF file...") -# subset_mutation_dataframe(sample_name, mut_file, out_maf, json_filters, req_plots) - -# if __name__ == '__main__': -# main() - - -sample_name_ = sys.argv[1] -mut_file = sys.argv[2] -o3d_seq_file_ = sys.argv[3] -sample_name_out_ = sys.argv[4] -# out_maf = sys.argv[3] -# json_filters = sys.argv[4] -# req_plots = sys.argv[5] - - - -if __name__ == '__main__': - # maf = subset_mutation_dataframe(mut_file, json_filters) - # plot_manager(sample_name, maf, req_plots) - manager(mut_file, o3d_seq_file_, sample_name_, sample_name_out_) - - -# Init -run_name = "all_samples" - -deepcsa_run_dir = "/workspace/nobackup/bladder_ts/results/2024-06-20_deepCSA" -maf_file = os.path.join(deepcsa_run_dir, f"writemaf/{run_name}.filtered.tsv.gz") -o3d_datasets = "/workspace/nobackup/scratch/oncodrive3d/datasets_240506" -o3d_annotations = "/workspace/nobackup/scratch/oncodrive3d/annotations_240506" -fig3_data = "/workspace/projects/bladder_ts/notebooks/manuscript_figures/Fig3/data" - -gene_order = ["KMT2D","EP300","ARID1A","CREBBP","NOTCH2","KMT2C","STAG2","RB1", - "RBM10","KDM6A","TP53","FGFR3","CDKN1A","FOXQ1", - # "PIK3CA","TERT" - ] - - - -# Oncodrive3D - -o3d_seq_df = pd.read_table(f"{o3d_datasets}/seq_for_mut_prob.tsv") -o3d_annot_df = pd.read_csv(f"{o3d_annotations}/uniprot_feat.tsv", sep="\t") - -o3d_prob = f"{deepcsa_run_dir}/oncodrive3d/run/{run_name}/{run_name}.miss_prob.processed.json" -o3d_prob = json.load(open(o3d_prob, encoding="utf-8")) - -o3d_score = f"{deepcsa_run_dir}/oncodrive3d/run/{run_name}/{run_name}.3d_clustering_pos.csv" -o3d_score = pd.read_csv(o3d_score)[["Gene", "Pos", "Score", "Score_obs_sim", "pval", "C", "C_ext"]] - - -# MAF - -maf_df = pd.read_csv(maf_file, sep = "\t", na_values = custom_na_values) -maf_df_f, snv_df, trunc_df, synon_df, miss_df = preprocess_maf(maf_df) - - - - - -for gene in gene_order: - - plot_pars = init_plot_pars() - gene_len = len(o3d_seq_df[o3d_seq_df["Gene"] == gene].Seq.values[0]) - if gene_len > 400: - # if gene_len > 2500: - # gene_len = 2500 - ref_ratio = 500 / gene_len - fsize_x = round(gene_len / 26) - else: - ref_ratio = 1 - fsize_x = 15 - - plot_pars["fsize"] = fsize_x, 10 - plot_pars["ofml_cbar_coord"] = (0.04, 0.36, 1.1*ref_ratio, 0.85) - - - generate_plot(gene, - o3d_seq_df, - o3d_annot_df, - o3d_prob, - o3d_score, - maf_df_f, - snv_df, - trunc_df, - synon_df, - miss_df, - plot_pars, - output_path="Fig3_plots/Fig3b") - - - - - - - -# Plot -# ==== - -def init_plot_pars(): - - plot_pars = {"fsize" : (20,10), - "hspace" : 0.1, # General space between all tracks - "ofml_cbar_coord" : (0.01, 0.35, 0.9, 1), # Box to anchor coordinate for OncodriveFML color bar - "track_title_x_coord" : 0.83, # x-coordinate (respect to protein len) for track txt title - "score_txt_x_coord" : 1.13, # as track title but for track score txt - "track_title_fontsize" : 14, - "ylabel_fontsize" : 13.5, - "xlabel_fontsize" : 13.5, - "ylabel_pad" : 38, - "ticksize" : 10.5, - "legend_fontsize" : 12, - "legend_frameon" : True, - "dpi" : 300 - } - plot_pars["colors"] = {"ofml" : "viridis_r", - "omega_trunc" : "#FA5E32", - "omega_synon" : "#89E4A2", - "omega_miss" : "#FABE4A", - "o3d_score" : "#6DBDCC", - "o3d_cluster" : "skyblue", - "o3d_prob" : "darkgray", - "frameshift" : "#E4ACF4", - "inframe" : "C5", - "hv_lines" : "lightgray" # General horizontal and vertical lines (e.g., needle plot vline) - } - plot_pars["h_ratios"] = {"omega_trunc" : 0.3, - "omega_synon" : 0.3, - "omega_miss" : 0.3, - "space1" : 0.015, - "o3d" : 0.5, - "space2" : 0.015, - "ofml" : 0.5, - "space3" : 0.015, - "indels" : 0.5, - "space4" : 0.015, - "domain" : 0.07 - } - - return plot_pars - - -def plot_count_track(pos_track, - axes, - ax=0, - neg_track=None, - gene_len=None, - col_pos_track="gray", - col_neg_track="gray", - label_pos_track=None, - label_neg_track=None, - ymargin=None, - hv_lines="lightgray"): - - axes[ax].vlines(pos_track["Pos"], ymin=0, ymax=pos_track["Count"], color=hv_lines, lw=1, zorder=1, alpha=0.5) - axes[ax].scatter(pos_track["Pos"], pos_track["Count"], color='white', zorder=3, lw=1, ec="white") # To cover the overlapping needle top part - axes[ax].scatter(pos_track["Pos"].values, pos_track["Count"].values, color=col_pos_track, zorder=4, - alpha=0.7, lw=0.1, ec="black", s=60, label=label_pos_track) - - if isinstance(neg_track, pd.DataFrame): - if isinstance(gene_len, int): - axes[ax].hlines(0, xmin=0, xmax=gene_len, color="gray", lw=0.6, zorder=1) - axes[ax].vlines(neg_track["Pos"], ymin=-neg_track["Count"], ymax=0, color=hv_lines, lw=1, zorder=1, alpha=0.5) - axes[ax].scatter(neg_track["Pos"], -neg_track["Count"], color='white', zorder=3, lw=1, ec="white") # To cover the overlapping needle top part - axes[ax].scatter(pos_track["Pos"].values, -neg_track["Count"].values, color=col_neg_track, zorder=4, - alpha=0.7, lw=0.1, ec="black", s=60, label=label_neg_track) - axes[ax].set_yticklabels(abs(axes[ax].get_yticks())) - - if ymargin is not None: - axes[ax].set_ylim(-np.max(neg_track["Count"])-ymargin, np.max(pos_track["Count"])+ymargin) - - -def add_ax_text(score, pvalue, y_text, x_text, ax, y_shift=0.2, y_adjust=0, equal_less=True): - - ax.text(x_text, y_text+(y_text*y_shift)+y_adjust, - fr'$\mathit{{Score}}$ = {np.round(score, 2)}', ha='center', va='center', fontsize=13.5, color="black") - - equal = "≤" if equal_less else "=" - ax.text(x_text, y_text-(y_text*y_shift)+y_adjust, - fr'$\mathit{{p}}$-value {equal} {pvalue}', ha='center', va='center', fontsize=13.5, color="black") - - -def get_transcript_ids(gene, maf_df_f, o3d_seq_df): - - canonical_tr = maf_df_f[maf_df_f["SYMBOL"] == gene].canonical_Feature.unique()[0] - o3d_tr = o3d_seq_df[o3d_seq_df["Gene"] == gene].Ens_Transcr_ID.values[0] - - return canonical_tr, o3d_tr - - -def plot_gene_selection_signals(maf_trunc_count_gene, - maf_miss_count_gene, - maf_synon_count_gene, - frameshift_indels_count_gene, - inframe_indels_count_gene, - ofml_muts_score_gene, - o3d_score_gene, - o3d_prob_gene, - plot_pars, - domain_df=None, - add_track_title_text=True, - add_score_text=False, - rm_spines=False, - light_spines=False, - output_path=None): - - - gene_len = len(o3d_prob_gene) - fig, axes = plt.subplots(len(plot_pars["h_ratios"]), 1, - figsize=plot_pars["fsize"], - sharex=True, - gridspec_kw={'hspace': plot_pars["hspace"], - 'height_ratios': plot_pars["h_ratios"].values()}) - colors = plot_pars["colors"] - - - # Omega - # ===== - - # Omega trunc - ax=0 - plot_count_track(maf_trunc_count_gene, axes, ax=ax, gene_len=gene_len, col_pos_track=colors["omega_trunc"]) - if add_score_text: - add_ax_text(score=np.round(omega_truncating_score, 2), - pvalue=omega_truncating_pvalue, - y_text=np.max(maf_synon_count_gene["Count"])/2, - x_text=gene_len*plot_pars["score_txt_x_coord"], - ax=axes[ax], - y_shift=0.9, - y_adjust=2.4) - if add_track_title_text: - axes[ax].text(gene_len*plot_pars["track_title_x_coord"], 4.7, - fr'$\mathbf{{Omega}}$ $\mathbf{{Truncating}}$', - ha='center', va='center', fontsize=plot_pars["track_title_fontsize"], color="black") - axes[ax].set_ylabel('Truncating\ncount', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - - n_max = np.max(maf_trunc_count_gene["Count"]) - j = 12 - i = 24 - axes[ax].set_ylim(-n_max/i, n_max + n_max/j) - - axes[ax].yaxis.set_major_locator(MaxNLocator(integer=True, nbins=3)) - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - if rm_spines: - axes[ax].spines['top'].set_visible(False) - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['top'].set_color(colors["hv_lines"]) - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - # Omega synonym - ax=1 - plot_count_track(maf_synon_count_gene, axes, ax=ax, gene_len=gene_len, col_pos_track=colors["omega_synon"], hv_lines=colors["hv_lines"]) - axes[ax].set_ylabel('Synonymous\ncount', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - - n_max = np.max(maf_synon_count_gene["Count"]) - j = 10 - i = 85.714 - axes[ax].set_ylim(-n_max/i, n_max + n_max/j) - - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - axes[ax].yaxis.set_major_locator(MaxNLocator(integer=True, nbins=3)) - axes[ax].spines['top'].set_visible(False) - if rm_spines: - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - # Omega miss - ax=2 - plot_count_track(maf_miss_count_gene, axes, ax=ax, gene_len=gene_len, col_pos_track=colors["omega_miss"], hv_lines=colors["hv_lines"]) - y_text=np.max(maf_miss_count_gene["Count"])/2 - if add_score_text: - add_ax_text(score=np.round(omega_misss_score, 2), - pvalue=omega_misss_pvalue, - y_text=y_text, - x_text=gene_len*plot_pars["score_txt_x_coord"], - y_shift=0.3, - ax=axes[ax], - y_adjust=4) - if add_track_title_text: - axes[ax].text(gene_len*plot_pars["track_title_x_coord"], y_text+y_text*0.6935, - fr'$\mathbf{{Omega}}$ $\mathbf{{Missense}}$', - ha='center', va='center', fontsize=plot_pars["track_title_fontsize"], color="black") - axes[ax].set_ylabel('Missense\ncount', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - - n_max = np.max(maf_miss_count_gene["Count"]) - j = 10.74 - i = 14.5 - axes[ax].set_ylim(-n_max/i, n_max + n_max/j) - - axes[ax].spines['top'].set_visible(False) - axes[ax].yaxis.set_major_locator(MaxNLocator(integer=True, nbins=3)) - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - if rm_spines: - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - - # Oncodrive3D - # =========== - - ax=4 - axes[ax].plot(range(1, gene_len+1), o3d_score_gene["O3D_score_norm"], zorder=2, color=colors["o3d_score"], lw=1, label="Clustering score") - axes[ax].plot(range(1, gene_len+1), -o3d_prob_gene, zorder=3, color=colors["o3d_prob"], lw=1, label="Missense mut probability") - axes[ax].fill_between(o3d_score_gene['Pos'], 0, o3d_score_gene["O3D_score_norm"], where=(o3d_score_gene['Cluster'] == 1), - color=colors["o3d_cluster"], alpha=0.3, label='Cluster', zorder=0, lw=2) - if add_score_text: - add_ax_text(score=np.round(o3d_score, 2), - pvalue=o3d_pvalue, - y_text=(np.max(o3d_score_gene["O3D_score_norm"]) + np.max(o3d_prob_gene))/2, - x_text=gene_len*plot_pars["score_txt_x_coord"], - ax=axes[ax], - y_shift=0.15, - y_adjust=-0.0075, - equal_less=True) - if add_track_title_text: - axes[ax].text(gene_len*plot_pars["track_title_x_coord"], 0.015, - fr'$\mathbf{{3D}}$-$\mathbf{{Clustering}}$', - ha='center', va='center', fontsize=plot_pars["track_title_fontsize"], color="black") - - axes[ax].set_ylabel('Score &\nprobability', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - axes[ax].yaxis.set_major_locator(MaxNLocator(integer=False, nbins=4)) - tick_labels = [np.round(abs(tick), 3) for tick in axes[ax].get_yticks()] - axes[ax].set_yticklabels(tick_labels) - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - - axes[ax].legend(loc="upper left", fontsize=plot_pars["legend_fontsize"], frameon=plot_pars["legend_frameon"]) - #axes[ax].legend(bbox_to_anchor=[0, 0.35, 1, 0]) - if rm_spines: - axes[ax].spines['top'].set_visible(False) - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['top'].set_color(colors["hv_lines"]) - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - - # OncodriveFML - # ============ - - ax=6 - axes[ax].vlines(ofml_muts_score_gene["Pos"].values, ymin=0, ymax=ofml_muts_score_gene["Count"], color=colors["hv_lines"], lw=1, zorder=1, alpha=0.5) - fml_scatter = axes[ax].scatter(ofml_muts_score_gene["Pos"].values, ofml_muts_score_gene["Count"].values, - c=ofml_muts_score_gene["CADD_score"].values, cmap=colors["ofml"], zorder=4, - alpha=0.7, lw=0.1, ec="black", s=60, label="SNV") - - y_text=np.round(np.max(ofml_muts_score_gene["Count"])/2) - if add_score_text: - add_ax_text(score=np.round(ofml_score, 2), - pvalue=ofml_pvalue, - y_text=y_text, - x_text=gene_len*plot_pars["score_txt_x_coord"], - ax=axes[ax], - y_adjust=3, - equal_less=True) - - if add_track_title_text: - axes[ax].text(gene_len*plot_pars["track_title_x_coord"], y_text+y_text*0.8, - fr'$\mathbf{{Functional}}$ $\mathbf{{Impact}}$', - ha='center', va='center', fontsize=plot_pars["track_title_fontsize"], color="black") - - inset_ax = inset_axes(axes[ax], width="10%", height="10%", loc='center left', - bbox_to_anchor=plot_pars["ofml_cbar_coord"], bbox_transform=axes[ax].transAxes, borderpad=0) - cbar = plt.colorbar(fml_scatter, cax=inset_ax, orientation='horizontal') - cbar.set_label('Impact score', fontsize=plot_pars["legend_fontsize"]) - cbar.ax.tick_params(labelsize=plot_pars["ticksize"]) - - n_max = np.max(ofml_muts_score_gene["Count"]) - j = 14.5 - i = 19.33 - axes[ax].set_ylim(-n_max/i, n_max + n_max/j) - - axes[ax].set_ylabel('SNV count', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - - if rm_spines: - axes[ax].spines['top'].set_visible(False) - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['top'].set_color(colors["hv_lines"]) - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - - # Frameshift enrichment - # ===================== - - ax=8 - plot_count_track(frameshift_indels_count_gene, axes, ax=ax, neg_track=inframe_indels_count_gene, gene_len=gene_len, - col_pos_track=colors["frameshift"], col_neg_track=colors["inframe"], label_pos_track="Frameshift", label_neg_track="Inframe", - ymargin=1, hv_lines=colors["hv_lines"]) - if add_score_text: - add_ax_text(score=np.round(indels_score, 2), - pvalue=np.round(indels_pvalue, 2), - y_text=(np.max(frameshift_indels_count_gene["Count"]) + np.max(inframe_indels_count_gene["Count"]))/2, - x_text=gene_len*plot_pars["score_txt_x_coord"], - ax=axes[ax], - y_shift=0.24, - y_adjust=0.4, - equal_less=False) - if add_track_title_text: - axes[ax].text(gene_len*plot_pars["track_title_x_coord"], 3.9, - fr'$\mathbf{{Frameshift}}$ $\mathbf{{enrichment}}$', - ha='center', va='center', fontsize=plot_pars["track_title_fontsize"], color="black") - axes[ax].legend(loc="upper left", fontsize=plot_pars["legend_fontsize"], frameon=plot_pars["legend_frameon"]) - axes[ax].set_ylabel('Indels count', fontsize=plot_pars["ylabel_fontsize"], rotation=45, labelpad=plot_pars["ylabel_pad"], va='center') - - axes[ax].yaxis.set_major_locator(MaxNLocator(integer=True, nbins=5)) - tick_labels = [abs(int(tick)) for tick in axes[ax].get_yticks()] - axes[ax].set_yticklabels(tick_labels) - - if rm_spines: - axes[ax].spines['top'].set_visible(False) - axes[ax].spines['right'].set_visible(False) - elif light_spines: - axes[ax].spines['top'].set_color(colors["hv_lines"]) - axes[ax].spines['right'].set_color(colors["hv_lines"]) - - - # Domain - # ====== - - if isinstance(domain_df, pd.DataFrame): - ax=10 - domain_color_dict = {} - - for n, name in enumerate(domain_df["Description"].unique()): - domain_color_dict[name] = f"C{n}" - - n = 0 - added_domain = [] - for i, row in domain_df.iterrows(): - if pd.Series([row["Description"], row["Begin"], row["End"]]).isnull().any(): - continue - - name = row["Description"] - start = int(row["Begin"]) - end = int(row["End"]) - axes[ax].fill_between(range(start, end+1), -0.45, 0.45, alpha=0.5, color=domain_color_dict[name]) - if name not in added_domain: - y = -0.04 - axes[ax].text(((start + end) / 2)+0.5, y, name, ha='center', va='center', fontsize=10, color="black") - added_domain.append(name) - axes[ax].set_yticks([]) - else: - axes[10].remove() - - - # Spaces - # ===== - axes[3].remove() - axes[5].remove() - axes[7].remove() - axes[9].remove() - - - # X axes - # ====== - axes[ax].tick_params(axis='y', labelsize=plot_pars["ticksize"]) - axes[ax].tick_params(axis='x', labelsize=plot_pars["ticksize"]) - axes[ax].set_xlabel('Protein position', fontsize=plot_pars["xlabel_fontsize"]) - - if output_path is not None: - print(f"Saving {output_path}") - plt.savefig(output_path, dpi=plot_pars["dpi"], bbox_inches='tight') - plt.show() - - -def generate_plot(gene, - o3d_seq_df, - o3d_annot_df, - o3d_prob, - o3d_score, - maf_df_f, - snv_df, - trunc_df, - synon_df, - miss_df, - plot_pars, - domain_df=None, - add_track_title_text=False, - add_score_text=False, - rm_spines=True, - light_spines=False, - output_path=None): - - # Prepare data - # ============ - - # O3D gene - o3d_score_df_gene, o3d_score_gene, o3d_prob_gene, o3d_top_pos_score_gene, uni_id, gene_pos = get_o3d_gene_data(gene, - o3d_seq_df, - o3d_prob, - o3d_score) - - # MAF gene - snv_count_gene = get_maf_pos_count(gene, snv_df, gene_pos) - trunc_count_gene = get_maf_pos_count(gene, trunc_df, gene_pos) - synon_count_gene = get_maf_pos_count(gene, synon_df, gene_pos) - miss_count_gene = get_maf_pos_count(gene, miss_df, gene_pos) - - # OFML gene - ofml_score_gene = get_ofml_score_gene(gene, fig3_data) - - # Indels - frameshift_indels_df, inframe_indels_df = get_frameshift_indels_maf(maf_df_f) - frameshift_indels_count_gene, inframe_indels_count_gene = get_frameshift_indels_gene(gene, - frameshift_indels_df, - inframe_indels_df, - gene_pos) - - # Domain - domain_gene = o3d_annot_df[(o3d_annot_df["Gene"] == gene)].reset_index(drop=True) - - # Transcripts and Uniprot ID - canonical_tr, o3d_tr = get_transcript_ids(gene, maf_df_f, o3d_seq_df) - print(f"> {gene} - {canonical_tr} - {o3d_tr} - {uni_id}") - - - # Plot - # ==== - - if output_path is not None: - output_path = f"{output_path}/{gene}.png" - plot_gene_selection_signals(trunc_count_gene, - miss_count_gene, - synon_count_gene, - frameshift_indels_count_gene, - inframe_indels_count_gene, - ofml_score_gene, - o3d_score_df_gene, - o3d_prob_gene, - plot_pars, - domain_gene, - add_track_title_text, - add_score_text, - rm_spines, - light_spines, - output_path) - - - - - -# Prepare data -# ============ - -def get_o3d_gene_data(gene, - seq_df, - prob_dict, - score_df): - - # Subset gene - seq_df_gene = seq_df[seq_df["Gene"] == gene] - gene_len = len(seq_df_gene.Seq.values[0]) - gene_pos = pd.DataFrame({"Pos" : range(1, gene_len+1)}) - uni_id, af_f = seq_df_gene[["Uniprot_ID", "F"]].values[0] - prob_gene = np.array(prob_dict[f"{uni_id}-F{af_f}"]) - score_gene_df = score_df[score_df["Gene"] == gene].reset_index(drop=True) - score_gene, top_pos_score_gene = score_df[["Score", "Score_obs_sim"]].values[0] - - ## O3D score vector - score_gene_df = gene_pos.merge(score_gene_df[["Pos", "Score_obs_sim", "C", "C_ext"]], how="left", on="Pos") - - # Don't include Extended clusters - score_gene_df["C"] = (score_gene_df["C"] == 1) & (score_gene_df["C_ext"] == 0) - score_gene_df["C"] = score_gene_df["C"].astype(int) - score_gene_df = score_gene_df.drop(columns=["C_ext"]) - - score_gene_df.columns = "Pos", "O3D_score", "Cluster" - score_gene_df["O3D_score"] = score_gene_df["O3D_score"].fillna(0) - score_gene_df["Cluster"] = score_gene_df["Cluster"].fillna(0) - score_gene_df["O3D_score_norm"] = score_gene_df["O3D_score"] / sum(score_gene_df["O3D_score"]) - - return score_gene_df, score_gene, prob_gene, top_pos_score_gene, uni_id, gene_pos - - -def preprocess_maf(maf_df): - - maf_df["CLEAN_SAMPLE_ID"] = maf_df["SAMPLE_ID"].apply(lambda x: "_".join(x.split("_")[1:3])) - - maf_df_f = maf_df.loc[(maf_df["VAF"] <= 0.35) & - (~maf_df["FILTER.not_in_exons"]) & - (~maf_df["FILTER.not_covered"]) & - (~maf_df["FILTER.no_pileup_support"]) & # avoid variants w/o VAF recomputed - (~maf_df["FILTER.n_rich"]) & - (~maf_df["FILTER.low_mappability"]) - ].reset_index(drop = True) - - # SNV - snvs_maf = maf_df_f[(maf_df_f["TYPE"] == 'SNV') & - (maf_df_f["canonical_SYMBOL"].isin(gene_order)) & - (maf_df_f["canonical_Protein_position"] != '-' ) - ].reset_index(drop = True) - snvs_maf["canonical_Protein_position"] = snvs_maf["canonical_Protein_position"].astype(int) - - # Omega - maf_df_trunc = snvs_maf[snvs_maf["canonical_Consequence_broader"].isin(['nonsense', "essential_splice"])] - maf_df_synon = snvs_maf[snvs_maf["canonical_Consequence_broader"] == 'synonymous'] - maf_df_miss = snvs_maf[snvs_maf["canonical_Consequence_broader"] == 'missense'] - - return maf_df_f, snvs_maf, maf_df_trunc, maf_df_synon, maf_df_miss - - -def get_maf_pos_count(gene, maf_df, gene_pos): - - # Subset gene and cols - cols = ["SYMBOL", - "TYPE", - "Feature", - "canonical_Feature", - "canonical_Consequence", - "canonical_Consequence_broader", - "canonical_Protein_position"] - - df_gene = maf_df.loc[maf_df["SYMBOL"] == gene, cols] - - # Get per-position mutations count - df_gene_count = df_gene.groupby("canonical_Protein_position").apply(lambda x: len(x)).reset_index().astype(int) - df_gene_count.columns = "Pos", "Count" - df_gene_count = gene_pos.merge(df_gene_count, how="left", on="Pos") - - return df_gene_count - - -def get_ofml_score_gene(gene, fig3_data_path): - - mutations_in_gene = pd.read_table(f"{fig3_data_path}/mutations_scored.{gene}.tsv", na_values = custom_na_values) - ofml_muts_score_gene = mutations_in_gene.groupby(by = "canonical_Protein_position").agg( { "MUT_ID" : 'count', - "CADDscore" : 'mean'}).reset_index() - ofml_muts_score_gene.columns = "Pos", "Count", "CADD_score" - - return ofml_muts_score_gene - - -def get_frameshift_indels_maf(maf_df_f): - - # Filter and somatic only - indel_maf_df = maf_df_f.loc[ - (maf_df_f["TYPE"].isin(["INSERTION", "DELETION"])) - ].reset_index(drop = True) - indel_maf_df["INDEL_LENGTH"] = (indel_maf_df["REF"].str.len() - indel_maf_df["ALT"].str.len()).abs() - indel_maf_df["INDEL_INFRAME"] = [ x % 3 == 0 for x in indel_maf_df["INDEL_LENGTH"] ] - indel_maf_df.loc[indel_maf_df["INDEL_LENGTH"] >= 15 , "INDEL_INFRAME"] = False - indel_maf_df.loc[indel_maf_df["canonical_Consequence_broader"] == 'nonsense' , "INDEL_INFRAME"] = False - indel_maf_df.loc[indel_maf_df["canonical_Consequence_broader"] == 'essential_splice' , "INDEL_INFRAME"] = False - - # Observed frameshift indels + inframes of length >= 5 AA - frameshift_indels = indel_maf_df[(~indel_maf_df["INDEL_INFRAME"]) & - (indel_maf_df["canonical_Protein_position"] != '-' )].reset_index(drop = True) - - # Inframe indels - inframe_indels = indel_maf_df[(indel_maf_df["INDEL_INFRAME"]) & - (indel_maf_df["canonical_Protein_position"] != '-' )].reset_index(drop = True) - - return frameshift_indels, inframe_indels - - -def get_frameshift_indels_gene(gene, frameshift_indels_df, inframe_indels_df, gene_pos): - - ## Frameshift indels - frameshift_indels_gene = frameshift_indels_df[frameshift_indels_df["canonical_SYMBOL"] == gene] - - # Count by pos considering first pos as protein pos - frameshift_indels_gene.canonical_Protein_position = frameshift_indels_gene.canonical_Protein_position.apply(lambda x: x.split("-")[0]) - frameshift_indels_gene.canonical_Protein_position = frameshift_indels_gene.canonical_Protein_position[ - frameshift_indels_gene.canonical_Protein_position.apply(lambda x: x.isdigit())] - frameshift_indels_count_gene = frameshift_indels_gene.groupby("canonical_Protein_position").apply(lambda x: len(x)).reset_index().astype(int) - frameshift_indels_count_gene.columns = "Pos", "Count" - frameshift_indels_count_gene = gene_pos.merge(frameshift_indels_count_gene, how="left", on="Pos") - - ## Inframe indels - inframe_indels_gene = inframe_indels_df[inframe_indels_df["canonical_SYMBOL"] == gene] - - # Count - inframe_indels_gene.canonical_Protein_position = inframe_indels_gene.canonical_Protein_position.apply(lambda x: x.split("-")[0]) - inframe_indels_gene.canonical_Protein_position = inframe_indels_gene.canonical_Protein_position[ - inframe_indels_gene.canonical_Protein_position.apply(lambda x: x.isdigit())] - inframe_indels_count_gene = inframe_indels_gene.groupby("canonical_Protein_position").apply(lambda x: len(x)).reset_index().astype(int) - inframe_indels_count_gene.columns = "Pos", "Count" - inframe_indels_count_gene = gene_pos.merge(inframe_indels_count_gene, how="left", on="Pos") - - return frameshift_indels_count_gene, inframe_indels_count_gene - -seq_regions_df = pd.read_table("/workspace/nobackup/scratch/oncodrive3d/datasets_240506/seq_for_mut_prob.tsv") -tb = tabix.open("/workspace/datasets/CADD/v1.6/hg38/whole_genome_SNVs.tsv.gz") - - -def get_scores_gene(gene_info): - chromosome = gene_info["Chr"].values[0] - reverse_strand = gene_info["Reverse_strand"].values[0] == 1 - exon_coords = eval(gene_info["Exons_coord"].values[0]) - if reverse_strand: - exon_coords = [(y, x) for x, y in exon_coords] - - - position_mutation = dict() - for start, stop in exon_coords: - for e in tb.query(str(chromosome), start, stop): - if int(e[1]) not in position_mutation: - position_mutation[int(e[1])] = dict() - position_mutation[int(e[1])][e[3]] = float(e[5]) - - start2end_coords = { x:y for x, y in enumerate(sorted(position_mutation.keys(), reverse = not reverse_strand))} - - return position_mutation, start2end_coords - - - - -def score_mutations_from_gene(snv_data): - scores_muts = [] - for ind, row in snv_data.iterrows(): - try: - scores_muts.append(scores_all_mutations_in_gene[row["POS"]].get(row["ALT"], 0)) - print("found") - except: - scores_muts.append(0) - print(row["Consequence"], "not found") - snv_data["CADDscore"] = scores_muts - return snv_data - - - -# for gene in ["TP53"]: -for gene in gene_order: - gene_info = seq_regions_df[seq_regions_df["Gene"] == gene] - scores_all_mutations_in_gene, scores_of_protein_position = get_scores_gene(gene_info) - - mutations_in_gene = snvs_maf[snvs_maf["canonical_SYMBOL"] == gene].copy() - scored_mutations_in_gene = score_mutations_from_gene(mutations_in_gene) - - scored_mutations_in_gene.to_csv(f"{data_dir}/mutations_scored.{gene}.tsv", sep = '\t', header= True, index = False) - - diff --git a/bin/check_contamination.py b/bin/check_contamination.py index a004360a..4fa836c7 100755 --- a/bin/check_contamination.py +++ b/bin/check_contamination.py @@ -6,435 +6,661 @@ import click -import pandas as pd import matplotlib.pyplot as plt +import pandas as pd +import logging import seaborn as sns +import operator from read_utils import custom_na_values - +from utils_filter import germline_mask, somatic_mask + +# Logging +logging.basicConfig( + format="%(asctime)s | %(levelname)s | %(name)s - %(message)s", + level=logging.INFO, + datefmt="%m/%d/%Y %I:%M:%S %p" +) + +LOG = logging.getLogger("check_contamination") + +# Constants +GERMLINE_LABEL = "Germline Samples" +SOMATIC_LABEL = "Somatic Samples" +CONTAMINATION_PROPORTION_THRESHOLD = 0.5 +# Deliberately more restrictive than the --somatic-vaf-boundary used for the between-samples +# analysis: this is the numerator of a QC metric, so only confidently low-VAF calls should count. +SNP_CONTAMINATION_VAF_THRESHOLD = 0.05 +VAF_COLUMNS = ["VAF", "vd_VAF", "VAF_AM"] +VARIANT_DETAIL_COLUMNS = [ + "SAMPLE_ID", + "MUT_ID", + "canonical_SYMBOL", + "ALT_DEPTH", + "DEPTH", + "VAF", + "canonical_Consequence_broader", + "FILTER", +] +OPS = { + ">": operator.gt, + "<": operator.lt, + ">=": operator.ge, + "<=": operator.le +} # Assuming somatic_variants and germline_variants are loaded as pandas DataFrames def compute_shared_variants(somatic_variants, germline_variants): - """ - # Example usage: - # shared_variants_matrix = compute_shared_variants(somatic_variants, germline_variants) + """Count mutations shared between each somatic sample and each germline sample. + + Parameters + ---------- + somatic_variants : pd.DataFrame + Variant table with at least ``SAMPLE_ID`` and ``MUT_ID`` columns, + providing the somatic side of the comparison. + germline_variants : pd.DataFrame + Variant table with at least ``SAMPLE_ID`` and ``MUT_ID`` columns, + providing the germline side of the comparison. + + Returns + ------- + pd.DataFrame + Integer matrix indexed by somatic ``SAMPLE_ID`` (rows) and germline + ``SAMPLE_ID`` (columns); each cell is the number of ``MUT_ID`` values + shared between that pair of samples. """ - unique_somatic_samples = sorted(somatic_variants['SAMPLE_ID'].unique()) - unique_germline_samples = sorted(germline_variants['SAMPLE_ID'].unique()) + unique_somatic_samples = sorted(somatic_variants["SAMPLE_ID"].unique()) + unique_germline_samples = sorted(germline_variants["SAMPLE_ID"].unique()) # Create a DataFrame to store counts (avoid .fillna downcasting warning) shared_counts = pd.DataFrame(0, index=unique_somatic_samples, columns=unique_germline_samples, dtype=int) # Iterate through each somatic sample for somatic_sample in unique_somatic_samples: - somatic_mutations = set(somatic_variants[somatic_variants['SAMPLE_ID'] == somatic_sample]['MUT_ID']) + somatic_mutations = set(somatic_variants[somatic_variants["SAMPLE_ID"] == somatic_sample]["MUT_ID"]) # Compare with germline mutations of all other samples for germline_sample in unique_germline_samples: - germline_mutations = set(germline_variants[germline_variants['SAMPLE_ID'] == germline_sample]['MUT_ID']) + germline_mutations = set(germline_variants[germline_variants["SAMPLE_ID"] == germline_sample]["MUT_ID"]) # Count shared mutations shared_counts.loc[somatic_sample, germline_sample] = len(somatic_mutations & germline_mutations) return shared_counts - - - -def contamination_detection_between_samples(maf_df, somatic_maf_df): - - # this is if we were to consider both unique and no-unique variants - vaf_threshold = 0.2 - germline_vars_all_samples = maf_df.loc[(maf_df["VAF"] > vaf_threshold) & (maf_df["vd_VAF"] > vaf_threshold) & (maf_df["VAF_AM"] > vaf_threshold), - ["SAMPLE_ID", "MUT_ID"]].drop_duplicates() - - print(germline_vars_all_samples["MUT_ID"].shape) - print(len(germline_vars_all_samples["MUT_ID"].unique())) - - - somatic_variants = somatic_maf_df[["SAMPLE_ID", "MUT_ID"]] - print(somatic_variants.shape) - - - all_variants = maf_df[["SAMPLE_ID", "MUT_ID"]] - print(all_variants.shape) - - - ## Somatic vs Germline - - shared_variants_somatic2germline_matrix = compute_shared_variants(somatic_variants, germline_vars_all_samples) - - plt.figure(figsize=(18, 15)) - - # Compute total number of germline mutations per sample - germline_counts = germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_variants_somatic2germline_matrix.columns) - - # Create custom column labels with germline mutation counts - col_labels = [f"(n={germline_counts[col]}) {col}" for col in shared_variants_somatic2germline_matrix.columns] - - # Build annotation DataFrame without using deprecated DataFrame.applymap - mask = shared_variants_somatic2germline_matrix > 30 - annot = shared_variants_somatic2germline_matrix.where(mask) - # convert selected values to nullable int then to string, replace missing with empty string - annot = annot.round(0).astype('Int64').astype(str).replace('', '').fillna('') - +def create_heatmap(variants_matrix: pd.DataFrame, + annot: pd.DataFrame, + col_labels: list | None, + xlabel: str, + ylabel: str, + title: str, + output_file: str, + annot_kws_color: str = "black", + size: tuple[int, int] = (18, 15)): + """Create heatmap for specified set of mutations.""" + plt.figure(figsize=size) sns.heatmap( - shared_variants_somatic2germline_matrix, + variants_matrix, annot=annot, fmt="", cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - xticklabels=col_labels, - yticklabels=shared_variants_somatic2germline_matrix.index, + cbar_kws={"label": "Shared Mutations"}, + xticklabels=col_labels if col_labels is not None else False, + yticklabels=variants_matrix.index, linewidths=0.5, - annot_kws={"color": "black", "fontsize": 10} + annot_kws={"color": annot_kws_color, "fontsize": 10}, ) - plt.xlabel("Germline Samples", fontsize=14) - plt.ylabel("Somatic Samples", fontsize=14) - plt.title("Somatic mutations that are germline in other samples", fontsize=16) - plt.savefig("somatic_vs_germline.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() - - + plt.xlabel(xlabel, fontsize=14) + plt.ylabel(ylabel, fontsize=14) + plt.title(title, fontsize=16) + plt.savefig(output_file, bbox_inches="tight", dpi=100) + plt.close() +def prepare_datasets(maf_df, somatic_maf_df, somatic_vaf_boundary): + """Prepare datasets for contamination analysis. + + Parameters + ---------- + maf_df : pd.DataFrame + Full mutation table for all samples (used to derive germline variants and the + all-variants set), with at least ``SAMPLE_ID``, ``MUT_ID``, ``VAF``, ``vd_VAF`` and + ``VAF_AM`` columns. + somatic_maf_df : pd.DataFrame + Filtered somatic mutation table, with at least ``SAMPLE_ID`` and ``MUT_ID`` columns. + somatic_vaf_boundary : float + VAF threshold passed to ``germline_mask`` to identify germline variants (a variant is + germline when all of ``VAF``/``vd_VAF``/``VAF_AM`` exceed it). + + Returns + ------- + tuple + Tuple containing: + - germline_vars_all_samples: DataFrame of germline variants across all samples. + - somatic_variants: DataFrame of somatic variants. + - all_variants: DataFrame of all variants. + """ + # Consider both unique and non-unique variants when collecting germline variants + germline_vars_all_samples = maf_df.loc[ + germline_mask(maf_df, somatic_vaf_boundary), ["SAMPLE_ID", "MUT_ID"] + ].drop_duplicates() + LOG.info(f"Total variants: {germline_vars_all_samples.shape}") + LOG.info(f"Unique germline variants: {len(germline_vars_all_samples['MUT_ID'].unique())}") + + somatic_variants = somatic_maf_df[["SAMPLE_ID", "MUT_ID"]] + LOG.info(f"Somatic variants: {somatic_variants.shape}") + all_variants = maf_df[["SAMPLE_ID", "MUT_ID"]] + LOG.info(f"All variants: {all_variants.shape}") + + return germline_vars_all_samples, somatic_variants, all_variants + + +def two_way_comparison(df_a: pd.DataFrame, df_b: pd.DataFrame, annotation_threshold: float, operator_str: str, normalize: bool) -> tuple[pd.DataFrame, pd.DataFrame, list[str]]: + """Detect cross-sample contamination, compare one set of mutations with another set of mutations. + + Parameters + ---------- + df_a : pd.DataFrame + First DataFrame of variants (e.g. somatic); its samples become the rows of the matrix. + df_b : pd.DataFrame + Second DataFrame of variants (e.g. germline); its samples become the columns of the + matrix and provide the per-sample counts used for the labels and the normalization. + annotation_threshold : float + Threshold for annotating shared variants. Ignored when ``normalize`` is True, where + cells are annotated when the proportion falls between 0.8 and 1. + operator_str : str + String representing the comparison operator (e.g., ">", "<", ">=", "<="). + normalize : bool + Whether to divide the shared counts by the number of variants of each df_b sample. + """ + shared_df = compute_shared_variants(df_a, df_b) - ## All vs Germline + # Compute total number of counts of mutations per sample -> the columns of shared_df are the samples of df_b + counts = ( + df_b["SAMPLE_ID"].value_counts().reindex(shared_df.columns) + ) - shared_all_vs_germline_variants_matrix = compute_shared_variants(all_variants, germline_vars_all_samples) + if normalize: + shared_df = shared_df.divide( + counts, axis=1 + ) - # Compute total number of germline mutations per sample - germline_counts = germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_all_vs_germline_variants_matrix.columns) + # Create custom column labels with mutation counts + col_labels = [f"(n={counts[col]}) {col}" for col in shared_df.columns] + # Build annotation DataFrame without using deprecated DataFrame.applymap + if normalize: + mask = (shared_df > 0.8) & ( + shared_df < 1) + else: + mask = OPS[operator_str](shared_df, annotation_threshold) + + annot = shared_df.where(mask) + if normalize: + annot = annot.round(2).astype("string").fillna("") + else: + # convert selected values to nullable int then to string, replace missing with empty string + annot = annot.round(0).astype("Int64").astype("string").replace("", "").fillna("") + + return shared_df, annot, col_labels + +def find_contaminated_pairs(proportion_matrix: pd.DataFrame, counts_matrix: pd.DataFrame, threshold: float) -> pd.DataFrame: + """Identify receiver samples carrying the germline variants of another sample. + + A sample is flagged as a receiver when its highest proportion of non-shared germline + variants coming from any single other sample exceeds ``threshold``. All sources tied at + that maximum are reported. + + Parameters + ---------- + proportion_matrix : pd.DataFrame + Matrix of the proportion of each source sample's non-shared germline variants that are + present as non-germline variants in each receiver sample (receivers as rows, sources as + columns). + counts_matrix : pd.DataFrame + Matrix of the corresponding raw shared variant counts, with the same shape and labels. + threshold : float + Minimum proportion above which a receiver sample is considered contaminated. + + Returns + ------- + pd.DataFrame + One row per contaminated receiver, with columns ``SAMPLE_ID``, + ``MAX_PROPORTION_GERMLINE_FROM_SOURCE`` and ``SOURCE_SAMPLEID_COUNTS``, the latter + holding the list of ``(shared_variant_count, source_sample_id)`` pairs tied at the + maximum proportion. + """ + max_prop_per_sample = proportion_matrix.max(axis="columns") - normalized_shared_all_vs_germline_variants_matrix = shared_all_vs_germline_variants_matrix.divide(germline_counts, axis=1) + receiver_source_pairs = [] + for sample, max_val in max_prop_per_sample[max_prop_per_sample > threshold].items(): + sample_vals = proportion_matrix.loc[sample, :] + sample_vals_count = counts_matrix.loc[sample, :] + source_sampleids = sample_vals[sample_vals == max_val].index.values + receiver_source_pairs.append( + ( + sample, + round(max_val, 3), + list(zip([sample_vals_count[x].item() for x in source_sampleids], source_sampleids)), + ) + ) + LOG.info( + f"{sample} has {max_val:.2f} proportion of the germline variants of {source_sampleids[0]} " + f"with a VAF not corresponding to germline variants " + f"(shared variants count: {sample_vals_count[source_sampleids[0]]})." + ) - # Count shared mutations between somatic and germline samples + return pd.DataFrame( + receiver_source_pairs, columns=["SAMPLE_ID", "MAX_PROPORTION_GERMLINE_FROM_SOURCE", "SOURCE_SAMPLEID_COUNTS"] + ) - plt.figure(figsize=(18, 15)) +def contaminated_pairs_to_long(contaminated_pairs: pd.DataFrame) -> pd.DataFrame: + """Expand the tied-source lists into one row per receiver/source pair. - # Create custom column labels with germline mutation counts - col_labels = [f"(n={germline_counts[col]}) {col}" for col in normalized_shared_all_vs_germline_variants_matrix.columns] + Parameters + ---------- + contaminated_pairs : pd.DataFrame + Output of ``find_contaminated_pairs``; must be non-empty. + Returns + ------- + pd.DataFrame + Long-format table with columns ``SAMPLE_ID``, ``MAX_PROPORTION_GERMLINE_FROM_SOURCE``, + ``SHARED_VARIANT_COUNT`` and ``SOURCE_SAMPLEID``. + """ + contaminated_pairs_long = contaminated_pairs.explode("SOURCE_SAMPLEID_COUNTS") + contaminated_pairs_long[["SHARED_VARIANT_COUNT", "SOURCE_SAMPLEID"]] = pd.DataFrame( + contaminated_pairs_long["SOURCE_SAMPLEID_COUNTS"].tolist(), index=contaminated_pairs_long.index + ) - # Annotation: show rounded values only when 0.8 < x < 1 - cond = (normalized_shared_all_vs_germline_variants_matrix > 0.8) & (normalized_shared_all_vs_germline_variants_matrix < 1) - annot = normalized_shared_all_vs_germline_variants_matrix.where(cond) - annot = annot.round(2).astype('string').fillna('') + return contaminated_pairs_long.drop(columns="SOURCE_SAMPLEID_COUNTS") + + +def export_germline_variants_in_receiver(maf_df: pd.DataFrame, + germline_vars_all_samples: pd.DataFrame, + receiver_sample: str, + source_sample: str): + """Write the source sample's germline variants alongside their VAF in the receiver sample. + + Parameters + ---------- + maf_df : pd.DataFrame + Full mutation table for all samples, with at least the ``VARIANT_DETAIL_COLUMNS``. + germline_vars_all_samples : pd.DataFrame + Germline variants across all samples, with at least ``SAMPLE_ID`` and ``MUT_ID`` columns. + receiver_sample : str + Sample suspected of being contaminated. + source_sample : str + Sample suspected of being the contamination source. + """ + variant_details = maf_df[VARIANT_DETAIL_COLUMNS] - sns.heatmap(normalized_shared_all_vs_germline_variants_matrix, - annot=annot, - fmt="", - cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - xticklabels=col_labels, yticklabels=normalized_shared_all_vs_germline_variants_matrix.index, - annot_kws={"color": "white", "fontsize": 10}, - linewidths=0.5) + receiver_variants = variant_details[variant_details["SAMPLE_ID"] == receiver_sample].drop("SAMPLE_ID", axis=1) - plt.xlabel("Germline Samples", fontsize = 14) - plt.ylabel("All mutations samples", fontsize = 14) - plt.title("All mutations that are germline in other samples", fontsize = 16) - plt.savefig("allmutations_vs_germline.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() + source_germline = germline_vars_all_samples[germline_vars_all_samples["SAMPLE_ID"] == source_sample] + source_variants = variant_details[ + (variant_details["SAMPLE_ID"] == source_sample) & (variant_details["MUT_ID"].isin(source_germline["MUT_ID"].values)) + ].drop("SAMPLE_ID", axis=1) + merged_samples = receiver_variants.merge( + source_variants, + on=["MUT_ID", "canonical_SYMBOL", "canonical_Consequence_broader"], + suffixes=("_dest", "_source"), + how="right", + ) + merged_samples.sort_values(by=["VAF_dest"], ascending=False).to_csv( + f"{source_sample}.germline_variants_in.{receiver_sample}.tsv", header=True, sep="\t", index=False + ) +def contamination_detection_between_samples(maf_df, somatic_maf_df, somatic_vaf_boundary): + """Detect cross-sample contamination by comparing somatic and germline mutations. + + Parameters + ---------- + maf_df : pd.DataFrame + Full mutation table for all samples (used to derive germline variants and the + all-variants set), with at least ``SAMPLE_ID``, ``MUT_ID``, ``VAF``, ``vd_VAF`` and + ``VAF_AM`` columns. + somatic_maf_df : pd.DataFrame + Filtered somatic mutation table, with at least ``SAMPLE_ID`` and ``MUT_ID`` columns. + somatic_vaf_boundary : float + VAF threshold passed to ``germline_mask`` to identify germline variants (a variant is + germline when all of ``VAF``/``vd_VAF``/``VAF_AM`` exceed it). + """ + # Prepare datasets + germline_vars_all_samples, somatic_variants, all_variants = prepare_datasets(maf_df, somatic_maf_df, somatic_vaf_boundary) + ## Somatic vs Germline + shared_variants_somatic2germline_matrix, somatic_vs_germline_annot, somatic_vs_germline_col_labels = two_way_comparison(somatic_variants, germline_vars_all_samples, annotation_threshold=30, operator_str=">", normalize=False) + create_heatmap( + variants_matrix=shared_variants_somatic2germline_matrix, + annot=somatic_vs_germline_annot, + col_labels=somatic_vs_germline_col_labels, + xlabel=GERMLINE_LABEL, + ylabel=SOMATIC_LABEL, + title="Somatic mutations that are germline in other samples", + output_file="somatic_vs_germline.pdf") ## Germline vs Germline + shared_germline_variants_matrix, germline_vs_germline_annot, germline_vs_germline_col_labels = two_way_comparison(germline_vars_all_samples, germline_vars_all_samples, annotation_threshold=0, operator_str="<", normalize=False) - shared_germline_variants_matrix = compute_shared_variants(germline_vars_all_samples, germline_vars_all_samples) - - plt.figure(figsize=(18, 15)) - - # Compute total number of germline mutations per sample - germline_counts = germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_germline_variants_matrix.columns) - - # Create custom column labels with germline mutation counts - col_labels = [f"(n={germline_counts[col]}) {col}" for col in shared_germline_variants_matrix.columns] - - - # Annotation: follow original logic (keep values where < 0, else blank) - mask = shared_germline_variants_matrix < 0 - annot = shared_germline_variants_matrix.where(mask) - annot = annot.astype('string').fillna('') - - sns.heatmap(shared_germline_variants_matrix, - annot=annot, - fmt="", - cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - xticklabels=col_labels, yticklabels=shared_germline_variants_matrix.index, - linewidths=0.5, - annot_kws={"fontsize": 8} - ) - - plt.xlabel("Germline Samples", fontsize = 14) - plt.ylabel("Germline Samples", fontsize = 14) - plt.title("Germline mutations that are germline in other samples", fontsize = 16) - plt.savefig("germline_vs_germline.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() - - - - - - # Compute total number of germline mutations per sample - germline_counts = germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_germline_variants_matrix.columns) - - normalized_share_germline_vs_germline_variants_matrix = shared_germline_variants_matrix.divide(germline_counts, axis=1) - - - plt.figure(figsize=(18, 15)) - - - # Create custom column labels with germline mutation counts - col_labels = [f"(n={germline_counts[col]}) {col}" for col in normalized_share_germline_vs_germline_variants_matrix.columns] - - - cond = (normalized_share_germline_vs_germline_variants_matrix > 0.8) & (normalized_share_germline_vs_germline_variants_matrix < 1) - annot = normalized_share_germline_vs_germline_variants_matrix.where(cond) - annot = annot.round(2).astype('string').fillna('') - - sns.heatmap(normalized_share_germline_vs_germline_variants_matrix, - annot=annot, - fmt="", - cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - xticklabels=col_labels, yticklabels=normalized_share_germline_vs_germline_variants_matrix.index, - annot_kws={"color": "white", "fontsize": 10}, - linewidths=0.5) - - plt.xlabel("Germline Samples", fontsize = 14) - plt.ylabel("Germline samples", fontsize = 14) - plt.savefig("normalized.germline_vs_germline.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() + create_heatmap( + variants_matrix=shared_germline_variants_matrix, + annot=germline_vs_germline_annot, + col_labels=germline_vs_germline_col_labels, + xlabel=GERMLINE_LABEL, + ylabel=GERMLINE_LABEL, + title="Germline mutations that are germline in other samples", + output_file="germline_vs_germline.pdf") + ## All vs Germline + normalized_shared_all_vs_germline_variants_matrix, all_vs_germline_annot, all_vs_germline_col_labels = two_way_comparison(all_variants, germline_vars_all_samples, annotation_threshold=0.3, operator_str=">", normalize=True) + + create_heatmap( + variants_matrix=normalized_shared_all_vs_germline_variants_matrix, + annot=all_vs_germline_annot, + col_labels=all_vs_germline_col_labels, + xlabel=GERMLINE_LABEL, + ylabel="All mutations samples", + title="All mutations that are germline in other samples", + output_file="allmutations_vs_germline.pdf", + annot_kws_color="white") + + # Germline mutations + normalized_germline_vs_germline_variants_matrix, norm_germline_vs_germline_annot, norm_germline_vs_germline_col_labels = two_way_comparison(germline_vars_all_samples, germline_vars_all_samples, annotation_threshold=0, operator_str="<", normalize=True) + create_heatmap( + variants_matrix=normalized_germline_vs_germline_variants_matrix, + annot=norm_germline_vs_germline_annot, + col_labels=norm_germline_vs_germline_col_labels, + xlabel=GERMLINE_LABEL, + ylabel=GERMLINE_LABEL, + title="Normalized germline mutations that are germline in other samples", + output_file="normalized.germline_vs_germline.pdf", + annot_kws_color="white") ## Somatic vs Remaining Germline + shared_all_vs_germline_variants_matrix = compute_shared_variants(all_variants, germline_vars_all_samples) shared_somatic_to_non_shared_germline = shared_all_vs_germline_variants_matrix - shared_germline_variants_matrix # Those cases where the number of mutations is smaller than 5 are set to 0 shared_somatic_to_non_shared_germline[shared_somatic_to_non_shared_germline < 5] = 0 # Compute total number of germline mutations per sample - germline_counts = germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_somatic_to_non_shared_germline.columns) - - - total_germline_available_per_sample = (germline_counts - shared_germline_variants_matrix) - - shared_somatic_to_non_shared_germline_proportion = (shared_somatic_to_non_shared_germline / total_germline_available_per_sample).fillna(0) - + germline_counts = ( + germline_vars_all_samples["SAMPLE_ID"].value_counts().reindex(shared_somatic_to_non_shared_germline.columns) + ) + total_germline_available_per_sample = germline_counts - shared_germline_variants_matrix - plt.figure(figsize=(22, 18)) + shared_somatic_to_non_shared_germline_proportion = ( + shared_somatic_to_non_shared_germline / total_germline_available_per_sample + ).fillna(0) cond = shared_somatic_to_non_shared_germline_proportion > 0.45 annot = shared_somatic_to_non_shared_germline_proportion.where(cond) - annot = annot.round(2).astype('string').fillna('') - - sns.heatmap(shared_somatic_to_non_shared_germline_proportion, - annot=annot, - fmt="", - cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - # xticklabels=col_labels, - yticklabels=shared_somatic_to_non_shared_germline_proportion.index, - annot_kws={"color": "black", "fontsize": 10}, - linewidths=0.5) - - plt.xlabel("Non-shared germline", fontsize = 14) - plt.ylabel("Somatic", fontsize = 14) - plt.title("Somatic mutations that are germline in other samples", fontsize = 16) - plt.savefig("contamination.somatic_vs_remaininggermline.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() - plt.close() - - - plt.figure(figsize=(22, 18)) + annot = annot.round(2).astype("string").fillna("") + create_heatmap( + variants_matrix=shared_somatic_to_non_shared_germline_proportion, + annot=annot, + col_labels=None, + xlabel="Non-shared germline", + ylabel="Somatic", + title="Somatic mutations that are germline in other samples", + output_file="contamination.somatic_vs_remaininggermline.pdf", + annot_kws_color="white", + size=(22, 18) + ) + cond = shared_somatic_to_non_shared_germline > 0 annot = shared_somatic_to_non_shared_germline.where(cond) # convert to nullable int then string, replace missing with empty string - annot = annot.round(0).astype('Int64').astype('string').replace('', '').fillna('') - - sns.heatmap(shared_somatic_to_non_shared_germline, - annot=annot, - fmt="", - cmap="Blues", - cbar_kws={'label': 'Shared Mutations'}, - # xticklabels=col_labels, - yticklabels=shared_somatic_to_non_shared_germline.index, - annot_kws={"color": "black", "fontsize": 10}, - linewidths=0.5) - - plt.xlabel("Non-shared germline", fontsize = 14) - plt.ylabel("Somatic", fontsize = 14) - plt.title("Somatic mutations that are germline in other samples (count)", fontsize = 16) - plt.savefig("contamination.somatic_vs_remaininggermline.numbers.pdf", bbox_inches = 'tight', dpi = 100) - plt.show() - plt.close() - + annot = annot.round(0).astype("Int64").astype("string").replace("", "").fillna("") - max_prop_per_sample = shared_somatic_to_non_shared_germline_proportion.max(axis = 'columns') + create_heatmap( + variants_matrix=shared_somatic_to_non_shared_germline, + annot=annot, + col_labels=None, + xlabel="Non-shared germline", + ylabel="Somatic", + title="Somatic mutations that are germline in other samples (count)", + output_file="contamination.somatic_vs_remaininggermline.numbers.pdf", + size=(22, 18) + ) ## Exploration of contaminated samples - receiver_source_pairs = [] - for sample, max_val in max_prop_per_sample[max_prop_per_sample>0.5].reset_index().values: - sample_vals = shared_somatic_to_non_shared_germline_proportion.loc[sample,:] - sample_vals_count = shared_somatic_to_non_shared_germline.loc[sample,:] + contaminated_pairs = find_contaminated_pairs( + shared_somatic_to_non_shared_germline_proportion, + shared_somatic_to_non_shared_germline, + CONTAMINATION_PROPORTION_THRESHOLD, + ) + contaminated_pairs.to_csv("contaminated_samples.detailed.tsv", header=True, sep="\t", index=False) - source_sampleids = sample_vals[sample_vals == max_val].index.values - source_sampleid = source_sampleids[0] - receiver_source_pairs.append((sample, round(max_val,3), - list(zip([sample_vals_count[x].item() for x in source_sampleids], source_sampleids)))) - - print(f'{sample} has {max_val:.2f} proportion of the germline variants of {source_sampleid} as with a VAF not corresponding to germline variants.') - print(f'Shared variants count: {sample_vals_count[source_sampleid]}') - print() - - - subseeeet = maf_df[["SAMPLE_ID", "MUT_ID", 'canonical_SYMBOL', "ALT_DEPTH", "DEPTH", "VAF", 'canonical_Consequence_broader', 'FILTER']] - p_dest = subseeeet[subseeeet["SAMPLE_ID"] == sample].drop("SAMPLE_ID", axis = 1) - - p_source_germ = germline_vars_all_samples[germline_vars_all_samples["SAMPLE_ID"] == source_sampleid] - p_source = subseeeet[(subseeeet["SAMPLE_ID"] == source_sampleid) - & (subseeeet["MUT_ID"].isin(p_source_germ["MUT_ID"].values)) - ].drop("SAMPLE_ID", axis = 1) - - merged_samples = p_dest.merge(p_source, - on = ["MUT_ID", 'canonical_SYMBOL', 'canonical_Consequence_broader'], - suffixes = ("_dest", "_source"), - how = 'right' - ) - - merged_samples.sort_values(by =["VAF_dest"], ascending=False - ).to_csv(f"{source_sampleid}.germline_variants_in.{sample}.tsv", - header = True, - sep = '\t', - index = False) - - # plt.figure(figsize=(8, 6)) - # plt.scatter(x = merged_samples["VAF_dest"].fillna(0), - # y = merged_samples["VAF_source"].fillna(0), - # # color = ['blue' if x == 0 else 'red' for x in merged_samples["VAF_dest"].fillna(0)] - # ) - - # plt.xscale('log') - # # plt.yscale('log') - # plt.xlabel("VAF_dest " + sample) - # plt.ylabel("VAF_source " + source_sampleid) - # plt.savefig(f"{source_sampleid}_germline_in_{sample}_VAF_scatter.pdf", bbox_inches = 'tight', dpi = 100) - # plt.show() - - # Store contamination results - contamination_detailed_df = pd.DataFrame(receiver_source_pairs, - columns=["SAMPLE_ID", "MAX_PROPORTION_GERMLINE_FROM_SOURCE", "SOURCE_SAMPLEID_COUNTS"]) - contamination_detailed_df.to_csv(f"contaminated_samples.detailed.tsv", - header = True, - sep = '\t', - index = False) - - if contamination_detailed_df.empty: - print("No contaminated samples detected.") + if contaminated_pairs.empty: + LOG.info("No contaminated samples detected.") return - contamination_detailed_df_long = contamination_detailed_df.explode("SOURCE_SAMPLEID_COUNTS") - expanded_df = pd.DataFrame(contamination_detailed_df_long["SOURCE_SAMPLEID_COUNTS"].tolist()) - expanded_df.columns = ["SHARED_VARIANT_COUNT", "SOURCE_SAMPLEID"] - contamination_detailed_df_long["SHARED_VARIANT_COUNT"] = expanded_df["SHARED_VARIANT_COUNT"].values - contamination_detailed_df_long["SOURCE_SAMPLEID"] = expanded_df["SOURCE_SAMPLEID"].values - contamination_detailed_df_long = contamination_detailed_df_long.drop("SOURCE_SAMPLEID_COUNTS", axis = 1) - contamination_detailed_df_long.to_csv(f"contaminated_samples.detailed.long.tsv", - header = True, - sep = '\t', - index = False) + for receiver in contaminated_pairs.itertuples(): + # Only the first of the sources tied at the maximum proportion is reported in detail + source_sampleid = receiver.SOURCE_SAMPLEID_COUNTS[0][1] + export_germline_variants_in_receiver(maf_df, germline_vars_all_samples, receiver.SAMPLE_ID, source_sampleid) + + contaminated_pairs_to_long(contaminated_pairs).to_csv( + "contaminated_samples.detailed.long.tsv", header=True, sep="\t", index=False + ) def data_loading(maf_path, somatic_maf_path): + """Load the full and somatic MAF tables, keeping only covered SNVs. + + Parameters + ---------- + maf_path : str + Path to the full MAF file; rows flagged ``FILTER.not_covered`` are dropped and only + ``TYPE == "SNV"`` rows are kept. + somatic_maf_path : str + Path to the filtered somatic MAF file; only ``TYPE == "SNV"`` rows are kept. + + Returns + ------- + tuple + ``(maf_df, somatic_maf_df)`` — the filtered full mutation table and the filtered + somatic mutation table, both as ``pd.DataFrame``. + """ maf_df = pd.read_table(maf_path, na_values=custom_na_values) print(maf_df.shape) - maf_df = maf_df[~(maf_df["FILTER.not_covered"]) - & (maf_df["TYPE"] == 'SNV') - ].reset_index() + maf_df = maf_df[~(maf_df["FILTER.not_covered"]) & (maf_df["TYPE"] == "SNV")].reset_index() print(maf_df.shape) somatic_maf_df = pd.read_table(somatic_maf_path, na_values=custom_na_values) print(somatic_maf_df.shape) - somatic_maf_df = somatic_maf_df[(somatic_maf_df["TYPE"] == 'SNV')] + somatic_maf_df = somatic_maf_df[(somatic_maf_df["TYPE"] == "SNV")] print(somatic_maf_df.shape) return maf_df, somatic_maf_df -def contamination_detection_in_snps(maf): - - snp_positions_maf = maf[maf["FILTER.gnomAD_SNP"]][ - ["SAMPLE_ID", "MUT_ID", "VAF"] - ].reset_index(drop = True) - - # being very restrictive in the VAF to count the occurrences of potentially contaminated mutations - somatic_snp_positions_maf = snp_positions_maf[snp_positions_maf["VAF"] < 0.05].reset_index(drop = True) - germline_snp_positions_maf = snp_positions_maf[snp_positions_maf["VAF"] >= 0.05].reset_index(drop = True) - - unique_SNP_positions = snp_positions_maf["MUT_ID"].unique() - number_unique_SNP_positions = len(unique_SNP_positions) - - sample_SNP_mutation_freq = [] - for sample in snp_positions_maf["SAMPLE_ID"].unique(): - germline_count = len(germline_snp_positions_maf[germline_snp_positions_maf["SAMPLE_ID"] == sample]) - somatic_count = len(somatic_snp_positions_maf[somatic_snp_positions_maf["SAMPLE_ID"] == sample]) - remaining_germline = number_unique_SNP_positions-germline_count - sample_SNP_mutation_freq.append([sample, - germline_count, - remaining_germline, - somatic_count, - somatic_count / remaining_germline if remaining_germline > 0 else 1 - ]) - sample_SNP_mutation_freq_df = pd.DataFrame(sample_SNP_mutation_freq) - sample_SNP_mutation_freq_df.columns = ["SAMPLE_ID", "germline_count", "remaining_germline", "somatic_count", "prop_somatic_SNPs"] - - # identify outliers in the "prop_somatic_SNPs" column - sample_SNP_mutation_freq_df = sample_SNP_mutation_freq_df.sort_values(by = "prop_somatic_SNPs", ascending = False) - sample_SNP_mutation_freq_df.to_csv("sample_SNP_mutation_freq.tsv", header = True, sep = '\t', index = False) - +def prepare_snp_datasets(maf: pd.DataFrame, vaf_threshold: float) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]: + """Split the known SNP positions into somatic-looking and germline-looking variants. + + Germline is defined as the complement of somatic rather than through ``germline_mask``, so + that the two sets partition the SNP rows. A variant whose VAF estimates disagree (some at or + below the threshold, some above) is therefore counted as germline: it is not a confident + somatic call, so it must not be left out of both sets and silently inflate the denominator of + ``prop_somatic_SNPs``. + + Rows with an undefined VAF cannot be classified either way and are dropped, since + ``somatic_mask`` is False for them and the germline complement would otherwise absorb them. + + Parameters + ---------- + maf : pd.DataFrame + Full mutation table with at least ``SAMPLE_ID``, ``MUT_ID``, the ``VAF_COLUMNS`` and the + boolean ``FILTER.gnomAD_SNP`` column. + vaf_threshold : float + Upper bound (inclusive) on all VAF estimates for a variant to be called somatic. + + Returns + ------- + tuple + Tuple containing: + - snp_positions_maf: DataFrame of all classifiable variants at known SNP positions. + - somatic_snp_positions_maf: DataFrame of the somatic ones. + - germline_snp_positions_maf: DataFrame of the remaining ones. + """ + snp_positions_maf = maf.loc[maf["FILTER.gnomAD_SNP"], ["SAMPLE_ID", "MUT_ID", *VAF_COLUMNS]] + + undefined_vaf = snp_positions_maf[VAF_COLUMNS].isna().any(axis="columns") + if undefined_vaf.any(): + LOG.warning(f"Discarding {undefined_vaf.sum()} SNP variants with an undefined {VAF_COLUMNS} VAF.") + snp_positions_maf = snp_positions_maf.loc[~undefined_vaf] + LOG.info(f"SNP variants: {snp_positions_maf.shape}") + + is_somatic = somatic_mask(snp_positions_maf, vaf_threshold) + somatic_snp_positions_maf = snp_positions_maf.loc[is_somatic] + germline_snp_positions_maf = snp_positions_maf.loc[~is_somatic] + LOG.info(f"Somatic SNP variants: {somatic_snp_positions_maf.shape}") + LOG.info(f"Germline SNP variants: {germline_snp_positions_maf.shape}") + + return snp_positions_maf, somatic_snp_positions_maf, germline_snp_positions_maf + + +def compute_snp_somatic_proportion(snp_positions_maf: pd.DataFrame, + somatic_snp_positions_maf: pd.DataFrame, + germline_snp_positions_maf: pd.DataFrame) -> pd.DataFrame: + """Compute, per sample, the proportion of non-germline SNP positions that look somatic. + + Parameters + ---------- + snp_positions_maf : pd.DataFrame + All classifiable variants at known SNP positions, as returned by ``prepare_snp_datasets``. + somatic_snp_positions_maf : pd.DataFrame + The somatic subset of ``snp_positions_maf``. + germline_snp_positions_maf : pd.DataFrame + The germline subset of ``snp_positions_maf``. + + Returns + ------- + pd.DataFrame + One row per sample with columns ``SAMPLE_ID``, ``germline_count``, ``remaining_germline``, + ``somatic_count`` and ``prop_somatic_SNPs``, sorted by the latter in descending order. + """ + samples = snp_positions_maf["SAMPLE_ID"].unique() + number_unique_snp_positions = snp_positions_maf["MUT_ID"].nunique() + + germline_count = germline_snp_positions_maf["SAMPLE_ID"].value_counts().reindex(samples, fill_value=0) + somatic_count = somatic_snp_positions_maf["SAMPLE_ID"].value_counts().reindex(samples, fill_value=0) + + # NOTE: the denominator counts every SNP position of the cohort that is not germline in this + # sample, including positions where the sample has no variant at all. Those can never end up in + # the numerator, so `prop_somatic_SNPs` is driven as much by how many SNP positions were called + # in a sample as by how many of them look somatic. It is comparable across samples of a cohort + # sequenced with the same panel, not an absolute contamination rate. + remaining_germline = number_unique_snp_positions - germline_count + + sample_snp_mutation_freq_df = pd.DataFrame( + { + "germline_count": germline_count, + "remaining_germline": remaining_germline, + "somatic_count": somatic_count, + "prop_somatic_SNPs": (somatic_count / remaining_germline).where(remaining_germline > 0, 1), + } + ).rename_axis("SAMPLE_ID").reset_index() + + return sample_snp_mutation_freq_df.sort_values(by="prop_somatic_SNPs", ascending=False) + + +def create_snp_proportion_plot(sample_snp_mutation_freq_df: pd.DataFrame, output_file: str): + """Plot the distribution of the per-sample proportion of somatic SNPs across the cohort. + + Parameters + ---------- + sample_snp_mutation_freq_df : pd.DataFrame + Per-sample table with a ``prop_somatic_SNPs`` column, as returned by + ``compute_snp_somatic_proportion``. + output_file : str + Path of the plot to write. + """ plt.figure(figsize=(6, 3)) - sns.violinplot(data=sample_SNP_mutation_freq_df, x="prop_somatic_SNPs", - fill= False, color="lightgray", inner=None) - sns.swarmplot(data=sample_SNP_mutation_freq_df, x="prop_somatic_SNPs", color="black", size=3) + sns.violinplot(data=sample_snp_mutation_freq_df, x="prop_somatic_SNPs", fill=False, color="lightgray", inner=None) + sns.swarmplot(data=sample_snp_mutation_freq_df, x="prop_somatic_SNPs", color="black", size=3) plt.title("Proportion of all SNPs across samples\ndetected as somatic") plt.xlabel("Proportion of somatic SNPs per sample") plt.ylabel("Density") - plt.savefig("sample_SNP_mutation_freq.pdf", dpi=300, bbox_inches="tight") + plt.savefig(output_file, dpi=300, bbox_inches="tight") plt.close() -@click.command() -@click.option('--maf_path', type=click.Path(exists=True), required=True, help='Path to the MAF file.') -@click.option('--somatic_maf', type=click.Path(exists=True), required=True, help='Path to the filtered somatic mutations file.') -def main(maf_path, somatic_maf): +def contamination_detection_in_snps(maf): + """Estimate per-sample contamination from the VAF distribution at known SNP positions. + + Restricts to gnomAD SNP positions, splits them into somatic-looking and germline-looking + sets by VAF, computes the per-sample proportion of SNP positions that look somatic, writes + the resulting table, and plots its distribution across samples. + + Parameters + ---------- + maf : pd.DataFrame + Full mutation table with at least ``SAMPLE_ID``, ``MUT_ID``, the ``VAF_COLUMNS`` and the + boolean ``FILTER.gnomAD_SNP`` column. """ - CLI entry point for assessing contamination between samples using germline and somatic mutations. + snp_positions_maf, somatic_snp_positions_maf, germline_snp_positions_maf = prepare_snp_datasets( + maf, SNP_CONTAMINATION_VAF_THRESHOLD + ) + + sample_snp_mutation_freq_df = compute_snp_somatic_proportion( + snp_positions_maf, somatic_snp_positions_maf, germline_snp_positions_maf + ) + sample_snp_mutation_freq_df.to_csv("sample_SNP_mutation_freq.tsv", header=True, sep="\t", index=False) + + create_snp_proportion_plot(sample_snp_mutation_freq_df, "sample_SNP_mutation_freq.pdf") + + +@click.command() +@click.option("--maf_path", type=click.Path(exists=True), required=True, help="Path to the MAF file.") +@click.option( + "--somatic_maf", type=click.Path(exists=True), required=True, help="Path to the filtered somatic mutations file." +) +@click.option( + "--somatic-vaf-boundary", + type=float, + default=0.3, + show_default=True, + help="VAF boundary for somatic variants; a variant with all of VAF/vd_VAF/VAF_AM above it is germline.", +) +def main(maf_path, somatic_maf, somatic_vaf_boundary): + """Assess cross-sample contamination using germline and somatic mutations. + + Loads the input tables and runs both the between-samples and the SNP-based contamination + analyses, writing their tables and plots to the current working directory. + + Parameters + ---------- + maf_path : str + Path to the full MAF file. + somatic_maf : str + Path to the filtered somatic mutations file. + somatic_vaf_boundary : float, optional + VAF boundary separating somatic from germline variants. Default is 0.3. """ - + maf_df, somatic_maf_df = data_loading(maf_path, somatic_maf) print("Running contamination analysis between samples") - contamination_detection_between_samples(maf_df, somatic_maf_df) + contamination_detection_between_samples(maf_df, somatic_maf_df, somatic_vaf_boundary) print("Running general contamination analysis") contamination_detection_in_snps(maf_df) - -if __name__ == '__main__': - +if __name__ == "__main__": main() diff --git a/bin/filter_cohort.py b/bin/filter_cohort.py index 755eee1d..c6dca9a9 100755 --- a/bin/filter_cohort.py +++ b/bin/filter_cohort.py @@ -41,13 +41,12 @@ """ import logging -from pathlib import Path import click import pandas as pd -from utils import add_filter from read_utils import custom_na_values -from utils_filter import expand_filter_column +from utils import add_filter +from utils_filter import expand_filter_column, germline_mask, somatic_mask # Logging logging.basicConfig( @@ -55,9 +54,10 @@ ) LOG = logging.getLogger("filter_cohort") -def flag_repetitive_variants(maf_df: pd.DataFrame, - repetitive_variant_threshold: int, - somatic_vaf_boundary: float) -> pd.DataFrame: + +def flag_repetitive_variants( + maf_df: pd.DataFrame, repetitive_variant_threshold: int, somatic_vaf_boundary: float +) -> pd.DataFrame: """ Flags filter column for repetitive variants from the MAF dataframe. A variant is considered repetitive if it appears in at least ``repetitive_variant_threshold`` samples. Additionally, variants that consistently appear at the same position in reads @@ -87,13 +87,17 @@ def flag_repetitive_variants(maf_df: pd.DataFrame, # Work with already filtered df + somatic only to explore potential artifacts # take only variant and sample info from the df - maf_df_f_somatic = maf_df.loc[maf_df["VAF"] <= somatic_vaf_boundary][["MUT_ID","SAMPLE_ID", "PMEAN", "PSTD"]].reset_index(drop = True) + maf_df_f_somatic = maf_df.loc[somatic_mask(maf_df, somatic_vaf_boundary)][ + ["MUT_ID", "SAMPLE_ID", "PMEAN", "PSTD"] + ].reset_index(drop=True) # Group by 'MUT_ID' and count occurrences maf_df_f_somatic_pivot = maf_df_f_somatic.groupby("MUT_ID").size().reset_index(name="count") # Store repetitive variants - repetitive_variants = maf_df_f_somatic_pivot[maf_df_f_somatic_pivot["count"] >= repetitive_variant_threshold]["MUT_ID"] + repetitive_variants = maf_df_f_somatic_pivot[maf_df_f_somatic_pivot["count"] >= repetitive_variant_threshold][ + "MUT_ID" + ] LOG.info("%s repetitive_variants", len(repetitive_variants)) # Flag repetitive variants in the original dataframe @@ -104,10 +108,10 @@ def flag_repetitive_variants(maf_df: pd.DataFrame, maf_df = maf_df.drop("repetitive_variant", axis=1) # Use the position in read information to filter repetitive variants with a fixed position (likely artifacts) - maf_df_f_somatic_pos_info = maf_df_f_somatic[~(maf_df_f_somatic["PMEAN"].isna()) & - (maf_df_f_somatic["PMEAN"] != -1) & - (maf_df_f_somatic["PSTD"] == 0)] - + maf_df_f_somatic_pos_info = maf_df_f_somatic[ + ~(maf_df_f_somatic["PMEAN"].isna()) & (maf_df_f_somatic["PMEAN"] != -1) & (maf_df_f_somatic["PSTD"] == 0) + ] + # Check if there are any repetitive variants with a fixed position if maf_df_f_somatic_pos_info.shape[0] == 0: LOG.info("No repetitive variants with fixed position found.") @@ -127,19 +131,20 @@ def flag_repetitive_variants(maf_df: pd.DataFrame, # Flag these variants in the maf dataframe maf_df["repetitive_mapping_variant"] = maf_df["MUT_ID"].isin(variants_with_rep_position) LOG.info("%s muts flagged as repetitive_mapping_variant", maf_df["repetitive_mapping_variant"].sum()) - - maf_df["FILTER"] = maf_df[["FILTER","repetitive_mapping_variant"]].apply(lambda x: add_filter(x["FILTER"], x["repetitive_mapping_variant"], "repetitive_mapping_variant"), - axis = 1 - ) - maf_df = maf_df.drop("repetitive_mapping_variant", axis = 1) + + maf_df["FILTER"] = maf_df[["FILTER", "repetitive_mapping_variant"]].apply( + lambda x: add_filter(x["FILTER"], x["repetitive_mapping_variant"], "repetitive_mapping_variant"), axis=1 + ) + maf_df = maf_df.drop("repetitive_mapping_variant", axis=1) return maf_df -def flag_cohort_n_rich(maf_df: pd.DataFrame, - n_rich_cohort_proportion: float, - somatic_vaf_boundary: float) -> pd.DataFrame: + +def flag_cohort_n_rich( + maf_df: pd.DataFrame, n_rich_cohort_proportion: float, somatic_vaf_boundary: float +) -> pd.DataFrame: """ - Flags FILTER column for cohort_n_rich variants from the MAF dataframe. + Flags FILTER column for cohort_n_rich variants from the MAF dataframe. Parameters ---------- @@ -161,62 +166,59 @@ def flag_cohort_n_rich(maf_df: pd.DataFrame, if max_samples < 2: LOG.warning("Not enough samples to identify cohort_n_rich mutations!") return maf_df - + number_of_samples = max(2, (max_samples * n_rich_cohort_proportion) // 1) LOG.info(f"Flagging mutations that are n_rich in at least: {number_of_samples} samples as cohort_n_rich") # Work with already filtered df to explore potential artifacts # take only variant and sample info from the df. - maf_df_f = maf_df[["MUT_ID", "SAMPLE_ID", "VAF_Ns", "FILTER"]].reset_index(drop = True) + maf_df_f = maf_df[["MUT_ID", "SAMPLE_ID", "VAF_Ns", "FILTER"]].reset_index(drop=True) # Aggregate n_rich variants n_rich_vars_df = ( maf_df_f[maf_df_f["FILTER"].str.contains("n_rich")] .groupby("MUT_ID") - .agg( - N_rich_frequency=('SAMPLE_ID', 'count'), - VAF_Ns_threshold=('VAF_Ns', 'min') - ) - ) - + .agg(N_rich_frequency=("SAMPLE_ID", "count"), VAF_Ns_threshold=("VAF_Ns", "min")) + ) + # Flag variants that are n_rich in at least number_of_samples samples -> cohort_n_rich - n_rich_vars = set(n_rich_vars_df[n_rich_vars_df['N_rich_frequency'] >= number_of_samples].index) + n_rich_vars = set(n_rich_vars_df[n_rich_vars_df["N_rich_frequency"] >= number_of_samples].index) maf_df["cohort_n_rich"] = maf_df["MUT_ID"].isin(n_rich_vars) LOG.info("%s muts flagged as cohort_n_rich", maf_df["cohort_n_rich"].sum()) - maf_df["FILTER"] = maf_df[["FILTER","cohort_n_rich"]].apply(lambda x: add_filter(x["FILTER"], x["cohort_n_rich"], "cohort_n_rich"), - axis = 1 - ) - + maf_df["FILTER"] = maf_df[["FILTER", "cohort_n_rich"]].apply( + lambda x: add_filter(x["FILTER"], x["cohort_n_rich"], "cohort_n_rich"), axis=1 + ) + # Flag variants that are n_rich in at least 1 sample -> cohort_n_rich_uni - n_rich_vars_uni = set(n_rich_vars_df[n_rich_vars_df['N_rich_frequency'] > 0].index) + n_rich_vars_uni = set(n_rich_vars_df[n_rich_vars_df["N_rich_frequency"] > 0].index) maf_df["cohort_n_rich_uni"] = maf_df["MUT_ID"].isin(n_rich_vars_uni) LOG.info("%s muts flagged as cohort_n_rich_uni", maf_df["cohort_n_rich_uni"].sum()) - maf_df["FILTER"] = maf_df[["FILTER","cohort_n_rich_uni"]].apply(lambda x: add_filter(x["FILTER"], x["cohort_n_rich_uni"], "cohort_n_rich_uni"), - axis = 1 - ) - + maf_df["FILTER"] = maf_df[["FILTER", "cohort_n_rich_uni"]].apply( + lambda x: add_filter(x["FILTER"], x["cohort_n_rich_uni"], "cohort_n_rich_uni"), axis=1 + ) + # Flag variants that exceed the VAF_Ns threshold -> cohort_n_rich_threshold - maf_df = maf_df.merge(n_rich_vars_df, on = 'MUT_ID', how = 'left') - maf_df['N_rich_frequency'] = maf_df['N_rich_frequency'].fillna(0) - maf_df['VAF_Ns_threshold'] = maf_df['VAF_Ns_threshold'].fillna(1.1) + maf_df = maf_df.merge(n_rich_vars_df, on="MUT_ID", how="left") + maf_df["N_rich_frequency"] = maf_df["N_rich_frequency"].fillna(0) + maf_df["VAF_Ns_threshold"] = maf_df["VAF_Ns_threshold"].fillna(1.1) - maf_df["cohort_n_rich_threshold"] = maf_df["VAF_Ns"] >= maf_df['VAF_Ns_threshold'] + maf_df["cohort_n_rich_threshold"] = maf_df["VAF_Ns"] >= maf_df["VAF_Ns_threshold"] LOG.info("%s muts flagged as cohort_n_rich_threshold", maf_df["cohort_n_rich_threshold"].sum()) - maf_df["FILTER"] = maf_df[["FILTER","cohort_n_rich_threshold"]].apply(lambda x: add_filter(x["FILTER"], x["cohort_n_rich_threshold"], "cohort_n_rich_threshold"), - axis = 1 - ) + maf_df["FILTER"] = maf_df[["FILTER", "cohort_n_rich_threshold"]].apply( + lambda x: add_filter(x["FILTER"], x["cohort_n_rich_threshold"], "cohort_n_rich_threshold"), axis=1 + ) # Drop temporary columns - maf_df = maf_df.drop(["cohort_n_rich", "cohort_n_rich_uni", "cohort_n_rich_threshold"], axis = 1) - + maf_df = maf_df.drop(["cohort_n_rich", "cohort_n_rich_uni", "cohort_n_rich_threshold"], axis=1) + return maf_df -def flag_other_samples_snp(maf_df, - somatic_vaf_boundary: float) -> pd.DataFrame: + +def flag_other_samples_snp(maf_df, somatic_vaf_boundary: float) -> pd.DataFrame: """ Filters out SNPs from other samples from the MAF dataframe @@ -234,28 +236,27 @@ def flag_other_samples_snp(maf_df, """ LOG.info("Flagging SNPs from other samples...") # Get all germline variants from all samples, consider both unique and non-unique variants - germline_vars_all_samples = maf_df.loc[(maf_df["VAF"] > somatic_vaf_boundary) & - (maf_df["VAF_AM"] > somatic_vaf_boundary) & - (maf_df["vd_VAF"] > somatic_vaf_boundary), - "MUT_ID"].unique() - + germline_vars_all_samples = maf_df.loc[germline_mask(maf_df, somatic_vaf_boundary), "MUT_ID"].unique() + LOG.info(f"Using all germline variants of all samples, total: {len(germline_vars_all_samples)} variants.") # Identify variants that are germline in other samples but somatic in the current sample maf_df["other_sample_SNP"] = False - maf_df.loc[(maf_df["MUT_ID"].isin(germline_vars_all_samples)) & - (maf_df["VAF"] <= somatic_vaf_boundary), "other_sample_SNP"] = True - LOG.info("%s muts flagged as other_sample_SNP", maf_df['other_sample_SNP'].sum()) + maf_df.loc[ + (maf_df["MUT_ID"].isin(germline_vars_all_samples)) & somatic_mask(maf_df, somatic_vaf_boundary), + "other_sample_SNP", + ] = True + LOG.info("%s muts flagged as other_sample_SNP", maf_df["other_sample_SNP"].sum()) # Flag variants that are germline in other samples but somatic in the current sample - maf_df["FILTER"] = maf_df[["FILTER","other_sample_SNP"]].apply( - lambda x: add_filter(x["FILTER"], x["other_sample_SNP"], "other_sample_SNP"), - axis = 1 - ) - maf_df = maf_df.drop("other_sample_SNP", axis = 1) + maf_df["FILTER"] = maf_df[["FILTER", "other_sample_SNP"]].apply( + lambda x: add_filter(x["FILTER"], x["other_sample_SNP"], "other_sample_SNP"), axis=1 + ) + maf_df = maf_df.drop("other_sample_SNP", axis=1) return maf_df + def flag_gnomad_snp(maf_df: pd.DataFrame) -> pd.DataFrame: """ Flags gnomAD SNPs in the MAF dataframe @@ -274,19 +275,20 @@ def flag_gnomad_snp(maf_df: pd.DataFrame) -> pd.DataFrame: # Flag gnomAD SNPs if "gnomAD_SNP" in maf_df.columns: - maf_df["gnomAD_SNP"] = maf_df["gnomAD_SNP"].replace({"True": True, "False": False, '-' : False}).fillna(False).astype(bool) + maf_df["gnomAD_SNP"] = ( + maf_df["gnomAD_SNP"].replace({"True": True, "False": False, "-": False}).fillna(False).astype(bool) + ) LOG.info("Out of %d positions, %d are gnomAD SNPs", maf_df["gnomAD_SNP"].shape[0], maf_df["gnomAD_SNP"].sum()) - - maf_df["FILTER"] = maf_df[["FILTER","gnomAD_SNP"]].apply( - lambda x: add_filter(x["FILTER"], x["gnomAD_SNP"], "gnomAD_SNP"), - axis = 1 - ) - maf_df = maf_df.drop("gnomAD_SNP", axis = 1) + + maf_df["FILTER"] = maf_df[["FILTER", "gnomAD_SNP"]].apply( + lambda x: add_filter(x["FILTER"], x["gnomAD_SNP"], "gnomAD_SNP"), axis=1 + ) + maf_df = maf_df.drop("gnomAD_SNP", axis=1) return maf_df -def flag_vaf_ns_threshold(maf_df: pd.DataFrame, vaf_ns_threshold: float) -> pd.DataFrame: +def flag_vaf_ns_threshold(maf_df: pd.DataFrame, vaf_ns_threshold: float) -> pd.DataFrame: """ Flag variants that have a proportion of Ns higher than vaf_ns_threshold @@ -307,14 +309,14 @@ def flag_vaf_ns_threshold(maf_df: pd.DataFrame, vaf_ns_threshold: float) -> pd.D maf_df["high_n_vaf"] = maf_df[["VAF_Ns", "VAF_Ns_AM"]].ge(vaf_ns_threshold).any(axis=1) LOG.info("%s muts flagged as high_n_vaf", maf_df["high_n_vaf"].sum()) - maf_df["FILTER"] = maf_df[["FILTER","high_n_vaf"]].apply( - lambda x: add_filter(x["FILTER"], x["high_n_vaf"], "high_n_vaf"), - axis = 1 - ) - maf_df = maf_df.drop("high_n_vaf", axis = 1) + maf_df["FILTER"] = maf_df[["FILTER", "high_n_vaf"]].apply( + lambda x: add_filter(x["FILTER"], x["high_n_vaf"], "high_n_vaf"), axis=1 + ) + maf_df = maf_df.drop("high_n_vaf", axis=1) return maf_df + def flag_distorted_expanded(maf_df: pd.DataFrame) -> pd.DataFrame: """ If there is a column named VAF_distorted_expanded_sq, add a filter flag for variants with distorted VAF distribution. @@ -331,19 +333,22 @@ def flag_distorted_expanded(maf_df: pd.DataFrame) -> pd.DataFrame: """ LOG.info("Flagging variants with distorted VAF distribution...") - if 'VAF_distorted_expanded_sq' in maf_df.columns: - maf_df["FILTER"] = maf_df[["FILTER","VAF_distorted_expanded_sq"]].apply( - lambda x: add_filter(x["FILTER"], x["VAF_distorted_expanded_sq"], "VAF_distorted_expanded_sq"), - axis = 1 - ) + if "VAF_distorted_expanded_sq" in maf_df.columns: + maf_df["FILTER"] = maf_df[["FILTER", "VAF_distorted_expanded_sq"]].apply( + lambda x: add_filter(x["FILTER"], x["VAF_distorted_expanded_sq"], "VAF_distorted_expanded_sq"), axis=1 + ) return maf_df -def flag_maf(maf_df: pd.DataFrame, sample_name: str, - repetitive_variant_threshold: int, - somatic_vaf_boundary: float, - n_rich_cohort_proportion: float, - vaf_ns_threshold: float) -> None: + +def flag_maf( + maf_df: pd.DataFrame, + sample_name: str, + repetitive_variant_threshold: int, + somatic_vaf_boundary: float, + n_rich_cohort_proportion: float, + vaf_ns_threshold: float, +) -> None: """ Script to process a MAF (Mutation Annotation Format) file. It filters out repetitive variants, cohort_n_rich variants, and SNPs from other samples. @@ -386,36 +391,44 @@ def flag_maf(maf_df: pd.DataFrame, sample_name: str, maf_df = expand_filter_column(maf_df) ## Save final filtered MAF - maf_df.to_csv(f"{sample_name}.cohort.filtered.tsv.gz", - sep = "\t", - header = True, - index = False) - + maf_df.to_csv(f"{sample_name}.cohort.filtered.tsv.gz", sep="\t", header=True, index=False) + LOG.info("Cohort flagging complete!") + @click.command() -@click.option('--maf-df-file', required=True, type=click.Path(exists=True), help='Input gzipped MAF file (TSV)') -@click.option('--sample-name', required=True, type=str, help='Sample name for output file') -@click.option('--repetitive-variant-threshold', required=True, type=int, help='Threshold for repetitive variants') -@click.option('--somatic-vaf-boundary', required=True, type=float, help='VAF boundary for somatic variants') -@click.option('--n-rich-cohort-proportion', required=True, type=float, help='Proportion for n-rich cohort filtering') -@click.option('--vaf-ns-threshold', required=False, type=float, default=0.1, help='VAF of Ns threshold for filtering variants') -def main(maf_df_file: str, sample_name: str, repetitive_variant_threshold: int, - somatic_vaf_boundary: float, n_rich_cohort_proportion: float, vaf_ns_threshold: float): +@click.option("--maf-df-file", required=True, type=click.Path(exists=True), help="Input gzipped MAF file (TSV)") +@click.option("--sample-name", required=True, type=str, help="Sample name for output file") +@click.option("--repetitive-variant-threshold", required=True, type=int, help="Threshold for repetitive variants") +@click.option("--somatic-vaf-boundary", required=True, type=float, help="VAF boundary for somatic variants") +@click.option("--n-rich-cohort-proportion", required=True, type=float, help="Proportion for n-rich cohort filtering") +@click.option( + "--vaf-ns-threshold", required=False, type=float, default=0.1, help="VAF of Ns threshold for filtering variants" +) +def main( + maf_df_file: str, + sample_name: str, + repetitive_variant_threshold: int, + somatic_vaf_boundary: float, + n_rich_cohort_proportion: float, + vaf_ns_threshold: float, +): """ CLI wrapper for flag_maf function. """ # Load MAF dataframe - maf_df = pd.read_csv(maf_df_file, compression='gzip', header=0, sep='\t', na_values=custom_na_values) + maf_df = pd.read_csv(maf_df_file, compression="gzip", header=0, sep="\t", na_values=custom_na_values) LOG.debug(f"{maf_df_file}") # Flag MAF file - flag_maf(maf_df, + flag_maf( + maf_df, sample_name, repetitive_variant_threshold, - somatic_vaf_boundary, + somatic_vaf_boundary, n_rich_cohort_proportion, - vaf_ns_threshold) - + vaf_ns_threshold, + ) + -if __name__ == '__main__': - main() \ No newline at end of file +if __name__ == "__main__": + main() diff --git a/bin/test/test_check_contamination.py b/bin/test/test_check_contamination.py new file mode 100644 index 00000000..ab0d365d --- /dev/null +++ b/bin/test/test_check_contamination.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +""" +Unit tests for check_contamination.py + +Tests the functionality of checking contamination in genomic data. +""" + +import sys +import unittest +from pathlib import Path +import pandas as pd +import numpy as np + +# Add the bin directory to the path to import the module +sys.path.insert(0, str(Path(__file__).parent.parent)) +from check_contamination import compute_shared_variants + +class TestComputeSharedVariants(unittest.TestCase): + """Test suite for compute_shared_variants function.""" + + def _somatic_df(self): + """Return a minimal somatic variants DataFrame.""" + return pd.DataFrame({ + 'SAMPLE_ID': ['s1', 's1', 's2'], + 'MUT_ID': ['mut1', 'mut2', 'mut3'], + }) + + def _germline_df(self): + """Return a minimal germline variants DataFrame.""" + return pd.DataFrame({ + 'SAMPLE_ID': ['s1', 's1', 's2'], + 'MUT_ID': ['mut1', 'mut4', 'mut5'], + }) + + # ---- basic behaviour ------------------------------------------------ + + def test_returns_dataframe(self): + result = compute_shared_variants(self._somatic_df(), self._germline_df()) + self.assertIsInstance(result, pd.DataFrame) + + def test_dtype_is_int(self): + result = compute_shared_variants(self._somatic_df(), self._germline_df()) + self.assertEqual(result.dtypes['s1'], np.int64) + + def test_shared_mutations_count(self): + """s1 shares mut1 with s1 → cell [s1, s1] == 1.""" + result = compute_shared_variants(self._somatic_df(), self._germline_df()) + self.assertEqual(result.loc['s1', 's1'], 1) + + def test_no_shared_mutations(self): + """s2 has mut3 which s1/s2 don't have → cell [s2, s1] == 0.""" + result = compute_shared_variants(self._somatic_df(), self._germline_df()) + self.assertEqual(result.loc['s2', 's1'], 0) + + def test_multiple_shared_mutations(self): + som = pd.DataFrame({'SAMPLE_ID': ['s1'], 'MUT_ID': ['mut1']}) + ger = pd.DataFrame({ + 'SAMPLE_ID': ['s1', 's1', 's1'], + 'MUT_ID': ['mut1', 'mut1', 'mut2'], # mut1 appears twice + }) + result = compute_shared_variants(som, ger) + # mut1 is in both sets → intersection size is 1 (set semantics) + self.assertEqual(result.loc['s1', 's1'], 1) + + def test_empty_somatic(self): + result = compute_shared_variants( + pd.DataFrame({'SAMPLE_ID': [], 'MUT_ID': []}), + self._germline_df(), + ) + self.assertEqual(len(result), 0) + + def test_empty_germline(self): + result = compute_shared_variants( + self._somatic_df(), + pd.DataFrame({'SAMPLE_ID': [], 'MUT_ID': []}), + ) + self.assertEqual(len(result.columns), 0) + + def test_both_empty(self): + result = compute_shared_variants( + pd.DataFrame({'SAMPLE_ID': [], 'MUT_ID': []}), + pd.DataFrame({'SAMPLE_ID': [], 'MUT_ID': []}), + ) + self.assertIsInstance(result, pd.DataFrame) + self.assertEqual(result.shape, (0, 0)) + + def test_samples_sorted(self): + """Sample IDs should appear in sorted order in index/columns.""" + som = pd.DataFrame({'SAMPLE_ID': ['z', 'a'], 'MUT_ID': ['m1', 'm2']}) + ger = pd.DataFrame({'SAMPLE_ID': ['y', 'b'], 'MUT_ID': ['m3', 'm4']}) + result = compute_shared_variants(som, ger) + self.assertEqual(list(result.index), ['a', 'z']) + self.assertEqual(list(result.columns), ['b', 'y']) diff --git a/bin/test/test_utils_filter.py b/bin/test/test_utils_filter.py new file mode 100644 index 00000000..cfb630b6 --- /dev/null +++ b/bin/test/test_utils_filter.py @@ -0,0 +1,554 @@ +#!/usr/bin/env python3 +""" +Unit tests for utils_filter.py. + +Covers: + - somatic_mask / germline_mask (each and combined) + - filter_maf (all criterion branches) + - load_filter_criteria (parsing, combining, prefix stripping) + - expand_filter_column (boolean column creation + required columns) + - extract_flagged_regions_bed (empty and non-empty BED output) +""" + +import os +import sys +import tempfile +import unittest +from pathlib import Path + +import pandas as pd + +# Add the bin directory to the path to import sibling modules +sys.path.insert(0, str(Path(__file__).parent.parent)) +from utils_filter import ( + expand_filter_column, + extract_flagged_regions_bed, + filter_maf, + germline_mask, + load_filter_criteria, + somatic_mask, +) + +THRESHOLD = 0.3 + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _make_df(rows: list[tuple[float, float, float]]) -> pd.DataFrame: + """Build a minimal MAF DataFrame with VAF, vd_VAF, and VAF_AM columns.""" + vafs, vd_vafs, vaf_ams = zip(*rows) + return pd.DataFrame({"VAF": list(vafs), "vd_VAF": list(vd_vafs), "VAF_AM": list(vaf_ams)}) + + +def _make_maf(rows: list[dict]) -> pd.DataFrame: + """Build a MAF DataFrame from a list of row dicts, preserving column order.""" + return pd.DataFrame(rows) + + +# --------------------------------------------------------------------------- +# somatic_mask +# --------------------------------------------------------------------------- + + +class TestSomaticMask(unittest.TestCase): + """Tests for somatic_mask(maf_df, threshold).""" + + def test_all_below_threshold_is_somatic(self): + """All three VAF columns strictly below threshold → somatic True.""" + df = _make_df([(0.1, 0.2, 0.05)]) + result = somatic_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [True]) + + def test_all_above_threshold_is_not_somatic(self): + """All three VAF columns strictly above threshold → somatic False.""" + df = _make_df([(0.5, 0.6, 0.4)]) + result = somatic_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False]) + + def test_boundary_equality_is_somatic(self): + """All three VAF columns exactly equal to threshold → somatic True (≤ is inclusive).""" + df = _make_df([(THRESHOLD, THRESHOLD, THRESHOLD)]) + result = somatic_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [True]) + + def test_asymmetric_is_not_somatic(self): + """Mixed columns (some ≤ threshold, some > threshold) → somatic False.""" + df = _make_df([(0.1, 0.1, 0.5)]) + result = somatic_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False]) + + def test_multiple_rows(self): + """Full four-case matrix in a single DataFrame.""" + df = _make_df( + [ + (0.1, 0.2, 0.05), # all-below → True + (0.5, 0.6, 0.4), # all-above → False + (THRESHOLD, THRESHOLD, THRESHOLD), # boundary → True + (0.1, 0.1, 0.5), # asymmetric → False + ] + ) + result = somatic_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [True, False, True, False]) + + +# --------------------------------------------------------------------------- +# germline_mask +# --------------------------------------------------------------------------- + + +class TestGermlineMask(unittest.TestCase): + """Tests for germline_mask(maf_df, threshold).""" + + def test_all_below_threshold_is_not_germline(self): + """All three VAF columns strictly below threshold → germline False.""" + df = _make_df([(0.1, 0.2, 0.05)]) + result = germline_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False]) + + def test_all_above_threshold_is_germline(self): + """All three VAF columns strictly above threshold → germline True.""" + df = _make_df([(0.5, 0.6, 0.4)]) + result = germline_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [True]) + + def test_boundary_equality_is_not_germline(self): + """All three VAF columns exactly equal to threshold → germline False (> is exclusive).""" + df = _make_df([(THRESHOLD, THRESHOLD, THRESHOLD)]) + result = germline_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False]) + + def test_asymmetric_is_not_germline(self): + """Mixed columns (some ≤ threshold, some > threshold) → germline False.""" + df = _make_df([(0.1, 0.1, 0.5)]) + result = germline_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False]) + + def test_multiple_rows(self): + """Full four-case matrix in a single DataFrame.""" + df = _make_df( + [ + (0.1, 0.2, 0.05), # all-below → False + (0.5, 0.6, 0.4), # all-above → True + (THRESHOLD, THRESHOLD, THRESHOLD), # boundary → False + (0.1, 0.1, 0.5), # asymmetric → False + ] + ) + result = germline_mask(df, THRESHOLD) + self.assertListEqual(result.tolist(), [False, True, False, False]) + + +# --------------------------------------------------------------------------- +# somatic + germline masks are never simultaneously True +# --------------------------------------------------------------------------- + + +class TestMasksAreNotComplements(unittest.TestCase): + """Verify that somatic and germline masks are never simultaneously True.""" + + def test_no_row_is_true_in_both_masks(self): + """No row should satisfy both somatic and germline conditions at once.""" + df = _make_df( + [ + (0.1, 0.2, 0.05), # all-below + (0.5, 0.6, 0.4), # all-above + (THRESHOLD, THRESHOLD, THRESHOLD), # boundary + (0.1, 0.1, 0.5), # asymmetric + ] + ) + somatic = somatic_mask(df, THRESHOLD) + germline = germline_mask(df, THRESHOLD) + both_true = (somatic & germline).tolist() + self.assertListEqual(both_true, [False, False, False, False]) + + def test_asymmetric_row_is_false_in_both_masks(self): + """An asymmetric row must be False in somatic AND False in germline (the 'neither' case).""" + df = _make_df([(0.1, 0.1, 0.5)]) + self.assertFalse(somatic_mask(df, THRESHOLD).iloc[0]) + self.assertFalse(germline_mask(df, THRESHOLD).iloc[0]) + + +# --------------------------------------------------------------------------- +# filter_maf +# --------------------------------------------------------------------------- + + +class TestFilterMaf(unittest.TestCase): + """Tests for filter_maf(maf_df, filter_criteria).""" + + def _base_maf(self) -> pd.DataFrame: + """Return a small MAF DataFrame exercising all criterion branches.""" + return _make_maf( + [ + {"MUT_ID": "M1", "VAF": 0.1, "DEPTH": 50, "FILTER": "PASS", "TYPE": "SNV", "FILTER.not_covered": False}, + { + "MUT_ID": "M2", + "VAF": 0.4, + "DEPTH": 30, + "FILTER": "n_rich;NM20", + "TYPE": "SNV", + "FILTER.not_covered": False, + }, + { + "MUT_ID": "M3", + "VAF": 0.2, + "DEPTH": 60, + "FILTER": "low_mappability", + "TYPE": "INDEL", + "FILTER.not_covered": True, + }, + { + "MUT_ID": "M4", + "VAF": 0.05, + "DEPTH": 80, + "FILTER": "PASS", + "TYPE": "SNV", + "FILTER.not_covered": False, + }, + ] + ) + + # --- numeric operator branch (len(operator)==2) --- + + def test_numeric_le_filters_correctly(self): + """('VAF', 'le 0.3') keeps only rows with VAF ≤ 0.3.""" + df = self._base_maf() + result = filter_maf(df, [("VAF", "le 0.3")]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M3", "M4"]) + + def test_numeric_ge_filters_correctly(self): + """('DEPTH', 'ge 50') keeps only rows with DEPTH ≥ 50.""" + df = self._base_maf() + result = filter_maf(df, [("DEPTH", "ge 50")]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M3", "M4"]) + + def test_numeric_lt_filters_correctly(self): + """('VAF', 'lt 0.2') keeps only rows with VAF < 0.2.""" + df = self._base_maf() + result = filter_maf(df, [("VAF", "lt 0.2")]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M4"]) + + def test_multiple_numeric_criteria_are_anded(self): + """Combining two numeric criteria narrows the result set.""" + df = self._base_maf() + result = filter_maf(df, [("VAF", "le 0.3"), ("DEPTH", "ge 50")]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M3", "M4"]) + + # --- notcontains / contains branch --- + + def test_notcontains_excludes_matching_filter_token(self): + """('FILTER', 'notcontains n_rich') removes rows whose FILTER cell contains 'n_rich'.""" + df = self._base_maf() + result = filter_maf(df, [("FILTER", "notcontains n_rich")]) + # M2 has 'n_rich' in its FILTER → removed + self.assertNotIn("M2", result["MUT_ID"].tolist()) + self.assertIn("M1", result["MUT_ID"].tolist()) + + def test_notcontains_respects_semicolon_split(self): + """notcontains splits on ';' so a token that is a prefix of another is handled correctly.""" + df = _make_maf( + [ + {"MUT_ID": "A", "FILTER": "NM20;PASS"}, + {"MUT_ID": "B", "FILTER": "NM200;PASS"}, # 'NM20' is NOT a token here + {"MUT_ID": "C", "FILTER": "PASS"}, + ] + ) + result = filter_maf(df, [("FILTER", "notcontains NM20")]) + self.assertNotIn("A", result["MUT_ID"].tolist()) + self.assertIn("B", result["MUT_ID"].tolist()) + self.assertIn("C", result["MUT_ID"].tolist()) + + def test_contains_keeps_only_matching_filter_token(self): + """('FILTER', 'contains n_rich') keeps only rows whose FILTER cell contains 'n_rich'.""" + df = self._base_maf() + result = filter_maf(df, [("FILTER", "contains n_rich")]) + self.assertListEqual(result["MUT_ID"].tolist(), ["M2"]) + + # --- boolean column branch --- + + def test_boolean_criterion_true_selects_matching_rows(self): + """('FILTER.not_covered', True) keeps only rows where FILTER.not_covered is True.""" + df = self._base_maf() + result = filter_maf(df, [("FILTER.not_covered", True)]) + self.assertListEqual(result["MUT_ID"].tolist(), ["M3"]) + + def test_boolean_criterion_false_selects_matching_rows(self): + """('FILTER.not_covered', False) keeps only rows where FILTER.not_covered is False.""" + df = self._base_maf() + result = filter_maf(df, [("FILTER.not_covered", False)]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M2", "M4"]) + + # --- plain-value (equality) branch --- + + def test_plain_value_equality_match(self): + """('TYPE', 'SNV') keeps only rows where TYPE == 'SNV'.""" + df = self._base_maf() + result = filter_maf(df, [("TYPE", "SNV")]) + self.assertListEqual(sorted(result["MUT_ID"].tolist()), ["M1", "M2", "M4"]) + + def test_plain_value_no_match_returns_empty(self): + """('TYPE', 'NONEXISTENT') returns an empty DataFrame.""" + df = self._base_maf() + result = filter_maf(df, [("TYPE", "NONEXISTENT")]) + self.assertEqual(len(result), 0) + + # --- no criteria leaves DataFrame unchanged --- + + def test_empty_criteria_returns_all_rows(self): + """An empty criteria list leaves the DataFrame unchanged.""" + df = self._base_maf() + result = filter_maf(df, []) + self.assertEqual(len(result), len(df)) + + +# --------------------------------------------------------------------------- +# load_filter_criteria +# --------------------------------------------------------------------------- + + +class TestLoadFilterCriteria(unittest.TestCase): + """Tests for load_filter_criteria(filters, somatic_filters).""" + + def test_extracts_notcontains_entries_from_filters(self): + """Items starting with 'notcontains ' in filters are returned with prefix stripped.""" + result = load_filter_criteria("notcontains n_rich,notcontains NM20", "") + self.assertListEqual(sorted(result), ["NM20", "n_rich"]) + + def test_extracts_notcontains_entries_from_somatic_filters(self): + """Items from somatic_filters starting with 'notcontains ' are included.""" + result = load_filter_criteria("", "notcontains low_mappability") + self.assertListEqual(result, ["low_mappability"]) + + def test_combines_both_lists(self): + """Entries from both arguments are merged before filtering.""" + result = load_filter_criteria("notcontains n_rich", "notcontains NM20") + self.assertListEqual(sorted(result), ["NM20", "n_rich"]) + + def test_non_notcontains_entries_are_excluded(self): + """Entries that do not start with 'notcontains ' are silently dropped.""" + result = load_filter_criteria("notcontains n_rich,VAF le 0.3,PASS", "") + self.assertListEqual(result, ["n_rich"]) + + def test_empty_strings_return_empty_list(self): + """Both arguments being empty strings yields an empty list.""" + result = load_filter_criteria("", "") + self.assertListEqual(result, []) + + def test_whitespace_is_trimmed_around_items(self): + """Leading/trailing whitespace around comma-separated items is stripped.""" + result = load_filter_criteria(" notcontains n_rich , notcontains NM20 ", "") + self.assertListEqual(sorted(result), ["NM20", "n_rich"]) + + def test_duplicate_entries_are_preserved(self): + """Duplicates across both arguments are preserved (no deduplication contract).""" + result = load_filter_criteria("notcontains n_rich", "notcontains n_rich") + self.assertEqual(result.count("n_rich"), 2) + + +# --------------------------------------------------------------------------- +# expand_filter_column +# --------------------------------------------------------------------------- + + +class TestExpandFilterColumn(unittest.TestCase): + """Tests for expand_filter_column(maf_df).""" + + def _make_filter_df(self, filter_values: list[str]) -> pd.DataFrame: + """Build a minimal MAF DataFrame with only a FILTER column.""" + return pd.DataFrame({"FILTER": filter_values}) + + def test_creates_boolean_column_for_each_token(self): + """Each unique ';'-delimited token gets its own FILTER. boolean column.""" + df = self._make_filter_df(["n_rich;NM20", "PASS", "NM20"]) + result = expand_filter_column(df) + self.assertIn("FILTER.n_rich", result.columns) + self.assertIn("FILTER.NM20", result.columns) + self.assertIn("FILTER.PASS", result.columns) + + def test_boolean_values_are_correct(self): + """True only where the token is present in that row's FILTER value.""" + df = self._make_filter_df(["n_rich;NM20", "PASS", "NM20"]) + result = expand_filter_column(df) + # Row 0: n_rich and NM20 present + self.assertTrue(result.loc[0, "FILTER.n_rich"]) + self.assertTrue(result.loc[0, "FILTER.NM20"]) + self.assertFalse(result.loc[0, "FILTER.PASS"]) + # Row 1: only PASS present + self.assertFalse(result.loc[1, "FILTER.n_rich"]) + self.assertTrue(result.loc[1, "FILTER.PASS"]) + # Row 2: only NM20 present + self.assertFalse(result.loc[2, "FILTER.n_rich"]) + self.assertTrue(result.loc[2, "FILTER.NM20"]) + + def test_required_columns_always_exist(self): + """FILTER.not_covered and FILTER.not_in_exons are always created even if absent in data.""" + df = self._make_filter_df(["PASS", "PASS"]) + result = expand_filter_column(df) + self.assertIn("FILTER.not_covered", result.columns) + self.assertIn("FILTER.not_in_exons", result.columns) + + def test_required_columns_are_false_when_token_absent(self): + """Required columns are all False when neither token appears in the data.""" + df = self._make_filter_df(["PASS", "n_rich"]) + result = expand_filter_column(df) + self.assertFalse(result["FILTER.not_covered"].any()) + self.assertFalse(result["FILTER.not_in_exons"].any()) + + def test_single_token_per_row(self): + """A FILTER column with no semicolons creates one boolean column per distinct value.""" + df = self._make_filter_df(["alpha", "beta", "alpha"]) + result = expand_filter_column(df) + self.assertIn("FILTER.alpha", result.columns) + self.assertIn("FILTER.beta", result.columns) + self.assertListEqual(result["FILTER.alpha"].tolist(), [True, False, True]) + self.assertListEqual(result["FILTER.beta"].tolist(), [False, True, False]) + + def test_all_required_columns_true_when_token_present(self): + """FILTER.not_covered is True exactly for the rows that contain 'not_covered'.""" + df = self._make_filter_df(["not_covered;n_rich", "PASS", "not_covered"]) + result = expand_filter_column(df) + self.assertListEqual(result["FILTER.not_covered"].tolist(), [True, False, True]) + + +# --------------------------------------------------------------------------- +# extract_flagged_regions_bed +# --------------------------------------------------------------------------- + + +class TestExtractFlaggedRegionsBed(unittest.TestCase): + """Tests for extract_flagged_regions_bed(maf_df, name, filters, specification).""" + + def setUp(self): + """Switch into a fresh temporary directory for each test; restore on teardown.""" + self._tmpdir = tempfile.mkdtemp() + self._orig_dir = os.getcwd() + os.chdir(self._tmpdir) + + def tearDown(self): + """Restore original working directory.""" + os.chdir(self._orig_dir) + + def _make_expanded_maf(self, rows: list[dict]) -> pd.DataFrame: + """Build a MAF with CHROM/POS columns then run expand_filter_column.""" + df = pd.DataFrame(rows) + return expand_filter_column(df) + + # --- empty case --- + + def test_empty_case_creates_empty_bed_file(self): + """No flagged rows → an empty .bed file is touched and the function returns None.""" + df = self._make_expanded_maf( + [ + {"CHROM": "chr1", "POS": 100, "FILTER": "PASS"}, + {"CHROM": "chr1", "POS": 200, "FILTER": "PASS"}, + ] + ) + result = extract_flagged_regions_bed(df, "sample1", ["n_rich"]) + self.assertIsNone(result) + bed_path = Path("sample1.flagged-pos.bed") + self.assertTrue(bed_path.exists()) + self.assertEqual(bed_path.stat().st_size, 0) + + def test_empty_case_with_specification_uses_correct_filename(self): + """specification parameter is included in the BED file name for the empty case.""" + df = self._make_expanded_maf([{"CHROM": "chr1", "POS": 100, "FILTER": "PASS"}]) + extract_flagged_regions_bed(df, "sample1", ["n_rich"], specification="cohort-") + self.assertTrue(Path("sample1.cohort-flagged-pos.bed").exists()) + + def test_empty_case_no_matching_filter_columns(self): + """When filter names have no corresponding FILTER.* columns, result is empty BED.""" + df = self._make_expanded_maf([{"CHROM": "chr1", "POS": 100, "FILTER": "PASS"}]) + result = extract_flagged_regions_bed(df, "sampleX", ["nonexistent_filter"]) + self.assertIsNone(result) + self.assertTrue(Path("sampleX.flagged-pos.bed").exists()) + + # --- non-empty case --- + + def test_nonempty_case_writes_bed_with_correct_columns(self): + """BED file has four tab-separated columns: CHROM, START, END, FILTERS.""" + df = self._make_expanded_maf( + [ + {"CHROM": "chr1", "POS": 500, "FILTER": "n_rich"}, + {"CHROM": "chr2", "POS": 1000, "FILTER": "PASS"}, + ] + ) + extract_flagged_regions_bed(df, "sample2", ["n_rich"]) + bed_path = Path("sample2.flagged-pos.bed") + self.assertTrue(bed_path.exists()) + bed = pd.read_csv(bed_path, sep="\t", header=None, names=["CHROM", "START", "END", "FILTERS"]) + self.assertEqual(len(bed), 1) + row = bed.iloc[0] + self.assertEqual(row["CHROM"], "chr1") + self.assertEqual(row["START"], 500) + self.assertEqual(row["END"], 500) + self.assertIn("FILTER.n_rich", row["FILTERS"]) + + def test_nonempty_case_multiple_filters_joined_with_comma(self): + """When a position has two active filter flags, FILTERS column is comma-joined.""" + df = self._make_expanded_maf( + [ + {"CHROM": "chr1", "POS": 300, "FILTER": "n_rich;NM20"}, + ] + ) + extract_flagged_regions_bed(df, "sample3", ["n_rich", "NM20"]) + bed = pd.read_csv( + Path("sample3.flagged-pos.bed"), sep="\t", header=None, names=["CHROM", "START", "END", "FILTERS"] + ) + self.assertEqual(len(bed), 1) + filters_value = bed.iloc[0]["FILTERS"] + # Both filter column names should appear, joined by comma + self.assertIn("FILTER.n_rich", filters_value) + self.assertIn("FILTER.NM20", filters_value) + self.assertIn(",", filters_value) + + def test_nonempty_case_multiple_rows_all_written(self): + """Multiple flagged positions each produce a row in the BED file.""" + df = self._make_expanded_maf( + [ + {"CHROM": "chr1", "POS": 100, "FILTER": "n_rich"}, + {"CHROM": "chr1", "POS": 200, "FILTER": "NM20"}, + {"CHROM": "chr2", "POS": 50, "FILTER": "PASS"}, + ] + ) + extract_flagged_regions_bed(df, "sample4", ["n_rich", "NM20"]) + bed = pd.read_csv( + Path("sample4.flagged-pos.bed"), sep="\t", header=None, names=["CHROM", "START", "END", "FILTERS"] + ) + self.assertEqual(len(bed), 2) + self.assertSetEqual(set(bed["START"].tolist()), {100, 200}) + + def test_nonempty_case_returns_none(self): + """Function has no explicit return in the non-empty path, so returns None.""" + df = self._make_expanded_maf([{"CHROM": "chr1", "POS": 100, "FILTER": "n_rich"}]) + result = extract_flagged_regions_bed(df, "sample5", ["n_rich"]) + self.assertIsNone(result) + + def test_nonempty_case_with_specification_uses_correct_filename(self): + """specification parameter is included in the BED file name for the non-empty case.""" + df = self._make_expanded_maf([{"CHROM": "chr1", "POS": 100, "FILTER": "n_rich"}]) + extract_flagged_regions_bed(df, "sample6", ["n_rich"], specification="cohort-") + self.assertTrue(Path("sample6.cohort-flagged-pos.bed").exists()) + + def test_nonempty_case_bed_is_sorted_by_chrom_and_pos(self): + """BED rows are sorted by CHROM then POS (ascending).""" + df = self._make_expanded_maf( + [ + {"CHROM": "chr2", "POS": 800, "FILTER": "n_rich"}, + {"CHROM": "chr1", "POS": 999, "FILTER": "n_rich"}, + {"CHROM": "chr1", "POS": 100, "FILTER": "n_rich"}, + ] + ) + extract_flagged_regions_bed(df, "sample7", ["n_rich"]) + bed = pd.read_csv( + Path("sample7.flagged-pos.bed"), sep="\t", header=None, names=["CHROM", "START", "END", "FILTERS"] + ) + self.assertEqual(len(bed), 3) + self.assertEqual(bed.iloc[0]["CHROM"], "chr1") + self.assertEqual(bed.iloc[0]["START"], 100) + self.assertEqual(bed.iloc[1]["START"], 999) + self.assertEqual(bed.iloc[2]["CHROM"], "chr2") + + +if __name__ == "__main__": + unittest.main() diff --git a/bin/utils_filter.py b/bin/utils_filter.py index bfd201be..3209ae36 100644 --- a/bin/utils_filter.py +++ b/bin/utils_filter.py @@ -1,15 +1,18 @@ #!/usr/bin/env python import logging -import pandas as pd from pathlib import Path + +import pandas as pd + """ Utility functions for extracting filters from a MAF DataFrame. """ LOG = logging.getLogger(__name__) + def filter_maf(maf_df, filter_criteria): - ''' + """ Filter a MAF dataframe with filtering information coming from a list of tuples. This can be either a dictionary transformed to list with the .items() method or by directly creating a list of tuples. [('VAF', 'le 0.3'), ('VAF_AM', 'le 0.3'), ('vd_VAF', 'le 0.3'), @@ -17,76 +20,127 @@ def filter_maf(maf_df, filter_criteria): ('FILTER', 'notcontains cohort_n_rich_uni'), ('FILTER', 'notcontains NM20'), ('FILTER', 'notcontains no_pileup_support'), ('FILTER', 'notcontains other_sample_SNP'), ('FILTER', 'notcontains low_mappability')] - ''' + """ # Define mappings for operators used in criteria operators = { - 'eq': lambda x, y: x == y, - 'ne': lambda x, y: x != y, - 'lt': lambda x, y: x < y, - 'le': lambda x, y: x <= y, - 'gt': lambda x, y: x > y, - 'ge': lambda x, y: x >= y, - 'not': lambda x, y: x != y, - 'notcontains': lambda x, y: x.apply(lambda z : y not in z.split(";")), # (~maf_df["FILTER"].str.contains("not_in_panel")) - 'contains': lambda x, y: x.apply(lambda z : y in z.split(";")) + "eq": lambda x, y: x == y, + "ne": lambda x, y: x != y, + "lt": lambda x, y: x < y, + "le": lambda x, y: x <= y, + "gt": lambda x, y: x > y, + "ge": lambda x, y: x >= y, + "not": lambda x, y: x != y, + "notcontains": lambda x, y: x.apply( + lambda z: y not in z.split(";") + ), # (~maf_df["FILTER"].str.contains("not_in_panel")) + "contains": lambda x, y: x.apply(lambda z: y in z.split(";")), } # Apply filters based on criteria from the JSON file for col, criterion in filter_criteria: - if isinstance(criterion, bool): pref_len = maf_df.shape[0] maf_df = maf_df[maf_df[col] == criterion] - print(f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations.") + print( + f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations." + ) - elif ' ' in criterion: + elif " " in criterion: operator, value = criterion.split(maxsplit=1) if len(operator) == 2 and operator in operators: # 'VAF' : 'le 0.35' pref_len = maf_df.shape[0] maf_df = maf_df[operators[operator](maf_df[col], float(value))] - print(f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations.") + print( + f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations." + ) elif operator in operators: # 'FILTER' : 'notcontains n_rich', pref_len = maf_df.shape[0] maf_df = maf_df[operators[operator](maf_df[col], value)] - print(f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations.") + print( + f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations." + ) else: print(f"We have no filtering criteria defined for {col}:{criterion} filter.") - else: # 'TYPE' : 'SNV' pref_len = maf_df.shape[0] maf_df = maf_df[maf_df[col] == criterion] - print(f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations.") + print( + f"Applying {col}:{criterion} filter implied going from {pref_len} mutations to {maf_df.shape[0]} mutations." + ) return maf_df + +def somatic_mask(maf_df: pd.DataFrame, threshold: float) -> pd.Series: + """ + Return a boolean mask identifying somatic variants. + + Parameters + ---------- + maf_df : pd.DataFrame + MAF dataframe containing at least the columns ``VAF``, ``vd_VAF``, and + ``VAF_AM``. + threshold : float + Upper bound (inclusive) on VAF for a variant to be called somatic. + + Returns + ------- + pd.Series + Boolean series with the same index as *maf_df*; ``True`` where the + variant is somatic. + """ + return (maf_df["VAF"] <= threshold) & (maf_df["vd_VAF"] <= threshold) & (maf_df["VAF_AM"] <= threshold) + + +def germline_mask(maf_df: pd.DataFrame, threshold: float) -> pd.Series: + """ + Return a boolean mask identifying germline variants. + + Parameters + ---------- + maf_df : pd.DataFrame + MAF dataframe containing at least the columns ``VAF``, ``vd_VAF``, and + ``VAF_AM``. + threshold : float + Lower bound (exclusive) on VAF for a variant to be called germline. + + Returns + ------- + pd.Series + Boolean series with the same index as *maf_df*; ``True`` where the + variant is germline. + """ + return (maf_df["VAF"] > threshold) & (maf_df["vd_VAF"] > threshold) & (maf_df["VAF_AM"] > threshold) + + def load_filter_criteria(filters: str, somatic_filters: str) -> list[str]: """ Parse filter criteria from comma-separated strings. - + Parameters ---------- filters : str Comma-separated list of filter criteria somatic_filters : str Comma-separated list of somatic filter criteria - + Returns ------- list[str] List of filter names to apply """ # Parse comma-separated strings into lists - filter_list = [f.strip() for f in filters.split(',') if f.strip()] - somatic_filter_list = [f.strip() for f in somatic_filters.split(',') if f.strip()] - + filter_list = [f.strip() for f in filters.split(",") if f.strip()] + somatic_filter_list = [f.strip() for f in somatic_filters.split(",") if f.strip()] + # Combine both lists all_filters = filter_list + somatic_filter_list @@ -95,21 +149,22 @@ def load_filter_criteria(filters: str, somatic_filters: str) -> list[str]: LOG.info(f"Loaded {len(result)} filter criteria: {result}") return result + def expand_filter_column(maf_df: pd.DataFrame) -> pd.DataFrame: """ Expands the FILTER column by creating new columns for each unique filter. Each new column indicates if the corresponding filter is present (True/False). """ # Split FILTER column once per row and convert to set for O(1) lookup - filter_sets = maf_df["FILTER"].str.split(";").apply(lambda x: set(x) if x != [''] else set()) - + filter_sets = maf_df["FILTER"].str.split(";").apply(lambda x: set(x) if x != [""] else set()) + # Get all unique filter values (excluding empty strings) all_filters = set( - filter_val - for filter_val in maf_df["FILTER"].str.split(";").explode().unique() - if filter_val and filter_val != '' + filter_val + for filter_val in maf_df["FILTER"].str.split(";").explode().unique() + if filter_val and filter_val != "" ) - + # Ensure "not_covered" and "not_in_exons" exist required_filters = {"not_covered", "not_in_exons"} all_filters.update(required_filters) @@ -120,7 +175,10 @@ def expand_filter_column(maf_df: pd.DataFrame) -> pd.DataFrame: return maf_df -def extract_flagged_regions_bed(maf_df: pd.DataFrame, name: str, FILTERS: list[str], specification: str = "") -> pd.DataFrame | None: + +def extract_flagged_regions_bed( + maf_df: pd.DataFrame, name: str, filters: list[str], specification: str = "" +) -> pd.DataFrame | None: """ Returns a BED file with the regions discarded, including the list of filters applied to each mutation. Creates a properly formatted BED file with 0-based coordinates and half-open intervals. @@ -131,7 +189,7 @@ def extract_flagged_regions_bed(maf_df: pd.DataFrame, name: str, FILTERS: list[s Input MAF dataframe with filter columns. POS column should contain 1-based coordinates. name : str Sample name to be used in the output BED file name. - FILTERS : list[str] + filters : list[str] List of filter criteria to check for in the MAF dataframe. specification : str, optional Additional string to include in the output BED file name (e.g., "cohort-"), by default "". @@ -143,7 +201,7 @@ def extract_flagged_regions_bed(maf_df: pd.DataFrame, name: str, FILTERS: list[s Output coordinates are 0-based with half-open intervals [start, end). """ # List of filter columns you want to check for - filter_columns = [f"FILTER.{f}" for f in FILTERS if f"FILTER.{f}" in maf_df.columns] + filter_columns = [f"FILTER.{f}" for f in filters if f"FILTER.{f}" in maf_df.columns] maf_df_filters = maf_df[maf_df[filter_columns].any(axis=1)] if filter_columns else pd.DataFrame() @@ -157,24 +215,20 @@ def extract_flagged_regions_bed(maf_df: pd.DataFrame, name: str, FILTERS: list[s bed_df = maf_df_filters[["CHROM", "POS"] + filter_columns] # Transform to long format - _bed_melt = (pd.melt(bed_df, - id_vars=["CHROM", "POS"], - value_vars=filter_columns, - var_name="FILTERS") - .query("value == True") - ) + _bed_melt = pd.melt(bed_df, id_vars=["CHROM", "POS"], value_vars=filter_columns, var_name="FILTERS").query( + "value == True" + ) LOG.info("Mutations flagged: %s", _bed_melt.shape[0]) # Aggregate filters per position bed_annotated = ( - _bed_melt - .drop_duplicates() - .sort_values(by=["CHROM", "POS"]) - .groupby(["CHROM","POS"])["FILTERS"] - .agg(','.join) - .reset_index() - .rename(columns={"POS": "START"}) + _bed_melt.drop_duplicates() + .sort_values(by=["CHROM", "POS"]) + .groupby(["CHROM", "POS"])["FILTERS"] + .agg(",".join) + .reset_index() + .rename(columns={"POS": "START"}) ) # The idea is to filter depth files at these positions, so make END = START (1-based) @@ -183,6 +237,8 @@ def extract_flagged_regions_bed(maf_df: pd.DataFrame, name: str, FILTERS: list[s LOG.info("Unique regions flagged: %s", bed_annotated.shape[0]) # Write BED file without header or index - (bed_annotated[["CHROM", "START", "END", "FILTERS"]] - .to_csv(f"{name}.{specification}flagged-pos.bed", sep="\t", header=False, index=False) - ) \ No newline at end of file + ( + bed_annotated[["CHROM", "START", "END", "FILTERS"]].to_csv( + f"{name}.{specification}flagged-pos.bed", sep="\t", header=False, index=False + ) + ) diff --git a/conf/general_files_IRB.config b/conf/general_files_IRB.config index fbeb2e69..7d973664 100644 --- a/conf/general_files_IRB.config +++ b/conf/general_files_IRB.config @@ -13,12 +13,12 @@ params { cadd_scores = "/data/bbg/datasets/CADD/v1.7/hg38/whole_genome_SNVs.tsv.gz" cadd_scores_ind = "/data/bbg/datasets/CADD/v1.7/hg38/whole_genome_SNVs.tsv.gz.tbi" - dnds_biomart_ref = "/data/bbg/datasets/pipelines/deepCSA/reference_datasets/homo_sapiens.v111.MANE.biomart.tsv" - dnds_covariates = "/data/bbg/projects/prominent/analysis/dNdScv/data/reference_files/covariates_hg19_hg38_epigenome_pcawg.rda" + dnds_biomart_ref = "/data/bbg/datasets/pipelines/deepCSA/reference_datasets/dNdScv_reference/homo_sapiens.v111.MANE.biomart.tsv" + dnds_covariates = "/data/bbg/datasets/pipelines/deepCSA/reference_datasets/dNdScv_covariates/covariates_hg19_hg38_epigenome_pcawg.rda" // oncodrive3d - datasets3d = "/data/bbg/datasets/oncodrive3d/datasets/datasets-260221" - annotations3d = "/data/bbg/datasets/oncodrive3d/annotations/annotations-260221" + datasets3d = "/data/bbg/datasets/oncodrive3d/datasets/datasets-260603" + annotations3d = "/data/bbg/datasets/oncodrive3d/annotations/annotations-260603" domains_file = "/data/bbg/projects/prominent/dev/internal_development/domains/bbgdomains.v2025.annotated_deepCSA.tsv" // Nanoseq masks diff --git a/conf/modules.config b/conf/modules.config index 9d3dbe0b..ea499691 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -357,6 +357,10 @@ process { ext.prop_samples_nrich = params.prop_samples_nrich } + withName: CONTAMINATION { + ext.germline_threshold = params.germline_threshold + } + withName: "TABLE2GROUP" { ext.unique_identifier = params.features_unique_identifier ext.feature_groups = params.features_groups_list diff --git a/modules/local/contamination/main.nf b/modules/local/contamination/main.nf index 7185b5d4..84dc1de0 100644 --- a/modules/local/contamination/main.nf +++ b/modules/local/contamination/main.nf @@ -16,10 +16,12 @@ process COMPUTE_CONTAMINATION { path "versions.yml" , topic: versions script: + def somatic_vaf_boundary = task.ext.germline_threshold ? "--somatic-vaf-boundary ${task.ext.germline_threshold}" : "" """ check_contamination.py \\ --maf_path ${maf} \\ - --somatic_maf ${somatic_maf} + --somatic_maf ${somatic_maf} \\ + ${somatic_vaf_boundary} cat <<-END_VERSIONS > versions.yml "${task.process}":