# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:scvi]
#     language: python
#     name: conda-env-scvi-py
# ---

# + [markdown] nbpresent={"id": "1f01a376-e4e2-4d3c-8723-0f0d91682b22"}
# # X-Inactivation Cell Type Differences DV Analysis
#
# Detection of X-inactivation via differential variance
#
# No detectable differences.
# -

# %env PYTHONPATH='/Users/mincheolkim/Github/scrna-parameter-estimation/simplesc'

# + nbpresent={"id": "9bc26761-c243-4e99-98ca-0cc6425e76f5"}
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')

# + nbpresent={"id": "928db49a-7d74-456c-bf37-9192ad38ba99"}
data_path = '/Users/mincheolkim/Google Drive/UCSF/research/parameter_estimation/x_inactivation_data/'
# -

# ### Read the the Lupus AnnData object

adata = sc.read(data_path + 'lupus_annotated_nonorm_V6_x_genes.h5ad')

# ### Read cell type and individual information

#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)

# ### Read DE and DV results
#
# The results are computed in the `compute_x_inactivation_statistics.py` script in the Wynton cluster.

# +
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}

# -

# ### Male vs female analysis

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
# -

# ### CT differences

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



StarcoderdataJupyterScriptsDedupFiltered