<jupyter_start><jupyter_text>Differential Expression analysis About this labOn this lab, we will try to understand the underlying mechanisms and concepts behind differential expression analysis and multiple hypothesis testing. You will be asked to reflect on some core concepts of the methods and to write down your interpretations. It is important that you understand what you are doing, and thoughtful answers are expected.To successfully complete the lab you have to answer all questions and submit them in **PDF format** in Canvas. You can, and should, discuss with your classmates.You will also be provided with some optional bonus questions that involve a little more research and programming on your part, but successful completion of those will award you extra points in the final examination.Let's begin! As in every notebook, we first begin by importing packages.<jupyter_code>import pandas as pd
import numpy as np
import scipy as sp
import seaborn as sns
import matplotlib.pyplot as plt
from ipywidgets import interact, interact_manual
from ipywidgets import IntSlider, FloatSlider, Dropdown, Text

interact_gen=interact_manual.options(manual_name="Generate data")
interact_plot=interact_manual.options(manual_name="Plot")
interact_calc=interact_manual.options(manual_name="Calculate tests")<jupyter_output><empty_output><jupyter_text>1 Differential expression analysis, an intuitionAs you have learned, differential expression analysis tries to answer the following question: Is there any **significant** difference between the expression of two (or more) different conditions?In this lab we will first be exploring the concepts using simulated data, which allows us to finely control what we are looking at and most importantly, know the true model behind the data. Later you will be asked to explore these concepts using real life data, on which you will have to look for the effects without knowing if there is really an undelying one.Our simulated data is very basic, and corresponds to two separately generated sets.On the first the data is generated by sampling from **two separate** normal distributions, with the data points colored according to which distribution they are generated from.On the second set, the data is generated by a **single** normal distribution, but with the same mean and standard deviation as the combined distribution from the previous set, and where the points are colored at random.In the language of differential expression, the colors represent the different conditions we want to test (e.g. treated/untreated, healthy/sick, etc.).The second data set then represents our **null hypothesis** $H_0$, and it says that there is no difference between the "expression" of the two conditions.Finally, the first set then represents the **alternative hypothesis** $H_1$ that says *there is* a difference between the two conditions.I hope the explanation of this was clear, but in case you have any questions, please ask the TA as this is very important for the rest of the exercise.Ok, let's try to get an intuition on differential expression fist. After executing the code below you will be presented with 2 plots and some sliders.The sliders control the parameters used for generating the samples from the two different sets, and plot the results.**The catch here is that you do not know which plot represents data under $H_0$ and which represents data under $H_1$.**Every time you move one of the sliders, the data is regenerated, and $H_0$ and $H_1$ **will be assigned randomly the left or right panes**, and you will have to reasses which one do you think represents the null and alternative hypothesis.Play around with the sliders as long as you feel necessary to answer these questions.An explanation of what the slides control:* Distance: the difference between the mean of each distribution in $H_1$* Dispersion: the standard deviation $\sigma$ of those distributions* Samples: The number of samples from each distribuition. Questions* 1.1 What is the influence of **each** variable in your ability to distinguish between $H_0$ and $H_1$?* 1.2 How do the different variables interact with each other? And how does it relate to your ability to distinguish the data?* 1.3 It is easy to see that often data under $H_1$ looks like data under $H_0$, but do you think you could ever misidentify data under $H_0$ as data under $H_1$?<jupyter_code>def generate_data(distance,dispersion, n_points):
    dispersion = [[dispersion,0],[0,dispersion]]
    cond1 = np.random.multivariate_normal([0,0], dispersion, n_points)
    cond2 = np.random.multivariate_normal([distance,0], dispersion, n_points)
    data1 = pd.DataFrame(cond1, columns=['x','y'])
    data1['Condition'] = 1
    data2 = pd.DataFrame(cond2, columns=['x','y'])
    data2['Condition'] = 2
    data = pd.concat([data1,data2])
    return data

def generate_h0(distance,dispersion, n_points):
    dispersion = [[dispersion + (distance/2)**2,0],[0,dispersion]]
    cond1 = np.random.multivariate_normal([distance/2,0], dispersion, n_points)
    cond2 = np.random.multivariate_normal([distance/2,0], dispersion, n_points)
    data1 = pd.DataFrame(cond1, columns=['x','y'])
    data1['Condition'] = 1
    data2 = pd.DataFrame(cond2, columns=['x','y'])
    data2['Condition'] = 2
    data = pd.concat([data1,data2])
    return data

def plot_random(data,data_s):
    f, ax = plt.subplots(1, 2)
    r1 = np.random.permutation(2)
    colors = np.random.permutation(['red','blue']).tolist()
    sns.scatterplot(x='x', y='y', data=data, hue = 'Condition', palette=colors, ax = ax[r1[0]], legend  = False)
    sns.scatterplot(x='x', y='y', data=data_s, hue = 'Condition', palette=colors, ax = ax[r1[1]], legend = False)

def dist_plot(Distance,Dispersion,Samples):
    plt.rcParams["figure.figsize"] = (20,10)
    data = generate_data(Distance,Dispersion, Samples)
    data_s = generate_h0(Distance,Dispersion, Samples)
    plot_random(data,data_s)
    plt.show()
interact(dist_plot, Distance = FloatSlider(min=0,max=10,value=1, continuous_update=False), Dispersion = FloatSlider(min=0,max=100,value=1, continuous_update=False), Samples=IntSlider(min=5,max=500, continuous_update=False))<jupyter_output><empty_output><jupyter_text>2 The t-testNow that you have a feel for what differential expression tries to achieve, and some obstacles on the way of doing that, lets introduce a statistical test to help us.We will use the **t-test**, which will give us p-value to assist us in answering the question.(As a technicality, you may have realized that the samples only differ on the $x$ axis, so we are only performing the test on that dimension)The setup here is nearly the same as above, just now you will be presented with the p-value resulting from the t-test as the title for each plot (rounded to 5 digits). Again, you should play around with the sliders and answer the following questions: Questions* 2.1 Look up the Student's t-test on your favourite online encyclopedia. What are the assumptions this test makes about the data? Do this assumptions hold here?* 2.2 The t-test gives us a p-value. How do you interpret this value on this setup?* 2.3 Does the p-value confirm your intuitions from the previous part? How would you augment your previous answers (1.1 - 1.3) on that light?* 2.4 Our initial question was "is there any significant difference between the two cases". As an answer, it is usual to reject $H_0$ for p-values lower than 0.05. Pick a very small distance between both distributions and see if you could make it significant using the other sliders, specially the sample size, and then comment on the relations of statistical significance and biological significance.<jupyter_code>def plot_random_ttest(data,data_s):
    plt.rcParams["figure.figsize"] = (20,10)
    cond1 = data[data['Condition'] == 1]
    cond2 = data[data['Condition'] == 2]
    cond1_s = data_s[data_s['Condition'] == 1]
    cond2_s = data_s[data_s['Condition'] == 2]
    ttestp = sp.stats.ttest_ind(cond1, cond2).pvalue[0]
    ttestp_s= sp.stats.ttest_ind(cond1_s, cond2_s).pvalue[0]
    f, ax = plt.subplots(1, 2)
    r1 = np.random.permutation(2)
    colors = np.random.permutation(['red','blue']).tolist()
    sns.scatterplot(x='x', y='y', data=data, hue = 'Condition', palette=colors, ax = ax[r1[0]], legend  = False)
    ax[r1[0]].set_title('p = ' + str(ttestp), fontsize=24)
    sns.scatterplot(x='x', y='y', data=data_s, hue = 'Condition', palette=colors, ax = ax[r1[1]], legend = False)
    ax[r1[1]].set_title('p = ' + str(ttestp_s), fontsize=24)
    

def dist_plot_ttest(Distance,Dispersion,Samples):
    plt.rcParams["figure.figsize"] = (20,10)
    data = generate_data(Distance,Dispersion, Samples)
    data_s = generate_h0(Distance,Dispersion, Samples)
    plot_random_ttest(data,data_s)
    plt.show()
interact(dist_plot_ttest, Distance = FloatSlider(min=0,max=10,value=1, continuous_update=False), Dispersion = FloatSlider(min=0,max=100,value=1, continuous_update=False), Samples=IntSlider(min=5,max=500, continuous_update=False))<jupyter_output><empty_output><jupyter_text>Exercise, differential gene expression analysis with TCGA dataNow it is time to test all this with real data. Below we will load TCGA data on breast cancer, both clinical and gene expression data taken from the tumours.The code below loads the data files, performs the tests and plots the result (including p-value from the t-test). All you need to do is input the name of the gene you want to test (using HUGO nomenclature), as well how you want to separate the groups of patients that you want to test for differential expression.To find valid separations, you will have to look at the original file *brca_clin.tsv*. Remember that the input has to be exactly present in the file, and is case-sensitive.Questions:* 2.5 Gene expression data often comes in the form of count data. Is this type of data compatible with the assumptions of the t-test? If not, how to go around it?* 2.6 Look at the clinical data file (on the data directory, **brca_clin.tsv**) using an application of your choice (e.g. Excel, LibreOffice, etc) and select a few clinical conditions that you think will affect gene expression. Comment on the number of samples.* 2.7 Perform the test on some (at least 5) genes that you know the function of, use **relevant** clinical conditions to separate the groups, and interpret the results.<jupyter_code>clinical_data = pd.read_csv('../data/brca_clin.tsv.gz', sep ='\t', index_col=2)
clinical_data = clinical_data.iloc[4:,1:]
expression_data = pd.read_csv('../data/brca.tsv.gz', sep ='\t', index_col=1)
expression_data = expression_data.iloc[:,2:].T
clinical_data.head()
def gene_ttest(clinical_df, expression_df, gene, separator, cond1, cond2):
    try:
        expression = expression_df[gene]
    except:
        print('ERROR: Gene not found in data')
    try:
        group1 = clinical_df[separator] == cond1
        index1 = clinical_df[group1].index
        group2 = clinical_df[separator] == cond2
        index2 = clinical_df[group2].index
    except:
        print('ERROR: Clinical condition wrong')
    expression1 = expression[index1].dropna()
    expression2 = expression[index2].dropna()
    cond = np.append(np.repeat('Group_1',len(expression1)),(np.repeat('Group_2',len(expression2))))
    comb_expression = np.append(expression1.values, expression2.values)
    plot_data = pd.DataFrame([comb_expression,cond]).T
    plot_data.columns = ['Expression', 'Condition']
    plot_data['Expression'] = plot_data['Expression'].astype(float)
    p_val = sp.stats.ttest_ind(expression1, expression2).pvalue
    plt.rcParams["figure.figsize"] = (20,10)
    sns.set(font_scale=2)
    sns.boxplot(x='Expression', y='Condition', data=plot_data)
    plt.title("p = " + str(p_val))


def interact_gene_ttest(Gene, Criteria, Group_1, Group_2):
    gene_ttest(clinical_data, expression_data, Gene, Criteria, Group_1, Group_2)
interact_plot(interact_gene_ttest, Gene = Text('BRCA1'), Criteria=Text('Surgical procedure first'), Group_1 = Text('Simple Mastectomy'), Group_2=Text('Lumpectomy'))<jupyter_output><empty_output><jupyter_text>3 Multiple hypothesis testingSo far we have been interested in each gene individually. We ask the differential expression question for each gene and answer it with a p-value from our statistical test.However, with the rich data generated from high-throughput experiments, we usually want to ask several different question (i.e. propose multiple hypothesis), and we thus move to the world of multiple hypothesis testing, where the p-value is loses some of its usefulness because it leads to high **family-wise error rate**.We will here use the same data generation scheme as before, but done many times over for each $H_0$ and $H_1$. This takes much more time, so you will have the option of generating the data once and visualizing it in multiple ways.A good way visualizing how a statistical test behaves with multiple hypothesis is to visualize the distribution of the p-values it generates, and we will explore why in the questions below.Below you will be able to generate data and visualize it using a histogram. Play around with it a while, use your best discretion on how many bins to use for the histogram, and answer the questions below.**Obs:** if you regenerate data, the plot will not automatically redraw.New sliders:* Tests: the number of times the data is generated and the t-test performed under each $H$* Bins: the number of bins in the histogram. Questions* 3.1 When we say "multiple hypothesis", which are the hypothesis we are referring to?* 3.2 Explain with your own words how we can have high **family-wise error rate** even when using very strict p-value thresholds.* 3.3 Do you think it is easier to identify $H_0$ and $H_1$ in this setup? Why do you think that is the case?* 3.4 Comment on the distribution of p-values under $H_0$, its shape and how it changes from different realizations and how each parameter affects the distribution.* 3.5 Take one realization of your data under $H_0$ and roughly count how many p-values under 0.05 you have. Is this expected? Comment also in context of your answer to question 1.3.* 3.6 Comment on the distribution of p-values under $H_1$, its shape and how it changes from different realizations and how each parameter affects the distribution.<jupyter_code>def multiple_tests(distance, dispersion, n_points, n_tests_h0, n_tests_h1):
    h0_stats = []
    for i in range(n_tests_h0):
        data_h0 = generate_h0(distance,dispersion, n_points)
        cond1_s = data_h0[data_h0['Condition'] == 1]
        cond2_s = data_h0[data_h0['Condition'] == 2]
        ttestp_s= sp.stats.ttest_ind(cond1_s, cond2_s).pvalue[0]
        h0_stats.append(ttestp_s)

    h1_stats = []
    for i in range(n_tests_h1):
        data = generate_data(distance,dispersion, n_points)
        cond1 = data[data['Condition'] == 1]
        cond2 = data[data['Condition'] == 2]
        ttestp = sp.stats.ttest_ind(cond1, cond2).pvalue[0]
        h1_stats.append(ttestp)
        
    return h0_stats, h1_stats

def plot_hist(h0_stats, h1_stats, bins):
    plt.rcParams["figure.figsize"] = (20,10)
    f, ax = plt.subplots(1, 2)
    sns.distplot(h0_stats, bins = bins, ax = ax[0], kde=False, hist_kws={"range": [0,1]})
    sns.distplot(h1_stats, bins = bins, ax = ax[1], kde=False, hist_kws={"range": [0,1]})
    ax[0].set_title('H0')
    ax[0].set_xlim([0,1])
    ax[1].set_title('H1')
    ax[1].set_xlim([0,1])
    plt.show()
    
def interactive_generate(Distance, Dispersion, Samples, Tests):
    global G_h0_results, G_h1_results
    G_h0_results, G_h1_results = multiple_tests(Distance,Dispersion, Samples, Tests, Tests)
    
def interactive_double_hist(Bins):
    plot_hist(G_h0_results, G_h1_results, Bins)
interactive_generate(1, 1, 3, 100)
interact_gen(interactive_generate, Distance = FloatSlider(min=0.0,max=10.0,value=1), Dispersion = FloatSlider(min=1.0,max=100.0,value=1), Samples = IntSlider(min=5,max=100), Tests = IntSlider(min=5,max=1000,value=100))
interact_plot(interactive_double_hist, Bins=IntSlider(min=10,max=50, step = 10, continuous_update=False))<jupyter_output><empty_output><jupyter_text>4 Calculating FDRWhile we can, in this context, analyse both the data generated under $H_0$ and $H_1$ separately, in reality they are often mixed together, and we use the statistical test precisely to **tell if the data comes from $H_0$ or not**.As you have seen in the exercises above, given enough tests, you will almost always have instances where you reject the null hypothesis when it is true no matter which threshold of significance you choose. In this context, we call these errors **false discoveries**.The false discovery rate (FDR), then, is the rate on which these errors occur.On the next plot, you will see a **stacked histogram** on the left, showing the distribution of p-values of both $H$ visualized together, and you will see a red vertical line that shows your significance threshold (left of the line = significant).On the right you have a colection of small squares where each represents a test. They are ranked by p-value with the lowest on the top-left corner and the highest on the lower-right, and coloured by which $H$ it originated from ($H_0$&rarr;red, $H_1$&rarr;blue). As you move the threshold, you will notice that some of those will get shaded and change color (blue&rarr;green, red&rarr;orange), and those represent the test that are below your threshold.So, by moving the threshold slider, you will select on which tests you will reject the null, or in the context of differential expression, the data you believe to be *expressed differentially* between conditions, and that will contain data that is both **truly not coming from $H_0$** and those that were **incorrectly rejected**.In our very special case here, we know from before which data comes from $H_0$ and from $H_1$, so we can calculate the FDR exactly, and that will be shown on top of the pane on the right.New sliders:* H_ratio: the ratio of tests performed under each $H$, calculated as $H_0/H_1$* Threshold: p-value threshold for rejecting $H_0$ Questions* 4.1 What relations do you see between the FDR and each of the sliders. You may refer to the answers you gave on previous questions if you wish.* 4.2 How does the FDR relate to the shape of the distribution of p-values of $H_0$ and $H_1$?* 4.3 It is usual to think that a lower p-value threshold will lead to a lower FDR. Is this always true?<jupyter_code>def plot_hist_comb(h0_stats, h1_stats, bins, threshold):
    plt.rcParams["figure.figsize"] = (20,10)
    f, ax = plt.subplots(1, 2)
    ax[0].hist([h0_stats, h1_stats], bins = bins, stacked = True)
    ax[0].vlines(threshold, 0, ax[0].get_ylim()[1], color='r', linestyles = 'dashed')
    stats = pd.DataFrame(h1_stats + h0_stats, columns = ['p'])
    stats['h'] = np.repeat(1,len(h1_stats)).tolist() + np.repeat(0,len(h0_stats)).tolist()
    stats['good'] = (stats['p'] <= threshold)*1
    discoveries = stats.loc[stats['p'] <= threshold]
    fdr_true = 1 - sum(discoveries['h'])/len(discoveries['h'])
    stats = stats.sort_values('p')
    stack1 = stats['h']
    stack2 = stats['good']
    side = np.ceil(np.sqrt(len(stack1))).astype(int)
    pad1 = np.zeros(side*side)+9
    pad2 = np.zeros(side*side)+9
    pad1[:len(stack1)]=stack1
    pad2[:len(stack1)]=stack2 * 5
    ax[1].imshow(pad1.reshape(side,side), cmap = 'Set1')
    ax[1].imshow(pad2.reshape(side,side), cmap = 'Set1', alpha = 0.4)
    ax[1].set_title('FDR: ' + str(np.round(fdr_true, 3)), fontsize=24)
    plt.show()
    
def interactive_generate_n(Distance, Dispersion, Samples, Tests, H_Ratio):
    global G_h0_results_n, G_h1_results_n
    r = H_Ratio
    n_tests_h0 = np.round(r*Tests*2).astype(int)
    n_tests_h1 = np.round((1-r)*Tests*2).astype(int)
    G_h0_results_n, G_h1_results_n = multiple_tests(Distance,Dispersion, Samples, n_tests_h0, n_tests_h1)

def interactive_plot_hist_comb(Bins, Threshold):
    plot_hist_comb(G_h0_results_n, G_h1_results_n, Bins, Threshold)
interactive_generate_n(1, 10, 30, 100, 0.5)
interact_gen(interactive_generate_n, Distance = FloatSlider(min=0.0,max=10.0,value=1), Dispersion = FloatSlider(min=0.0,max=100.0,value=1), Samples = IntSlider(min=5,max=100), Tests = IntSlider(min=5,max=1000,value=100), H_Ratio = FloatSlider(min=0,max=1, step = 0.01, value = 0.5, continuous_update=False))
interact_plot(interactive_plot_hist_comb, Threshold = FloatSlider(min=0,max=0.5, step = 0.01,value=0.05, continuous_update=False), Bins=IntSlider(min=10,max=100, value=20, step = 10, continuous_update=False))<jupyter_output><empty_output><jupyter_text>5 Estimating the FDR, using $\pi_0$You probably realized by now that the p-values under $H_0$ follow a uniform distribution. This comes from the very definition of p-values itself (!), more on why on the bonus exercise. Therefore a very good way to see if your statistical test isn't testing for the wrong thing, is to see if the distribution of p-values from that test under $H_0$ follows those characteristics.Knowing that fact, if we know the number of hypothesis coming from $H_0$, we can estimate the FDR very accurately. Unfortunately, that is never the case, and we have to also estimate this number.We do so by estimating the proportion of data under $H_0$ from the total number of test we made, and we call this proportion $\pi_0$.Below, you will have the chance to estimate $\pi_0$ and then see if your estimation is correct.When you input your estimated $\pi_0$, there will be a horizontal line on the plot representing how a perfectly unifor distribuition would look like given that $\pi_0$ value.Every time you generate new data, the true $\pi_0$ will be selected at random from the range $[0,1]$, and you can try to estimate what it is by using the horizontal line and the distribution.See if you can find a good way to estimate it, and then answer some questions.New sliders:* Pi0: Your estimative for $\pi_0$, used to estimate the FDR* Separate: Separate or not $H_0$ and $H_1$ by color Questions* 5.1 Describe one way you came up with for estimating $\pi_0$, when it works and when it doesn't work. (For this question, use your intuition and the experience you had now playing with the plot, there is no need to refer to any literature on calculating $\pi_0$)* 5.2 Is it easier to estimate low or high $\pi_0$? Why?* 5.3 What does it mean to set $\pi_0$ to 1. Can we do this and still get reasonable results? Why?<jupyter_code>def plot_hist_pi0(h0_stats, h1_stats, bins, threshold, pi0, sep):
    plt.rcParams["figure.figsize"] = (20,10)
    stats = pd.DataFrame(h1_stats + h0_stats, columns = ['p'])
    stats['h'] = np.repeat(1,len(h1_stats)).tolist() + np.repeat(0,len(h0_stats)).tolist()
    stats = stats.loc[stats['p'] <= threshold]
    fdr_true = 1 - sum(stats['h'])/len(stats['h'])
    m = len(h0_stats + h1_stats)
    fdr_est = min((pi0*m*threshold)/len(stats['h']),1)
    h = pi0*m/bins
    if sep:
        plt.hist([h0_stats, h1_stats], bins = bins, stacked = True, color = ['tab:blue', 'tab:orange'])
        plt.title('Estimated FDR: ' + str(np.round(fdr_est, 3)) + ',   True FDR: ' + str(np.round(fdr_true, 3)), fontsize=24)
    else:
        plt.hist([h0_stats, h1_stats], bins = bins, stacked = True, color = ['tab:blue', 'tab:blue'])
        plt.title('Estimated FDR: ' + str(np.round(fdr_est, 3)), fontsize=24)
    plt.vlines(threshold, 0, plt.gca().get_ylim()[1], color='r', linestyles = 'dashed')
    plt.hlines(h, 0, 1, color='black', linestyles = 'dashed', linewidth=3)
    plt.show()
    
def interactive_plot_hist_pi0(Bins, Threshold, Pi0, Separate):
    plot_hist_pi0(G_h0_results_pi0, G_h1_results_pi0, Bins, Threshold, Pi0, Separate)
    
def interactive_generate_n_rand(Distance, Dispersion, Samples, Tests):
    global G_h0_results_pi0, G_h1_results_pi0
    r = np.random.uniform(0,1)
    n_tests_h0 = np.round(r*Tests*2).astype(int)
    n_tests_h1 = np.round((1-r)*Tests*2).astype(int)
    G_h0_results_pi0, G_h1_results_pi0 = multiple_tests(Distance,Dispersion, Samples, n_tests_h0, n_tests_h1)
interactive_generate_n_rand(1, 10, 30, 100)
interact_gen(interactive_generate_n_rand, Distance = FloatSlider(min=0.0,max=10.0,value=1), Dispersion = FloatSlider(min=0.0,max=100.0,value=1), Samples = IntSlider(min=5,max=100,value=10), Tests = IntSlider(min=5,max=1000,value=100))
interact_plot(interactive_plot_hist_pi0, Threshold = FloatSlider(min=0.01,max=0.5,value=0.05, step = 0.01, continuous_update=False), Bins=IntSlider(min=10,max=100, step = 10,value=20, continuous_update=False), Pi0 = FloatSlider(min=0,max=1, step = 0.01, continuous_update=False), Separate=False)<jupyter_output><empty_output><jupyter_text>Multiple hypothesis on TCGA dataNow all that is left to do the multiple hypothesis testing on our data. Here you will be able to choose two clinical groups and test all the genes for differential expression on those two groups.* 5.4 Using relevant clinical conditions (at least 5), perform the multiple tests.Report the resulting histograms and the estimated $\pi_0$, and give your interpretation of the results.<jupyter_code>def multiple_gene_ttest(clinical_df, expression_df, separator, cond1, cond2):
    p_vec = []
    try:
        group1 = clinical_df[separator] == cond1
        index1 = clinical_df[group1].index
        group2 = clinical_df[separator] == cond2
        index2 = clinical_df[group2].index
    except:
        print('Clinical condition wrong')
    for gene in expression_df.columns:
        try:
            expression = expression_df[gene]
        except:
            print('Gene not found in data')
        expression1 = expression[index1].dropna()
        expression2 = expression[index2].dropna()
        p_val = sp.stats.ttest_ind(expression1, expression2).pvalue
        p_vec.append(p_val)
    return p_vec

def plot_hist_pi0_TCGA(stats, bins, threshold, pi0):
    plt.rcParams["figure.figsize"] = (20,10)
    stats = np.array(stats)
    stats_accepted = stats[stats <= threshold]
    m = len(stats)
    fdr_est = min((pi0*m*threshold)/len(stats_accepted),1)
    h = pi0*m/bins
    plt.hist(stats, bins = bins)
    plt.title('Estimated FDR: ' + str(np.round(fdr_est, 3)), fontsize=24)
    plt.vlines(threshold, 0, plt.gca().get_ylim()[1], color='r', linestyles = 'dashed')
    plt.hlines(h, 0, 1, color='black', linestyles = 'dashed', linewidth=3)
    plt.show()


def interact_multiple_gene_ttest(Criteria, Group_1, Group_2):
    global TCGA_multiple_t
    TCGA_multiple_t = multiple_gene_ttest(clinical_data, expression_data, Criteria, Group_1, Group_2)
    
def interactive_plot_hist_TCGA(Bins, Threshold, Pi0):
    plot_hist_pi0_TCGA(TCGA_multiple_t, Bins, Threshold, Pi0)
interact_calc(interact_multiple_gene_ttest, Criteria=Text('Surgical procedure first'), Group_1 = Text('Simple Mastectomy'), Group_2=Text('Lumpectomy'))
interact_plot(interactive_plot_hist_TCGA, Threshold = FloatSlider(min=0.01,max=0.5,value=0.05, step = 0.01, continuous_update=False), Bins=IntSlider(min=10,max=100, step = 10,value=20, continuous_update=False), Pi0 = FloatSlider(min=0,max=1, step = 0.01, continuous_update=False))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10213094843767252, 'fraction_numerical': 0.025394722314232086, 'mean_word_length': 4.9288675540039275, 'pattern_counts': {'":': 2, '<': 37, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '319301', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 8673, 'n_tokens_neox': 8282, 'n_words': 3025, 'path': 'lab/Lab1.ipynb', 'repo_name': 'statisticalbiotechnology/cb2030', 'size': '27171'}