<jupyter_start><jupyter_text>Cell-cell communication analysis Performing differential expression (DE) analysis of the NanoString ROIs with `limma`<jupyter_code>%matplotlib inline

import numpy as np
import pandas as pd
import seaborn as sns
import scanpy as sc
import anndata
import matplotlib as mpl
import scipy
import matplotlib.pyplot as plt

sc.settings.verbosity = 3  # verbosity: errors (0), warnings (1), info (2), hints (3)
#sc.logging.print_versions()
sc.settings.set_figure_params(dpi=80)  # low dpi (dots per inch) yields small inline figures
nanostring_data_path = '../KidneyDataset/'<jupyter_output><empty_output><jupyter_text>Read in and prepare NanoString data<jupyter_code># target QC’d, filtered, and quantile normalized (Q3) target level count data
# going to feed this into limma
Kidney_Q3Norm_TargetCountMatrix = pd.read_csv(nanostring_data_path + 'Kidney_Q3Norm_TargetCountMatrix.txt',
                                             sep='\t', index_col=0)
# ROI metadata
Kidney_Sample_Annotations = pd.read_csv(nanostring_data_path + 'Kidney_Sample_Annotations.txt',
                                             sep='\t')
# here in Kidney_Sample_Annotations column 'region' there is no indication of whether tubule is proximal or distal
# need to add this depending on the SegmentLabel - neg means proximal tubules,  PanCK - distal

def add_exact_region(idx, table):
    curr_SegmentLabel = table.loc[idx,'SegmentLabel']
    if curr_SegmentLabel == 'Geometric Segment':
        return('glomerulus')
    elif curr_SegmentLabel == 'neg':
        return('proximal tubule')
    elif curr_SegmentLabel == 'PanCK':
        return('distal tubule')
    else:
        return('error')
Kidney_Sample_Annotations['idx'] = Kidney_Sample_Annotations.index
Kidney_Sample_Annotations['region_exact'] = Kidney_Sample_Annotations['idx'].apply(lambda x: add_exact_region(x,Kidney_Sample_Annotations))
Kidney_Sample_Annotations.set_index('SegmentDisplayName', inplace=True)
# getting rid of NaNs
Kidney_Sample_Annotations['pathology'] = [elem if str(elem) != 'nan' else 'NA' for elem in Kidney_Sample_Annotations['pathology']]
# We will be comparing all abnormal ROIs to all healthy ROIs (acc to pathology)
# and then this 1 list of DE I will artificially copy for each of the single cell clusters
comparison_cluster = 'pathology'
Kidney_Q3Norm_TargetCountMatrix = Kidney_Q3Norm_TargetCountMatrix.T
# make sure indices match
Kidney_Sample_Annotations = Kidney_Sample_Annotations.loc[Kidney_Q3Norm_TargetCountMatrix.index]
# object with all features that passed QC (.raw.X)
adata_full = anndata.AnnData(X = Kidney_Q3Norm_TargetCountMatrix,
                             obs = Kidney_Sample_Annotations,
                             var = pd.DataFrame(index = Kidney_Q3Norm_TargetCountMatrix.columns)
                            )<jupyter_output><empty_output><jupyter_text>_____________________________________________________________________________________________________________________________________________________________  DE with limma- all abnormal ROIs vs all healthy ROIs according to pathology annotation<jupyter_code># discard the ROIs that are not glomeruli (aka tubules)
adata_full = adata_full[adata_full.obs['pathology'].isin(['abnormal', 'healthy'])].copy()
# marker calling

t = adata_full.X.T
df = pd.DataFrame(data=t, columns= adata_full.obs.index, index=adata_full.var_names)

meta_df = pd.DataFrame(data={'Cell':list(adata_full.obs.index),
                             'cell_type':[ str(i) for i in adata_full.obs['pathology']],
                             #'sample':[ str(i) for i in adata_full.obs['sample']]
                            })
meta_df.set_index('Cell', inplace=True)
meta_df.reset_index(inplace=True)
%load_ext rpy2.ipython
outpath = './'
%%R
library(limma)
library(edgeR)
# because R replaces things
meta_df['Cell'] = [elem.replace(' | ','...') for elem in meta_df['Cell']]
meta_df['Cell'] = [elem.replace(' ','.') for elem in meta_df['Cell']]
#meta_df.head()

df.columns = [elem.replace(' | ','...') for elem in df.columns]
np.unique(meta_df['cell_type'], return_counts=True)
case = 'abnormal'
ctrl = 'healthy'
%%R -i df -i meta_df -i outpath -i ctrl -i case 

library(limma)
library(edgeR)

# Format
ex_mat=as.matrix(df)
rownames(meta_df) = meta_df$Cell

# subset meta
meta_df = subset(meta_df, cell_type %in% unlist(c(ctrl, case)) )

print(unique(meta_df$cell_type))

# Shared cells
shared_cells = intersect(rownames(meta_df), colnames(ex_mat))
message(length(shared_cells), ' shared cells')
ex_mat = ex_mat[, shared_cells]
meta_df = meta_df[shared_cells,]

# Filter lowly expressed genes
keep = rowSums(ex_mat, na.rm=T) > 0.1
ex_mat = ex_mat[ keep, ]
keep = aveLogCPM(ex_mat) > 0.1
ex_mat = ex_mat[ keep, ]

# Extract celltypes
cells = rownames(meta_df)
celltypes = unique(meta_df$cell_type)
covariates = meta_df$covariate

# Extract cells in cluster and rest
cells_case = rownames(subset(meta_df, cell_type == case))
cells_ctrl = rownames(subset(meta_df, cell_type == ctrl)) # changed from control to ctrl

# build cluster_type vector
cluster_type = rep(0, length(cells))
names(cluster_type) = cells
cluster_type[ cells_case ] = 'case'
cluster_type[ cells_ctrl ] = 'ctrl'

print(unique(cluster_type))

#design.matrix <- model.matrix(~ 0 + cluster_type + covariates)
design.matrix <- model.matrix(~ 0 + cluster_type)

# Now tell limma how do you want to compare (i.e. case vs control)
contrast.matrix <- makeContrasts(caseVScontrol = cluster_typecase - cluster_typectrl, levels = design.matrix)

# Make model and run contrasts
fit <- lmFit(ex_mat, design.matrix)
fit <- contrasts.fit(fit, contrast.matrix)
fit <- eBayes(fit)

# Make a dataframe containing the important data
results = topTable(fit, adjust="fdr", number = nrow(ex_mat), coef = 'caseVScontrol')

# Add and filter needed data
results$Gene = rownames(results)
results = results[ , c('Gene', 'logFC', 'P.Value', 'adj.P.Val')]
results$AveExpr_cluster = apply(ex_mat[ results$Gene, cells_case], 1, mean)
results$AveExpr_rest = apply(ex_mat[ results$Gene, cells_ctrl], 1, mean)
results$percentExpr_cluster = apply(ex_mat[ results$Gene, cells_case], 1, function(x) sum(c(x > 0)+0) ) / length(cells_case)
results$percentExpr_rest = apply(ex_mat[ results$Gene, cells_ctrl], 1, function(x) sum(c(x > 0)+0) ) / length(cells_ctrl)

results$AveExpr_cluster = round(results$AveExpr_cluster, 6)
results$AveExpr_rest = round(results$AveExpr_rest, 6)
results$percentExpr_cluster = round(results$percentExpr_cluster, 6)
results$percentExpr_rest = round(results$percentExpr_rest, 6)
# and store it as csv file
write.csv(results, file = paste0(outpath, case, '_vs_', ctrl, '_limma_DEGs.csv'), row.names = F, col.names = T, quote = F)
# let's have a look at the DEGs
DE_table = pd.read_csv('./abnormal_vs_healthy_limma_DEGs.csv', index_col=0)
DE_table
# how many significant DEGs?
DE_table_significant = DE_table[DE_table['adj.P.Val'] < 0.05]
DE_table_significant<jupyter_output><empty_output><jupyter_text>Creating a joint DE table out of limma outputs from NanoString data Because we don't have cell type resolution here, let's copy this DE table for every cell type artificially to then intersect with the expressed L/R from Young et al. data<jupyter_code># from Young et al.
meta = pd.read_csv('./cellphonedb_meta.tsv', sep='\t')
cell_types = np.unique(meta['cell_type'])
# get rid of 'celltype_'
cell_types = [elem[9:] for elem in cell_types]
cell_types
# let's have a look at the DEGs
DE_table = pd.read_csv('./abnormal_vs_healthy_limma_DEGs.csv')
DE_table
# artificially copying it for all the cell types
DE_tables = {}

for ct in cell_types:
    #print(ct)
    DE_tables[ct] = DE_table
DE_table
cluster_list = []

for ct in cell_types:
    #print(ct)
    cluster_list.append([ct]*len(DE_table))

cluster_list = [item for sublist in cluster_list for item in sublist]
# without any filtering
joint_DE_table = pd.concat(DE_tables.values())
joint_DE_table['cluster'] = cluster_list
joint_DE_table
joint_DE_table.to_csv('./joint_DEGs_list_all_cell_types_for_cellphone.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11057573511242896, 'fraction_numerical': 0.0053125772176921176, 'mean_word_length': 5.123298033282905, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '64816', 'license': 'MIT', 'n_tokens_mistral': 2757, 'n_tokens_neox': 2643, 'n_words': 681, 'path': 'Analysis/3_Cell-cell_comm_analysis/S4_DE_analysis_limma_ROIs.ipynb', 'repo_name': 'BioinfoTongLI/NanoHack21-1', 'size': '8094'}