<jupyter_start><jupyter_text><jupyter_code>%matplotlib inline
import os
os.chdir("/content/drive/Shareddrives/Documentation/Data/Yale_Perturb/")
#!unzip processed.zip
os.chdir("./processed")
!pip install scanpy
import scanpy
import pandas as pd
sc_object = scanpy.read_h5ad("crispr_shallow.h5ad")
sc_object<jupyter_output><empty_output><jupyter_text>Explain the meta file* infected: idicating if the cell is effected by SARS CoVid virus (Not indicating if the cell has been perturbed)* labels: Only 6435 cells are with labels, the labels distrubution are as below<jupyter_code>x = pd.DataFrame(sc_object.obs['labels'].value_counts())
x.plot.bar(title = "All cells")
x.drop(index = "nan").plot.bar(title = "Cells with labels")<jupyter_output><empty_output><jupyter_text>Meta dataframe<jupyter_code>sc_object.obs<jupyter_output><empty_output><jupyter_text>Information about genes<jupyter_code>sc_object.var.loc['ACE2',]
sc_object.var
sc_object.obs['Infected'].value_counts()
sc_object.obs['batch'].value_counts()<jupyter_output><empty_output><jupyter_text>Expression matrix<jupyter_code>sc_object.X<jupyter_output><empty_output><jupyter_text>EM algorithm to fit the importance of each perturbationAdapted code based on https://github.com/asncd/MIMOSCA/blob/master/contrived-em_example.ipynb<jupyter_code>import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
matplotlib.rcParams['pdf.fonttype'] = 42
import sklearn
from sklearn import linear_model
import seaborn as sns
sns.set_context('talk')
sns.set_style('white')
sns.set_style('ticks')<jupyter_output><empty_output><jupyter_text>Bayes_cov_col<jupyter_code>def bayes_cov_col(Y,X,cols,lm):
    """
    @Y    = Expression matrix, cells x x genes, expecting pandas dataframe
    @X    = Covariate matrix, cells x covariates, expecting pandas dataframe
    @cols = The subset of columns that the EM should be performed over, expecting list
    @lm   = linear model object
    """
    #EM iterateit
    Yhat=pd.DataFrame(lm.predict(X))
    Yhat.index=Y.index
    Yhat.columns=Y.columns
    SSE_all=np.square(Y.subtract(Yhat))
    X_adjust=X.copy()


    df_SSE   = []
    df_logit = []

    for curcov in cols:

        curcells=X[X[curcov]>0].index
        #print("These are currcells")
        #print(curcells)
        # Curcells are the index with perturbation in the last assignment
        
        # If we have more than 2 cells perturbed
        if len(curcells)>2:
            #print("Did I do the correction?")
            X_notcur=X.copy()
            # Set the column all zeros (Create X_0 in the paper)
            X_notcur[curcov]=[0]*len(X_notcur)
            #print(f"This is the numerators {sum(X_notcur[curcov])}")

            # Extract the covariates and expression of perturbed cells in the last assignment
            X_sub=X_notcur.loc[curcells]
            Y_sub=pd.DataFrame(Y.loc[curcells])
            
            #print(f"This is Y_sub type: {type(Y_sub)}")
            #print(f"This is Y_sub shape: {Y_sub.shape}")
            #print(f"This is Y_sub : {Y_sub}")
            ynp = np.array(Y_sub)
            dino = 2.0 * ynp.var(axis=0)
            GENE_var = pd.DataFrame(2.0 * ynp.var(axis=0))
            GENE_var.index = Y_sub.columns
            # Get the variance as the denominators here the variace is across one column 
            # (a gene's variance across all cells perturbed)
            #GENE_var= 2.0 * Y_sub.var(axis=0)
            #print(f"This is variance of sub :{Y_sub.var(axis=0)}")
            #print(f"This is GENE_var: {GENE_var}")

            # The variance could be negative????
            vargenes= GENE_var[GENE_var>0].index
            #print(f"This is the vargenes {vargenes}")

            # X_0 * beta
            Yhat_notcur=pd.DataFrame(lm.predict(X_sub))
            Yhat_notcur.index=Y_sub.index
            Yhat_notcur.columns=Y_sub.columns

            SSE_notcur=np.square(Y_sub.subtract(Yhat_notcur))
            SSE=SSE_all.loc[curcells].subtract(SSE_notcur)
            SSE_sum=SSE.sum(axis=1)

            #print(f"This is the SSE_sum: {SSE_sum}")
            #print(f"This is the SSE shape: {SSE.shape}")
            #print(f"This is before sum : {SSE.div(dino+0.5)[vargenes]}")
            SSE_transform = SSE.div(dino+0.5)[vargenes].sum(axis=1)

            #print(f"This is the SSE_transform: {SSE_transform}")
            logitify=np.divide(1.0,1.0+np.exp(SSE_transform))#sum))
            
            #print(f"This is the logitify: {logitify}")

            df_SSE.append(SSE_sum)
            df_logit.append(logitify)

            X_adjust[curcov].loc[curcells]=logitify

    return X_adjust<jupyter_output><empty_output><jupyter_text>Refined function (Bug free with new pandas package version)<jupyter_code>def run_model_bycol(Y,X,EM_cols=None,modalpha=0.005,verbose=0):   
    """
    @Y    = Expression matrix, expecting pandas dataframe, cells x genes
    @X    = Covariate matrix, expecting pandas dataframe, cells x covariates
    @EM_cols = The subset of columns that the EM should be performed over, list
    """
    # Use the initial expression matrix and covariates assignment (perturbation) fit the inital beta
    enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=modalpha,max_iter=10000)
    enet.fit(X,Y)

    if verbose==1:
        print(enet.score(X,Y))


    # Be is initiated as the initial beta
    Be=pd.DataFrame(enet.coef_)
    Be.columns=X.columns
    Be.index=Y.columns
    
    # Get the predicted expression matrix Yhat from X
    Yhat=pd.DataFrame(enet.predict(X))
    Yhat.index=Y.index
    Yhat.columns=Y.columns

    if EM_cols is not None:

        # Get the new covariate matrix
        X_adjust=bayes_cov_col(Y,X,EM_cols,enet)

        print('done with EM')

        #refit model
        
        # Fit the new beta
        enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=0.0004,max_iter=10000)
     
        enet.fit(X_adjust,Y)
        Yhat=pd.DataFrame(enet.predict(X_adjust))
        Yhat.index=Y.index
        Yhat.columns=Y.columns

        if verbose==1:
            print(enet.score(X_adjust,Y))

        Be=pd.DataFrame(enet.coef_)
        Be.columns=X_adjust.columns
        Be.index=Y.columns
    else:
        X_adjust=X.copy()

    RES_out=Y.subtract(Yhat)  

    return(Be,X_adjust,RES_out)

Y_sub<jupyter_output><empty_output><jupyter_text>Simulation poplulation with 2 classes cells mixed1000 cells in the population:Half (500) are class 1Half (500) are class 2 All cells with 200 genesBut only partial of these 200 genes are infrmative genes (?) (y1 and y2)The other 180 genes are all noisy genes (y3)<jupyter_code>#simulate two multivariate normal distributions for two classes
#initial parameters
ncells=1000
ngenes=200
halfcells=int(ncells/2.0)
tenthcells=int(ncells/10.0)
informativegenes=0.1
effectsize=1.0
#Simulate Class 0 cells
y1=pd.DataFrame(np.random.normal(loc=0,scale=1,size=(halfcells,int(informativegenes*ngenes))))

#Simulate Class 1 cells
y2=pd.DataFrame(np.random.normal(loc=effectsize,scale=1,size=(halfcells,int(informativegenes*ngenes))))

#Noisy genes
y3=pd.DataFrame(np.random.normal(loc=0,scale=1,size=(ncells,int(ngenes*(1.0-informativegenes)))))

#Concatenate to form cells x genes expression matrix
Y=pd.concat([y1,y2]).reset_index(drop=True)
Y=pd.concat([Y,y3],axis=1)
Y.columns=range(np.shape(Y)[1])
Y.head()
y3
X=pd.DataFrame()
classvec=[0]*halfcells
classvec.extend([1]*halfcells)
X['class']=classvec

X_noise=X.copy()

#X_noise.ix[0:tenthcells]=1   Here pandas abandoned .ix attribut of pandas so we need to change .ix to .iloc
X_noise.iloc[0:tenthcells]=1   # Here we use wrong labels of 1/10 cells
X_noise
#Kernel Density estimates of the two distributions across all genes
sns.kdeplot(Y.iloc[range(halfcells),:].values.flatten())
sns.kdeplot(Y.iloc[halfcells+1:,:].values.flatten())  # This is the true label (didn't use the psuedo wrong labels we've generated in the last step)
plt.xlabel('Expression Value')
plt.ylabel('Density Estimate')
#plot a single gene as a function of class membership
plt.scatter(X['class'],Y[0])
plt.xlabel('class label')
plt.ylabel('Expression')<jupyter_output><empty_output><jupyter_text>Here uses a simple linear regression to approximate the label (class)<jupyter_code>#Fit regression model
lm=sklearn.linear_model.Ridge()
lm.fit(X,Y)
B=pd.DataFrame(lm.coef_)
B
#coefficients for informative subset of genes
plt.plot(B[0])
plt.axvline(ngenes*informativegenes,c='red')
plt.axhline(effectsize,c='gray')
plt.xlabel('Gene Index')
plt.ylabel('Coefficient')
#Adjust covariates based on fitted coefficients
X_adjust=bayes_cov_col(Y,X,['class'],lm)
X_adjust
#Plot original and adjusted covariates
plt.plot(X,label='Original')
plt.plot(X_adjust,label='Adjusted')
plt.legend(loc='lower right')
plt.xlabel('Cell Index')
plt.ylabel('Class Probability')
#Fit regression model on noisy covariates (note how the first ten percent of cells are mostly reclassified as class 0)
lm_noise=sklearn.linear_model.Ridge()
lm_noise.fit(X_noise,Y)
B_noise=pd.DataFrame(lm_noise.coef_)
X_adjust_noise=bayes_cov_col(Y,X_noise,['class'],lm_noise)
plt.plot(X_noise,label='Original Noisy')
plt.plot(X_adjust_noise,label='Adjusted')
plt.legend(loc='lower right')
plt.xlabel('Cell Index')
plt.ylabel('Class Probability')
lm_noise_adjust=sklearn.linear_model.Ridge()
lm_noise_adjust.fit(X_adjust_noise,Y)
B_noise_adjust=pd.DataFrame(lm_noise_adjust.coef_)
lm_noise_adjust=sklearn.linear_model.Ridge()
lm_noise_adjust.fit(X_adjust_noise,Y)
B_noise_adjust=pd.DataFrame(lm_noise_adjust.coef_)
#plot coefficients for informative subset of genes
plt.plot(B_noise[0],label='Original')
plt.plot(B_noise_adjust[0],label='Adjusted')
plt.axvline(ngenes*informativegenes,c='red')
plt.axhline(effectsize,c='gray')
plt.legend(loc='upper right')
plt.xlabel('Gene Index')
plt.ylabel('Coefficient')<jupyter_output><empty_output><jupyter_text>Apply this to our real perturb data  Prepare data in the correct formatThe bayesian function expect the expression matrix in such format    @Y    = Expression matrix, cells x x genes, expecting pandas dataframe    @X    = Covariate matrix, cells x covariates, expecting pandas dataframe    @cols = The subset of columns that the EM should be performed over, expecting list    @lm   = linear model object  Prepare expression matrix Y The original sparce matrix is cells by genes, so we don't need to transpose it.Just convert sparse matrix to a pandas DataFrame with header and index.<jupyter_code>sc_object.X
import copy
# Convert sparse matrix to pandas dataframe and add index and column names to it
Y = pd.DataFrame.sparse.from_spmatrix(copy.deepcopy(sc_object.X))
Y.index = copy.deepcopy(sc_object.obs.index)
Y.columns = copy.deepcopy(sc_object.var.index)
Y.head()
!pwd
Y.to_csv("./stupid.csv")
import numpy as np
ynp = np.array(Y)
pd.DataFrame(ynp.var(axis = 0))
Y_loaded = pd.read_csv("./stupid.csv")
Y.var(axis = 0)<jupyter_output><empty_output><jupyter_text>Prepare covariate matrix Just use the meta file as this covariate matrix, of course we won't use all of the variables in it. The used variables will be specified in the next parameterAlso we need to convert the label information into a one-hot encoding matrix<jupyter_code>X = sc_object.obs.copy()
X
# Get one-hot encoding
y = pd.get_dummies(X.labels)
cols = list(y.columns)
X = pd.merge(X, y, left_index = True, right_index= True)
X['batch'].value_counts()
X['Infected'].value_counts()
"""
map_class_labels = pd.DataFrame(X['labels'].value_counts())#.to_dict()
map_class_labels["labels"] = np.ones(len(X['labels'].value_counts()))#np.arange(len(X['labels'].value_counts()))
map_class_labels = map_class_labels.to_dict()
map_class_labels = map_class_labels['labels']
map_class_labels['nan'] = 0
"""
# Convert the string levels into integer levels
map_class_batch = {"SARS2" : 0, "MOCK" : 1}
map_class_infected = {"Bystander" : 0, "Mock" : 1, "Infected": 3}

X['batch'] = X['batch'].map(map_class_batch)
X['Infected'] = X['Infected'].map(map_class_infected)
#X['labels'] = X['labels'].map(map_class_labels)
X['batch'] = X['batch'].astype('int64')#.astype('category', ordered=True)
X['Infected'] = X['Infected'].astype('int64')#.astype('category', ordered=True)
#X['labels'] = X['labels'].astype('int64')#.astype('category', ordered=True)
X<jupyter_output><empty_output><jupyter_text>Prepare colsThe subset of columns that the EM should be performed over, a list.<jupyter_code>cols = cols + ["batch", "Infected",]<jupyter_output><empty_output><jupyter_text>Prepare the linear model object<jupyter_code>X.loc[:,cols]
enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=modalpha,max_iter=10000)

#Fit regression model
import sklearn
lm=sklearn.linear_model.Ridge()
lm.fit(X.loc[:,cols],Y)
B=pd.DataFrame(lm.coef_)<jupyter_output><empty_output><jupyter_text>Visualize the coefficient<jupyter_code>#coefficients for informative subset of genes
plt.plot(B[0])
#plt.axvline(ngenes*informativegenes,c='red')
#plt.axhline(effectsize,c='gray')
plt.xlabel('Gene Index')
plt.ylabel('Coefficient')
plt.show()
Y.shape<jupyter_output><empty_output><jupyter_text>Run model<jupyter_code>Be,X_adjust,RES_out = run_model_bycol(Y, X.loc[:,cols] , cols)<jupyter_output>/usr/local/lib/python3.7/dist-packages/pandas/core/indexing.py:670: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  iloc._setitem_with_indexer(indexer, value)
/usr/local/lib/python3.7/dist-packages/pandas/core/indexing.py:670: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  iloc._setitem_with_indexer(indexer, value)
/usr/local/lib/python3.7/dist-packages/pandas/core/indexing.py:670: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  iloc._setit[...]<jupyter_text>Get the original $\beta$<jupyter_code>enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=0.005,max_iter=10000)
enet.fit(X.loc[:,cols],Y)
Be_original=pd.DataFrame(enet.coef_)
Be_original.columns=X.loc[:,cols].columns
Be_original.index=Y.columns
Be_original.loc['ACE2',:].plot()
Be_original['ACE2'].plot(title = "The extent a gene perturbed by ACE2", xlabel = "Genes")<jupyter_output><empty_output><jupyter_text>Visualize  Compare the initial covariates and the adjusted covariates<jupyter_code>X_adjust['ACE2']
original_cells = X[X['ACE2']>0].index
other_cells = X[X['ACE2']==0].index
reorder_indices = original_cells.append(other_cells)
First_one_thouthand = X.iloc[0:1000,:].index
plt.plot(X.loc[reorder_indices,'ACE2'],label='Original Covariates')
plt.plot(X_adjust.loc[reorder_indices,'ACE2'], label='Adjusted')
plt.legend(loc='lower right')
plt.xlabel('Cell Index')
plt.ylabel('Class Probability')
plt.title("ACE2 targeted probability")
plt.xticks([])
plt.show()<jupyter_output><empty_output><jupyter_text>Compare the original $\beta$ and the refitted $\beta$ $\beta$ is a Gene by Covariates matrix (Different from the paper), So we can inteprete the entry of $\beta$ as the effectness of each perturbation on gene<jupyter_code># This is adjusted Beta
Be<jupyter_output><empty_output><jupyter_text>For example, for ACE2 perturbation<jupyter_code>Be['ACE2'].plot()
plt.plot(Be_original['ACE2'],label='Original', alpha = 0.5)
plt.plot(Be['ACE2'], label='Adjusted', alpha = 0.5)
plt.legend(loc='upper left')
plt.xlabel('Gene Index')
plt.ylabel('Coef')
plt.title("ACE2 Effectness")
plt.xticks([])
plt.show()<jupyter_output><empty_output><jupyter_text>For ACE2 gene<jupyter_code>Be[Be['ACE2'] > 0.05].index
Be.loc['ACE2','ACE2']
X_adjust=bayes_cov_col(Y, X.loc[:,cols] , cols, lm)
X_adjust
X_adjust
X_adjust['ACE2'].value_counts()
#X['ACE2'].value_counts()
X['ACE2'].value_counts()<jupyter_output><empty_output><jupyter_text>**Based on the results in the above two cells, we see that the adjusted label just gave the 100% accurate prediction, which make me feel the current label is the version that has already corrected by the same algorithm**  EM iteration<jupyter_code>def run_model_bycol(Y,X,EM_cols=None,modalpha=0.005,verbose=0):   
    """
    @Y    = Expression matrix, expecting pandas dataframe, cells x genes
    @X    = Covariate matrix, expecting pandas dataframe, cells x covariates
    @EM_cols = The subset of columns that the EM should be performed over, list
    """

    enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=modalpha,max_iter=10000)
    enet.fit(X,Y)
    if verbose==1:
        print(enet.score(X,Y))

    Be=pd.DataFrame(enet.coef_)
    Be.columns=X.columns
    Be.index=Y.columns

    Yhat=pd.DataFrame(enet.predict(X))
    Yhat.index=Y.index
    Yhat.columns=Y.columns

    if EM_cols is not None:

        X_adjust=bayes_cov_col(Y,X,EM_cols,enet)

        #print('done with EM')

        #refit model

        enet=sklearn.linear_model.ElasticNet(precompute=True,l1_ratio=0.5,alpha=0.0004,max_iter=10000)
     
        enet.fit(X_adjust,Y)
        Yhat=pd.DataFrame(enet.predict(X_adjust))
        Yhat.index=Y.index
        Yhat.columns=Y.columns

        if verbose==1:
            print(enet.score(X_adjust,Y))

        Be=pd.DataFrame(enet.coef_)
        Be.columns=X_adjust.columns
        Be.index=Y.columns
    else:
        X_adjust=X.copy()

    RES_out=Y.subtract(Yhat)  

    return(Be,X_adjust,RES_out)
X_adjust['labels'].value_counts()
#X['labels'].value_counts()
pd.DataFrame(X.loc[:,'labels'] - X_adjust.loc[:,'labels']).value_counts()
#Plot original and adjusted covariates
plt.plot(X.loc[:,'labels'],label='Original')
plt.plot(X_adjust.loc[:,'labels'],label='Adjusted')
plt.legend(loc='lower right')
plt.xlabel('Cell Index')
plt.ylabel('Class Probability')

B[1]
from numpy import array
from scipy.sparse import coo_matrix
row  = array([0, 0, 1, 3, 1, 0, 0])
col  = array([0, 2, 1, 3, 1, 0, 0])
data = array([1, 1, 1, 1, 1, 1, 1])
A = coo_matrix((data, (row, col)), shape=(4, 4)).tocsc()
A<jupyter_output><empty_output><jupyter_text>Expression matrix visulization<jupyter_code>ACE2_indices = X[X['ACE2'] == 1].index
sum(Y.loc[ACE2_indices,'ACE2'])/len(ACE2_indices)
sum(Y['ACE2'].drop(ACE2_indices, axis = 0))/(len(Y) - len(ACE2_indices))<jupyter_output><empty_output><jupyter_text>Examples  pandas.DataFrame.var()if axis = 0 calculate the variance accross one columnif axis = 1 calculate the variance accross one row<jupyter_code>import pandas as pd

data = [[10, 10], [10, 15], [10, 20]]

df = pd.DataFrame(data)

print(df[0].var(axis = 0))
print(df.div([2,2]))<jupyter_output>0     1
0  5.0   5.0
1  5.0   7.5
2  5.0  10.0
{'fraction_non_alphanumeric': 0.13048966496607584, 'fraction_numerical': 0.017093567559038552, 'mean_word_length': 4.848661950169179, 'pattern_counts': {'":': 1, '<': 95, '<?xml version=': 0, '>': 100, 'https://': 4, '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': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '413996', 'license': 'MIT', 'n_tokens_mistral': 7173, 'n_tokens_neox': 6726, 'n_words': 1413, 'path': 'Perturb_dataset.ipynb', 'repo_name': 'JiaxinLi-lipluszn/Organized-Notebook', 'size': '19013'}