<jupyter_start><jupyter_text>Lab 4: EM Algorithm and Single-Cell RNA-seq Data  Name: <NAME> Here (Your netid here)  Due April 2, 2021 11:59 PM  Preamble (Don't change this)  Important Instructions - 1. Please implement all the *graded functions* in main.py file. Do not change function names in main.py.2. Please read the description of every graded function very carefully. The description clearly states what is the expectation of each graded function. 3. After some graded functions, there is a cell which you can run and see if the expected output matches the output you are getting. 4. The expected output provided is just a way for you to assess the correctness of your code. The code will be tested on several other cases as well.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%run main.py
module = Lab4()<jupyter_output><empty_output><jupyter_text>Part 1 : Expectation-Maximization (EM) algorithm for transcript quantification  IntroductionThe EM algorithm is a very helpful tool to compute maximum likelihood estimates of parameters in models that have some latent (hidden) variables.In the case of the transcript quantification problem, the model parameters we want to estimate are the transcript relative abundances $\rho_1,...,\rho_K$.The latent variables are the read-to-transcript indicator variables $Z_{ik}$, which indicate whether the $i$th read comes from the $k$th transcript (in which case $Z_{ik}=1$.In this part of the lab, you will be given the read alignment data.For each read and transcript pair, it tells you whether the read can be mapped (i.e., aligned) to that transcript.Using the EM algorithm, you will estimate the relative abundances of the trascripts.  Reading read transcript data - We have 30000 reads and 30 transcripts<jupyter_code>n_reads=30000
n_transcripts=30
read_mapping=[]
with open("read_mapping_data.txt",'r') as file :
    lines_reads=file.readlines()
for line in lines_reads :
    read_mapping.append([int(x) for x in line.split(",")])
read_mapping[:10]<jupyter_output><empty_output><jupyter_text>Rather than giving you a giant binary matrix, we encoded the read mapping data in a more concise way. read_mapping is a list of lists. The $i$th list contains the indices of the transcripts that the $i$th read maps to.  Reading true abundances and transcript lengths<jupyter_code>with open("transcript_true_abundances.txt",'r') as file :
    lines_gt=file.readlines()
ground_truth=[float(x) for x in lines_gt[0].split(",")]

with open("transcript_lengths.txt",'r') as file :
    lines_gt=file.readlines()
tr_lengths=[float(x) for x in lines_gt[0].split(",")]
ground_truth[:5]
tr_lengths[:5]<jupyter_output><empty_output><jupyter_text>Graded Function 1 : expectation_maximization (10 marks) Purpose : To implement the EM algorithm to obtain abundance estimates for each transcript.E-step :  In this step, we calculate the fraction of read that is assigned to each transcript (i.e., the estimate of $Z_{ik}$). For read $i$ and transicript $k$, this is calculated by dividing the current abundance estimate of transcript $k$ by the sum of abundance estimates of all transcripts that read $i$ maps to.M-step :  In this step, we update the abundance estimate of each transcript based on the fraction of all reads that is currently assigned to the transcript. First we compute the average fraction of all reads assigned to the transcript. Then, (if transcripts are of different lengths) we divide the result by the transcript length.Finally, we normalize all abundance estimates so that they add up to 1.Inputs - read_mapping (which is a list of lists where each sublist contains the transcripts to which a particular read belongs to. The length of this list is equal to the number of reads, i.e. 30000; tr_lengths (a list containing the length of the 30 transcripts, in order); n_iterations (the number of EM iterations to be performed)Output - a list of lists where each sublist contains the abundance estimates for a transcript across all iterations. The length of each sublist should be equal to the number of iterations plus one (for the initialization) and the total number of sublists should be equal to the number of transcripts.<jupyter_code>history=module.expectation_maximization(read_mapping,tr_lengths,20)
print(len(history))
print(len(history[0]))
print(history[0][-5:])
print(history[1][-5:])
print(history[2][-5:])<jupyter_output>30
21
[0.033769639494636614, 0.03381298624783303, 0.03384568373972949, 0.0338703482393148, 0.03388895326082054]
[0.0020082674603036053, 0.0019649207071071456, 0.0019322232152109925, 0.0019075587156241912, 0.0018889536941198502]
[0.0660581789629968, 0.06606927656035864, 0.0660765012689558, 0.06608120466668756, 0.0660842666518177]<jupyter_text>Expected Output - 3021[0.033769639494636614, 0.03381298624783303, 0.03384568373972948, 0.0338703482393148, 0.03388895326082054][0.0020082674603036053, 0.0019649207071071456, 0.0019322232152109925, 0.0019075587156241912, 0.0018889536941198502][0.0660581789629968, 0.06606927656035864, 0.06607650126895578, 0.06608120466668756, 0.0660842666518177] You can use the following function to visualize how the estimated relative abundances are converging with the number of iterations of the algorithm.<jupyter_code>def visualize_em(history,n_iterations) :
    #start code here
    fig, ax = plt.subplots(figsize=(8,6))
    for j in range(n_transcripts):   
        ax.plot([i for i in range(n_iterations+1)],[history[j][i] - ground_truth[j] for i in range(n_iterations+1)],marker='o')
    #end code here
visualize_em(history,20)<jupyter_output><empty_output><jupyter_text>Part 2 : Exploring Single-Cell RNA-seq data In a study published in 2015, Zeisel et al. used single-cell RNA-seq data to explore the cell diversity in the mouse brain. We will explore the data used for their study.You can read more about it [here](https://science.sciencemag.org/content/347/6226/1138).<jupyter_code>#reading single-cell RNA-seq data
lines_genes=[]
with open("Zeisel_expr.txt",'r') as file :
    lines_genes=file.readlines()
lines_genes[0][:300]<jupyter_output><empty_output><jupyter_text>Each line in the file Zeisel_expr.txt corresponds to one gene.The columns correspond to different cells (notice that this is the opposite of how we looked at this matrix in class).The entries of this matrix correspond to the number of reads mapping to a given gene in the corresponding cell.<jupyter_code># reading true labels for each cell
with open("Zeisel_labels.txt",'r') as file :
    true_labels = file.read().splitlines()<jupyter_output><empty_output><jupyter_text>The study also provides us with true labels for each of the cells.For each of the cells, the vector true_labels contains the name of the cell type.There are nine different cell types in this dataset.<jupyter_code>set(true_labels)<jupyter_output><empty_output><jupyter_text>Graded Function 2 : prepare_data (10 marks) :Purpose - To create a dataframe where each row corresponds to a specific cell and each column corresponds to the expressions levels of a particular gene across all cells. You should name the columns as "Gene_1", "Gene_2", and so on.We will iterate through all the lines in lines_genes list created above, add 1 to each value and take log.Each line will correspond to 1 column in the dataframeOutput - gene expression dataframe Note - All the values in the output dataframe should be rounded off to 5 digits after the decimal<jupyter_code>data_df=module.prepare_data(lines_genes)
print(data_df.shape)
print(data_df.iloc[0:3,:5])
print(data_df.columns)<jupyter_output>Index(['Gene_0', 'Gene_1', 'Gene_2', 'Gene_3', 'Gene_4', 'Gene_5', 'Gene_6',
       'Gene_7', 'Gene_8', 'Gene_9',
       ...
       'Gene_19962', 'Gene_19963', 'Gene_19964', 'Gene_19965', 'Gene_19966',
       'Gene_19967', 'Gene_19968', 'Gene_19969', 'Gene_19970', 'Gene_19971'],
      dtype='object', length=19972)<jupyter_text>Expected Output :``(3005, 19972)````   Gene_0    Gene_1    Gene_2   Gene_3    Gene_4``   ``0     0.0  1.38629  1.38629     0.0  0.69315````1     0.0  0.69315  0.69315     0.0  0.69315````2     0.0  0.00000  1.94591     0.0  0.69315``  Graded Function 3 : identify_less_expressive_genes (10 marks)Purpose : To identify genes (columns) that are expressed in less than 25 cells. We will create a list of all gene columns that have values greater than 0 for less than 25 cells.Input - gene expression dataframeOutput - list of column names which are expressed in less than 25 cells<jupyter_code>drop_columns = module.identify_less_expressive_genes(data_df)
print(len(drop_columns))
print(drop_columns[:10])<jupyter_output>5120
Index(['Gene_28', 'Gene_126', 'Gene_145', 'Gene_146', 'Gene_151', 'Gene_152',
       'Gene_167', 'Gene_168', 'Gene_170', 'Gene_173'],
      dtype='object')<jupyter_text>Expected Output : ``5120`` ``['Gene_28', 'Gene_126', 'Gene_145', 'Gene_146', 'Gene_151', 'Gene_152', 'Gene_167', 'Gene_168', 'Gene_170', 'Gene_173']``  Filtering less expressive genesWe will now create a new dataframe in which genes which are expressed in less than 25 cells will not be present<jupyter_code>df_new = data_df.drop(drop_columns, axis=1)
df_new.head()<jupyter_output><empty_output><jupyter_text>Graded Function 4 :  perform_pca (10 marks)Pupose - Perform Principal Component Analysis on the new dataframe and take the top 50 principal componentsInput - df_newOutput - numpy array containing the top 50 principal components of the data. Note - All the values in the output should be rounded off to 5 digits after the decimal Note - Please use random_state=365 for the PCA object you will create<jupyter_code>pca_data=module.perform_pca(df_new)
print(pca_data.shape)
print(type(pca_data))
print(pca_data[0:3,:5])<jupyter_output>(3005, 50)
<class 'numpy.ndarray'>
[[26.97148 -2.7244   0.62163 25.90148 -6.24736]
 [26.49135 -1.58774 -4.79315 24.01094 -7.25618]
 [47.82664  5.06799  2.15177 30.24367 -3.38878]]<jupyter_text>Expected Output : ``(3005, 50)````````[[26.97148  -2.7244  0.62163 25.90148 -6.24736]```` [26.49135 -1.58774 -4.79315 24.01094 -7.25618]`` `` [47.82664  5.06799  2.15177 30.24367 -3.38878]]``  (Non-graded) Function 5 :  perform_tsnePupose - Perform t-SNE on the pca_data and obtain 2 t-SNE componentsWe will use TSNE class of the sklearn.manifold package. Use random_state=1000 and perplexity=50Documenation can be found here - https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.htmlInput - pca_dataOutput - numpy array containing the top 2 tsne components of the data.**Note: This function will not be graded because of the random nature of t-SNE.**<jupyter_code>tsne_data50 = module.perform_tsne(pca_data)
print(tsne_data50.shape)
print(tsne_data50[:3,:])<jupyter_output>(3005, 2)
[[ 19.031317 -45.3434  ]
 [ 19.188553 -44.945473]
 [ 17.369982 -47.997364]]<jupyter_text>Expected Output :(These numbers can deviate a bit depending on your sklearn)``(3005, 2)````[[ 15.069608 -47.535984]```` [ 15.251476 -47.172073]`` `` [ 13.3932   -49.909657]]``<jupyter_code>fig, ax = plt.subplots(figsize=(12,8))
sns.scatterplot(tsne_data50[:,0], tsne_data50[:,1], hue=true_labels)
plt.show()<jupyter_output>/usr/local/lib/python3.9/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(
{'fraction_non_alphanumeric': 0.09594160076475189, 'fraction_numerical': 0.1074997827409403, 'mean_word_length': 5.375623268698061, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 2, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '113', 'license': 'MIT', 'n_tokens_mistral': 4475, 'n_tokens_neox': 3616, 'n_words': 1290, 'path': 'ECE365/genomics/Genomics_Lab4/ECE365-Genomics-Lab4-Spring21.ipynb', 'repo_name': 'debugevent90901/courseArchive', 'size': '11510'}