<jupyter_start><jupyter_text>Pathway activity analysis Modeled after Xiao et al. Nat Commun (2019) Fig. 2For the i-th metabolic gene, the mean expression level across cells of the j-th cell type is$E_{i,j} = \frac{\sum_{k=1}^{n_j} g_{i,k}}{n_j}, i \in 1 ... M, j \in 1 ... N $in which $n_j$ is the number of cells in the j-th cell type, $g_{i,k}$ is the expression level of the i-th gene in the k-th cell in this cell type, M is the number of metabolic genes, and N is the number of cell types.Relative expression level of gene i in cell type j is$r_{i,j} = \frac{E_{i,j}}{1/N \sum_j^N E_{i,j}}$Exclude genes with low expression level or high dropout rate:Exclude outliers in each pathway defined by genes with relative expression levels >3 x 75th percentile expression or <1/3 x 25th percentilePathway activity score for the t-th pathway and the j-th cell type is the weighted average of $r_{i,j}$ over all genes included in this pathway:$p_{t,j} = \frac{\sum_{i=1}^{m_t} w_i \times r_{i,j}}{\sum_{i=1}^{m_t} w_i}$where $m_t$ is the number of genes in the pathway t, $w_i$ is the weighting factor equal to the reciprocal of the number of pathways that include the i-th gene.<jupyter_code>pathways = [list(genes) for genes in list(metab_genes_dict.values())]

# wi is the reciprocal of the number of pathways that include the ith gene,
# summarized in array w
w = [[pathway.count(gene) / list(metab_gene_list).count(gene) for gene 
      in adata_metab.var_names.values] for pathway in pathways]
# Here, cell type is defined as individual clusters ('louvain')
# within a given group of immune cells (B cells, T cells, etc.)
# Clusters calculated from all genes rather than just metabolic genes
# in order to associate the metabolic gene set enrichment with the non-metabolic analyses
def pathway_activity(data, data_metab, cluster_label, print_labels=False, nanval=0, exclude_clusters=[]):
    cluster_ids = list(data.obs[cluster_label].unique())
    cluster_ids.sort()
    
    for cluster_id in exclude_clusters:
        cluster_ids.remove(cluster_id)
        
    mel = np.zeros([len(cluster_ids), len(data_metab.X[0])])
    for i, cluster_id in enumerate(cluster_ids):
        if print_labels:
            print(cluster_id)
        cluster = data_metab[data.obs[cluster_label]==cluster_id]
        g = cluster.X # Matrix of gene expression of ith gene in the kth cell of this cluster (cell type)
        n = len(cluster) # Number of cells in the cluster (cell type)
        # Calculate mean expression level
        mel[i] = np.sum(g, axis=0) / n # Sum gene expressions per gene across all cells in this cluster (cell type),
            # divided by number of cells in the cluster (cell type)        

    # Remove genes that are 0 in any cluster to avoid extreme ratio value
    mel = np.array(mel)
    m = min(mel[mel > 0])
    mel[mel==0] = m
            
    # Calculate relative expression level of each gene
    rel = mel / (np.sum(mel, axis=0) / len(mel))

    # Filtering
    lowExpr, highExpr = np.percentile(rel, [25, 75], axis=None) 
    
    tooHigh = (rel > 3 * highExpr) * rel
    tooLow = (rel < lowExpr / 3) * rel
    
    rel -= tooHigh
    rel -= tooLow
    rel = np.nan_to_num(rel, nan=nanval)
    
    # Calculate the pathway activity score
    # Pathway activity matrix
    # Each row is a pathway, each column is a cluster
    p = np.array([np.sum(w * rel[j], axis=1) / np.sum(w, axis=1) for j in range(len(rel))]).transpose()
    
    # Remove nan values
    p = np.nan_to_num(p, nan=1)
    
    return p
adata.obs['disease_state'].unique()
adata.obs['cell_type'].unique()
adata.obs['cell_type_per_dz'] = adata.obs['cell_type'].astype(str) + ' of ' + adata.obs['disease_state'].astype(str)
adata_metab.obs['cell_type_per_dz'] = adata_metab.obs['cell_type'].astype(str) + \
                ' of ' + adata_metab.obs['disease_state'].astype(str)
adata.obs['cell_type_per_dz'].unique()
adata.obs['cell_type'].hist()
adata.obs['disease_state'].hist()<jupyter_output><empty_output><jupyter_text>Pathway activity analysis between cell types<jupyter_code>p = pathway_activity(adata, adata_metab, 'cell_type_per_dz', print_labels=True,
                    exclude_clusters=['Non PE T-cells of Control', 'Non T-cells of Control',
                    'Non T-cells of COVID-19', 'Unconfident PE T-cells of COVID-19',
                    'PE T-cells of Control', 'Unconfident PE T-cells of Hyperacute HIV',
                    'Unconfident PE T-cells of Control', 'Unconfident PE T-cells of Bacterial Sepsis'])<jupyter_output>Non PE T-cells of Bacterial Sepsis
Non PE T-cells of COVID-19
Non PE T-cells of Hyperacute HIV
PE T-cells of Bacterial Sepsis
PE T-cells of COVID-19
PE T-cells of Hyperacute HIV<jupyter_text>Calculate statistical significance: random permutation test1. Randomly shuffle cell type labels 5,000 times to simulate null distribution of pathway activity scores2. Compare to the pathway activity scores in the original, non-shuffled dataset3. p-value = fraction of random pathway activity scores     - $> p_{t,j}$ if $p_{t,j} > 1$, meaning that pathway activity is significantly higher     - $< p_{t,j}$ if $p_{t,j} < 1$, meaning that pathway activity is significantly lower<jupyter_code># Statistical significance of higher or lower pathway activity in a specific cell type
# Random permutation test
enriched = (p > 1)
decreased = (p < 1)

# Shuffle cell cluster labels 1000x
n = 1000
compares = np.zeros((len(p), len(p[0])))

# Declare relevant variables before for loop to prevent memory leaks
obs_name = 'cell_type_per_dz'

# Shuffle cell cluster labels 1000x
data = adata[~adata.obs[obs_name].isin(['Non PE T-cells of Control', 'Non T-cells of Control',
                    'Non T-cells of COVID-19', 'Unconfident PE T-cells of COVID-19',
                    'PE T-cells of Control', 'Unconfident PE T-cells of Hyperacute HIV',
                    'Unconfident PE T-cells of Control', 'Unconfident PE T-cells of Bacterial Sepsis'])]
data_metab = adata_metab[data.obs[obs_name].index]

data.obs[obs_name + '_shuffle'] = pd.Series()
data_metab.obs[obs_name + '_shuffle'] = pd.Series()
clusters_shuffled = np.zeros(len(data.obs))
p_shuffle = np.zeros([len(w), len(data.obs[obs_name].unique())])

for i in range(n):
    if i % 100 == 0:
        print(i)
    clusters_shuffled = np.random.permutation(data.obs[obs_name])
    data.obs[obs_name + '_shuffle'] = clusters_shuffled
    data_metab.obs[obs_name + '_shuffle'] = clusters_shuffled
    
    # Compare new pathway activity scores to original dataset pathway activity scores
    p_shuffle = pathway_activity(data, data_metab, obs_name + '_shuffle')
    compares += enriched * (p_shuffle > p) + decreased * (p_shuffle < p)

compares /= n
# Non-significantly changed activity scores are reverted to 1 (no activity change)
significant = (compares < 0.01)
p_sig = p ** significant
# Save pathway activity scores
df_p = pd.DataFrame(p_sig)
df_p.to_csv(master_dir + 'output/disease_comparisons/' + obs_name + '-metab-pathway_activity.tabular', 
                       sep='\t', header=None, index=None)
# Read saved pathway activity scores
df_p = pd.read_csv(master_dir + 'output/disease_comparisons/' + obs_name + '-metab-pathway_activity.tabular', 
                   sep='\t', header=None)
df_p.head()
fig = sns.clustermap(df_p.loc[:, [1, 4, 0, 3, 2, 5]], cmap='coolwarm', 
                     xticklabels=pd.Series(['Non-proliferative-active T-cells, bacterial sepsis',
        'Non-proliferative-active T-cells, COVID-19', 'Non-proliferative-active T-cells, hyperacute HIV', 
                                            'Proliferative-active T-cells, bacterial sepsis',
            'Proliferative-active T-cells, COVID-19', 
            'Proliferative-active T-cells, hyperacute HIV']).iloc[[1, 4, 0, 3, 2, 5]], 
                     yticklabels=[], #yticklabels=metab_genes.index, 
                     norm=DivergingNorm(1.0), row_cluster=False, col_cluster=False, figsize=(3, 10))
plt.savefig(master_dir + 'output/disease_comparisons/metab-pathways-per-dz_heatmap.pdf')<jupyter_output><ipython-input-248-20223066ccc6>:8: MatplotlibDeprecationWarning: 
The DivergingNorm class was deprecated in Matplotlib 3.2 and will be removed two minor releases later. Use TwoSlopeNorm instead.
  norm=DivergingNorm(1.0), row_cluster=False, col_cluster=False, figsize=(3, 10))<jupyter_text>Venn diagrams of differential metabolic gene expression or metabolic pathway activities Gene expression-based<jupyter_code>metab_genes.index
from matplotlib_venn import venn3, venn3_circles
# Get list of differentially expressed metabolic pathways and genes among the PE groups in each disease
# PE groups only
data = adata[adata.obs[obs_name].isin(['PE T-cells of COVID-19', 'PE T-cells of Bacterial Sepsis',
                                      'PE T-cells of Hyperacute HIV'])]
data_metab = adata_metab[data.obs[obs_name].index]

sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        key_added='rank_genes_groups')
sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='PE T-cells of COVID-19', key_added='rank_genes_groups_vs_COVID')
sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='PE T-cells of Hyperacute HIV', key_added='rank_genes_groups_vs_HIV')
sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='PE T-cells of Bacterial Sepsis', key_added='rank_genes_groups_vs_sepsis')
a = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of COVID-19', 
                                pval_cutoff=0.05, key='rank_genes_groups')['names']
b = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of Bacterial Sepsis', 
                                pval_cutoff=0.05, key='rank_genes_groups')['names']
c = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of Hyperacute HIV', 
                                pval_cutoff=0.05, key='rank_genes_groups')['names']

d = set(a).intersection(set(b))
e = set(b).intersection(set(c))
f = set(c).intersection(set(a))
g = d.intersection(set(c))

venn3(subsets = [set(a), set(b), set(c)], set_labels = ('COVID-19 PE T cells', 
                                                        'Bacterial sepsis PE T cells', 'Hyperacute HIV PE T cells'))
# Another Venn diagram with overlaps of upregulated genes in each PE group vs. non-PE group of the same disease

# Get list of differentially expressed metabolic pathways and genes among the PE groups in each disease
# PE groups only
data = adata[adata.obs[obs_name].isin(['Non PE T-cells of Bacterial Sepsis',
        'Non PE T-cells of COVID-19', 'Non PE T-cells of Hyperacute HIV', 'PE T-cells of Bacterial Sepsis',
        'PE T-cells of COVID-19', 'PE T-cells of Hyperacute HIV'])]
data_metab = adata_metab[data.obs[obs_name].index]

sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='Non PE T-cells of COVID-19', key_added='rank_genes_groups_COVID')
sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='Non PE T-cells of Hyperacute HIV', key_added='rank_genes_groups_HIV')
sc.tl.rank_genes_groups(data_metab, groupby=obs_name, n_genes=len(data_metab.var), 
                        reference='Non PE T-cells of Bacterial Sepsis', key_added='rank_genes_groups_sepsis')
a = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of COVID-19', 
                                pval_cutoff=0.05, key='rank_genes_groups_COVID')['names']
b = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of Bacterial Sepsis', 
                                pval_cutoff=0.05, key='rank_genes_groups_sepsis')['names']
c = sc.get.rank_genes_groups_df(data_metab, 'PE T-cells of Hyperacute HIV', 
                                pval_cutoff=0.05, key='rank_genes_groups_HIV')['names']

d = set(a).intersection(set(b))
e = set(b).intersection(set(c))
f = set(c).intersection(set(a))
g = d.intersection(set(c))

venn3(subsets = [set(a), set(b), set(c)], set_labels = ('COVID-19 PE T cells', 
                                                        'Bacterial sepsis PE T cells', 
                                                        'Hyperacute HIV PE T cells'))
plt.savefig(master_dir + 'output/disease_comparisons/gex_high_vs_non-PE/venn.pdf')
diseases = ['COVID', 'Sepsis', 'HIV', 'COVID and sepsis', 'Sepsis and HIV', 'COVID and HIV', 'All']
ai = set(a).difference(set(b)).difference(set(c))
bi = set(b).difference(set(a)).difference(set(c))
ci = set(c).difference(set(b)).difference(set(a))
di = d.difference(g)
ei = e.difference(g)
fi = f.difference(g)

# Save the lists of these genes
for i, gene_list in enumerate([ai, bi, ci, di, ei, fi, g]):
    f = open(master_dir + 'output/disease_comparisons/gex_high_vs_non-PE/' + diseases[i] + '.txt', 'w')
    for j in gene_list:
        f.write(j + '\n')
    f.close()
len(fi)<jupyter_output><empty_output><jupyter_text>Metabolic pathway-based<jupyter_code>df_metab_pathways = df_p.loc[:, [1, 4, 0, 3, 2, 5]]
df_metab_pathways.index = metab_genes.index
df_metab_pathways.columns = pd.Series(['Non PE T-cells of Bacterial Sepsis',
'Non PE T-cells of COVID-19', 'Non PE T-cells of Hyperacute HIV', 'PE T-cells of Bacterial Sepsis',
'PE T-cells of COVID-19', 'PE T-cells of Hyperacute HIV']).iloc[[1, 4, 0, 3, 2, 5]]

adata_metab_pathways = sc.AnnData(df_p.loc[:, [1, 4, 0, 3, 2, 5]].transpose())
adata_metab_pathways.var = pd.DataFrame(metab_genes.index.values).set_index(metab_genes.index.values)
adata_metab_pathways.obs = pd.DataFrame({obs_name : pd.Series(['Non PE T-cells of Bacterial Sepsis',
'Non PE T-cells of COVID-19', 'Non PE T-cells of Hyperacute HIV', 'PE T-cells of Bacterial Sepsis',
'PE T-cells of COVID-19', 'PE T-cells of Hyperacute HIV']).iloc[[1, 4, 0, 3, 2, 5]]})
adata_metab_pathways.obs[obs_name] = adata_metab_pathways.obs[obs_name].astype('category')
a = df_metab_pathways[df_metab_pathways['PE T-cells of COVID-19'] > 1.0].index.values
b = df_metab_pathways[df_metab_pathways['PE T-cells of Bacterial Sepsis'] > 1.0].index.values
c = df_metab_pathways[df_metab_pathways['PE T-cells of Hyperacute HIV'] > 1.0].index.values

venn3(subsets = [set(a), set(b), set(c)], set_labels = ('COVID-19 PE T cells', 
                                                        'Bacterial sepsis PE T cells', 
                                                        'Hyperacute HIV PE T cells'))
a = df_metab_pathways[df_metab_pathways['PE T-cells of COVID-19'] > df_metab_pathways['Non PE T-cells of COVID-19']].index.values
b = df_metab_pathways[df_metab_pathways['PE T-cells of Bacterial Sepsis'] > df_metab_pathways['Non PE T-cells of Bacterial Sepsis']].index.values
c = df_metab_pathways[df_metab_pathways['PE T-cells of Hyperacute HIV'] > df_metab_pathways['PE T-cells of Hyperacute HIV']].index.values

venn3(subsets = [set(a), set(b), set(c)], set_labels = ('COVID-19 PE T cells', 
                                                        'Bacterial sepsis PE T cells', 
                                                        'Hyperacute HIV PE T cells'))<jupyter_output>/home/jwlee/.local/lib/python3.8/site-packages/matplotlib_venn/_venn3.py:61: UserWarning: Circle C has zero area
  warnings.warn("Circle C has zero area")<jupyter_text>Box plot of metabolic pathway activities per cell group<jupyter_code>f = plt.figure(figsize=(4, 2))
f = sns.boxplot(x=np.repeat(df_metab_pathways.columns.values, 90), 
            y=np.array(df_metab_pathways.transpose()).flatten())
f.set_xticklabels(labels=df_metab_pathways.columns.values, rotation=90)
f.set_ylim([0, 2.5])
plt.savefig(master_dir + 'output/disease_comparisons/metab-pathways-per-dz.pdf')
temp = df_metab_pathways.columns.copy()
for i in temp:
    temp = temp.drop(i)
    for j in temp:
        print('\n', i, 'vs', j, ':', wilcoxon(df_metab_pathways[i], df_metab_pathways[j]).pvalue)
df_metab_pathways.columns.drop('PE T-cells of COVID-19')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12128543954339599, 'fraction_numerical': 0.013958682300390842, 'mean_word_length': 3.4815123714206284, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '327676', 'license': 'MIT', 'n_tokens_mistral': 5823, 'n_tokens_neox': 5316, 'n_words': 1440, 'path': 'disease_comparisons/.ipynb_checkpoints/enrichment_analysis-disease-comparison-cd8-checkpoint.ipynb', 'repo_name': 'jihoonlee0/SARS-CoV-2_metab', 'size': '16119'}