<jupyter_start><jupyter_text>GSEA analysis on leukemia dataset<jupyter_code>%load_ext autoreload
%autoreload 2
from gsea import *
import numpy as np
%pylab
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>genes, D, C = read_expression_file("data/leukemia.txt")
gene_sets, gene_set_names = read_genesets_file("data/pathways.txt", genes)
gene_set_hash = {}
for i in range(len(gene_sets)):
    gene_set_hash[gene_set_names[i][0]] = {'indexes':gene_sets[i],'desc':gene_set_names[i][1]}
# verify that the dimensions make sense
len(genes),D.shape,len(C)<jupyter_output><empty_output><jupyter_text>Enrichment score calculationsWe graphically present the calculation of ES.<jupyter_code>L,r = rank_genes(D,C)<jupyter_output><empty_output><jupyter_text>See if the first genes in *L* are indeed correlated with *C*<jupyter_code>scatter(D[L[1],:],C)
scatter(D[L[-1],:],C)
scatter(D[L[1000],:],C)<jupyter_output><empty_output><jupyter_text>Graphical ilustration of ES calculations<jupyter_code>p_exp = 1
def plot_es_calculations(name, L, r):
    S = gene_set_hash[name]['indexes']
    N = len(L)
    S_mask = np.zeros(N)
    S_mask[S] = 1
    # reorder gene set mask
    S_mask = S_mask[L]
    N_R = sum(abs(r*S_mask)**p_exp)
    P_hit = np.cumsum(abs(r*S_mask)**p_exp)/N_R if N_R!=0 else np.zeros_like(S_mask)
    N_H = len(S)
    P_mis = np.cumsum((1-S_mask))/(N-N_H) if N!=N_H else np.zeros_like(S_mask)
    idx = np.argmax(abs(P_hit - P_mis))
    print("ES =", P_hit[idx]-P_mis[idx])
    f, axarr = plt.subplots(3, sharex=True)
    axarr[0].plot(S_mask)
    axarr[0].set_title('gene set %s' % name)
    axarr[1].plot(r)
    axarr[1].set_title('correlation with phenotype')
    axarr[2].plot(P_hit-P_mis)
    axarr[2].set_title('random walk')
L,r = rank_genes(D,C)
plot_es_calculations('CBF_LEUKEMIA_DOWNING_AML', L, r)<jupyter_output>ES = 0.214604376897<jupyter_text>Random phenotype labelsNow let's assign phenotype labels randomly. Is the ES much different?<jupyter_code>N, k = D.shape
pi = np.array([np.random.randint(0,2) for i in range(k)])
L, r = rank_genes(D,pi)
print(pi)
plot_es_calculations('CBF_LEUKEMIA_DOWNING_AML', L, r)<jupyter_output>[0 1 0 1 1 0 1 1 0 1 1 0 0 0 0 1 1 1 0 0 1 0 1 1 1 0 1 0 1 1 0 1 1 0 0 1 0
 0 0 1 0 0 0 1 0 1 1 0 1 1 0 0 1 0 1 0 1 1 1 1 0 1 1 1 0 1 1 1 0 0 1 0]
ES = -0.446675459428<jupyter_text>GSEA analysis<jupyter_code># use `n_jobs=-1` to use all cores
%time order, NES, p_values = gsea(D, C, gene_sets, n_jobs=-1)
from IPython.display import display, Markdown
s = "| geneset | NES | p-value | number of genes in geneset |\n |-------|---|---|---|\n "
for i in range(len(order)):
    s = s + "| **%s** | %.3f | %.7f | %d |\n" % (gene_set_names[order[i]][0], NES[i], p_values[i], len(gene_sets[order[i]]))
display(Markdown(s))<jupyter_output><empty_output><jupyter_text>Multiple Hypotesis testingWe present two example gene sets. One with a high *NES* and low *p-value* and one with a low *NES* and a high *p-value*. We plot the histograms of null distribution for ES.<jupyter_code>name = 'DNA_DAMAGE_SIGNALLING'
L,r = rank_genes(D,C)
plot_es_calculations(name, L, r)
n = 1000
S =  gene_set_hash[name]['indexes']
L, r = rank_genes(D,C)
ES = enrichment_score(L,r,S)
ES_pi = np.zeros(n)
for i in range(n):
    pi = np.array([np.random.randint(0,2) for i in range(k)])
    L, r = rank_genes(D,pi)
    ES_pi[i] = enrichment_score(L,r,S)
hist(ES_pi,bins=100)
plot([ES,ES],[0,20],'r-',label="ES(S)")
title("Histogram of ES vlues for random phenotype labels.\nRed line is ES for the selected gene set.")

name = 'tcrPathway'
L,r = rank_genes(D,C)
plot_es_calculations(name, L, r)
n = 1000
S =  gene_set_hash[name]['indexes']
L, r = rank_genes(D,C)
ES = enrichment_score(L,r,S)
ES_pi = np.zeros(n)
for i in range(n):
    pi = np.array([np.random.randint(0,2) for i in range(k)])
    L, r = rank_genes(D,pi)
    ES_pi[i] = enrichment_score(L,r,S)
hist(ES_pi,bins=100)
plot([ES,ES],[0,20],'r-',label="ES(S)")
title("Histogram of ES vlues for random phenotype labels.\nRed line is ES for the selected gene set.")<jupyter_output><empty_output><jupyter_text>Performance optimizations<jupyter_code>%timeit L,R = rank_genes(D,C)
%timeit ES = enrichment_score(L,r,S)
%prun order, NES, p_values = gsea(D, C, gene_sets)<jupyter_output>
{'fraction_non_alphanumeric': 0.15596544478169508, 'fraction_numerical': 0.03478869950968947, 'mean_word_length': 5.272327964860907, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '14562', 'license': 'MIT', 'n_tokens_mistral': 1979, 'n_tokens_neox': 1800, 'n_words': 423, 'path': 'Leukemia.ipynb', 'repo_name': 'mrcinv/GSEA.py', 'size': '4283'}