<jupyter_start><jupyter_text>HDACInvestigate the differences in transcriptome variability between drugs that target the same system<jupyter_code>import scanpy as sc
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import itertools
import scipy.stats as stats
import sys
sys.path.append('/home/ssm-user/Github/scrna-parameter-estimation/dist/memento-0.0.6-py3.8.egg')
sys.path.append('/home/ssm-user/Github/misc-seq/miscseq/')
import encode
import memento
import gseapy as gp
from gseapy.plot import gseaplot

pd.set_option('display.max_columns', None)
data_path = '/data_volume/memento/sciplex/'
ct = 'A549'
adata = sc.read(data_path + 'h5ad/{}_HDAC.h5ad'.format(ct))<jupyter_output><empty_output><jupyter_text>Setup memento<jupyter_code>adata.obs['q'] = 0.05
memento.setup_memento(adata, q_column='q', filter_mean_thresh=0.07)<jupyter_output>Version 0.0.6<jupyter_text>One sample test<jupyter_code>subset = adata[(adata.obs.target == 'HDAC') & (adata.obs.dose == 10000)].copy().copy()

memento.create_groups(subset, label_columns=['product_name'])
memento.compute_1d_moments(subset, min_perc_group=.9)
memento.ht_1d_moments(
    subset, 
    formula_like='1',
    cov_column='Inter', 
    num_boot=10000, 
    verbose=1,
    num_cpus=90)

# subset.write(data_path + 'hdac/{}_vs_vehicle.h5ad'.format(drug))
df = memento.get_1d_ht_result(subset)
df['dv_fdr'] = memento.util._fdrcorrect(df['dv_pval'])
df['de_fdr'] = memento.util._fdrcorrect(df['de_pval'])<jupyter_output><empty_output><jupyter_text>GSEA<jupyter_code>import gseapy as gp
from gseapy.plot import gseaplot

pre_res = gp.prerank(
    rnk=df[['gene','dv_coef']].sort_values('dv_coef'), 
    gene_sets='GO_Biological_Process_2018',
    processes=14,
    permutation_num=100, # reduce number to speed up testing
    outdir=None, seed=6)<jupyter_output><empty_output><jupyter_text>Comparison to control, using dosage as treatmentUse different drugs as strata<jupyter_code>drug_counts = adata.obs.query('~vehicle').product_name.value_counts()
drug_list = drug_counts[drug_counts > 0].index.tolist()
drug_list = [
    'Pracinostat (SB939)',
    'Quisinostat (JNJ-26481585) 2HCl',
    'Abexinostat (PCI-24781)',
    'Mocetinostat (MGCD0103)',
    'Divalproex Sodium',
    'M344',
    'Resminostat',
    'Panobinostat (LBH589)'
]
drug = drug_list[0]
# subset = adata[((adata.obs.target == target))].copy().copy()
subset = adata.copy().copy()
# subset.obs['dose_indicator'] = (subset.obs.product_name == drug).astype(int)
subset.obs['dose_indicator'] = subset.obs.dose.apply(lambda x: {0:0.0 ,10:1.0, 100:2.0, 1000:3.0, 10000:4.0}[x])
subset.obs['dose_level'] = 'dose_' + subset.obs['dose'].astype(str)

memento.create_groups(subset, label_columns=['dose_indicator', 'product_name'])
memento.compute_1d_moments(subset, min_perc_group=.9)
memento.ht_1d_moments(
    subset, 
    formula_like='1 + dose_indicator + product_name',
    treatment_col='dose_indicator', 
    num_boot=20000, 
    verbose=1,
    num_cpus=92,
    resampling='permutation',
    approx=False)
df = memento.get_1d_ht_result(subset)
df['dv_fdr'] = memento.util._fdrcorrect(df['dv_pval'])
df['de_fdr'] = memento.util._fdrcorrect(df['de_pval'])
plt.scatter(df['de_coef'], df['dv_coef'], s=1)
df.query('gene == "STAT3"')
pre_res = gp.prerank(rnk=df.sort_values('de_coef', ascending=False)[['gene', 'de_coef']], 
                     gene_sets='WikiPathways_2019_Human',
                     processes=4,
                     permutation_num=100, # reduce number to speed up testing
                     outdir=None, seed=6)
pre_res.res2d.sort_index().query('fdr < 0.1')

enr = gp.enrichr(gene_list=df.query('de_fdr < 0.01 & de_coef > 0').sort_values('dv_coef', ascending=False).gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
df.query('de_fdr < 0.01 & de_coef > 0 & dv_fdr < 0.1 & dv_coef < 0').sort_values('dv_coef', ascending=False).gene.shape
enr = gp.enrichr(gene_list=df.query('de_fdr < 0.1 & de_coef > 0').sort_values('dv_coef', ascending=False).gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.sort_values('Adjusted P-value').head(10)
enr = gp.enrichr(gene_list=df.query('dv_fdr < 0.1').sort_values('dv_coef', ascending=False).gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.sort_values('Adjusted P-value').head(10)
enr = gp.enrichr(gene_list=df.query('dv_fdr < 0.1 & dv_coef > 0').sort_values('dv_coef', ascending=False).gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.sort_values('Adjusted P-value').head(10)
enr = gp.enrichr(gene_list=df.query('dv_fdr < 0.1 & dv_coef < 0').sort_values('dv_coef', ascending=False).gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.sort_values('Adjusted P-value').head(10)<jupyter_output><empty_output><jupyter_text>Multiclass comparison<jupyter_code>subset = adata[~adata.obs.vehicle].copy().copy()
subset.obs['dose_level'] = 'dose_' + subset.obs['dose'].astype(str)

memento.create_groups(subset, label_columns=['product_name', 'dose_level'])
memento.compute_1d_moments(subset, min_perc_group=.9)
memento.ht_1d_moments(
    subset, 
    formula_like='1 + product_name + dose_level',
    treatment_col='product_name', 
    num_boot=10000, 
    verbose=1,
    num_cpus=94,
    resampling='permutation',
    approx=False)<jupyter_output><empty_output><jupyter_text>Drug-drug comparisonsInvestigate if including the dosage changes anything<jupyter_code>drug_counts = adata.obs.query('~vehicle').product_name.value_counts()
drug_list = drug_counts[drug_counts > 0].index.tolist()
drug_list = drug_counts[drug_counts > 0].index.tolist()
drug_list = [
    'Pracinostat (SB939)',
    'Quisinostat (JNJ-26481585) 2HCl',
    'Abexinostat (PCI-24781)',
    'Mocetinostat (MGCD0103)',
    'Divalproex Sodium',
    'M344',
    'Resminostat',
    'Panobinostat (LBH589)'
]
len(list(itertools.combinations(drug_list,2)))
# for drug1, drug2 in itertools.combinations(drug_list,2):
    
#     print(drug1, drug2)

#     subset = adata[adata.obs.product_name.isin([drug1, drug2]) & (adata.obs.dose == 10000)].copy().copy()
#     subset.obs['is_drug1'] = (subset.obs.product_name==drug1).astype(int)
    
#     memento.create_groups(subset, label_columns=['is_drug1'])
#     memento.compute_1d_moments(subset, min_perc_group=.9)
#     memento.ht_1d_moments(
#         subset, 
#         formula_like='1 + is_drug1',
#         cov_column='is_drug1', 
#         num_boot=10000, 
#         verbose=0,
#         num_cpus=94)
    
#     subset.write(data_path + 'hdac/{}_vs_{}_10000.h5ad'.format(drug1, drug2))
for drug1, drug2 in itertools.combinations(drug_list,2):
    
    print(drug1, drug2)

    subset = adata[adata.obs.product_name.isin([drug1, drug2])].copy().copy()
    subset.obs['is_drug1'] = (subset.obs.product_name==drug1).astype(int)
    subset.obs['dose_level'] = 'dose_' + subset.obs['dose'].astype(str)
    
    memento.create_groups(subset, label_columns=['is_drug1', 'dose_level'])
    memento.compute_1d_moments(subset, min_perc_group=.9)
    memento.ht_1d_moments(
        subset, 
        formula_like='1 + is_drug1 + dose_level',
        treatment_col='is_drug1', 
        num_boot=10000, 
        verbose=0,
        num_cpus=94,
        resampling='permutation',
        approx=False)
    
    break
    
#     subset.write(data_path + 'hdac/{}_vs_{}_stratified.h5ad'.format(drug1, drug2))
df = memento.get_1d_ht_result(subset)
df.head(5)
df['dv_fdr'] = memento.util._fdrcorrect(df['dv_pval'])
df['de_fdr'] = memento.util._fdrcorrect(df['de_pval'])

plt.hist(df.dv_pval)
df.shape
np.log(5)
df.query('de_fdr < 0.01 & dv_fdr < 0.1')
plt.scatter(df['de_coef'], df['dv_coef'], s=1)
plt.scatter(df.query('dv_fdr < 0.1')['de_coef'], df.query('dv_fdr < 0.1')['dv_coef'], s=10)<jupyter_output><empty_output><jupyter_text>Combine P-values<jupyter_code>dv_pvals = pd.DataFrame()
start = True
for drug1, drug2 in itertools.combinations(drug_list,2):
    
    subset = sc.read(data_path + 'hdac/{}_vs_{}_stratified.h5ad'.format(drug1, drug2))
    df = memento.get_1d_ht_result(subset)
    if start:
        
        dv_pvals['gene'] = df['gene']
        dv_pvals['{}_{}'.format(drug1, drug2)] = df['dv_pval']
        start = False
    else:
        dv_pvals = dv_pvals.merge(df[['gene', 'dv_pval']].rename(columns={'dv_pval':'{}_{}'.format(drug1, drug2)}), how='outer')
        
# dv_pvals = pd.DataFrame()
# start = True
# for drug1, drug2 in itertools.combinations(drug_list,2):
    
#     subset = sc.read(data_path + 'hdac/{}_vs_{}_10000.h5ad'.format(drug1, drug2))
#     df = memento.get_1d_ht_result(subset)
#     if start:
        
#         dv_pvals['gene'] = df['gene']
#         dv_pvals['{}_{}'.format(drug1, drug2)] = df['de_pval']
#         start = False
#     else:
#         dv_pvals = dv_pvals.merge(df[['gene', 'de_pval']].rename(columns={'de_pval':'{}_{}'.format(drug1, drug2)}), how='outer')
combined_pvals = []
for g in dv_pvals.gene:
    pvalues = dv_pvals.query('gene == "{}"'.format(g)).iloc[0, 1:].values.astype(float)
    chisq, pv = stats.combine_pvalues(pvalues[np.isfinite(pvalues)])
    combined_pvals.append((g, chisq, pv))
combined_pvals = pd.DataFrame(combined_pvals, columns=['gene', 'chisq', 'pv'])
combined_pvals['fdr'] = memento.util._fdrcorrect(combined_pvals['pv'])
g = 'IDI1;ACLY;INSIG1;DHCR24;FDFT1'.split(';')<jupyter_output><empty_output><jupyter_text>GSEA + EnrichrIDI1;ACLY;INSIG1;DHCR24;FDFT1<jupyter_code># for gene in combined_pvals.query('fdr < 0.01').sort_values('fdr').gene:
#     print(gene)
import gseapy as gp
from gseapy.plot import gseaplot

# pre_res = gp.prerank(rnk=combined_pvals.sort_values('chisq', ascending=False)[['gene', 'chisq']], gene_sets='GO_Biological_Process_2021',
#                      processes=90,
#                      permutation_num=100, # reduce number to speed up testing
#                      outdir=None,seed=6)
# pre_res.res2d.sort_index().head(30)

import gseapy as gp
from gseapy.plot import gseaplot

enr = gp.enrichr(gene_list=combined_pvals.query('fdr < 0.01').sort_values('fdr').gene,
                 gene_sets=['GO_Biological_Process_2021'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.sort_values('Adjusted P-value').head(30)
enr.results.sort_values('Adjusted P-value').head(30)
enr = gp.enrichr(gene_list=combined_pvals.query('fdr < 0.01').sort_values('fdr').gene,
                 gene_sets=['WikiPathways_2019_Human'],
                 organism='Human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None,
                 cutoff=0.5 # test dataset, use lower value from range(0,1)
                )
enr.results.head(50)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15169660678642716, 'fraction_numerical': 0.03351630073186959, 'mean_word_length': 3.8003992015968064, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 43, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '161669', 'license': 'MIT', 'n_tokens_mistral': 4915, 'n_tokens_neox': 4505, 'n_words': 687, 'path': 'analysis/sciplex/hdac.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '12024'}