<jupyter_start><jupyter_text>X-Inactivation Cell Type Differences DV AnalysisDetection of X-inactivation via differential varianceNo detectable differences.<jupyter_code>%env PYTHONPATH='/Users/mincheolkim/Github/scrna-parameter-estimation/simplesc'
import pandas as pd
import matplotlib.pyplot as plt
import scanpy.api as sc
import scipy as sp
import itertools
import numpy as np
import scipy.stats as stats
from scipy.integrate import dblquad
import scipy.sparse as sparse
import seaborn as sns
import imp
import time
from statsmodels.stats.multitest import fdrcorrection
import sys
sys.path.append('/Users/mincheolkim/Github/scrna-parameter-estimation/simplesc')
import simplesc
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
import warnings
warnings.filterwarnings('ignore')
data_path = '/Users/mincheolkim/Google Drive/UCSF/research/parameter_estimation/x_inactivation_data/'<jupyter_output><empty_output><jupyter_text>Read the the Lupus AnnData object<jupyter_code>adata = sc.read(data_path + 'lupus_annotated_nonorm_V6_x_genes.h5ad')<jupyter_output><empty_output><jupyter_text>Read cell type and individual information<jupyter_code>#ct_list = pd.read_csv(data_path + 'lupus_ct_list.csv')['ct_cov'].tolist()
ct_list = ['Tc', 'Th', 'cM', 'B', 'NK']
ind_list = pd.read_csv(data_path + 'lupus_ind_list.csv')
males = set(ind_list.query('Female == 0.0').ind_cov)
females = set(ind_list.query('Female == 1.0').ind_cov)<jupyter_output><empty_output><jupyter_text>Read DE and DV resultsThe results are computed in the `compute_x_inactivation_statistics.py` script in the Wynton cluster.<jupyter_code>de_t_stats = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_de_t_stats.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}
dv_t_stats = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_dv_t_stats.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}

de_pval = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_de_pvals.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}
dv_pval = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_dv_pvals.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}

mean = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_ct_mean.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}
var = {ct:pd.read_csv(data_path + 'ind_ct_combined_statistics/{}_ct_var.csv'.format(ct), index_col=0).T.dropna(how='all') for ct in ct_list}<jupyter_output><empty_output><jupyter_text>Male vs female analysis<jupyter_code>for ct in ct_list:

    mean_pvals = np.zeros(mean[ct].shape[0])
    var_pvals = np.zeros(var[ct].shape[0])
    
    var_diff = np.zeros(var[ct].shape[0])
    
    # Compute difference in means
    idx = 0
    for gene, row in mean[ct].iterrows():
        _, mean_pvals[idx] = stats.mannwhitneyu(row[males], row[females])
        idx += 1
        
    # Compute difference in vars
    idx = 0
    for gene, row in var[ct].iterrows():
        var_diff[idx] = row[females].mean() - row[males].mean()
        _, var_pvals[idx] = stats.mannwhitneyu(row[males], row[females])
        idx += 1
        
        
    mean[ct]['sex_diff_pval'] = mean_pvals
    var[ct]['sex_diff_pval'] = var_pvals
    var[ct]['var_diff'] = var_diff
    
    _, mean[ct]['sex_diff_fdr'] = fdrcorrection(mean_pvals)
    _, var[ct]['sex_diff_fdr'] = fdrcorrection(var_pvals)
    
    print(ct)
    print('mean', mean[ct].query('sex_diff_fdr < 0.1').index.tolist())
    print('var', var[ct].query('sex_diff_fdr < 0.1 & var_diff > 0').index.tolist())
    print()
mean_T = {ct:mean[ct].T for ct in ct_list}
var_T = {ct:var[ct].T for ct in ct_list}

for ct in ct_list:
    mean_T[ct]['sex'] = mean_T[ct].index.map(lambda x: 'male' if x in males else 'female')
    var_T[ct]['sex'] = var_T[ct].index.map(lambda x: 'male' if x in males else 'female')
sns.stripplot(x='sex', y='RPS4X', data=mean_T['Tc'], jitter=True)
#sns.boxplot(x='sex', y='RPL36A', data=mean_T['Tc'])
sns.stripplot(x='sex', y='RPS4X', data=var_T['Tc'], jitter=True)
#sns.boxplot(x='sex', y='XIST', data=var_T['Tc'])
sns.stripplot([0, 1], [
            var[ct].loc[gene][males].dropna(),
            var[ct].loc[gene][females].dropna()])
gene = 'CD99'
plt.figure(figsize=(20, 5))
plt.subplots_adjust(wspace=0.7, hspace=0.5)

for idx, ct in enumerate(ct_list):
    try:
        plt.subplot(2, len(ct_list), idx+1);
        plt.title(ct + ' Mean')
        plt.boxplot([
            mean[ct].loc[gene][males].dropna(),
            mean[ct].loc[gene][females].dropna()])
    except: 
        continue

for idx, ct in enumerate(ct_list):
    try:
        plt.subplot(2, len(ct_list), len(ct_list) + idx+1);
        plt.title(ct + ' Var')
        plt.boxplot()
    except: 
        continue
gene = 'XIST'
plt.figure(figsize=(20, 5))
plt.subplots_adjust(wspace=0.7, hspace=0.5)

for idx, ct in enumerate(ct_list):
    try:
        plt.subplot(2, len(ct_list), idx+1);
        plt.title(ct + ' Mean')
        plt.boxplot([
            mean[ct].loc[gene][males].dropna(),
            mean[ct].loc[gene][females].dropna()])
    except: 
        continue

for idx, ct in enumerate(ct_list):
    try:
        plt.subplot(2, len(ct_list), len(ct_list) + idx+1);
        plt.title(ct + ' Var')
        plt.boxplot([
            var[ct].loc[gene][males].dropna(),
            var[ct].loc[gene][females].dropna()])
    except: 
        continue<jupyter_output><empty_output><jupyter_text>CT differences<jupyter_code>female_dv = pd.read_csv(data_path + 'female_specific_ct_dv.csv', index_col=0)
def plot_ct_diff(gene):
    
    plt.figure(figsize=(20, 5));
    plt.subplots_adjust(wspace=0.1, hspace=0.5)
    
    plt.subplot(2, 2, 1);
    plt.boxplot([mean[ct].loc[gene][females].dropna().values for ct in ct_list]);
    plt.xticks(np.arange(len(ct_list))+1, ct_list)
    plt.title('Mean - female');

    plt.subplot(2, 2, 2);
    plt.boxplot([mean[ct].loc[gene][males].dropna().values for ct in ct_list]);
    plt.xticks(np.arange(len(ct_list))+1, ct_list)
    plt.title('Mean - male');
    
    plt.subplot(2, 2, 3);
    plt.boxplot([var[ct].loc[gene][females].dropna().values for ct in ct_list]);
    plt.xticks(np.arange(len(ct_list))+1, ct_list)
    plt.title('Var - female');

    plt.subplot(2, 2, 4);
    plt.boxplot([var[ct].loc[gene][males].dropna().values for ct in ct_list]);
    plt.xticks(np.arange(len(ct_list))+1, ct_list)
    plt.title('Var - male');
mean['cM'].shape
temp = np.diag(np.load(data_path + 'IGTB986_B_ct_cov.npy'))
temp = pd.read_csv(data_path + 'ind_ct_combined_statistics/B_ct_mean.csv', index_col=0).T
temp
mean['B'].shape
temp.shape
temp.shape<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1638058597218112, 'fraction_numerical': 0.012577685705830127, 'mean_word_length': 4.402877697841727, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 28, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '529970', 'license': 'MIT', 'n_tokens_mistral': 2759, 'n_tokens_neox': 2637, 'n_words': 438, 'path': 'analysis/x_inactivation/x_inactivation_ind_ct_dff_analysis.ipynb', 'repo_name': 'yelabucsf/scrna-parameter-estimation', 'size': '6758'}