# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Setup

# +
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc

plt.rcParams['pdf.fonttype'] = 42
plt.rcParams['ps.fonttype'] = 42
# -

# # Load host data
#

# Read the raw counts file

bulk_data = pd.read_csv('../data/swab_gene_counts.csv', index_col=0)
display(bulk_data.shape)
bulk_data.head()

# And convert the gene identifiers to gene symbol

import mygene
mg = mygene.MyGeneInfo()

mygene_converter = mg.querymany(list(bulk_data.index),scopes='all', species='human', as_dataframe=True)

bulk_data['symbol'] = mygene_converter['symbol']
bulk_data = bulk_data.set_index('symbol')
bulk_data.head()

# Finally write the csv to use in CibersortX.

bulk_data.to_csv('../data/swab_gene_counts_with_gene_symbol.csv')

# # Single cell datasets

# The next step is to prepare the single cell dataset. We are using **Travaglini et al, Human Lung Cell Atlas** (bioRxiv 2019)

# Note: If you don't have the data you can easily download it from synapse

# +
# import synapseclient 
# import synapseutils 
 
# syn = synapseclient.Synapse() 
# syn.login('username','userpassword') 
# krasnow_facs = synapseutils.syncFromSynapse(syn, 'syn21560510') 
# krasnow_droplet = synapseutils.syncFromSynapse(syn, 'syn21560511')
# -

# But here we're using the local copies

datasets_dir = '/mnt/ibm_lg/angela/sfaira/data/human/lung/'

datasets = ['droplet_normal_lung_blood_scanpy.20200205.RC4.h5ad',
 'facs_normal_lung_blood_scanpy.20200205.RC4.h5ad']

# We grouped the cells in the dataset according to their overall category.

# +
cell_dict_facs = {
 'B':'B cell',
 'Basal':'Basal',
 'Basophil/Mast 1':'Basophil/Mast',
 'CD4+ Memory/Effector T':'T cell',
 'CD4+ Naive T':'T cell',
 'CD8+ Memory/Effector T':'T cell',
 'CD8+ Naive T':'T cell',
 'Ciliated':'Ciliated',
 'Classical Monocyte':'monocytes/macrophages',
 'Club':'Club',
 'Dendritic':'Dendritic',
 'Differentiating Basal':'Basal',
 'Goblet':'Goblet',
 'IGSF21+ Dendritic':'Dendritic',
 'Intermediate Monocyte':'monocytes/macrophages',
 'Macrophage':'monocytes/macrophages',
 'Myeloid Dendritic Type 2':'Dendritic',
 'Natural Killer':'T cell',
 'Natural Killer T':'T cell',
 'Neutrophil':'Neutrophil',
 'Nonclassical Monocyte':'monocytes/macrophages',
 'Plasmacytoid Dendritic':'Dendritic',
 'Proliferating NK/T':'T cell'
}

cell_dict_droplet = {
 'B':'B cell',
 'Basal':'Basal',
 'Basophil/Mast 1':'Basophil/Mast',
 'Basophil/Mast 2':'Basophil/Mast',
 'CD4+ Memory/Effector T':'T cell',
 'CD4+ Naive T':'T cell',
 'CD8+ Memory/Effector T':'T cell',
 'CD8+ Naive T':'T cell',
 'Ciliated':'Ciliated',
 'Classical Monocyte':'monocytes/macrophages',
 'Club':'Club',
 'Dendritic':'Dendritic',
 'Differentiating Basal':'Basal',
 'EREG+ Dendritic':'Dendritic',
 'Goblet':'Goblet',
 'IGSF21+ Dendritic':'Dendritic',
 'Intermediate Monocyte':'monocytes/macrophages',
 'Ionocyte':'Ionocyte',
 'Macrophage':'monocytes/macrophages',
 'Myeloid Dendritic Type 1':'Dendritic',
 'Myeloid Dendritic Type 2':'Dendritic',
 'Natural Killer':'T cell',
 'Natural Killer T':'T cell',
 'Neutrophil':'Neutrophil',
 'Nonclassical Monocyte':'monocytes/macrophages',
 'OLR1+ Classical Monocyte':'monocytes/macrophages',
 'Plasmacytoid Dendritic':'Dendritic',
 'Proliferating Basal':'Basal',
 'Proliferating Macrophage':'monocytes/macrophages',
 'Proliferating NK/T':'T cell',
 'Proximal Basal':'Basal',
 'Proximal Ciliated':'Ciliated',
 'TREM2+ Dendritic':'Dendritic'
}
   

# -

for ds in datasets:
    print(ds)
    adata = sc.read_h5ad(datasets_dir+ds)
    
    adata.obs['free_annotation'] = [f.split('_')[0] for f in adata.obs['free_annotation']]
    
    if ds == 'droplet_normal_lung_blood_scanpy.20200205.RC4.h5ad':
        adata = adata[adata.obs['free_annotation'].isin(cell_dict_droplet.keys())].copy()
    else:    
        adata = adata[adata.obs['free_annotation'].isin(cell_dict_facs.keys())].copy()
        
    adata.obs['free_annotation'] = adata.obs['free_annotation'].astype(str)
    adata.obs['free_annotation'] = adata.obs['free_annotation'].astype('category')
    sc.tl.rank_genes_groups(adata, 'free_annotation', method='wilcoxon')
    
    adata_rank_genes = pd.DataFrame(adata.uns['rank_genes_groups']['names'])
    
    adata_rank_genes_list = []
    for i in adata_rank_genes.columns:
        adata_rank_genes_list = adata_rank_genes_list+list(adata_rank_genes[i])
    
    adata_rank_genes_list = list(set(adata_rank_genes_list))
    try:
        adata.X = adata.raw.X
    except:
        print("no raw slot available")
    adata = adata[:,adata_rank_genes_list]

    adatapd = pd.DataFrame(adata.X.transpose(),
                          columns=adata.obs['free_annotation'],
                         index = adata_rank_genes_list)

    if ds == 'droplet_normal_lung_blood_scanpy.20200205.RC4.h5ad':
        adatapd.columns = [cell_dict_droplet[c] for c in adatapd.columns]
    else:
        adatapd.columns = [cell_dict_facs[c] for c in adatapd.columns]
        
    path = 'cell-type-signatures-for-deconvolution_'+ds.split(".")[0]+'.txt'
    adatapd.to_csv(path, sep='\t')

adatapd.head()

# This produces to txt files ready to use with CiberSortX.

# # CiberSort matrices after cell selection

# ## Load metadata

host_metadata = pd.read_csv(
    '../data/metatable_with_viral_status.csv')
host_metadata.head()

# ## FACS

# Read the CiberSortX results using the smartseq2 dataset

cibersort_facs_abs = pd.read_excel('../results/deconvolution/CIBERSORTx_Krasnow_facs_small_list_Results.xlsx',index_col=0)
cibersort_facs_abs.head()

cibersort_facs_abs.shape

facs_deconv = cibersort_facs_abs
facs_deconv = facs_deconv[facs_deconv.columns[:-3]]
facs_deconv = facs_deconv.merge(host_metadata[['CZB_ID','viral_status']],left_on='Mixture',right_on='CZB_ID')
facs_deconv = facs_deconv.set_index('CZB_ID')
facs_deconv.head()

df_facs = pd.pivot_table(facs_deconv.reset_index(), 
                    values=['Ciliated', 'Club', 'Goblet', 'Basal', 'T cell','monocytes/macrophages',
                            'Neutrophil', 'Basophil/Mast', 'Dendritic','B cell'],
               index=['viral_status','CZB_ID']).reset_index()
df_facs = df_facs.sort_values(by='viral_status')
df_facs['Viral Load'] = df_facs['viral_status']
df_facs.loc[df_facs['viral_status']=='no_virus','Viral Load'] = 'Non-Viral ARI'
df_facs.loc[df_facs['viral_status']=='other_virus','Viral Load'] = 'Other Viral ARI'
df_facs.loc[df_facs['viral_status']=='SC2','Viral Load'] = 'COVID-19'
df_facs.head()


# ## Droplet

# Read the CiberSortX results using the 10X dataset

cibersort_droplet_abs = pd.read_excel(
    '../results/deconvolution/CIBERSORTx_Krasnow_droplet_small_list_Results.xlsx',index_col=0)
cibersort_droplet_abs.head()

droplet_deconv = cibersort_droplet_abs
droplet_deconv = droplet_deconv[droplet_deconv.columns[:-3]]
droplet_deconv = droplet_deconv.merge(host_metadata[['CZB_ID','viral_status']],left_on='Mixture',right_on='CZB_ID')
droplet_deconv = droplet_deconv.set_index('CZB_ID')
droplet_deconv.head()

df = pd.pivot_table(droplet_deconv.reset_index(), 
                    values=['Dendritic', 'B cell', 'monocytes/macrophages', 'T cell',
                            'Club','Basal', 'Basophil/Mast', 'Ciliated', 'Ionocyte','Goblet'],
               index=['viral_status','CZB_ID']).reset_index()
df = df.sort_values(by='viral_status')
df['Viral Load'] = df['viral_status']
df.loc[df['viral_status']=='no_virus','Viral Load'] = 'Non-Viral ARI'
df.loc[df['viral_status']=='other_virus','Viral Load'] = 'Other Viral ARI'
df.loc[df['viral_status']=='SC2','Viral Load'] = 'COVID-19'
df.head()


# ## FACS + droplet

# Combine both deconvolutions to get the complete picture

# +
df_facs_droplet = pd.DataFrame()
df_facs_droplet['B cell'] = (df['B cell'] + df_facs['B cell'])/2
df_facs_droplet['Basal'] = (df['Basal'] + df_facs['Basal'])/2
df_facs_droplet['Basophil/Mast'] = (df['Basophil/Mast'] + df_facs['Basophil/Mast'])/2
df_facs_droplet['Ciliated'] = (df['Ciliated'] + df_facs['Ciliated'])/2
df_facs_droplet['Club'] = (df['Club'] + df_facs['Club'])/2
df_facs_droplet['Dendritic'] = (df['Dendritic'] + df_facs['Dendritic'])/2
df_facs_droplet['Goblet'] = (df['Goblet'] + df_facs['Goblet'])/2
df_facs_droplet['T cell'] = (df['T cell'] + df_facs['T cell'])/2
df_facs_droplet['Monocytes/macrophages'] = (df['monocytes/macrophages'] + df_facs['monocytes/macrophages'])/2

df_facs_droplet['Neutrophil'] = df_facs['Neutrophil']
df_facs_droplet['Ionocyte'] = df['Ionocyte']

df_facs_droplet['czb_id'] = df_facs['CZB_ID']
df_facs_droplet['viral_load'] = df_facs['Viral Load']

df_facs_droplet = df_facs_droplet.sort_index(axis=1)

df_facs_droplet['sum'] = df_facs_droplet.iloc[:, :11].sum(axis=1)
df_facs_droplet.iloc[:, :11] = df_facs_droplet.iloc[:, :11].div(df_facs_droplet['sum'],axis=0)
df_facs_droplet['sum'] = df_facs_droplet.iloc[:, :11].sum(axis=1)

df_facs_droplet.head()
# -

# And write the results to a csv file.

df_facs_droplet.to_csv('../results/deconvolution/CIBERSORTx_Results_Krasnow_facs_droplet.csv')

# # Visualizations and statistical test

# +
from statannot import add_stat_annotation

cells = ['B cell', 'Basal', 'Basophil/Mast', 'Ciliated', 
       'Dendritic', 'Goblet', 'Ionocyte', 'T cell',
       'Monocytes/macrophages', 'Neutrophil']
i = 0
j = 0

fig, axes = plt.subplots(nrows=5, ncols=2,figsize=(10,18))
fig.subplots_adjust(hspace=.25)

covid_color_dict = {'Non-Viral ARI':'#111111','Other Viral ARI':'#E84646','COVID-19':'#65ADC2'}

for ax,c in zip(axes.flatten(),cells):
    
    sns.boxplot(data = df_facs_droplet, order = ['Non-Viral ARI','Other Viral ARI','COVID-19'],
               x = 'viral_load', y=c, palette=covid_color_dict,ax=ax,showfliers=False)
    try:
        add_stat_annotation(ax,data=df_facs_droplet, x = 'viral_load', y=c, 
                            order = ['Non-Viral ARI','Other Viral ARI','COVID-19'],
                                box_pairs=[('Non-Viral ARI','Other Viral ARI'),
                                           ('Other Viral ARI','COVID-19'),
                                           ('Non-Viral ARI','COVID-19')],
                                test='Mann-Whitney', loc='outside', verbose=0, text_format='full')# use text_format="star") for different p-value format
        
    except:
        print(c + ' didnt work')
    sns.swarmplot(data = df_facs_droplet, order = ['Non-Viral ARI','Other Viral ARI','COVID-19'],
               x = 'viral_load', y=c, palette=covid_color_dict,ax=ax)
    for patch in ax.artists:
        r, g, b, a = patch.get_facecolor()
        patch.set_facecolor((r, g, b, .1))
    if j<1:
        j+=1
    else:
        i+=1
        j = 0
    
    ax.set_xlabel('')

plt.tight_layout()
plt.savefig('../figures/deconvolution/Host_deconvolution_Krasnow_facs_droplet_MWW_Bonferroni.pdf')


# +
f = plt.figure()
df_facs_droplet.groupby(['viral_load'])['B cell','Basal','Basophil/Mast','Ciliated','Club','Dendritic',
                     'Goblet','Ionocyte','T cell','Monocytes/macrophages', 'Neutrophil'].agg(np.mean).reindex(
    ['Non-Viral ARI','Other Viral ARI','COVID-19']).plot(
    kind='bar',stacked=True, ax=f.gca(),colormap='tab20')

plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))

plt.tight_layout()
plt.savefig('../figures/deconvolution/Host_deconvolution_Krasnow_facs_droplet_cell_fractions.pdf')
# -



StarcoderdataJupyterScriptsDedupFiltered