<jupyter_start><jupyter_text>Statistical background<jupyter_code>def score(dataset, cutoffs):
    cols = [col for col in dataset.columns if (col not in ["name"] and "score" not in col)]
    for c in cutoffs:
        name = "score"+str(c)
        dataset[name] = dataset[dataset[cols] >= c][cols].fillna(0).sum(axis=1)
def gumbel(x, mi, beta):
    z = (x - mi) / beta
    return np.exp(-(z + np.exp(-z))) / beta
cutoffs = [i/10.0 for i in xrange(1,10)]
colnames = ["score"+str(i) for i in cutoffs]
cnum = len(cutoffs)
all_data = pd.DataFrame(columns=["mol", "sample", "size", "cutoff", "score"])
k = 0
for i in glob("../data/statistical_model/results/" + method + "/*tab.bz2"):
    sample = i.split('/')[-1].split('.')[0]
    if k%10 == 0:
        print "processed", k, "samples"
    k += 1
    size = sample.split("_")[-1]
    dataset = pd.read_table(i, sep="\t")
    score(dataset, cutoffs)
    for j in xrange(cnum):
        tmp = dataset[["name", colnames[j]]]
        tmp.columns = ["mol", "score"]
        tmp["cutoff"] = cutoffs[j]
        tmp["size"] = size
        tmp["sample"] = sample
        all_data = pd.concat([all_data, tmp])
print "processed", k, "samples"
all_data["size"] = all_data["size"].astype(float)
all_data["score"] = all_data["score"].astype(float)
all_data["cutoff"] = all_data["cutoff"].astype(float)

print all_data.shape
all_data.head()
grid = sns.FacetGrid(all_data, col="cutoff", size=5, col_wrap=3, sharex=False, sharey=False, legend_out=True)
grid = (grid.map(plt.scatter, "size", "score").add_legend())
grid.savefig("../figures/supplementary/" + method + "_random_scores.pdf", bbox_inches="tight");
cutoff_groups = all_data.groupby(["cutoff", "size"])
means = cutoff_groups.mean()
means.reset_index(inplace=True)
means.head()
stdev = cutoff_groups.std()
stdev.reset_index(inplace=True)
stdev.head()
def raw_mean(x, a):
    return a*x

def raw_std(x, a, b):
    return a*(x**b)
params = {}

for c in cutoffs:
    params[c] = {}
    tmp = means[means["cutoff"] == c]
    x = tmp["size"]
    y = tmp["score"]
    try:
        best_vals, covar = curve_fit(raw_mean, x, y, p0=[1.0])
        params[c]["mean"] = list(best_vals)
    except:
        print "cannot find avg params for c == %s" % c
        del params[c]
        continue

    tmp = stdev[stdev["cutoff"] == c]
    x = tmp["size"]
    y = tmp["score"]
    try:
        best_vals, covar = curve_fit(raw_std, x, y, p0=[1.0, 1.0])
        params[c]["std"] = list(best_vals)
    except:
        print "cannot find std params for c == %s" % c
        del params[c]
params
f = open(("../data/statistical_model/" + method + ".params"), "w")
f.write(str(params))
f.close()
j = 0
grid_mean = sns.FacetGrid(means, col="cutoff", size=5, col_wrap=3, sharex=False, sharey=False, legend_out=True)
grid_mean = (grid_mean.map(plt.scatter, "size", "score").add_legend())
for ax in grid_mean.axes.flat:
    cutoff = cutoffs[j]
    if cutoff in params:
        x = np.linspace(0,1000)
        ax.plot(x, raw_mean(x, *params[cutoff]["mean"]), c="red", ls="--")
    j += 1
grid_mean.add_legend();
grid_mean.savefig("../figures/supplementary/" + method + "_random_mean.pdf", bbox_inches="tight")
j = 0
grid_std = sns.FacetGrid(stdev, col="cutoff", size=5, col_wrap=3, sharex=False, sharey=False, legend_out=True)
grid_std = grid_std.map(plt.scatter, "size", "score")
for ax in grid_std.axes.flat:
    cutoff = cutoffs[j]
    if cutoff in params:
        x = np.linspace(0,1000)
        ax.plot(x, raw_std(x, *params[cutoff]["std"]), c="red", ls="--")
    j += 1
grid_std.add_legend();
grid_std.savefig("../figures/supplementary/" + method + "_random_sd.pdf", bbox_inches="tight")
def normalize(row):
    cutoff = row["cutoff"]
    if cutoff in params:
        size = row["size"]
        score = row["score"]
        z_score = ((row["score"] - raw_mean(size, *params[cutoff]["mean"])) / 
                    raw_std(size, *params[cutoff]["std"]))
    else:
        z_score = np.nan
    return z_score
all_data["z_score"] = all_data.apply(normalize, axis=1)
all_data = all_data.dropna()
np.mean(all_data["z_score"]), np.std(all_data["z_score"])
beta = math.sqrt(6.0) / math.pi
mi = -0.5772156649*beta
values = all_data["z_score"]
bins = np.arange(np.min(values), np.max(values), 1)
y_real, bin_edges = np.histogram(values, bins=bins)
y_real = y_real / float(len(values))
x = bin_edges[:-1] + (bin_edges[1] - bin_edges[0]) / 2
f, ax = plt.subplots()
ax.plot(x, y_real, label="Real")
ax.plot(x, gumbel(x, mi, beta), label="Theoretical")
ax.legend()
f.savefig("../figures/supplementary/" + method + "_fit.pdf", bbox_inches="tight");
f, ax = plt.subplots()
ax.plot(x[:10], y_real[:10], 'o', label="Real")
x_smooth = np.linspace(-5, 10)
ax.plot(x_smooth, gumbel(x_smooth, mi, beta), label="Theoretical")
ax.legend()
f.savefig("../figures/supplementary/" + method + "_fit_notail.pdf", bbox_inches="tight");<jupyter_output><empty_output><jupyter_text>Results<jupyter_code>def pvalue(z):
    z = np.array(z)
    x = -np.exp(-z * math.pi / math.sqrt(6) - 0.5772156649)
    idx = np.where(z < 28)[0]
    res = np.empty_like(z)
    res = -(x + (x**2)/2.0 + (x**3)/6.0)
    res[idx] = 1.0 - np.exp(x[idx])
    return res
targets = []
datasets = {}

for i in glob("../data/results/"+method+"/*tab.bz2"):
    #print i
    name = i.split('/')[-1].split('.')[0]
    targets.append(name)
    datasets[name] = pd.read_table(i)
print targets
drug_targets = pd.read_csv("../data/drug_targets.csv", index_col=0)
drug_targets.head()
np.sum(drug_targets["pred"] == "new")
np.sum(drug_targets["pred"] == "old")
np.sum(drug_targets["pred"] == "drugbank")
trivial = pd.read_csv("../data/trivial_hits.csv", index_col=0)
chemb_low_activity = pd.read_csv("../data/chembl_low_activity_hits.csv", index_col=0)<jupyter_output><empty_output><jupyter_text>Raw scores<jupyter_code>all_data = pd.DataFrame(columns=["target", "drug", "cutoff", "score", "size"])
cutoffs.remove(0.9)
cnum = len(cutoffs)
cutoffs, cnum, method, len(targets)
for t in targets:
    print t
    score(datasets[t], cutoffs)
    size = datasets[t].shape[1] - cnum - 1
    for i in xrange(cnum):
        col = colnames[i]
        cutoff = cutoffs[i]
        tmp = datasets[t][["name", col]]
        tmp.columns = ["drug", "score"]
        tmp["target"] = t
        tmp["cutoff"] = cutoff
        tmp["size"] = datasets[t].shape[1] - cnum - 1
        all_data = pd.concat([all_data, tmp])

print all_data.shape
all_data.head()<jupyter_output>SCN5A
CRHR2
MC3R
ADORA3
TBXA2R
BDKRB2
PDE3A
CACNA1C
EDNRB
PDE4D
BDKRB1
GPR109A
ADRB1
HRH2
CACNA1B
AVPR1A
PTGS1
NPY2R
MLN
NR1I2
SLC29A1
AVPR2
NTSR1
MAOA
GIPR
NPY1R
(136448, 5)<jupyter_text>Z-scores, p-values, e-values<jupyter_code>all_data["z_score"] = all_data.apply(normalize, axis=1)
all_data = all_data.dropna()
all_data.head()
grid_z = sns.FacetGrid(all_data, col="cutoff", size=5, col_wrap=3, sharex=False, sharey=False,
                                      legend_out=True)
grid_z = (grid_z.map(sns.boxplot, "z_score").add_legend())
grid_z;
all_data["p_value"] = pvalue(all_data["z_score"])
all_data.head()
N = 73 * 656
N
all_data["e_value"] = all_data["p_value"] * N
all_data.head()
all_data.to_csv("../data/results/"+method+"_all_results.csv")<jupyter_output><empty_output><jupyter_text>Summary<jupyter_code>e_range = np.arange(-10, 1)
e_range
predictions = {e: {t: set() for t in targets} for e in e_range}

for e in e_range:
    pred = all_data[(all_data["e_value"] <= 10**e)].groupby("target")
    for t, group in pred:
        predictions[e][t] = predictions[e][t] | set(list(group["drug"]))

low_act = {}
for t in targets:
    low_act[t] = set(list(chemb_low_activity[(chemb_low_activity["target"] == t)]["drug"].unique()))

triv = {}
for t in targets:
    triv[t] = set(list(trivial[(trivial["target"] == t)]["drug"].unique()))

known = {}
for t in targets:
    known[t] = (set(list(drug_targets[(drug_targets["target"] == t)]["drug"].unique())) | triv[t] | low_act[t])
summary = pd.DataFrame(columns=["e_value", "target", "pred", "true", "not trivial", "low activity", "false",
                                "not found"])

i = 0
for e in e_range:
    for t in targets:
        pr = predictions[e][t]
        tr = len(pr & known[t])
        nt = len((pr & known[t]) - triv[t])
        la = len(pr & low_act[t])
        f = len(pr - known[t])
        nf = len((known[t] - pr))        
        summary.loc[i] = [10**e, t, len(pr), tr, nt, la, f, nf]
        i += 1

summary.head()
summary.to_csv("../data/results/"+method+"_summary.csv")
summary_e = summary.groupby("e_value").sum()[["pred", "true", "not trivial", "low activity", "false", "not found"]]
summary_e
summary_e.to_csv("../data/results/"+method+"_summary_e.csv")<jupyter_output><empty_output><jupyter_text>---  Summary for DeCAF<jupyter_code>tmp = pd.read_csv("../data/results/DeCAF_all_results.csv", index_col=0)
tmp = tmp.loc[tmp.apply(lambda x: (x["target"] in targets) and (x["cutoff"] != 0.9), axis=1)]
tmp.head()
predictions = {e: {t: set() for t in targets} for e in e_range}

for e in e_range:
    pred = tmp[(tmp["e_value"] <= 10**e)].groupby("target")
    for t, group in pred:
        predictions[e][t] = predictions[e][t] | set(list(group["drug"]))

low_act = {}
for t in targets:
    low_act[t] = set(list(chemb_low_activity[(chemb_low_activity["target"] == t)]["drug"].unique()))

triv = {}
for t in targets:
    triv[t] = set(list(trivial[(trivial["target"] == t)]["drug"].unique()))

known = {}
for t in targets:
    known[t] = (set(list(drug_targets[(drug_targets["target"] == t)]["drug"].unique())) | triv[t] | low_act[t])
summary = pd.DataFrame(columns=["e_value", "target", "pred", "true", "not trivial", "low activity", "false",
                                "not found"])

i = 0
for e in e_range:
    for t in targets:
        pr = predictions[e][t]
        tr = len(pr & known[t])
        nt = len((pr & known[t]) - triv[t])
        la = len(pr & low_act[t])
        f = len(pr - known[t])
        nf = len((known[t] - pr))        
        summary.loc[i] = [10**e, t, len(pr), tr, nt, la, f, nf]
        i += 1

summary.head()
summary_e = summary.groupby("e_value").sum()[["pred", "true", "not trivial", "low activity", "false", "not found"]]
summary_e
summary_e.to_csv("../data/results/DeCAF_subset_summary_e.csv")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19508403769552124, 'fraction_numerical': 0.016224618672884484, 'mean_word_length': 3.781235485369252, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '549107', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4186, 'n_tokens_neox': 3998, 'n_words': 814, 'path': 'supplementary/notebooks/5_usrcat.ipynb', 'repo_name': 'uibcdf/DeCAF', 'size': '10293'}