<jupyter_start><jupyter_text>Creating DataFrames for good and badly predicted drugs We are going to check for differences in drug pairs that are highly predictable or that are not.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy 
from scipy import stats<jupyter_output><empty_output><jupyter_text>We import our dataset and the DatFrame with the Pearson score of the drug pairs. drug_comb_commutative has the training set commutative but not the test set.<jupyter_code>path = "/content/drive/MyDrive/archivos_tfm/drug_comb_commutative.csv"

df_drug_comb = pd.read_csv(path, index_col = 0)
df_drug_comb.columns
drug_pair_pc = pd.read_csv("/content/drive/MyDrive/archivos_tfm/Drug_comb_Pearson_test_set_lb.csv", index_col = 0)
drug_pair_pc
good_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] >= 0.5]
good_drugs
bad_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] <= 0]
bad_drugs<jupyter_output><empty_output><jupyter_text>Now we want to create DataFrames of monotherapy information of the drugs that yield synergy to later check the p values of each feature and see if the null hypothesis is true.<jupyter_code>dataset_size = df_drug_comb.shape[0]
complete_list = np.arange(dataset_size, dtype = int)
mask_train = df_drug_comb["Dataset"]=="train"

test_idx = complete_list[~mask_train]

test_df = df_drug_comb.iloc[test_idx]

test_df
valid_xs_perm = pd.DataFrame()
valid_xs_perm[['Cell line name', 'Compound A', 
               'Compound B', 'Max. conc. A',
       'Max. conc. B', 'IC50 A', 
       'H A', 'Einf A', 
       'IC50 B', 'H B', 'Einf B',
       'Synergy score', 'QA', 
       'Combination ID', 'Challenge', 'Dataset']] = test_df[['Cell line name', 'Compound B', 
               'Compound A', 'Max. conc. B',
       'Max. conc. A', 'IC50 B', 
       'H B', 'Einf B', 
       'IC50 A', 'H A', 'Einf A','Synergy score', 'QA', 
       'Combination ID', 'Challenge', 'Dataset']]

      
commutative_df = pd.concat([test_df, valid_xs_perm])
commutative_df<jupyter_output><empty_output><jupyter_text>Now we have our commutative test set which has account for the fact that the drug pair AB is the same as BA. Now we only have to check for the p values of the features of drug A for example. mono_df is conmutative<jupyter_code>mono_df = commutative_df[['Max. conc. A','IC50 A', 'H A', 'Einf A','Combination ID']]
mono_df

good_drugs_mono = mono_df[mono_df["Combination ID"].isin(good_drugs["Combination ID"])]
good_drugs_mono
bad_drugs_mono = mono_df[mono_df["Combination ID"].isin(bad_drugs["Combination ID"])]
bad_drugs_mono<jupyter_output><empty_output><jupyter_text>Analysis of good and bad drugs<jupyter_code>fig1, ax1 = plt.subplots()
ax1.set_title('Basic Plot')
data = [good_drugs_mono["Max. conc. A"],bad_drugs_mono["Max. conc. A"]]
labels = ["Good","Bad"]
ax1.boxplot(data, labels =labels)
def boxplot_feature(feature, good_drugs_mono, bad_drugs_mono):
  fig1, ax1 = plt.subplots()
  ax1.set_title(feature)
  data = [good_drugs_mono[feature],bad_drugs_mono[feature]]
  labels = ["Good","Bad"]
  ax1.boxplot(data, labels =labels)

for feature in ['Max. conc. A', 'IC50 A', 'H A', 'Einf A']:
  boxplot_feature(feature, good_drugs_mono, bad_drugs_mono)
feature = 'Max. conc. A'
xA, xB = np.random.normal(1, 0.1, len(good_drugs_mono[feature])), np.random.normal(3, 0.1, len(bad_drugs_mono[feature]))

plt.scatter(xA, good_drugs_mono[feature])
plt.scatter(xB, bad_drugs_mono[feature])
plt.show()
def scatterplot_feature(feature, good_drugs_mono, bad_drugs_mono):
  xA, xB = np.random.normal(1, 0.1, len(good_drugs_mono[feature])), np.random.normal(3, 0.1, len(bad_drugs_mono[feature]))

  plt.scatter(xA, good_drugs_mono[feature])
  plt.scatter(xB, bad_drugs_mono[feature])
  plt.legend(["Good", "Bad"])
  plt.title(feature)
  plt.show()

for feature in ['Max. conc. A', 'IC50 A', 'H A', 'Einf A']:
  scatterplot_feature(feature, good_drugs_mono, bad_drugs_mono)<jupyter_output><empty_output><jupyter_text>Automating thresholds<jupyter_code>"up_thres is the threshold for highly predictable drugs"
"low_thres for the lowly predictable drugs" 
def good_bad_split(drug_pair_pc, up_thres, low_thres):
  good_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] >= up_thres]
  bad_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] <= low_thres]
  good_drugs_mono = mono_df[mono_df["Combination ID"].isin(good_drugs["Combination ID"])]
  bad_drugs_mono = mono_df[mono_df["Combination ID"].isin(bad_drugs["Combination ID"])]

  return good_drugs_mono, bad_drugs_mono

good_drugs_mono, bad_drugs_mono = good_bad_split(drug_pair_pc, up_thres = 0.9, low_thres = -0.8)


for feature in ['Max. conc. A', 'IC50 A', 'H A', 'Einf A']:
  
  boxplot_feature(feature, good_drugs_mono, bad_drugs_mono)
plt.show() 
for feature in ['Max. conc. A', 'IC50 A', 'H A', 'Einf A']:
  scatterplot_feature(feature, good_drugs_mono, bad_drugs_mono)<jupyter_output>/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  X = np.atleast_1d(X.T if isinstance(X, np.ndarray) else np.asarray(X))
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  X = np.atleast_1d(X.T if isinstance(X, np.ndarray) else np.asarray(X))
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tupl[...]<jupyter_text>We don't see much difference after playing with the threshold, meaning that it is not only a matter of monotherapy features, of each drugs, interaction has to be taken into account.  Checking abs(diff)We were unable to see significant differences between the values of monotherapy features of good and bad predicted drugs. We will now try with more complex features such as the absolute value of the difference of each monotherapy value between the drug pair and its mean.For that matter we have now to create a diff_df and av_df that are not commutative, because we don't want to overcount.<jupyter_code>test_df.columns
diff_df = pd.DataFrame()
diff_df["diff_Max conc"] = np.abs(test_df['Max. conc. A'] -test_df['Max. conc. B'] )
diff_df["diff_IC50"] = np.abs(test_df['IC50 A'] -test_df['IC50 B'] )
diff_df["diff_'H"] = np.abs(test_df['H A'] -test_df['H B'] )
diff_df["diff_Einf"] = np.abs(test_df['Einf A'] -test_df['Einf B'] )
diff_df['Combination ID'] = test_df['Combination ID']<jupyter_output><empty_output><jupyter_text>We modify good_bad_split so that we can generalize for every DataFrame of feautres<jupyter_code>"up_thres is the threshold for highly predictable drugs"
"low_thres for the lowly predictable drugs" 

def good_bad_split(drug_pair_pc, up_thres, low_thres, feature_df):
  good_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] >= up_thres]
  bad_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] <= low_thres]
  good_drugs_feat = feature_df[feature_df["Combination ID"].isin(good_drugs["Combination ID"])]
  bad_drugs_feat = feature_df[feature_df["Combination ID"].isin(bad_drugs["Combination ID"])]

  return good_drugs_feat, bad_drugs_feat

good_drugs_diff, bad_drugs_diff = good_bad_split(drug_pair_pc,
                                                 up_thres = 0.5, 
                                                 low_thres = 0, 
                                                 feature_df = diff_df)

for feature in ["diff_Max conc", "diff_IC50", "diff_'H", "diff_Einf"]:
  
  boxplot_feature(feature, good_drugs_diff, bad_drugs_diff)
plt.show() 
for feature in ["diff_Max conc", "diff_IC50", "diff_'H", "diff_Einf"]:
  scatterplot_feature(feature, good_drugs_diff, bad_drugs_diff)
test_df.head(10)
diff_df.head(10)<jupyter_output><empty_output><jupyter_text>We don't see anything, again.We try averages next:<jupyter_code>av_df = pd.DataFrame()
av_df["av_Max conc"] = (test_df['Max. conc. A'] +test_df['Max. conc. B'])/2
av_df["av_IC50"] = (test_df['IC50 A'] +test_df['IC50 B'] )/2
av_df["av_'H"] = (test_df['H A'] +test_df['H B'] )/2
av_df["av_Einf"] = (test_df['Einf A'] +test_df['Einf B'] )/2
av_df['Combination ID'] = test_df['Combination ID']
"up_thres is the threshold for highly predictable drugs"
"low_thres for the lowly predictable drugs" 

def good_bad_split(drug_pair_pc, up_thres, low_thres, feature_df):
  good_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] >= up_thres]
  bad_drugs = drug_pair_pc[drug_pair_pc["Pearson coefficient"] <= low_thres]
  good_drugs_feat = feature_df[feature_df["Combination ID"].isin(good_drugs["Combination ID"])]
  bad_drugs_feat = feature_df[feature_df["Combination ID"].isin(bad_drugs["Combination ID"])]

  return good_drugs_feat, bad_drugs_feat

good_drugs_av, bad_drugs_av = good_bad_split(drug_pair_pc,
                                                 up_thres = 0.5, 
                                                 low_thres = 0, 
                                                 feature_df = av_df)

for feature in ["av_Max conc", "av_IC50", "av_'H", "av_Einf"]:
  
  boxplot_feature(feature, good_drugs_av, bad_drugs_av)
plt.show() 
for feature in ["av_Max conc", "av_IC50", "av_'H", "av_Einf"]:
  scatterplot_feature(feature, good_drugs_av, bad_drugs_av)<jupyter_output>/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  X = np.atleast_1d(X.T if isinstance(X, np.ndarray) else np.asarray(X))
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  X = np.atleast_1d(X.T if isinstance(X, np.ndarray) else np.asarray(X))
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1376: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tupl[...]<jupyter_text>Distribution of WPC and plotsWe are going to do some statistics, we can do:1.   Histogram of WPC with and without low number of cell lines. Check if there is a difference and a suitable threshold for good and bad combinations2.   Notched Boxplot for WPC depending on n_cl. 3.   Rank plot where we rank predicted and true synergy score. Then throw a theoretical dotted line y = x and then calculate Spearman rank correlation. Or check if that's what we should check.4.   We can do some statistics with the predicted and expected synergy using Notched Boxplot, we can check the WPC for cell lines depending of tissue of origin, that's interesting.<jupyter_code>"Without restriction of cell line"
sns.histplot(drug_pair_pc["Pearson coefficient"])
sns.histplot(drug_pair_pc["n_cl"], bins=4)
sns.set(style="darkgrid")
notched = sns.boxplot(x=drug_pair_pc["n_cl"], y=drug_pair_pc["Pearson coefficient"], notch=True)
sns.stripplot(x=drug_pair_pc["n_cl"], y=drug_pair_pc["Pearson coefficient"],
              size=4, color=".3", linewidth=0)
notched.set(xlabel = "Number of cell lines", ylabel="Pearson Correlation")
fig = notched.get_figure()
fig.savefig("/content/drive/MyDrive/archivos_tfm/notched_box_wpc_n_cl.pdf")<jupyter_output><empty_output><jupyter_text>Now that we have added the LB dataset we see that we have a minimum of 5 cell lines per drug pair. Not bad.Now we do the rank plot for predicted and expected synergy on the test set. We use the test set that included the LB split.<jupyter_code>y_pred = np.loadtxt("/content/drive/MyDrive/archivos_tfm/y_pred_current_best_lb.txt")<jupyter_output><empty_output><jupyter_text>We add the observed synergy to a copy of the test_df dataframe<jupyter_code>pred_obs = test_df.copy()
pred_obs.columns = ['Cell line name', 'Compound A', 'Compound B', 'Max. conc. A',
       'Max. conc. B', 'IC50 A', 'H A', 'Einf A', 'IC50 B', 'H B', 'Einf B',
       'Experimental synergy score', 'QA', 'Combination ID', 'Challenge', 'Dataset']
pred_obs["Predicted synergy score"] = y_pred
pred_obs.columns
pred_obs.sort_values(by = "Experimental synergy score", inplace = True)
plt.plot(pred_obs["Experimental synergy score"], 
                pred_obs['Predicted synergy score'], "bo")
plt.show()
print(stats.spearmanr(pred_obs["Experimental synergy score"], 
                                               pred_obs['Predicted synergy score']))

print(stats.pearsonr(pred_obs["Experimental synergy score"], 
                                               pred_obs['Predicted synergy score']))<jupyter_output><empty_output><jupyter_text>Nothing to see, if the results were perfect we would see the points falling in the y=x line. These experiments have a lot of noise and we already knew that this was a bad metric altogether.We can also see the WPC with the rank<jupyter_code>plt.plot(np.arange(len(drug_pair_pc)), 
         drug_pair_pc.sort_values(by = "Pearson coefficient")["Pearson coefficient"], "bo")

plt.xlabel("Rank", fontsize = 13)
plt.ylabel("Pearson Correlation", fontsize = 13)
plt.savefig("/content/drive/MyDrive/archivos_tfm/pearson_correlation_rank.pdf")
drug_pair_pc["Pearson coefficient"]<jupyter_output><empty_output><jupyter_text>GDSC Tissue descriptor and synergyWe are going to compare the predicted and observed synergy depending on things such as gdsc descriptor, pathway, function, etc. For that we use the dictionaries obtained for visualizing the embeddings<jupyter_code>path = "/content/drive/MyDrive/archivos_tfm/saezr_nc19_DataS1_modificado.xlsx"
sheet_to_df_map = pd.read_excel(path, sheet_name=None)
df_cell_lines = sheet_to_df_map["Cell lines"]

df_drug_port = sheet_to_df_map["Drug portfolio"]

gdsc_dict = dict(zip(df_cell_lines["Cell line name"], df_cell_lines["GDSC tissue descriptor 1"]))
function_dict = dict(zip(df_drug_port["Challenge drug name"], df_drug_port["Function"]))
pathway_dict = dict(zip(df_drug_port["Challenge drug name"], df_drug_port["Pathway"]))<jupyter_output>/usr/local/lib/python3.7/dist-packages/openpyxl/worksheet/_reader.py:312: UserWarning: Unknown extension is not supported and will be removed
  warn(msg)<jupyter_text>We create the new columns that have the additional information of drugs and cell lines<jupyter_code>pred_obs["GDSC tissue descriptor 1"] = pred_obs['Cell line name'].map(gdsc_dict)

pred_obs["Compound A Function"] = pred_obs['Compound A'].map(function_dict)
pred_obs["Compound B Function"] = pred_obs['Compound B'].map(function_dict)

pred_obs["Compound A Pathway"] = pred_obs['Compound A'].map(pathway_dict)
pred_obs["Compound B Pathway"] = pred_obs['Compound B'].map(pathway_dict)
pred_obs.to_csv("/content/drive/MyDrive/archivos_tfm/pred_obs_lb.csv")
pred_obs.columns
sns.set(rc={'figure.figsize':(8.7,6.27)})

notched = sns.boxplot(x=pred_obs["GDSC tissue descriptor 1"], y=pred_obs['Experimental synergy score'], notch=True)

fig = notched.get_figure()<jupyter_output><empty_output><jupyter_text>We cap the experimental synergy to 100 for better visualization<jupyter_code>def capped(x):
  if np.abs(x) <= 100:
    return x
  else:
    return np.sign(x)*100
exp_syn_cap = [capped(x) for x in pred_obs['Experimental synergy score']]
pred_obs['Experimental synergy score capped'] = exp_syn_cap
sns.set_theme(style="ticks")
fig1 = plt.figure(figsize=(9, 11))


sns.boxplot(x='Experimental synergy score capped', 
            y="GDSC tissue descriptor 1", data=pred_obs,
            whis=[0, 100], width=.6, palette="vlag", notch=True)
sns.stripplot(x='Experimental synergy score capped', y="GDSC tissue descriptor 1", data=pred_obs,
              size=4, color=".3", linewidth=0)
ax.set(xlabel ='Experimental synergy score' ,ylabel="GDSC tissue descriptor 1")
#sns.despine(trim=True, left=True)
plt.savefig("/content/drive/MyDrive/archivos_tfm/experimental_synergy_tissue.pdf",
            bbox_inches = 'tight')
sns.set_theme(style="ticks")
fig2 = plt.figure(figsize=(9, 11))


sns.boxplot(x='Predicted synergy score', y="GDSC tissue descriptor 1", data=pred_obs,
            whis=[0, 100], width=.6, palette="vlag",notch=True)
sns.stripplot(x='Predicted synergy score', y="GDSC tissue descriptor 1", data=pred_obs,
              size=4, color=".3", linewidth=0)
ax.set(xlabel ='Predicted synergy score' ,ylabel="GDSC tissue descriptor 1")
#sns.despine(trim=True, left=True)
plt.savefig("/content/drive/MyDrive/archivos_tfm/predicted_synergy_tissue.pdf",
            bbox_inches = 'tight')
plt.subplot(1, 2, 1)
ax1 =sns.boxplot(x='Experimental synergy score capped', 
            y="GDSC tissue descriptor 1", data=pred_obs,
            whis=[0, 100], width=.6, palette="vlag", notch=True)
sns.stripplot(x='Experimental synergy score capped', y="GDSC tissue descriptor 1", data=pred_obs,
              size=4, color=".3", linewidth=0)
ax1.set(xlabel ='Experimental synergy score' ,ylabel="GDSC tissue descriptor 1")

plt.subplot(1, 2, 2)
ax2 = sns.boxplot(x='Predicted synergy score', y="GDSC tissue descriptor 1", data=pred_obs,
            whis=[0, 100], width=.6, palette="vlag",notch=True)
sns.stripplot(x='Predicted synergy score', y="GDSC tissue descriptor 1", data=pred_obs,
              size=4, color=".3", linewidth=0)
ax2.set(xlabel ='Predicted synergy score', ylabel = "")
ax2.set(yticklabels=[])

plt.show()
plt.savefig("/content/drive/MyDrive/archivos_tfm/predicted_observed_synergy_tissue.pdf",
            bbox_inches = 'tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11961669787421522, 'fraction_numerical': 0.012446304659103425, 'mean_word_length': 5.077309236947791, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '200047', 'license': 'Apache-2.0', 'n_tokens_mistral': 6585, 'n_tokens_neox': 6211, 'n_words': 1600, 'path': 'Drug_pair_analysis.ipynb', 'repo_name': 'edmanft/Drug_Synergy', 'size': '18158'}