Plate_bias

0.1 Parameters

0.2 functions

source_from_github(repositoy = "DEG_functions",version = "0.2.3")

ℹ SHA-1 hash of file is 042341a0866ca219ce0c03aa0c61b7ca4b1dc928

0.3 Data

acc1_cancer_cells = readRDS(file = "./Data/acc1_cancer_cells_15KnCount_V4.RDS")

1 Title

2 cnmf programs are by plate

from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_' + nfeatures+ 'vargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}
print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),combine = T),
          title = "metagenes expression")

metagenes expression

NA

print_tab(plt = DimPlot(acc1_cancer_cells,pt.size = 1,group.by = "orig.ident")
          ,title = "umap by plate")

umap by plate

NA

nmf_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("metagene.1","metagene.2", "orig.ident"))
# myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )

plt = ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.1",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))+

ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.2",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))

print_tab(plt = plt, title = "plates diff")

plates diff

NA

acc1_cancer_cells = SetIdent(object = acc1_cancer_cells, value = "orig.ident")
print_tab(plt = 
            DotPlot(object = acc1_cancer_cells, features = c("metagene.1","metagene.2","metagene.3","metagene.4"),scale = F)
          ,title = "dot plot")

dot plot

NA

3 regress nFeature_RNA

#regress nFeature_RNA
acc1_cancer_cells <- ScaleData(acc1_cancer_cells, vars.to.regress = c("nFeature_RNA","percent.mt","nCount_RNA"))
acc1_cancer_cells <- RunPCA(acc1_cancer_cells, features = VariableFeatures(object = acc1_cancer_cells))

ElbowPlot(acc1_cancer_cells, ndims = 50) # checking the dimensionality 

pc2use = 1:10

acc1_cancer_cells <- FindNeighbors(acc1_cancer_cells, dims = pc2use)
acc1_cancer_cells <- FindClusters(acc1_cancer_cells, resolution = 0.5)
acc1_cancer_cells <- RunUMAP(acc1_cancer_cells, dims = pc2use)

DimPlot(acc1_cancer_cells,pt.size = 1,group.by = "orig.ident")

3.1 cNMF

library(reticulate)

#write expression
nfeatures = 2000
nfeatures_name = (nfeatures/1000) %>% as.character()
acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = nfeatures)
vargenes = VariableFeatures(object = acc1_cancer_cells)
hmsc_scaled_expression = FetchData(object = acc1_cancer_cells,vars = vargenes, slot = "scale.data")
hmsc_scaled_expression[hmsc_scaled_expression<0]= 0

write.table(x = hmsc_scaled_expression ,file = paste0('./Data/cNMF/hmsc_scaled_expression',nfeatures_name,'Kvargenes.txt'),sep = "\t")

from cnmf import cNMF
import numpy as np
nfeatures_name = r.nfeatures_name
name = 'HMSC_cNMF_scaled_'+nfeatures_name+'Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/hmsc_scaled_expression'+nfeatures_name+'Kvargenes.txt'
tpm_fn = counts_fn ## This is a weird case where because this dataset is not 3' end umi sequencing, we opted to use the TPM matrix as the input matrix rather than the count matrix

cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, seed=14,tpm_fn=tpm_fn)

cnmf_obj.factorize(worker_i=0, total_workers=1)

cnmf_obj.combine()
cnmf_obj.k_selection_plot()

3.2 Save object

# import pickle
# f = open('./Data/cNMF/HMSC_cNMF_scaled_'+nfeatures_name+'Kvargenes/cnmf_obj.pckl', 'wb')
# pickle.dump(cnmf_obj, f)
# f.close()

3.3 Load object

from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_scaled_' + nfeatures+ 'vargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

knitr::include_graphics('./Data/cNMF/HMSC_cNMF_scaled_2Kvargenes/HMSC_cNMF_scaled_2Kvargenes.k_selection.png')

selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm

4 regress nFeature_RNA programs

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}
print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),combine = T),
          title = "metagenes expression")

metagenes expression

NA

print_tab(plt = DimPlot(acc1_cancer_cells,pt.size = 1,group.by = "orig.ident")
          ,title = "umap by plate")

umap by plate

NA

nmf_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("metagene.1","metagene.2", "orig.ident"))
# myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )

plt = ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.1",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))+

ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.2",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))

print_tab(plt = plt, title = "plates diff")

plates diff

NA

acc1_cancer_cells = SetIdent(object = acc1_cancer_cells, value = "orig.ident")
print_tab(plt = 
            DotPlot(object = acc1_cancer_cells, features = c("metagene.1","metagene.2","metagene.3","metagene.4"),scale = F)
          ,title = "dot plot")

dot plot

NA

5 Harmony

acc1_cancer_cells = readRDS(file = "./Data/acc1_cancer_cells_15KnCount_V4.RDS")

sc <- import('scanpy', convert = FALSE)
gene_expression = t(as.matrix(GetAssayData(acc1_cancer_cells,slot='data'))) 

gene_expression = 2**gene_expression #convert log2(TPM+1) to TPM+1
gene_expression = gene_expression-1 #convert TPM+1 to TPM

adata <- sc$AnnData(
  X   = gene_expression, 
  obs = acc1_cancer_cells[[]],
  var = GetAssay(acc1_cancer_cells)[[]]
)
rm(gene_expression)

# %matplotlib inline
# %load_ext autoreload
# %autoreload 2
import scipy
import scanpy as sc
import matplotlib.pyplot as plt
import harmonypy
from harmonypy import run_harmony
import sys
from cnmf import cNMF

import numpy as np
import pandas as pd
from sklearn.decomposition import NMF
import seaborn as sns
import yaml
## copied from harmonypy.py
def moe_correct_ridge(Z_orig, Z_cos, Z_corr, R, W, K, Phi_Rk, Phi_moe, lamb):
    Z_corr = Z_orig.copy()
    for i in range(K):
        Phi_Rk = np.multiply(Phi_moe, R[i,:])
        x = np.dot(Phi_Rk, Phi_moe.T) + lamb
        W = np.dot(np.dot(np.linalg.inv(x), Phi_Rk), Z_orig.T)
        W[0,:] = 0 # do not remove the intercept
        Z_corr -= np.dot(W.T, Phi_Rk)
    Z_cos = Z_corr / np.linalg.norm(Z_corr, ord=2, axis=0)
    return Z_cos, Z_corr, W, Phi_Rk

5.1 Scale high var genes+ PCA

adata = r.adata


# sc.pp.normalize_per_cell(adata) #Data is already normlized
adata.var['highly_variable'] = adata.var['vst.variable']
adata = adata[:,adata.var['highly_variable']].copy()

sc.pp.scale(adata, zero_center=False, max_value=50)
sc.tl.pca(adata, use_highly_variable=True)
sc.pl.pca_variance_ratio(adata, log=True)

#Run harmony

harmony_res = harmonypy.run_harmony(adata.obsm['X_pca'], adata.obs, 'orig.ident',plot_convergence = True, theta=3.5)

5.2 Adujst correction

X = np.array(adata[:,adata.var['highly_variable']].X)


_, X_corr, _, _ = moe_correct_ridge(X.T, None, None, harmony_res.R, None, harmony_res.K,
                                            None, harmony_res.Phi_moe, harmony_res.lamb)
X_corr = X_corr.T

5.3 non negative correction

adata_corr = adata[:,adata.var['highly_variable']].copy()
  
X_corr_nonneg = X_corr.copy()
X_corr_nonneg[X_corr_nonneg<0]= 0
adata_corr_nonneg = adata_corr.copy()
adata_corr_nonneg.X = X_corr_nonneg

5.4 import to seurat for plotting UMAP

exprs <- t(py$adata_corr_nonneg$X)
colnames(exprs) <- py$adata$obs_names$to_list()
rownames(exprs) <- py$adata$var_names$to_list()

# Create the Seurat object
lung_corr_nonneg <- CreateSeuratObject(counts = exprs)
lung_corr_nonneg$orig.ident <- py$adata$obs["orig.ident"]

5.5 PCA after correction

lung_corr_nonneg@assays$RNA@scale.data = exprs
lung_corr_nonneg <- RunPCA(lung_corr_nonneg, features = rownames(lung_corr_nonneg),verbose = F)

ElbowPlot(lung_corr_nonneg)

5.6 Run dim reduction

5.7 UMAP after correction

DimPlot(lung_corr_nonneg, reduction = "umap",group.by = "orig.ident")

Warning in grSoftVersion() :
  unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
  libXt.so.6: cannot open shared object file: No such file or directory

DimPlot(lung_corr_nonneg, reduction = "umap")

sc.write('./Data/cNMF/acc1CancerCells_Harmony_NoNeg.h5ad', adata_corr_nonneg)

cNMF:

name = 'HMSC_cNMF_harmony_2Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/acc1CancerCells_Harmony_NoNeg.h5ad'
tpm_fn = counts_fn


cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, n_iter=10, seed=14,tpm_fn=None,densify=True)

cnmf_obj.factorize(worker_i=0, total_workers=1)

cnmf_obj.combine()
cnmf_obj.k_selection_plot()

knitr::include_graphics("./Data/cNMF/HMSC_cNMF_harmony_2Kvargenes/HMSC_cNMF_harmony_2Kvargenes.k_selection.png")

5.8 Save object

#import pickle
#f = open('./Data/cNMF/HMSC_cNMF_harmony_2Kvargenes/cnmf_obj.pckl', 'wb')
#pickle.dump(cnmf_obj, f)
#f.close()

5.9 Load object

from cnmf import cNMF
import pickle
nfeatures = "2K"
f = open('./Data/cNMF/HMSC_cNMF_harmony_2Kvargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

selected_k = 3
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

gep_scores = py$gep_scores
gep_tpm = py$gep_tpm
all_metagenes= py$usage_norm

6 Harmony results

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}
print_tab(plt = 
            FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes),combine = T),
          title = "metagenes expression")

metagenes expression

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  lung_corr_nonneg = AddMetaData(object = lung_corr_nonneg,metadata = metage_metadata)
}
print_tab(plt = 
            FeaturePlot(object = lung_corr_nonneg,features = colnames(all_metagenes),combine = T),
          title = "metagenes expression")

metagenes expression

NA

all_markers = FindAllMarkers(object = lung_corr_nonneg,features = rownames(lung_corr_nonneg),logfc.threshold = 0,min.pct = 0,densify = T)

print_tab(plt = DimPlot(lung_corr_nonneg, reduction = "umap"),title = "UMAP")

UMAP

all_markers = all_markers %>% 
  mutate(fdr = p.adjust(p_val,method = "fdr"))%>%  #add fdr 
  dplyr::filter(fdr<0.05) %>%   dplyr::filter(abs(avg_log2FC)>1)


canonical_pathways = msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::filter(gs_subcat != "CGP") %>%  dplyr::distinct(gs_name, gene_symbol) 

plt_list = list()

for (cluster_num in unique(all_markers$cluster)) {
   deg = all_markers %>% dplyr::filter(cluster == cluster_num)
  
  p = enrichment_analysis(markers = deg, background = rownames(lung_corr_nonneg), ident.1 =  paste("cluster ", 
                                                                                                   cluster_num),
                          ident.2 =paste("cluster ", cluster_num),fdr_Cutoff = 0.05,custom_pathways = canonical_pathways)
   plt_list[[cluster_num]] = p
}

print_tab(plt = ggarrange(plotlist = plt_list),title = " clusters enrichment")

clusters enrichment

NA

print_tab(plt = DimPlot(acc1_cancer_cells,pt.size = 1,group.by = "orig.ident")
          ,title = "umap by plate")

nmf_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("metagene.1","metagene.2", "orig.ident"))
# myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )

plt = ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.1",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))+

ggboxplot(nmf_vs_plate, x = "orig.ident", y = "metagene.2",
          palette = "jco",
          add = "jitter")+
  stat_compare_means(method = "wilcox.test",comparisons = list(c("ACC.plate2","ACC1.P3")))

print_tab(plt = plt, title = "plates diff")

plates diff

NA

acc1_cancer_cells = SetIdent(object = acc1_cancer_cells, value = "orig.ident")
print_tab(plt = 
            DotPlot(object = acc1_cancer_cells, features = c("metagene.1","metagene.2","metagene.3","metagene.4"),scale = F)
          ,title = "dot plot")

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),100) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T)
   
  plt_list[[i]] = res$plt
}
print_tab(plt = ggarrange(plotlist = plt_list),title = " geps enrichment")