<jupyter_start><jupyter_text>Linear Covariance TestingUse linear mixed effect models to test the effect of the guide on the covariance of IRF4 and BATF and target genes.  Import<jupyter_code>from IPython.core.display import display, HTML
import warnings
warnings.filterwarnings('ignore')
display(HTML("<style>.container { width:100% !important; }</style>"))
%matplotlib inline
repo_path = '/Users/mincheolkim/Github/'
data_path = '/Users/mincheolkim/Documents/'
import sys
sys.path.append(repo_path + 'scVI')
sys.path.append(repo_path + 'scVI-extensions')
import statsmodels.api as sm
import statsmodels.formula.api as smf
import imp
import numpy as np
import pandas as pd
import scvi_extensions.dataset.supervised_data_loader as sdl
import scvi_extensions.dataset.cropseq as cs
import scvi_extensions.inference.supervised_variational_inference as svi
import scvi_extensions.hypothesis_testing.mean as mn
import scvi_extensions.hypothesis_testing.variance as vr
import scvi_extensions.dataset.label_data_loader as ldl
import scanpy.api as sc<jupyter_output><empty_output><jupyter_text>Load a dataset<jupyter_code>h5_filename = data_path + 'raw_gene_bc_matrices_h5.h5'
metadata_filename = data_path + 'nsnp20.raw.sng.km_vb1_default.norm.meta.txt'
imp.reload(cs)
# Load the dataset
gene_dataset = cs.CropseqDataset(
    filename=h5_filename,
    metadata_filename=metadata_filename,
    batch='wells',
    use_labels='gene',
    save_path='')<jupyter_output>Preprocessing CROP-seq dataset
Number of cells kept after filtering with metadata: 283634
Number of cells kept after removing all zero cells: 283634
Finished preprocessing CROP-seq dataset<jupyter_text>Read into scanpy<jupyter_code>embedded_adata = sc.read('/Users/mincheolkim/Documents/nsnp20.raw.sng.km_vb1_default.pc60.norm.h5ad')
pd.Series(gene_dataset.ko_gene.reshape(-1)).value_counts()<jupyter_output><empty_output><jupyter_text>Create a dataframe with relevant cells<jupyter_code>ko_genes = ['NO_GUIDE', 'BATF', 'IRF4', 'JUNB', 'STAT1']
genes_of_interest = [
    'BATF', 
    'IRF4', 
    'JUNB', 
    'RORC',
    'BCL6', 'MAF', 'IL10', 'IL17A', 'IL17B', 'CYCS', 'IL2', 'CMIP', 'IFNG', 'IL4', 'IL5', 'STAT1',
    'IL21', 'IL23R']#,'IL1RB1',]
goi_indices = [np.where(gene_dataset.gene_names == gene)[0][0] for gene in genes_of_interest]
gene_dataset.ko_gene_lookup
dfs=[]
for ko_gene in ko_genes:
    indices = (gene_dataset.labels == np.where(gene_dataset.ko_gene_lookup == ko_gene)[0][0]).reshape(-1)
    print(ko_gene, indices.sum())
    expr = gene_dataset.X[indices, :].astype(float)
    umi_counts = expr.sum(axis=1)
    df = pd.DataFrame(expr[:, goi_indices].todense()/umi_counts, columns=genes_of_interest)
    
    for gene in genes_of_interest:
        df[gene + '_std'] = (df[gene] - df[gene].mean())/np.sqrt(df[gene].var())
    df['ko_gene'] = ko_gene
    df['donor'] = pd.Series(gene_dataset.donor_batches[indices].reshape(-1)).astype(int)
    df['louvain'] = gene_dataset.louvain[indices].reshape(-1)
    df['well'] = gene_dataset.wells[indices].reshape(-1)
    dfs.append(df)
df = pd.get_dummies(pd.concat(dfs), columns=['ko_gene'])<jupyter_output>NO_GUIDE 87336
BATF 915
IRF4 961
JUNB 856
STAT1 1458<jupyter_text>Compute covariances of interest and decide the group<jupyter_code>def compute_point_cov(s1, s2):
   return df[s1 + '_std']*df[s2 + '_std']
groups = 'donor'<jupyter_output><empty_output><jupyter_text>Get some general sense of how sparse these genes areIRF4 isnt detected for most cells, so I'm going to test BATF using IRF4 KO cells.<jupyter_code>df['cov_BATF_RORC'] = compute_point_cov('BATF', 'RORC')
df['cov_IRF4_RORC'] = compute_point_cov('IRF4', 'RORC')
df['cov_BATF_IFNG'] = compute_point_cov('BATF', 'IFNG')
df['cov_IRF4_IFNG'] = compute_point_cov('IRF4', 'IFNG')
df['cov_BATF_RORC'] = compute_point_cov('BATF', 'RORC')
df['cov_IRF4_RORC'] = compute_point_cov('IRF4', 'RORC')
df['cov_STAT1_IFNG'] = compute_point_cov('STAT1', 'IFNG')
df['cov_IRF4_IFNG'] = compute_point_cov('IRF4', 'IFNG')<jupyter_output><empty_output><jupyter_text>Test differential BATF/RORC covariance between IRF4 KO cells and NO_GUIDE cells<jupyter_code>print(smf.mixedlm(
    'cov_STAT1_IFNG ~ ko_gene_NO_GUIDE',
    df.query('ko_gene_STAT1 > 0 | ko_gene_NO_GUIDE > 0'),
    groups=df.query('ko_gene_STAT1 > 0 | ko_gene_NO_GUIDE > 0')[groups]).fit(maxiter=2000, method='nm').summary())<jupyter_output>Mixed Linear Model Regression Results
============================================================
Model:            MixedLM Dependent Variable: cov_STAT1_IFNG
No. Observations: 88794   Method:             REML          
No. Groups:       9       Scale:              0.6086        
Min. group size:  7571    Likelihood:         -103952.9886  
Max. group size:  11927   Converged:          Yes           
Mean group size:  9866.0                                    
------------------------------------------------------------
                  Coef.  Std.Err.   z    P>|z| [0.025 0.975]
------------------------------------------------------------
Intercept         -0.037    0.020 -1.806 0.071 -0.077  0.003
ko_gene_NO_GUIDE   0.021    0.021  1.041 0.298 -0.019  0.062
Group Var          0.000    0.000                           
============================================================<jupyter_text>Scratch<jupyter_code>temp = (gene_dataset.X > 0).sum(axis=0)
(temp > gene_dataset.X.shape[0]*0.05).sum()
min(1, 300/400)
700/250000<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16789260757631483, 'fraction_numerical': 0.042111070246414126, 'mean_word_length': 4.2398843930635834, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 39, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '201844', 'license': 'MIT', 'n_tokens_mistral': 2141, 'n_tokens_neox': 1968, 'n_words': 377, 'path': 'notebooks/linear_covariance_testing.ipynb', 'repo_name': 'mincheoly/scVI-extensions', 'size': '5438'}