<jupyter_start><jupyter_text>Comparative Analysis with MELD<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scprep
import meld
import sklearn
import pickle
import cna
import time

# making sure plots & clusters are reproducible
np.random.seed(42)<jupyter_output><empty_output><jupyter_text>Test on Sepsis<jupyter_code># Import sepsis expression data
d = cna.read('/data/srlab/lrumker/MCSC_Project/sepsis_data/pbmc.h5ad')
sepsis_data = d.to_df() # Cells x genes
sepsis_metadata = d.obs
## Define phenotype of interest as ANY sepsis to match CNA
any_sepsis = np.zeros(sepsis_metadata.shape[0])
any_sepsis[np.where(sepsis_metadata['pheno']=='Bac-SEP')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='ICU-SEP')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='URO')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='Int-URO')[0]] = 1

sepsis_metadata['AnySepsis'] = any_sepsis
# PCA for dimensionality reduction with default parameters (100 components)
data_pca = scprep.reduce.pca(sepsis_data)
# Used default MELD parameters
time_a = time.time()
meld_op = meld.MELD()
sample_densities = meld_op.fit_transform(data_pca, sample_labels=sepsis_metadata['patient'])
time_b = time.time()
print(time_b-time_a)
# Save sample densities object
pickle.dump(sample_densities, open( "MELD_sample_densities_sepsis.p", "wb" ) )
# Load saved object
sample_densities = pickle.load( open( "MELD_sample_densities_sepsis.p", "rb" ) )
# Apply row-wise L1 normalization
sample_likelihoods = sample_densities/sample_densities.sum(axis=1)[:,None]
# Identify case sample colums
sepsis_samples = np.unique(sepsis_metadata['patient'][np.where(sepsis_metadata['AnySepsis']==1)[0]])
sepsis_metadata['likelihood_sepsis'] = sample_likelihoods[sepsis_samples].mean(axis=1).values<jupyter_output><empty_output><jupyter_text>Compare per-cell values from MELD and CNA<jupyter_code># Load CNA results
cna_res = pd.read_csv("/data/srlab/lrumker/MCSC_Project/mcsc_scratch/sepsisres.csv")
# Correlation between per-cell values from MELD and CNA
np.corrcoef(cna_res['ncorrs'], sepsis_metadata['likelihood_sepsis'])[0,1]<jupyter_output><empty_output><jupyter_text>Analysis of statistical properties of MELD<jupyter_code># read sepsis data
import cna
d = cna.read('/data/srlab/lrumker/MCSC_Project/sepsis_data/pbmc.h5ad')

sepsis_data = d.to_df() # Cells x genes
sepsis_metadata = d.obs

## Define phenotype of interest as ANY sepsis
any_sepsis = np.zeros(sepsis_metadata.shape[0])
any_sepsis[np.where(sepsis_metadata['pheno']=='Bac-SEP')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='ICU-SEP')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='URO')[0]] = 1
any_sepsis[np.where(sepsis_metadata['pheno']=='Int-URO')[0]] = 1
sepsis_metadata['AnySepsis'] = any_sepsis
# initiatlize sampleXmeta for sepsis
sampleXmeta = sepsis_metadata[['patient', 'AnySepsis']].drop_duplicates().set_index('patient', drop=True)
true_pheno = sampleXmeta.AnySepsis.values.astype(bool)
# obtain memory requirement for storing NAM
cna.tl.nam(d)
d.uns['NAM.T'].values.nbytes / np.power(2, 20)
res = cna.tl.association(d, true_pheno)
# read MELD sample densities
# path = '/data/srlab/lrumker/MCSC_Project/sepsis_data/sepsis_sample_densities.p'
path = "MELD_sample_densities_sepsis.p"
sample_densities = pickle.load(open(path, 'rb'))
sample_densities /= sample_densities.sum(axis=1)[:,None] # This is termed sample "likelihoods" in the publication
color_scale_threshold = 0.0025
# perform non-null analysis
val = sample_densities[sampleXmeta.index[true_pheno]].mean(axis=1) - sample_densities.mean(axis=1)
plt.scatter(d.obsm['X_tsne'][:,0], d.obsm['X_tsne'][:,1],
           alpha=0.5, c = val, cmap = "seismic", s=0.2, vmin=-color_scale_threshold, vmax=color_scale_threshold)
# create null phenotype visualizations
np.random.seed(0)
for i in range(5):
    ix = np.argsort(np.random.randn(len(true_pheno)))
    null = true_pheno[ix]
    nullval = sample_densities[sampleXmeta.index[null]].mean(axis=1) - sample_densities.mean(axis=1)
    plt.scatter(d.obsm['X_tsne'][:,0], d.obsm['X_tsne'][:,1],
               alpha=0.5, c=nullval, cmap = "seismic", s=0.2, vmin=-color_scale_threshold, vmax=color_scale_threshold)
    plt.show()
    nullres = cna.tl.association(d, null, local_test=False)
    print('CNA p =', nullres.p)
# compute null distribution for MELD scores
T = 500
null_scores = np.zeros((T, len(sample_densities)))
for i in range(500):
    print('.', end='')
    ix = np.argsort(np.random.randn(len(true_pheno)))
    null = true_pheno[ix]
    null_scores[i] = sample_densities[sampleXmeta.index[null]].mean(axis=1) - sample_densities.mean(axis=1)
null_scores = null_scores.T
# compute FDRs
fdrs = []
thresholds = np.arange(50, 100, 1)
for thresh in np.percentile(np.abs(val), thresholds):
    print('.', end='')
    discoveries = (np.abs(val) >= thresh).sum()
    Efd = (np.abs(null_scores) >= thresh).sum(axis=0).mean()
    fdrs.append(Efd/discoveries)
CNA_thresholds = np.array([(np.abs(res.ncorrs) <= t).mean()*100 for t in res.fdrs.threshold])
CNA_fdrs = res.fdrs.fdr.values
plt.plot(thresholds, fdrs, label='MELD')
plt.plot(CNA_thresholds, CNA_fdrs, label='CNA')
plt.xlabel('Per-Cell Value Threshold (%ile)')
plt.ylim(0,0.2)
plt.axhline(y=0.05, color='k', linestyle='--')
plt.ylabel('FDR')
plot_thresholds = [0.001, 0.01, 0.05, 0.1, 0.15, 0.2]
total_cells = d.obs.shape[0]
MELD_cells_passing = []
for t in plot_thresholds:
    if np.min(fdrs)>t:
        MELD_cells_passing.append(0)
    else:
        i_first_below = np.min(np.where(np.array(fdrs)<t)[0])
        numcells = total_cells*(1-thresholds[i_first_below]/100)
        MELD_cells_passing.append(numcells)
MELD_cells_passing = np.ceil(MELD_cells_passing).astype(int)        
total_cells = d.obs.shape[0]
CNA_cells_passing = []
for t in plot_thresholds:
    if np.min(CNA_fdrs)>t:
        CNA_cells_passing.append(0)
    else:
        i_first_below = np.min(np.where(np.array(CNA_fdrs)<t)[0])
        numcells = res.fdrs['num_detected'].iloc[i_first_below]
        CNA_cells_passing.append(numcells)
CNA_cells_passing = np.array(CNA_cells_passing).astype(int)   
x = np.arange(len(plot_thresholds))  # the label locations
width = 0.3  # the width of the bars

fig, ax = plt.subplots()
rects1 = ax.bar(x - width/2, CNA_cells_passing, width, label='CNA')
rects2 = ax.bar(x + width/2, MELD_cells_passing, width, label='MELD')

# Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_ylabel('# Cells Passing FDR Threshold')
ax.set_xlabel('FDR Threshold')
ax.set_xticks(x)
ax.set_xticklabels(plot_thresholds)
ax.legend()

fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Create Figure<jupyter_code>fig = plt.figure(figsize = (10,5))
gs = fig.add_gridspec(nrows=2, ncols=15, hspace = 0.55, wspace=0.5)
ax1 = fig.add_subplot(gs[0,0:4])
ax2 = fig.add_subplot(gs[0,5:9])
ax3 = fig.add_subplot(gs[0,10:15])
ax4 = fig.add_subplot(gs[1,0:3])
ax5 = fig.add_subplot(gs[1,3:6])
ax6 = fig.add_subplot(gs[1,6:9])
ax7 = fig.add_subplot(gs[1,9:12])
ax8 = fig.add_subplot(gs[1,12:15])


ax = ax1 # Scores on true sepsis values
val = sample_densities[sampleXmeta.index[true_pheno]].mean(axis=1) - sample_densities.mean(axis=1)
ax.scatter(d.obsm['X_tsne'][:,0], d.obsm['X_tsne'][:,1],
           alpha=0.5, c = val, cmap = "seismic", s=0.2, vmin=-color_scale_threshold, vmax=color_scale_threshold)
ax.set_title('Sepsis MELD Scores')
ax.axis('off')

ax = ax2 # Correlation between per-cell scores 
ax.scatter(cna_res['ncorrs'], sample_densities[sampleXmeta.index[true_pheno]].mean(axis=1), s=0.2)
ax.set_xticks([-0.4, 0, 0.4])
ax.set_xlabel('CNA Neighborhood Coefficient')
ax.set_ylabel('MELD Score')

ax = ax3
x = np.arange(len(plot_thresholds))  # the label locations
width = 0.3  # the width of the bars
rects1 = ax.bar(x - width/2, CNA_cells_passing/1000, width, label='CNA')
rects2 = ax.bar(x + width/2, MELD_cells_passing/1000, width, label='MELD')

# Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_ylabel('# Sig. Cells (Thousands)')
ax.set_xlabel('FDR Threshold')
ax.set_xticks(x)
ax.set_xticklabels(plot_thresholds)
ax.legend()

np.random.seed(0)
for ax in [ax4, ax5, ax6, ax7, ax8]: 
    ix = np.argsort(np.random.randn(len(true_pheno)))
    null = true_pheno[ix]
    nullval = sample_densities[sampleXmeta.index[null]].mean(axis=1) - sample_densities.mean(axis=1)
    ax.scatter(d.obsm['X_tsne'][:,0], d.obsm['X_tsne'][:,1],
               alpha=0.5, c=nullval, cmap = "seismic", s=0.2, vmin=-color_scale_threshold, vmax=color_scale_threshold)
    nullres = cna.tl.association(d, null, local_test=False)
    ax.text(0.13, -0.05, 'CNA p = '+str(np.round(nullres.p,2)),
        transform=ax.transAxes, fontsize=10, color="black")
    ax.axis('off')
    
ax=ax6
ax.set_title('MELD Scores on Randomly-Shuffled Per-Sample Values') 
    
plt.tight_layout()
plt.savefig('../_figs/suppfig.meld.pdf')<jupyter_output>performing association test
performing association test
performing association test
performing association test
{'fraction_non_alphanumeric': 0.15438983986747654, 'fraction_numerical': 0.024406405300938707, 'mean_word_length': 5.7632561613144135, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 24, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '17696', 'license': 'MIT', 'n_tokens_mistral': 3861, 'n_tokens_neox': 3580, 'n_words': 616, 'path': 'meld/MELD_comparative_analysis.ipynb', 'repo_name': 'yakirr/cna-display', 'size': '9055'}