Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8+ T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.

Code ↔ Paper

12 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 12 matches

[1] § Methods › Pathway enrichment and gene module score analysis ↔ DESeq2_prot_clin.trial_bulk.R, lines 580–665 · score 0.99 · ATP6V0D1, ATP6V1B2, ATP6V1G1, FCGR1A, P2RX7, CORO1A
[2] § Methods › Pathway enrichment and gene module score analysis ↔ V5_prot_clin_trial.R, lines 1514–1594 · score 0.99 · ATP6V0D1, ATP6V1B2, ATP6V1G1, FCGR1A, P2RX7, CORO1A
[3] § Results › Nasal Protollin reduces the classical monocyte pro-inflammatory/apoptotic signature in AD subjects ↔ figs_protollin_clin_trial_paper.R, lines 61–104 · score 0.84 · IFI44L, S100A12, S100A9, IL1B, ISG15, CCL5
[4] § Methods › Single-cell RNA-seq analysis ↔ V5_prot_clin_trial.R, lines 559–595 · score 0.83 · FindClusters, FindNeighbors, dimensional reduction, UMAP, Seurat, resolution
[5] § Results › Monocytes from Protollin-treated subjects have an increased phagocytic signature ↔ DESeq2_prot_clin.trial_bulk.R, lines 580–665 · score 0.77 · FCGR1A, P2RX7, CORO1A, ROS, CD36, CD93
[6] § Results › Monocytes from Protollin-treated subjects have an increased phagocytic signature ↔ V5_prot_clin_trial.R, lines 1514–1594 · score 0.77 · FCGR1A, P2RX7, CORO1A, ROS, CD36, CD93
[7] § Results › Nasal Protollin downregulates CD8 + T cell cytotoxicity-related genes ↔ figs_protollin_clin_trial_paper.R, lines 271–333 · score 0.73 · MAP3K8, effector CD8, CD81, CEBPB, TNFAIP3, CD74
[8] § Results › Protollin reduces the expression of costimulatory molecules in classical monocytes, which are related to T cell-myeloid interactions ↔ figs_protollin_clin_trial_paper.R, lines 271–333 · score 0.68 · CLEC2B, nonclassical monocytes, KLRF1, CD58, CD28, CD55
[9] § Results › In vitro treatment of AD monocytes with Protollin reverses their pro-inflammatory and apoptotic-related profile ↔ DESeq2_prot_pien.R, lines 538–602 · score 0.68 · PPP2CA, FADD, ROCK1, THBD, TNFRSF21, TRIM4
[10] § Results › In vitro treatment of AD monocytes with Protollin reverses their pro-inflammatory and apoptotic-related profile ↔ DESeq2_prot_pien.R, lines 372–441 · score 0.66 · G0S2, IL1B, HC monocytes, NFKBID, JUNB, endocytic
[11] § Methods › Differential gene expression analysis using bulk RNA-seq ↔ DESeq2_prot_pien.R, lines 768–854 · score 0.57 · variable genes, DESeq2, outliers, distance, Transcript, clustering
[12] § Methods › Cell–cell communication inference using single-cell data ↔ V5_prot_clin_trial.R, lines 2351–2414 · score 0.56 · ligand receptor, expressed genes, interactions, cell

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 2,957 lines · 116 KB · no license · 4 matches

# Calling libraries
# Seurat v5
library(openxlsx)
library(Seurat)
library(harmony)
library(reticulate)
library(BPCells)
library(dplyr)
#library(SeuratData)
library(SeuratWrappers)
library(Azimuth)
library(Matrix)
library(sctransform)
library(car)
library(scater)
library(ggplot2)
library(patchwork)
library(ggrepel)
options(future.globals.maxSize = 3e+09)
options(Seurat.object.assay.version = "v5")
library(BiocParallel)
register(MulticoreParam(14))
#library(CellChat)
#library(SeuratDisk)
######################### BPCells ###################
################### First try ----->>> The one that worked
##################################################### PBMC
# # Set working dir
# setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/PBMC/")
#
# ## Loop through h5 files and output BPCells matrices on-disk
# file.dir <- "../../../cellranger/"
#
# files.set <- c(
# "GEM5",
# "GEM6",
# "GEM7",
# "GEM8",
# "GEM13",
# "GEM14",
# "GEM15",
# "GEM16",
# "GEM20",
# "GEM21",
# "GEM22"
# )
#
# data.list <- c()
#
# for (i in 1:length(files.set)) {
# ## Create binary matrices
# path <- paste0(file.dir, files.set[i], "/outs/filtered_feature_bc_matrix.h5")
# data <- open_matrix_10x_hdf5(path = path)
# write_matrix_dir(mat = data, dir = paste0(files.set[i], "_BP2"))
# ## Load in BP matrices
# mat <- open_matrix_dir(dir = paste0(files.set[i], "_BP2"))
# mat <- Azimuth:::ConvertEnsembleToSymbol(mat = mat, species = "human")
# dataset_name <- files.set[i]
# mat2 <- CreateSeuratObject(counts = mat, min.features = 500) %>% PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
# subset(subset = nFeature_RNA < 5000 & percent.mt < 20 #& ("FCGR3A" > 0 | "CD14" > 0)
# ) %>% NormalizeData(normalization.method = "LogNormalize") %>%
# FindVariableFeatures(selection.method = "vst")
# mat2$GEM <- dataset_name
# data.list[[i]] <- mat2
# rm(mat)
# rm(mat2)
# rm(data)
# }
# # Name layers
# names(data.list) <- files.set
#
# # Merge layers and create seurat obj during merging
# features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 2000)
# prot.combined <- merge(data.list[[1]], y = data.list[2:length(data.list)],
# add.cell.ids = files.set, merge.data = T)
# #DefaultAssay(prot.combined) <- "RNA"
# VariableFeatures(prot.combined) <- features
# #prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$GEM)
#
# # Normalize and scale merged obj
# prot.combined <- NormalizeData(prot.combined)
# prot.combined <- FindVariableFeatures(prot.combined)
# prot.combined <- ScaleData(prot.combined)
#
# ## Dimensionality reduction and integration
# prot.combined <- RunPCA(prot.combined)
# gc()
# ElbowPlot(prot.combined, ndims = 50)
# prot.combined <- FindNeighbors(prot.combined, dims = 1:30, reduction = "pca")
# prot.combined <- FindClusters(prot.combined, resolution = 0.5, cluster.name = "unintegrated_clusters")
#
# prot.combined <- RunUMAP(prot.combined, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
# gc()
#
#
#
#
#
#
#
#
#
#
#
# #################################################### Monocytes
# # Set working dir
# setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/Mono/")
#
# ## Loop through h5 files and output BPCells matrices on-disk
# file.dir <- "../../../cellranger/"
#
# files.set <- c(
# "GEM1",
# "GEM2",
# "GEM3",
# "GEM4",
# "GEM9",
# "GEM10",
# "GEM11",
# "GEM12",
# "GEM17",
# "GEM18",
# "GEM19"
# )
#
# data.list <- c()
#
# for (i in 1:length(files.set)) {
# ## Create binary matrices
# path <- paste0(file.dir, files.set[i], "/outs/filtered_feature_bc_matrix.h5")
# data <- open_matrix_10x_hdf5(path = path)
# write_matrix_dir(mat = data, dir = paste0(files.set[i], "_BP2"))
# ## Load in BP matrices
# mat <- open_matrix_dir(dir = paste0(files.set[i], "_BP2"))
# mat <- Azimuth:::ConvertEnsembleToSymbol(mat = mat, species = "human")
# dataset_name <- files.set[i]
# mat2 <- CreateSeuratObject(counts = mat, min.features = 500) %>% PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
# subset(subset = nFeature_RNA < 5000 & percent.mt < 20 #& ("FCGR3A" > 0 | "CD14" > 0)
# ) %>% NormalizeData(normalization.method = "LogNormalize") %>%
# FindVariableFeatures(selection.method = "vst")
# mat2$GEM <- dataset_name
# data.list[[i]] <- mat2
# rm(mat)
# rm(mat2)
# rm(data)
# }
# # Name layers
# names(data.list) <- files.set
#
# # Merge layers and create seurat obj during merging
# features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 2000)
# prot.combined <- merge(data.list[[1]], y = data.list[2:length(data.list)],
# add.cell.ids = files.set, merge.data = T)
# #DefaultAssay(prot.combined) <- "RNA"
# VariableFeatures(prot.combined) <- features
# #prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$GEM)
#
# # Normalize and scale merged obj
# prot.combined <- NormalizeData(prot.combined)
# prot.combined <- FindVariableFeatures(prot.combined)
# prot.combined <- ScaleData(prot.combined)
#
# ## Dimensionality reduction and integration
# prot.combined <- RunPCA(prot.combined)
# gc()
# ElbowPlot(prot.combined, ndims = 50)
# prot.combined <- FindNeighbors(prot.combined, dims = 1:30, reduction = "pca")
# prot.combined <- FindClusters(prot.combined, resolution = 0.9, cluster.name = "unintegrated_clusters")
#
# prot.combined <- RunUMAP(prot.combined, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
# gc()
#
# ## add metadata
# meta.data <- read.csv("meta.txt", header = T)
# prot.combined$patient <- prot.combined$GEM
# prot.combined$day <- prot.combined$GEM
# prot.combined$group <- prot.combined$GEM
#
# for (i in 1:length(meta.data$GEM)) {
# prot.combined$patient <- recode(prot.combined$patient, "meta.data$GEM[i] = meta.data$Patient[i]")
# prot.combined$day <- recode(prot.combined$day, "meta.data$GEM[i] = meta.data$Day[i]")
# prot.combined$group <- recode(prot.combined$group, "meta.data$GEM[i] = meta.data$Group[i]")
# }
# prot.combined$patient_day <- paste(prot.combined$patient, prot.combined$day, sep = "_")
#
#
#
# ## Save and Load data
# saveRDS(object = prot.combined, file = "obj_BP2_unintegrated_mono.Rds")
# prot.combined <- readRDS("./obj_BP2_unintegrated_mono.Rds")
#
# DimPlot(prot.combined, reduction = "umap.unintegrated", raster = F,
# ncol = 3,
# label = T,
# #group.by = "GEM"
# split.by = "patient_day"
# )#, combine = F)
#
# FeaturePlot(prot.combined, features = "IL1B",#c("CD14","FCGR3A"),
# pt.size = 0.1,
# #ncol = 3,
# reduction = "umap.unintegrated", raster = F#,
# #split.by = "GEM"
# )
#
#
# ################################################### PBMC and Monocytes Prot treated samples day 0 and 15
# # Set working dir
# setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/prot_treated/")
#
# ## Loop through h5 files and output BPCells matrices on-disk
# file.dir <- "../../../cellranger/"
#
# files.set <- c(
# "GEM1",
# "GEM10",
# "GEM11",
# "GEM12",
# "GEM13",
# "GEM14",
# "GEM15",
# "GEM16",
# "GEM2",
# "GEM3",
# "GEM4",
# "GEM5",
# "GEM6",
# "GEM7",
# "GEM8",
# "GEM9"
# )
#
# data.list <- c()
#
# for (i in 1:length(files.set)) {
# ## Create binary matrices
# # path <- paste0(file.dir, files.set[i], "/outs/filtered_feature_bc_matrix.h5")
# # data <- open_matrix_10x_hdf5(path = path)
# # write_matrix_dir(mat = data, dir = paste0(files.set[i], "_BP2"))
# ## Load in BP matrices
# mat <- open_matrix_dir(dir = paste0(files.set[i], "_BP2"))
# mat <- Azimuth:::ConvertEnsembleToSymbol(mat = mat, species = "human")
# dataset_name <- files.set[i]
# mat2 <- CreateSeuratObject(counts = mat, min.features = 500) %>% PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
# subset(subset = nFeature_RNA < 5000 & percent.mt < 20 #& ("FCGR3A" > 0 | "CD14" > 0)
# ) %>% NormalizeData(normalization.method = "LogNormalize") %>%
# FindVariableFeatures(selection.method = "vst")
# mat2$GEM <- dataset_name
# data.list[[i]] <- mat2
# rm(mat)
# rm(mat2)
# rm(data)
# }
# # Name layers
# names(data.list) <- files.set
#
# # Merge layers and create seurat obj during merging
# features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 3000)
# prot.combined <- merge(data.list[[1]], y = data.list[2:length(data.list)],
# add.cell.ids = files.set, merge.data = T)
# #DefaultAssay(prot.combined) <- "RNA"
# VariableFeatures(prot.combined) <- features
# #prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$GEM)
#
# prot.combined <- NormalizeData(prot.combined, normalization.method = "LogNormalize")
# prot.combined <- FindVariableFeatures(prot.combined, selection.method = "vst", nfeatures = 3000)
# prot.combined <- ScaleData(prot.combined, vars.to.regress = c("percent.mt","nCount_RNA"), #latent.data = "nFeature_RNA",
# model.use = "linear")#, features = all.genes)
# gc()
#
# ## Dimensionality reduction and integration
# prot.combined <- RunPCA(prot.combined)
# gc()
# ElbowPlot(prot.combined, ndims = 50)
# prot.combined <- FindNeighbors(prot.combined, dims = 1:40, reduction = "pca")
# prot.combined <- FindClusters(prot.combined, resolution = 1.2, cluster.name = "unintegrated_clusters")
#
# prot.combined <- RunUMAP(prot.combined, dims = 1:40, reduction = "pca", reduction.name = "umap.unintegrated")
# gc()
#
# ## add metadata
# meta.data <- read.csv("meta.txt", header = T)
# prot.combined$patient <- prot.combined$GEM
# prot.combined$day <- prot.combined$GEM
# prot.combined$group <- prot.combined$GEM
# prot.combined$dose <- prot.combined$GEM
#
# for (i in 1:length(meta.data$GEM)) {
# prot.combined$patient <- recode(prot.combined$patient, "meta.data$GEM[i] = meta.data$Patient[i]")
# prot.combined$day <- recode(prot.combined$day, "meta.data$GEM[i] = meta.data$Day[i]")
# prot.combined$group <- recode(prot.combined$group, "meta.data$GEM[i] = meta.data$Group[i]")
# prot.combined$dose <- recode(prot.combined$dose, "meta.data$GEM[i] = meta.data$Dose[i]")
# }
# prot.combined$patient_day <- paste(prot.combined$patient, prot.combined$day, sep = "_")
#
# DimPlot(prot.combined, reduction = "umap.unintegrated", raster = F,
# ncol = 1,
# label = T,
# repel = T,
# #group.by = "cohort2"
# #split.by = "cohort2"
# )#, combine = F)
#
# ## Save and Load data
# saveRDS(object = prot.combined, file = "obj_BP2_unintegrated_prot.Rds")
# prot.combined <- readRDS("./obj_BP2_unintegrated_prot.Rds")
#
# DimPlot(prot.combined, reduction = "umap.unintegrated", raster = F,
# #ncol = 3,
# label = T,
# #group.by = "GEM"
# #split.by = "patient_day"
# )#, combine = F)
#
# FeaturePlot(prot.combined, features = "IL1B",#c("CD14","FCGR3A"),
# pt.size = 0.1,
# #ncol = 3,
# reduction = "umap.unintegrated", raster = F#,
# #split.by = "GEM"
# )
#
#
#
#
# ##################################### PBMC and Monocytes -------- Without P11
# # Set working dir
# setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/noP11/")
#
# ## Loop through h5 files and output BPCells matrices on-disk
# file.dir <- "../../../cellranger/"
#
# files.set <- c(
# "GEM1",
# "GEM10",
# "GEM11",
# "GEM12",
# "GEM13",
# "GEM14",
# "GEM15",
# "GEM16",
# # "GEM17",
# "GEM18",
# "GEM19",
# "GEM2",
# # "GEM20",
# "GEM21",
# "GEM22",
# "GEM3",
# "GEM4",
# "GEM5",
# "GEM6",
# "GEM7",
# "GEM8",
# "GEM9"
# )
#
# data.list <- c()
#
# for (i in 1:length(files.set)) {
# ## Create binary matrices
# path <- paste0(file.dir, files.set[i], "/outs/filtered_feature_bc_matrix.h5")
# data <- open_matrix_10x_hdf5(path = path)
# write_matrix_dir(mat = data, dir = paste0(files.set[i], "_BP2"))
# ## Load in BP matrices
# mat <- open_matrix_dir(dir = paste0(files.set[i], "_BP2"))
# mat <- Azimuth:::ConvertEnsembleToSymbol(mat = mat, species = "human")
# dataset_name <- files.set[i]
# mat2 <- CreateSeuratObject(counts = mat, min.features = 500) %>% PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
# subset(subset = nFeature_RNA < 5000 & percent.mt < 20 #& ("FCGR3A" > 0 | "CD14" > 0)
# ) %>% NormalizeData(normalization.method = "LogNormalize") %>%
# FindVariableFeatures(selection.method = "vst")
# mat2$GEM <- dataset_name
# data.list[[i]] <- mat2
# rm(mat)
# rm(mat2)
# rm(data)
# }
# # Name layers
# names(data.list) <- files.set
#
# # Merge layers and create seurat obj during merging
# features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 2000)
# prot.combined <- merge(data.list[[1]], y = data.list[2:length(data.list)],
# add.cell.ids = files.set, merge.data = T)
# #DefaultAssay(prot.combined) <- "RNA"
# VariableFeatures(prot.combined) <- features
# #prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$GEM)
#
# # Normalize and scale merged obj
# prot.combined <- NormalizeData(prot.combined)
# prot.combined <- FindVariableFeatures(prot.combined)
# prot.combined <- ScaleData(prot.combined)
#
# ## Dimensionality reduction and integration
# prot.combined <- RunPCA(prot.combined)
# gc()
# ElbowPlot(prot.combined, ndims = 50)
# prot.combined <- FindNeighbors(prot.combined, dims = 1:40, reduction = "pca")
# prot.combined <- FindClusters(prot.combined, resolution = 1, cluster.name = "unintegrated_clusters")
#
# prot.combined <- RunUMAP(prot.combined, dims = 1:40, reduction = "pca", reduction.name = "umap.unintegrated")
# gc()
#
# ## add metadata
# meta.data <- read.csv("meta.txt", header = T)
# prot.combined$patient <- prot.combined$GEM
# prot.combined$day <- prot.combined$GEM
# prot.combined$group <- prot.combined$GEM
# prot.combined$dose <- prot.combined$GEM
# prot.combined$treatment <- prot.combined$GEM
#
# for (i in 1:length(meta.data$GEM)) {
# prot.combined$patient <- recode(prot.combined$patient, "meta.data$GEM[i] = meta.data$Patient[i]")
# prot.combined$day <- recode(prot.combined$day, "meta.data$GEM[i] = meta.data$Day[i]")
# prot.combined$group <- recode(prot.combined$group, "meta.data$GEM[i] = meta.data$Group[i]")
# prot.combined$dose <- recode(prot.combined$dose, "meta.data$GEM[i] = meta.data$Dose[i]")
# prot.combined$treatment <- recode(prot.combined$treatment, "meta.data$GEM[i] = meta.data$Treatment[i]")
# }
# prot.combined$patient_day <- paste(prot.combined$patient, prot.combined$day, sep = "_")
#
# ## Save and Load data
# saveRDS(object = prot.combined, file = "obj_BP2_unintegrated_all.Rds")
# prot.combined <- readRDS("./obj_BP2_unintegrated_all.Rds")
#
# DimPlot(prot.combined, reduction = "umap.unintegrated", raster = F,
# #ncol = 3,
# label = T,
# #group.by = "GEM"
# #split.by = "treatment"
# )#, combine = F)
#
# FeaturePlot(prot.combined, features = c("CD14","FCGR3A"),
# pt.size = 0.1,
# #ncol = 3,
# reduction = "umap.unintegrated", raster = F#,
# #split.by = "GEM"
# )
# classic.mono <- c("CD14","CCR2","CCR5","SELL")
# nonclassic.mono <- c("FCGR3A","CX3CR1","HLA-DRA")
# inter.mono <- c("CD14","HLA-DRA","ITGAX","CD68")
# Macrophages <- c("CD14", "FCGR3A", "CD64", "CD68", "CD71", "CCR5")
# bcells <- c("CD19","IGKC","IGHM","CD27", "CD1D", "CD22","CD86", "MS4A1", "IGLC2", "IGLL5", "IGLC7", "IGLC3")
# dc <- c("ITGAX","FCER1A", "CCR7","CD1C","NRP1")
# modc <- c("FCER1A","ZBTB46", "IRF4","CD1C","KLF4","ITGAM","SIRPA","MRC1")
# baso <- c("ITGB2","PECAM1","IL3RA","LAMP1")
# platelets <- c("CD41","CD42b","CD61","CD31","PPBP","PF4","GNG11","SDPR","CLU","CD41","CD110")
# all.markers <- c("CD14","CCR2","CCR5","SELL","FCGR3A","CX3CR1","HLA-DRA","ITGAX","CD68","FCER1A","CD1C")
#
# VlnPlot(prot.combined, features = "CD48"#Macrophages
# , pt.size = 0.1, ncol = 1, raster = F, split.by = "day"#, idents = "6"
# #, combine = F
# )
#
# prot.combined <- JoinLayers(prot.combined)
# prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$patient)
# prot.combined <- RunHarmony(prot.combined, group.by.vars = "patient")
# gc()
#
#
# zk.response0 <- FindMarkers(myeloids, ident.1 = "5_Progressor",
# ident.2 = c("5_RRMS","5_Nonprogressor","5_HC"),#NULL,
# slot = "data",
# assay = "RNA",
# features = NULL,
# logfc.threshold = 0,
# test.use = "wilcox",
# min.pct = 0.5,
# min.diff.pct = -Inf,
# verbose = TRUE,
# only.pos = F,
# max.cells.per.ident = Inf,
# random.seed = 1,
# latent.vars = NULL,
# min.cells.feature = 3,
# min.cells.group = 3,
# pseudocount.use = 1,
# mean.fxn = NULL,
# fc.name = NULL,
# base = 2,
# densify = FALSE,
# recorrect_umi = TRUE
# )
# write.xlsx(as.data.frame(zk.response0), rowNames = T, file="wilcox_5_progx5_RRMS-HC_Nonprog_DEGs.xlsx")
# rm(zk.response0)
#
#
# prot.combined <- JoinLayers(prot.combined)
# Idents(prot.combined) <- "seurat_clusters"
# #Idents(prot.combined) <- "disease.state"
#
# #prot.markers <- FindAllMarkers(prot.combined, assay = "RNA", slot = "data", only.pos = T, min.pct = 0.3) %>% group_by(cluster)
# #write.csv(as.data.frame(prot.markers), file="All_clus_markers.csv")
################### PBMC and Monocytes
# Set working dir
setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/")
## Loop through h5 files and output BPCells matrices on-disk
file.dir <- "../../cellranger/"
files.set <- c(
"GEM1",
"GEM10",
"GEM11",
"GEM12",
"GEM13",
"GEM14",
"GEM15",
"GEM16",
"GEM17",
"GEM18",
"GEM19",
"GEM2",
"GEM20",
"GEM21",
"GEM22",
"GEM3",
"GEM4",
"GEM5",
"GEM6",
"GEM7",
"GEM8",
"GEM9"
)
data.list <- c()
for (i in 1:length(files.set)) {
## Create binary matrices
# path <- paste0(file.dir, files.set[i], "/outs/filtered_feature_bc_matrix.h5")
# data <- open_matrix_10x_hdf5(path = path)
# write_matrix_dir(mat = data, dir = paste0(files.set[i], "_BP2"))
## Load in BP matrices
mat <- open_matrix_dir(dir = paste0(files.set[i], "_BP2"))
mat <- Azimuth:::ConvertEnsembleToSymbol(mat = mat, species = "human")
dataset_name <- files.set[i]
mat2 <- CreateSeuratObject(counts = mat, min.features = 500) %>% PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
subset(subset = nFeature_RNA < 5000 & percent.mt < 20 #& ("FCGR3A" > 0 | "CD14" > 0)
) %>% NormalizeData(normalization.method = "LogNormalize") %>%
FindVariableFeatures(selection.method = "vst")
mat2$GEM <- dataset_name
data.list[[i]] <- mat2
rm(mat)
rm(mat2)
rm(data)
}
# Name layers
names(data.list) <- files.set
# Merge layers and create seurat obj during merging
features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 2000)
prot.combined <- merge(data.list[[1]], y = data.list[2:length(data.list)],
add.cell.ids = files.set, merge.data = T)
#DefaultAssay(prot.combined) <- "RNA"
VariableFeatures(prot.combined) <- features
#prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$GEM)
# # Visualize QC metrics as a violin plot
# VlnPlot(prot.combined, pt.size = 0.0, #group.by = "cohort",
# raster = F,
# features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# plot1 <- FeatureScatter(prot.combined, feature1 = "nCount_RNA", feature2 = "percent.mt",raster = F, pt.size = 0.05)
# plot2 <- FeatureScatter(prot.combined, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",raster = F, pt.size = 0.05)
# plot1 + plot2
# Normalize and scale merged obj
prot.combined <- NormalizeData(prot.combined)
prot.combined <- FindVariableFeatures(prot.combined)
prot.combined <- ScaleData(prot.combined)
#### Dimensionality reduction and integration
prot.combined <- RunPCA(prot.combined)
gc()
ElbowPlot(prot.combined, ndims = 50)
prot.combined <- FindNeighbors(prot.combined, dims = 1:40, reduction = "pca")
prot.combined <- FindClusters(prot.combined, resolution = 1, cluster.name = "unintegrated_clusters")
prot.combined <- RunUMAP(prot.combined, dims = 1:40, reduction = "pca", reduction.name = "umap.unintegrated")
gc()
# Visualization on UMAP
DimPlot(prot.combined, reduction = "umap.unintegrated", label = TRUE, repel = TRUE,
raster = F
#group.by = 'customclassif',
#split.by = "condition"
)
## add metadata
meta.data <- read.csv("meta.txt", header = T)
prot.combined$patient <- prot.combined$GEM
prot.combined$day <- prot.combined$GEM
prot.combined$group <- prot.combined$GEM
prot.combined$dose <- prot.combined$GEM
prot.combined$treatment <- prot.combined$GEM
prot.combined$condition <- prot.combined$GEM
for (i in 1:length(meta.data$GEM)) {
prot.combined$condition <- recode(prot.combined$condition, "meta.data$GEM[i] = meta.data$condition[i]")
# prot.combined$patient <- recode(prot.combined$patient, "meta.data$GEM[i] = meta.data$Patient[i]")
# prot.combined$day <- recode(prot.combined$day, "meta.data$GEM[i] = meta.data$Day[i]")
# prot.combined$group <- recode(prot.combined$group, "meta.data$GEM[i] = meta.data$Group[i]")
# prot.combined$dose <- recode(prot.combined$dose, "meta.data$GEM[i] = meta.data$Dose[i]")
# prot.combined$treatment <- recode(prot.combined$treatment, "meta.data$GEM[i] = meta.data$Treatment[i]")
}
prot.combined$patient_day <- paste(prot.combined$patient, prot.combined$day, sep = "_")
## Save and Load data
saveRDS(object = prot.combined, file = "obj_BP2_unintegrated_all.Rds")
prot.combined <- readRDS("./obj_BP2_unintegrated_all.Rds")
prot.combined <- JoinLayers(prot.combined)
#prot.combined[["RNA"]] <- split(prot.combined[["RNA"]], f = prot.combined$patient)
prot.combined[["RNA"]]$scale.data <- NULL
prot.combined[["RNA3"]] <- as(object = prot.combined[["RNA"]], Class = "Assay")
############## Identify cell types with scType
# lapply(c("dplyr","Seurat","HGNChelper"), library, character.only = T)
# source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/gene_sets_prepare.R");
# source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/sctype_score_.R")
#
# # DB file
# db_ = "https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/ScTypeDB_full.xlsx";
# tissue = "Immune system" # e.g. Immune system,Pancreas,Liver,Eye,Kidney,Brain,Lung,Adrenal,Heart,Intestine,Muscle,Placenta,Spleen,Stomach,Thymus
#
# # prepare gene sets
# gs_list = gene_sets_prepare(db_, tissue)
#
# # get cell-type by cell matrix
# es.max = sctype_score(scRNAseqData = prot.combined[["SCT"]]@scale.data, scaled = TRUE, gs = gs_list$gs_positive, gs2 = gs_list$gs_negative)
# # NOTE: scRNAseqData parameter should correspond to your input scRNA-seq matrix.
# # In case Seurat is used, it is either pbmc[["RNA"]]@scale.data (default), pbmc[["SCT"]]@scale.data, in case sctransform is used for normalization,
# # or pbmc[["integrated"]]@scale.data, in case a joint analysis of multiple single-cell datasets is performed.
#
# # merge by cluster
# cL_results = do.call("rbind", lapply(unique([email hidden]$seurat_clusters), function(cl){
# es.max.cl = sort(rowSums(es.max[ ,rownames([email hidden][[email hidden]$seurat_clusters==cl, ])]), decreasing = !0)
# head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum([email hidden]$seurat_clusters==cl)), 10)
# }))
# sctype_scores = cL_results %>% group_by(cluster) %>% top_n(n = 1, wt = scores)
#
# # set low-confident (low ScType score) clusters to "unknown"
# sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown"
# print(sctype_scores[,1:3])
#
# [email hidden]$customclassif = ""
# for(j in unique(sctype_scores$cluster)){
# cl_type = sctype_scores[sctype_scores$cluster==j,];
# [email hidden]$customclassif[[email hidden]$seurat_clusters == j] = as.character(cl_type$type[1])
# }
#
# # Visualization on UMAP
# DimPlot(prot.combined, reduction = "umap", label = TRUE, repel = TRUE,
# #group.by = 'customclassif',
# #split.by = "condition"
# )
######### celltypes records
library(ggstatsplot)
Idents(prot.combined) <- "celltypes"
myeloids <- subset(prot.combined, idents = c("0","4","8","13","27","11","7","22","12"))
#prot.combined2$disease.patient <- paste(prot.combined2$disease.state, prot.combined2$GEM, sep = "_")
num.cells <- as.data.frame(table(prot.combined$patient_day, prot.combined$patient_day))
num.cells <- num.cells[num.cells[,3] !=0,][,2:3]
num.cells.celltype <- as.data.frame.matrix(table(prot.combined$patient_day, prot.combined$seurat_clusters))
freq.num.cells.celltype <- num.cells.celltype / num.cells[,2]*100
tfreq.num.cells.celltype <- t(freq.num.cells.celltype)
tfreq.num.cells.celltype <- tfreq.num.cells.celltype[rowSums(tfreq.num.cells.celltype) > 0,]
write.csv(as.data.frame(freq.num.cells.celltype), file="freq.num.clusters_harmony_no.unk_patient_day.csv")
#tfreq.num.cells.celltype <- read.csv("treated_tfreq.num.cells.cluster_myeloids_patient_day.csv", header = T, row.names = 1)
f_freqs <- read.csv("freq.num.cell.types_harmony_no.unk_patient_day.csv",header = T)
cell.types <- c(
Classical.monocytes,
Intermediate.monocytes,
Nonclassical.monocytes,
iNK,
mNK,
Treg,
Memory.B.cells,
Naive.B.cells,
mature.B.cells,
Plasma.B.cells,
Platelets,
mo.DC,
pDC,
CD8..NKT.like,
Memory.T.CD8.,
Naive.CD8..T.cells,
Effector.CD8..T.cells,
Memory.T.CD4.,
Naive.CD4..T.cells
)
plt <- ggbetweenstats(data = f_freqs,
x = dose,
y = Classical.monocytes,
p.adjust.method = "none",
type = "np"
)
ggsave(filename = "Classical.monocytes_np_freq_dose.pdf",
plot = plt,
width = 4,
height = 6,
device = "pdf")
plt <- ggbetweenstats(data = f_freqs,
x = dose,
y = CD8..NKT.like,
p.adjust.method = "none",
type = "p"
)
ggsave(filename = "CD8..NKT.like_p_freq_dose.pdf",
plot = plt,
width = 4,
height = 6,
device = "pdf")
library(RColorBrewer)
paletteLength <- 40
colors <- colorRampPalette( rev(brewer.pal(11, "Set3")))(paletteLength)
n <- 60
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
#pie(rep(1,n), col=sample(col_vector, n))
barplot(tfreq.num.cells.celltype, col = sample(col_vector, n), legend.text = rownames(tfreq.num.cells.celltype),
xlim = c(0,12), main = "Cell types distribution")
## Dimensionality reduction of integrated data
prot.combined <- RunHarmony(prot.combined, group.by.vars = "patient")
gc()
ElbowPlot(prot.combined, ndims = 50, reduction = "harmony")
prot.combined <- RunUMAP(prot.combined, dims = 1:40, reduction = "harmony", reduction.name = "umap")
prot.combined <- FindNeighbors(prot.combined, reduction = "harmony", dims = 1:40)
prot.combined <- FindClusters(prot.combined, resolution = 2, cluster.name = "harmony_clusters")
gc()
no.placebo <- subset(prot.combined, subset = patient %in% c("P12","P13","P15","P16"))
DimPlot(prot.combined, reduction = "umap", raster = F,
#ncol = 4,
label = T,
group.by = "seurat_clusters"
#split.by = "patient_day"
)#, combine = F)
## Save and Load data
saveRDS(object = prot.combined, file = "obj_BP2_harmony_no.unk.Rds")
prot.combined <- readRDS("./obj_BP2_harmony_no.unk.Rds")
FeaturePlot(prot.combined, features = c("IL1B","PTGS2","CXCL8","NLRP3"),#Macrophages,#c("CD8A","CD8B"),
pt.size = 0.1,
ncol = 2,
reduction = "umap",
raster = F,
#min.cutoff = 1,
#max.cutoff = 1.5,
#split.by = "treatment"
)
FeaturePlot(prot.combined, features = c("G0S2","EGR1","SDC2","GPR183"),#Macrophages,#c("CD8A","CD8B"),
pt.size = 0.1,
ncol = 2,
reduction = "umap",
raster = F,
#min.cutoff = 1,
#max.cutoff = 1.5,
#split.by = "treatment"
)
FeaturePlot(prot.combined, features = c("RGCC","PPIF","MIR222HG","TRIB1"),#Macrophages,#c("CD8A","CD8B"),
pt.size = 0.1,
ncol = 2,
reduction = "umap",
raster = F,
#min.cutoff = 1,
#max.cutoff = 1.5,
#split.by = "treatment"
)
VlnPlot(prot.combined, features = "NCAM1"
, pt.size = 0.1, ncol = 1, raster = F#, split.by = "day"#, idents = "6"
#, combine = F
)
FeaturePlot(prot.combined, features = "IL16",
pt.size = 0.1,
#ncol = 2,
reduction = "umap", raster = F,
#split.by = "treatment"
)
cell.types <- c(
"Classical monocytes",
"Intermediate monocytes",
"Nonclassical monocytes",
"iNK",
"mNK",
"Treg",
"Memory B cells",
"Naive B cells",
#"mature B cells",
"Plasma B cells",
"Platelets",
"mo-DC",
"pDC",
"CD8+ NKT-like",
"Memory T CD8+",
"Naive CD8+ T cells",
"Effector CD8+ T cells",
"Memory T CD4+",
"Naive CD4+ T cells"
)
p2 <- VlnPlot(prot.combined, features = "TOX",#c("JUN","STAT1", "CCL3", "CCL3L1"),#c("IRF1","IFNG","IFNGR1","IFNGR2"),#c("CD8A","CD4","CD19"),#c("TMEM176A","TMEM176B"),
#split.by = "disease",
group.by = "patient_day",
#group.by = "celltypes",
pt.size = 0.05,
raster = F,
#ncol = 5,
#slot = "counts",
#add.noise = F,
#log = T,
#sort = "increasing",
idents = #c("5"),
"Effector CD8+ T cells"
#"8", # moDC
#"2", # Nonclassical
#"6", # Intermediate
# c("0","1","3","4","5","9","17"),#,"13","25"), # Classical monocytes
) + scale_y_continuous(limits = c(0.000, 4)) +
stat_summary(fun = mean, geom = "point",size = 35, colour = "black", shape = 95)
p2$layers[[2]]$aes_params$alpha <- 0.3
p2
DotPlot(prot.combined, features = c("TOX","CD8A","IRF1","GZMA","GZMH","GZMM","IFNG","TNF","PRF1","BATF","TYROBP","GSDMB"
#"FOXP3","ITGA4","HLA-DRA","TIGIT"
#"ITGAE","CCR2","TNFRSF18","GPA33","ICOS","CD4","IL2RA","SELL"
#"CCR7","RELA","CXCR3","ENTPD1","IL10","PDCD1","RELB","PECAM1","IL2RB"
),
#cols = c("blue","blue"),#"blue"),#"red"),#"green","yellow","gray","pink","brown","lightblue"),
col.max = 20, #idents = #c("Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C"),#"Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"),
#c("Classical Mono_1","Classical Mono_2","Intermediate Mono","Nonclassical Mono",
#"pDC_AD","pDC_C",
#"mo-DC_AD","mo-DC_C"
#),
idents = "Effector CD8+ T cells",
# c("NK_4_C","NK_4_AD","NK_8_C","NK_8_AD","NK_21_C","NK_21_AD",
# "CD8+ NKT-like_C", "CD8+ NKT-like_AD"#, "NK_AD", "NK_C"
# c("NK_AD","NK_C","Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C","Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"
#c("62","67","49","47"),
dot.scale = 10,
cluster.idents = F, group.by = "dose",
#scale = F,
#split.by = "disease"
) + RotatedAxis()
classic.mono <- c("CD14","CCR2","CCR5","SELL")
nonclassic.mono <- c("FCGR3A","CX3CR1","HLA-DRA")
inter.mono <- c("CD14","HLA-DRA","ITGAX","CD68")
Macrophages <- c("CD14", "FCGR3A", "CD64", "CD68", "CD71", "CCR5")
bcells <- c("CD19","IGKC","IGHM","CD27", "CD1D", "CD22","CD86", "MS4A1", "IGLC2", "IGLL5", "IGLC7", "IGLC3")
dc <- c("ITGAX","FCER1A", "CCR7","CD1C","NRP1")
modc <- c("FCER1A","ZBTB46", "IRF4","CD1C","KLF4","ITGAM","SIRPA","MRC1")
baso <- c("ITGB2","PECAM1","IL3RA","LAMP1")
platelets <- c("CD41","CD42b","CD61","CD31","PPBP","PF4","GNG11","SDPR","CLU","CD41","CD110")
all.markers <- c("CD14","CCR2","CCR5","SELL","FCGR3A","CX3CR1","HLA-DRA","ITGAX","CD68","FCER1A","CD1C")
DotPlot(prot.combined, features = c("TNF","IFNG","KLRC1","NCAM1","IL2RB","IL7R","TBX21","EOMES",# iNKs
"PRF1","GZMM","GZMH","GZMA","GZMB", # mNKs
"LILRB1", "KLRB1", "ZBTB16", # NKT-like
"CD3E","CD3D", # T cells
"CD8A","CD8B","PTPRC","CCL5", #"CD244", # T cells
"CD4","S100A4","SELL", # T cells
"FOXP3", "IL2RA", # Tregs
"TRDC","TRDV1","TRDV2","TRGC2","TRGV9","TRGC1", #gamma delta
"ITGB2","PECAM1","IL3RA","LAMP1",#, #baso
"CD64", "CD68", "CD71", "CCR5","ITGAM", # Macrophages
"CD1C","ITGAX","FCER1A","CCR7","NRP1", # DC
"CD14","FCGR3A", #monocytes
"PF4", # platelets
"CD19","IGKC","IGHM","CD27","CD1D","CD22","CD86","MS4A1","IGLC2","IGLC3","IGHD","CD79A","CD79B","AIM2", "BANK1","RALGPS2","TNFRSF13B", # B cells
"IL4R","CXCR4", "BTG1", "TCL1A", "YBX3", # Naive B cells
"COCH", "SSPN", "TEX9", "TNFRSF13C", "LINC01781", # Memory B cells
"LINC01857", # mature B cells
"IGHA2","MZB1","TNFRSF17","DERL3","TXNDC5","POU2AF1","CPNE5","NT5DC2"# plasma cells
# "IL10"#,"IL1B","IL15","IL7" #"IL2","IL3","IL6","IL4","IL22"
),
#cols = c("blue","blue"),#"blue"),#"red"),#"green","yellow","gray","pink","brown","lightblue"),
col.max = 20, #idents = #c("Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C"),#"Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"),
#c("Classical Mono_1","Classical Mono_2","Intermediate Mono","Nonclassical Mono",
#"pDC_AD","pDC_C",
#"mo-DC_AD","mo-DC_C"
#),
#idents = "CD8+ TEM",
# c("NK_4_C","NK_4_AD","NK_8_C","NK_8_AD","NK_21_C","NK_21_AD",
# "CD8+ NKT-like_C", "CD8+ NKT-like_AD"#, "NK_AD", "NK_C"
# c("NK_AD","NK_C","Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C","Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"
#c("62","67","49","47"),
dot.scale = 10,
cluster.idents = T, #group.by = "patho",
#scale = F,
#split.by = "disease"
) + RotatedAxis()
prot.combined <- RenameIdents(prot.combined,
`0` = "Classical monocytes_1",
`1` = "Naive CD4+ T cells_1",
`2` = "Memory T CD4+_1",
`3` = "mNK_1",
`4` = "Classical monocytes_2",
`5` = "Effector CD8+ T cells_1",
`6` = "Naive CD8+ T cells_1",
`7` = "Nonclassical monocytes_1",
`8` = "Classical monocytes_3",
`9` = "Naive B cells_1",
`10` = "Naive B cells_2",
`11` = "Intermediate monocytes_1",
`12` = "mo-DC_1",
`13` = "Classical monocytes_4",
#`14` = "Unk_1",
`15` = "Memory B cells_1",
`16` = "CD8+ NKT-like_1",
`17` = "Memory T CD8+_1",
`18` = "Treg_1",
#`19` = "Unk_2",
`20` = "iNKT_1",
`21` = "pDC_1",
`22` = "Nonclassical monocytes_3",
`23` = "Naive B cells_3",
`24` = "Naive B cells_4",
`25` = "mature B cells_1",
#`26` = "Unk_3",
`27` = "Classical monocytes_5",
`28` = "iNK_1",
`29` = "Memory T CD4+_1",
`30` = "CD8+ NKT-like_1",
`31` = "mNK_2",
`32` = "mNK_3",
#`33` = "Unk_4",
`34` = "Nonclassical monocytes_2",
#`35` = "Unk_5",
`36` = "Platelets_1",
`37` = "Naive B cells_5",
`38` = "mature B cells_2",
`39` = "mo-DC_2",
`40` = "Plasma B cells_1",
`41` = "pDC_2"
)
## Final annotations
Idents(prot.combined) <- "seurat_clusters"
prot.combined$cluster.dose <- paste(Idents(prot.combined), prot.combined$dose, sep = "_")
Idents(prot.combined) <- "cluster.dose"
prot.combined$clus.celltypes <- Idents(prot.combined)
prot.combined$celltypes <- sub("(.*)_.*", "\\1", Idents(prot.combined))
Idents(prot.combined) <- "celltypes"
prot.combined$celltypes.dose <- paste(prot.combined$celltypes, prot.combined$dose, sep = "_")
Idents(prot.combined) <- "celltypes.dose"
prot.combined$celltypes.patientday <- paste(prot.combined$celltypes, prot.combined$patient_day, sep = "_")
Idents(prot.combined) <- "celltypes.patientday"
prot.combined$celltypes.condition <- paste(prot.combined$celltypes, prot.combined$treatment, sep = "_")
Idents(prot.combined) <- "celltypes.condition"
prot.combined$cluster.patientday <- paste(prot.combined$new_clusters, prot.combined$patient_day, sep = "_")
Idents(prot.combined) <- "cluster.patientday"
## Save and Load data
saveRDS(object = prot.combined, file = "obj_BP2_integrated_annotated_all.Rds")
prot.combined <- readRDS("./obj_BP2_integrated_annotated_all.Rds")
prot.combined2 <- subset(prot.combined, idents = c("Classical monocytes",
"CD4+ Tcm",
"CD4+ Naive",
"mNK",
"CD8+ Teff",
"CD8+ Tcm",
"Nonclassical monocytes",
"B cells",
"Intermediate monocytes",
"mo-DC",
"CD8+ NKT-like",
"CD8+ Naive",
"Treg",
"pDC",
"iNK",
"Platelets"))
saveRDS(object = prot.combined2, file = "obj_BP2_integrated_annotated_no.unk.Rds")
prot.combined <- readRDS("./obj_BP2_integrated_annotated_no.unk.Rds")
DimPlot(prot.combined, reduction = "umap", raster = F,
#ncol = 4,
repel = T,
label = T,
group.by = "celltypes",
#split.by = "dose"
)#, combine = F)
###### Diff exp analysis
cell.types <- c("0","4","7","8","11","12","13","20")
cell.types <- c(
"Classical monocytes",
"Intermediate monocytes",
"Nonclassical monocytes",
#"iNK",
#"iNKT",
"mNK",
"Treg",
"Memory B cells",
"Naive B cells",
#"mature B cells",
#"Plasma B cells",
#"Platelets",
"mo-DC",
"pDC",
"CD8+ NKT-like",
"Memory T CD8+",
"Naive CD8+ T cells",
"Effector CD8+ T cells",
"Memory T CD4+",
"Naive CD4+ T cells"
)
up.list <- c()
down.list <- c()
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_P12_D15")
,ident.2 =
paste0(cell.types[i], "_P12_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.0,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_P12_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response1 <- zk.response0[zk.response0$avg_log2FC > 0,]
zk.response2 <- zk.response0[zk.response0$avg_log2FC < 0,]
zk.response0b <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_P13_D15")
,ident.2 =
paste0(cell.types[i], "_P13_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0b <- zk.response0b[zk.response0b$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0b), rowNames = T,file=paste0("wilcox_P13_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response1b <- zk.response0b[zk.response0b$avg_log2FC > 0,]
zk.response2b <- zk.response0b[zk.response0b$avg_log2FC < 0,]
zk.response0c <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_P15_D15")
,ident.2 =
paste0(cell.types[i], "_P15_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0c <- zk.response0c[zk.response0c$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0c), rowNames = T,file=paste0("wilcox_P15_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response1c <- zk.response0c[zk.response0c$avg_log2FC > 0,]
zk.response2c <- zk.response0c[zk.response0c$avg_log2FC < 0,]
zk.response0d <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_P16_D15")
,ident.2 =
paste0(cell.types[i], "_P16_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0d <- zk.response0d[zk.response0d$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0d), rowNames = T,file=paste0("wilcox_P16_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response1d <- zk.response0d[zk.response0d$avg_log2FC > 0,]
zk.response2d <- zk.response0d[zk.response0d$avg_log2FC < 0,]
zk.response1f <- zk.response1[!(row.names(zk.response1) %in% row.names(zk.response1b)),]
zk.response1f <- zk.response1f[!(row.names(zk.response1f) %in% row.names(zk.response1c)),]
zk.response1f <- zk.response1f[!(row.names(zk.response1f) %in% row.names(zk.response1d)),]
write.xlsx(as.data.frame(zk.response1f), rowNames = T,file=paste0("wilcox_exclusive_up_P12_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response2f <- zk.response2[!(row.names(zk.response2) %in% row.names(zk.response2b)),]
zk.response2f <- zk.response2f[!(row.names(zk.response2f) %in% row.names(zk.response2c)),]
zk.response2f <- zk.response2f[!(row.names(zk.response2f) %in% row.names(zk.response2d)),]
write.xlsx(as.data.frame(zk.response2f), rowNames = T,file=paste0("wilcox_exclusive_down_P12_D15xD0_",cell.types[i],"_DEGs.xlsx"))
zk.response1e <- zk.response1b[row.names(zk.response1b) %in% row.names(zk.response1),]
zk.response1e <- zk.response1e[row.names(zk.response1e) %in% row.names(zk.response1c),]
zk.response1e <- zk.response1e[row.names(zk.response1e) %in% row.names(zk.response1d),]
zk.response2e <- zk.response2b[row.names(zk.response2b) %in% row.names(zk.response2),]
zk.response2e <- zk.response2e[row.names(zk.response2e) %in% row.names(zk.response2c),]
zk.response2e <- zk.response2e[row.names(zk.response2e) %in% row.names(zk.response2d),]
up.list[[i]] <- row.names(zk.response1e)
down.list[[i]] <- row.names(zk.response2e)
rm(zk.response1)
rm(zk.response2)
rm(zk.response0)
rm(zk.response1b)
rm(zk.response2b)
rm(zk.response0b)
rm(zk.response1c)
rm(zk.response2c)
rm(zk.response0c)
rm(zk.response1d)
rm(zk.response2d)
rm(zk.response0d)
rm(zk.response1e)
rm(zk.response2e)
rm(zk.response1f)
rm(zk.response2f)
gc()
}
names(up.list) <- cell.types
up.list2 <- t(plyr::ldply(up.list, rbind))
colnames(up.list2) <- up.list2[1,]
up.list2 <- up.list2[-c(1), ]
write.xlsx(as.data.frame(up.list2), rowNames = F,file="wilcox_celltype_D15xD0_converged_up_per.patient_DEGs.xlsx")
names(down.list) <- cell.types
down.list2 <- t(plyr::ldply(down.list, rbind))
colnames(down.list2) <- down.list2[1,]
down.list2 <- down.list2[-c(1), ]
write.xlsx(as.data.frame(down.list2), rowNames = F,file="wilcox_celltype_D15xD0_converged_down_per.patient_DEGs.xlsx")
cell.types <- c(
"Classical monocytes",
"Nonclassical monocytes",
"Intermediate monocytes",
"mo-DC"
)
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = c(paste0(cell.types[i], "_P12_D15"),
paste0(cell.types[i], "_P13_D15")
)
,ident.2 =
c(#paste0(cell.types[i], "_P04_D15"),
#paste0(cell.types[i], "_P07_D15"),
paste0(cell.types[i], "_P12_D0"),
paste0(cell.types[i], "_P13_D0"),
#paste0(cell.types[i], "_P15_D0"),
#paste0(cell.types[i], "_P16_D0")
)
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
zk.response0 <- zk.response0[order(zk.response0$avg_log2FC),]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_P13_D15xD0_", cell.types[i], "_DEGs.xlsx"))
rm(zk.response0)
}
#### Wilcox
cell.types <- c(
"Classical monocytes",
"Intermediate monocytes",
"Nonclassical monocytes",
#"iNK",
#"iNKT",
"mNK",
"Treg",
"Memory B cells",
"Naive B cells",
#"mature B cells",
#"Plasma B cells",
"Platelets",
"mo-DC",
"pDC",
"CD8+ NKT-like",
"Memory T CD8+",
"Naive CD8+ T cells",
"Effector CD8+ T cells",
"Memory T CD4+",
"Naive CD4+ T cells"
)
Idents(prot.combined) <- "celltypes.patientday"
up.list <- c()
down.list <- c()
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_1")
,ident.2 =
paste0(cell.types[i], "_0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
zk.response1 <- zk.response0[zk.response0$avg_log2FC > 0.1,]
zk.response2 <- zk.response0[zk.response0$avg_log2FC < -0.1,]
zk.response0b <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_1.5")
,ident.2 =
paste0(cell.types[i], "_0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0b <- zk.response0b[zk.response0b$p_val_adj < 0.1,]
zk.response1b <- zk.response0b[zk.response0b$avg_log2FC > 0.1,]
zk.response2b <- zk.response0b[zk.response0b$avg_log2FC < -0.1,]
zk.response1c <- zk.response1b[row.names(zk.response1b) %in% row.names(zk.response1),]
zk.response2c <- zk.response2b[row.names(zk.response2b) %in% row.names(zk.response2),]
up.list[[i]] <- row.names(zk.response1c)
down.list[[i]] <- row.names(zk.response2c)
rm(zk.response1)
rm(zk.response2)
rm(zk.response0)
rm(zk.response1b)
rm(zk.response2b)
rm(zk.response0b)
rm(zk.response1c)
rm(zk.response2c)
}
names(up.list) <- cell.types
up.list2 <- t(plyr::ldply(up.list, rbind))
colnames(up.list2) <- up.list2[1,]
up.list2 <- up.list2[-c(1), ]
write.xlsx(as.data.frame(up.list2), rowNames = F,file="wilcox_1.5-1x0_converged_up_FC0.1_DEGs.xlsx")
names(down.list) <- cell.types
down.list2 <- t(plyr::ldply(down.list, rbind))
colnames(down.list2) <- down.list2[1,]
down.list2 <- down.list2[-c(1), ]
write.xlsx(as.data.frame(down.list2), rowNames = F,file="wilcox_1.5-1x0_converged_down_FC0.1_DEGs.xlsx")
DotPlot(no.placebo, features = up.list[[15]],
col.max = 20,
idents = "Naive CD8+ T cells",
dot.scale = 10,
cluster.idents = F,
group.by = "patient_day",
#group.by = "dose",
#scale = F,
#split.by = "disease"
) + RotatedAxis()
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_1.5")
,ident.2 =
paste0(cell.types[i], "_0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
zk.response1 <- zk.response0[zk.response0$avg_log2FC > 0,]
zk.response2 <- zk.response0[zk.response0$avg_log2FC < 0,]
write.xlsx(as.data.frame(zk.response1), rowNames = T,file=paste0("wilcox_1.5x0_", cell.types[i], "_up_DEGs.xlsx"))
rm(zk.response1)
write.xlsx(as.data.frame(zk.response2), rowNames = T,file=paste0("wilcox_1.5x0_", cell.types[i], "_down_DEGs.xlsx"))
rm(zk.response2)
rm(zk.response0)
}
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = paste0(cell.types[i], "_1")
,ident.2 =
paste0(cell.types[i], "_0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_Prot_1x0ug_", cell.types[i], "_DEGs.xlsx"))
rm(zk.response0)
}
numbers <- 0:41
numbers <- as.character(numbers)
numbers <- numbers[-grep("29",numbers)]
for (i in 1:length(numbers)) {
zk.response0 <- FindMarkers(prot.combined, ident.1 = paste0(numbers[i], "_1")
,ident.2 =
paste0(numbers[i], "_0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.1,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_Prot_1x0ug_clus_", numbers[i], "_DEGs.xlsx"))
rm(zk.response0)
}
FeaturePlot(prot.combined, features = "PDCD1",#c("CD14","FCGR3A"),
pt.size = 0.1,
#ncol = 3,
reduction = "umap", raster = F,
split.by = "dose"
)
VlnPlot(prot.combined, features = "TLR2",
#c("IL18","IL1B","NFKB1","NLRP3","PTGS2","CXCL8"),
pt.size = 0.1,
ncol = 1,
raster = F,
idents = "Classical monocytes",#c("0","8","4","13","11","7","12"),#c("CD4+ Tcm","CD4+ Naive"),
#c("mNK","CD8+ NKT-like","CD8+ Teff","CD4+ Tcm","Treg"),#c("0","5"),
split.by = "patient_day",
#group.by = "celltypes"
#, combine = F
)
phagocytosis <- c(
"FCN1",
"CD93",
"LYST",
"NCF2",
"ANXA1",
"BIN2",
"ITGAM",
"RAB14",
"CLCN3",
"PRKCD",
"CORO1A",
"RAC1",
"SPG11",
"ITGAL",
"P2RX7",
"IRF8",
"PECAM1",
"TICAM2",
"FCGR1A",
"VAV1",
"CCR2",
"NCF4",
"CD14",
"ICAM3",
"PAK1",
"ARHGAP25",
"ITGB2",
"HCK",
"CD36",
"ABL1"
)
phagocytosis2 <- c(
"CD93",
"LYST",
"CLCN3",
"CORO1A",
"ITGAL",
"P2RX7",
"PECAM1",
"FCGR1A",
"CCR2",
"ITGB2",
"HCK",
"CD36"
)
endocytosis <- c(
"CD93",
"LYST",
"SNX17",
"CLCN3",
"CORO1A",
"ITGAL",
"P2RX7",
"PECAM1",
"LIPA",
"DNAJC13",
"FCGR1A",
"CAP1",
"CCR2",
"HEATR5A",
"ASGR1",
"PYCARD",
"DENND1A",
"FKBP15",
"ITGB2",
"HCK",
"LRRK2",
"CD36",
"RIN2",
"SNX10"
)
ROS <- c(
"HVCN1",
"ATP6V0D1",
"NCF1",
"ATP6V1G1",
"ATP6V1B2",
"CYBB",
"RAC2"
)
classic.mono.up <- read.delim("classic.mono.up.txt",header = T)
classic.mono.up1 <- t(classic.mono.up[,1])
classic.mono.up1.5 <- t(classic.mono.up[,2])
classic.mono.up.f <- classic.mono.up1[classic.mono.up1 %in% classic.mono.up1.5]
write.xlsx(as.data.frame(classic.mono.up.f), rowNames = T, file="classic.mono_up_1n1.5x0_DEGs.xlsx")
toll.likes <- c("TLR1","TLR2","TLR4","TLR6","TLR7","TLR8","TLR9")
targets <- c("GZMA","GZMM","GZMH","PRF1","TYROBP","IRF1","GSDMB","IFNG","CD8A","BATF","TOX")
irf8_targets <- t(read.delim("irf8_targets.txt", header = F))
foxos <- c("FOXO1","FOXO3","FOXO4","FOXO6","FOXP2")
list2 <- t(read.delim("down-biomarkers_CD8.txt", header = F))
list1 <- t(read.csv("IL11.csv", header = F))
DotPlot(prot.combined, features = list2,#"IRF8",
#cols = c("blue", "red","yellow","green","pink"),
#col.max = 20,
#col.min = 5,
#dot.min = 3,
dot.scale = 12,
#cluster.idents = T,
group.by = "patient_day",
#scale = F,
#split.by = "disease.state",
idents = c(#"CD8+ Teff",
"CD8+ Tcm")#"mo-DC",
# c("mo-DC_HC","mo-DC_RRMS","mo-DC_Progressor","mo-DC_Nonprogressor"
# ,"mo-DC_PPMS"
# )
# c("Classical monocytes_P12_D15", "Classical monocytes_P12_D0", "Classical monocytes_P04_D15",
# "Classical monocytes_P07_D15", "Classical monocytes_P13_D0", "Classical monocytes_P15_D15",
# "Classical monocytes_P15_D0", "Classical monocytes_P16_D15", "Classical monocytes_P11_D15",
# "Classical monocytes_P16_D0", "Classical monocytes_P13_D15"
#"iNK_0", "iNK_1", "iNK_1.5"#,
#"mNK_0", "mNK_1", "mNK_1.5"#,
#"CD8+ NKT-like_0", "CD8+ NKT-like_1", "CD8+ NKT-like_1.5"#,
#"CD8+ Teff_0", "CD8+ Teff_1", "CD8+ Teff_1.5",
#"Classical monocytes_0", "Classical monocytes_1", "Classical monocytes_1.5"
#"Nonclassical monocytes_0", "Nonclassical monocytes_1", "Nonclassical monocytes_1.5"
#"CD4+ Tcm_0", "CD4+ Tcm_1", "CD4+ Tcm_1.5"#,
#"Treg_0", "Treg_1", "Treg_1.5"
#)
) + RotatedAxis()
all.markers <- c("CD14","CCR2",#"CCR5",
"SELL","FCGR3A","CX3CR1","HLA-DRA",#"ITGAX","CD68",
"FCER1A","CD1C")
DoHeatmap(prot.combined, features = targets, group.by = "",
slot = "data",
cells = 1:30000,
size = 5,
disp.max = 2.5,
disp.min = 0
)
cell.types <- c(
"Classical monocytes",
"CD4+ Tcm",
"CD4+ Naive",
"mNK",
"CD8+ Teff",
"CD8+ Tcm",
"Nonclassical monocytes",
"B cells",
"Intermediate monocytes",
"mo-DC",
"CD8+ NKT-like",
"CD8+ Naive",
"Treg",
"pDC",
"iNK",
"Platelets"
)
############################### OLeg P01
prot.combined <- readRDS("./obj_BP2_harmony_no.unk.Rds")
Idents(prot.combined) <- "patient"
prot.combined2 <- subset(prot.combined, idents = c("P12","P13","P15","P16"))
prot.combined3 <- subset(prot.combined, idents = c("P12","P13","P11","P04","P07","P16"))
prot.combined4 <- subset(prot.combined, idents = c("P12","P13"))
prot.combined5 <- subset(prot.combined, idents = c("P15","P16"))
prot.combined6 <- subset(prot.combined, idents = c("P12","P13","P11","P04","P07"))
prot.combined7 <- subset(prot.combined, idents = c("P15","P16","P04","P07"))
Idents(prot.combined) <- "celltypes"
Idents(prot.combined2) <- "celltypes"
Idents(prot.combined3) <- "celltypes"
Idents(prot.combined4) <- "celltypes"
Idents(prot.combined5) <- "celltypes"
Idents(prot.combined6) <- "celltypes"
Idents(prot.combined7) <- "celltypes"
# genes1 <- c("HLA-DRB1","CD33","GRN","MS4A6A","FCER1G","BIN1","TREM1","CCR2","PSEN1","HAVCR2","MAPT","APP","GSAP","BACE1" #,"RORA","APOE","TREM2"
# ) # classic Monocytes
genes1 <- c("HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB5","HLA-DQA1","FCER1G","CTSH","CTSB","CD33","GRN","ABI3","BCKDK","LILRB2","SPI1","PILRA","ZYX",
"NDUFS2","RIN3","MS4A6A","SCARB2","CCR2",
"PSEN1","APP","HAVCR2","SCIMP","SNX1","PICALM"
#"BIN1","TREM1","MAPT","GSAP","BACE1"#,"RORA","APOE","TREM2"
)
# genes1 <- c("HLA-DRA","HLA-DRB1","HLA-DRB5","FCER1G","CTSB","CD33","GRN","ABI3","BCKDK","LILRB2","SPI1","PILRA","ZYX",
# "NDUFS2","RIN3","MS4A6A","SCARB2","CCR2",
# "PSEN1","APP","HAVCR2","SCIMP","SNX1","PICALM"
# #"BIN1","TREM1","MAPT","GSAP","BACE1"#,"RORA","APOE","TREM2""HLA-DQB1","HLA-DQA1","CTSH",
# )
DotPlot(prot.combined6, features = rev(genes1),
cols = "RdBu",
col.max = 20,
dot.scale = 15,
# cluster.idents = T,
#scale = F,
group.by = "dose",
idents = "Classical monocytes",
) + RotatedAxis() + coord_flip()
DotPlot(prot.combined6, features = phagocytosis,
#cols = "RdBu",
col.max = 20,
dot.scale = 15,
# cluster.idents = T,
# scale = F,
group.by = "dose",
idents = "Classical monocytes",
) + RotatedAxis() + coord_flip()
my_comparisons <- c("P04_D15","P07_D15","P12_D0","P13_D0","P15_D0","P16_D0",
"P11_D15","P12_D15","P13_D15","P15_D15","P16_D15")
prot.combined$patient_day <- factor(prot.combined$patient_day, levels = my_comparisons)
DotPlot(prot.combined2, features = rev(genes1),
cols = "RdBu",
col.max = 20,
dot.scale = 15,
# cluster.idents = T,
#scale = F,
group.by = "patient_day",
idents = "Classical monocytes",
) + RotatedAxis() + coord_flip()
my_comparisons <- c("P12_D0","P13_D0","P15_D0","P16_D0",
"P12_D15","P13_D15","P15_D15","P16_D15")
my_comparisons <- c("P12_D0","P12_D15","P13_D0","P13_D15",
"P15_D0","P15_D15","P16_D0","P16_D15")
prot.combined2$patient_day <- factor(prot.combined2$patient_day, levels = my_comparisons)
DotPlot(prot.combined2, features = rev(genes1),
cols = "RdBu",
col.max = 20,
dot.scale = 15,
# cluster.idents = T,
#scale = F,
group.by = "patient_day",
idents = "Intermediate monocytes",
) + RotatedAxis() + coord_flip()
genes1 <- c("HLA-DRB1","CD33","GRN","MS4A6A","FCER1G","BIN1","TREM1","CCR2","PSEN1","HAVCR2","MAPT","APP","GSAP","BACE1" #,"RORA","APOE","TREM2"
) # classic Monocytes
FeaturePlot(prot.combined2, features = #c("HLA-DRB1","CD33","GRN"),
c(#"MS4A6A","FCER1G","BIN1"#,
#"TREM1","CCR2","PSEN1"#,
#"HAVCR2","MAPT","APP"#,
#"GSAP","BACE1" #,"RORA","APOE","TREM2"
"CCR2"
),
pt.size = 2,
#min.cutoff = 0.75,
#max.cutoff = 1,
#ncol = 3,
reduction = "umap", raster = T,
split.by = "dose"
) + theme(legend.position = "right")
# zk.combined2$celltypes.condition <- paste(zk.combined2$celltypes, zk.combined2$condition, sep = "_")
# zk.combined2$celltypes.condition <- paste(zk.combined2$customclassif, zk.combined2$condition, sep = "_")
# T.cells2$cluster.condition <- paste(Idents(T.cells2), T.cells2$condition, sep = "_")
# zk.combined2$cluster.condition <- paste(zk.combined2$seurat_clusters, zk.combined2$condition, sep = "_")
# Idents(T.cells2) <- "seurat_clusters"
# Idents(zk.combined2) <- "celltypes"
# Idents(zk.combined2) <- "celltypes.condition"
# Idents(T.cells2) <- "cluster.condition"
# Idents(zk.combined2) <- "customclassif"
# zk.response0 <- FindMarkers(prot.combined, ident.1 = "26",
# ident.2 = NULL,#c("5_RRMS","5_Nonprogressor","5_HC"),#NULL,
# slot = "data",
# assay = "RNA",
# features = NULL,
# logfc.threshold = 0.4,
# test.use = "wilcox",
# min.pct = 0.5,
# min.diff.pct = -Inf,
# verbose = TRUE,
# only.pos = F,
# max.cells.per.ident = Inf,
# random.seed = 1,
# latent.vars = NULL,
# min.cells.feature = 3,
# min.cells.group = 3,
# pseudocount.use = 1,
# mean.fxn = NULL,
# fc.name = NULL,
# base = 2,
# densify = FALSE,
# recorrect_umi = TRUE
# )
# write.xlsx(as.data.frame(zk.response0), rowNames = T, file="wilcox_clus.26_markers_DEGs.xlsx")
# rm(zk.response0)
#
# prot.combined <- JoinLayers(prot.combined)
#
# prot.markers <- FindAllMarkers(prot.combined, assay = "RNA", slot = "data", only.pos = T,
# min.pct = 0.5,
# logfc.threshold = 0.4,
# ) %>% group_by(cluster)
# write.xlsx(as.data.frame(prot.markers), rowNames = T, file="All_clus_pos_markers.xlsx")
#
#
#
# mo.dcs <- subset(prot.combined, idents = "7")
# monocytes <- subset(prot.combined, idents = c("0","1","2","3","4","5","6"#,"12","22"
# )
# #subset = (CD14 > 2 | FCGR3A > 2)
# )
# monocytes <- RenameIdents(monocytes, `0` = "Classical", `1` = "Classical", `2` = "Classical",
# `3` = "Nonclassical", `4` = "Intermediate", `5` = "Classical",
# `6` = "Classical")
#
#
# myeloids <- subset(prot.combined, idents = c("0","1","2","3","4","5","6","7"))
# myeloids <- RenameIdents(myeloids, `0` = "Classical", `1` = "Classical", `2` = "Classical",
# `3` = "Nonclassical", `4` = "Intermediate", `5` = "Classical",
# `6` = "Classical", `7` = "mo-DC")
# myeloids$celltypes <- Idents(myeloids)
#
# ## Save and Load data
# saveRDS(object = monocytes, file = "mono_BP2_unintegrated_annotated.Rds")
# saveRDS(object = mo.dcs, file = "mo-dc_BP2_unintegrated_annotated.Rds")
# saveRDS(object = myeloids, file = "myeloid_BP2_unintegrated_annotated.Rds")
#
# monocytes <- readRDS("./mono_BP2_unintegrated_annotated.Rds")
# mo.dcs <- readRDS("./mo-dc_BP2_unintegrated_annotated.Rds")
# myeloids <- readRDS("./myeloid_BP2_unintegrated_annotated.Rds")
#
# myeloids2 <- subset(myeloids, subset = disease.state %in% c("RRMS","HC","Nonprogressor","Progressor"))
#
# DimPlot(myeloids, reduction = "umap.unintegrated", raster = F,
# #ncol = 2,
# label = T,
# #group.by = "celltypes"
# split.by = "disease.state"
# )#, combine = F)
#
# FeaturePlot(myeloids, features = "IL15",#c("CD14","FCGR3A"),
# pt.size = 0.1,
# #ncol = 2,
# reduction = "umap.unintegrated", raster = F,
# #split.by = "disease_sex"
# )
#
#
# VlnPlot(myeloids, features = #ptgs,
# "CD48",
# pt.size = 0.1,
# #ncol = 5,
# raster = F,
# idents = "Classical monocytes",#c("0","5"),
# split.by = "treatment",
# #group.by = "celltypes"
# #, combine = F
# )
#
# p2 <- VlnPlot(prot.combined, features = "IL15",#c("JUN","STAT1", "CCL3", "CCL3L1"),#c("IRF1","IFNG","IFNGR1","IFNGR2"),#c("CD8A","CD4","CD19"),#c("TMEM176A","TMEM176B"),
# #split.by = "disease",
# group.by = "dose",
# #group.by = "celltypes",
# pt.size = 0.05,
# raster = F,
# #ncol = 5,
# #slot = "counts",
# #add.noise = F,
# #log = T,
# #sort = "increasing",
# idents = c("Classical monocytes"),
# #"Intermediate monocytes"
# #"8", # moDC
# #"2", # Nonclassical
# #"6", # Intermediate
# # c("0","1","3","4","5","9","17"),#,"13","25"), # Classical monocytes
# ) + scale_y_continuous(limits = c(0.000, 5.5)) +
# stat_summary(fun = mean, geom = "point",size = 30, colour = "black", shape = 95)
# p2$layers[[2]]$aes_params$alpha <- 0.05
# p2
#
# irf8_targets <- t(read.delim("irf8_targets.txt", header = F))
# foxos <- c("FOXO1","FOXO3","FOXO4","FOXO6","FOXP2")
# list2 <- t(read.delim("5progxrrms-hc-nonprog.txt", header = F))
# list1 <- t(read.csv("IL11.csv", header = F))
# DotPlot(myeloids, features = modc.up,#"IRF8",
# #cols = c("blue", "red","yellow","green","pink"),
# #col.max = 20,
# #col.min = 5,
# #dot.min = 3,
# dot.scale = 12,
# cluster.idents = T,
# #scale = F,
# #split.by = "disease.state",
# idents = #"mo-DC",
# c("mo-DC_HC","mo-DC_RRMS","mo-DC_Progressor","mo-DC_Nonprogressor"
# # ,"mo-DC_PPMS"
# )
# ) + RotatedAxis()
#
# all.markers <- c("CD14","CCR2",#"CCR5",
# "SELL","FCGR3A","CX3CR1","HLA-DRA",#"ITGAX","CD68",
# "FCER1A","CD1C")
# DoHeatmap(myeloids, features = cell.markers,
# slot = "data",
# cells = 1:30000,
# size = 5,
# disp.max = 2.5,
# disp.min = 0
# )
# cell.markers <- c("CD14","FCGR3A","FCER1A","CD1C"#,"CD68"
# )
# DoHeatmap(myeloids, features = markers[,2][!markers[,2]==""], group.by = "patho", slot = "")
#
# #Idents(prot.combined) <- "seurat_clusters"
# #Idents(prot.combined) <- "disease.state"
# myeloids$celltypes <- Idents(myeloids)
# Idents(myeloids) <- "celltypes"
# Idents(myeloids) <- "seurat_clusters"
# myeloids$celltypes.disease <- paste(Idents(myeloids), myeloids$disease, sep = "_")
# Idents(myeloids) <- "celltypes.disease"
# myeloids$celltypes.prog <- paste(Idents(myeloids), myeloids$prog, sep = "_")
# Idents(myeloids) <- "celltypes.prog"
# myeloids$celltypes.disease.state <- paste(Idents(myeloids), myeloids$disease.state, sep = "_")
# Idents(myeloids) <- "celltypes.disease.state"
# myeloids$celltypes.disease.state_sex <- paste(myeloids$celltypes, myeloids$disease.state_sex, sep = "_")
# Idents(myeloids) <- "celltypes.disease.state_sex"
#myeloids$celltypes <- sub("(.*)_.*", "\\1", Idents(myeloids))
ptgs <- c("PTGS2","PTGES2","PTGES3")
prostaglandin <- c("ALOX5",
"PTGS2","PTGIS","EDN2",
"PTGES3","PTGES","EDN1",
"DAGLB","PLA2G4A","CBR1","PTGS1","PLA2G10",
"MIF","PLA2G4F", "AKR1C3",
"TBXAS1","PTGDS","PRXL2B",
"PNPLA8","PTGES2","HPGDS",
"CD74"
)
irfs <- c("IRF1","IRF2","IRF3", "IRF4","IRF5","IRF6","IRF7","IRF8","IRF9")
atfs <- c("ATF1","ATF2","ATF3","ATF4","ATF5","ATF6","ATF7")
foxo1 <- c(
"ADIPOR1",
"CDKN1B",
"CXCR4",
"FABP4",
"IGFBP1",
"IL6",
"TNFSF10",
"AR",
"EGR1",
"FSHB",
"TXNIP",
"ANGPT2",
"EDN1",
"HYOU1",
"IRS2",
"KLF2",
"LHB",
"NEUROG3",
"NKX6-1",
"NLK",
"PDGFA",
"PDGFB",
"PRL",
"RAG1"
)
fos.act <-c(
"CASP9",
"CCK",
"CREM",
"DDIT3",
"ERCC4",
"EZR",
"FAS",
"FMO4",
"HSPH1",
"MMP1",
"MMP9",
"NEFL",
"NOS2",
"PTGS2",
"SMAD7",
"SOX7",
"SRR"
)
fos.rep <- c(
"ACTA1",
"BATF3",
"BCL2L1",
"CCK",
"CD40LG",
"CRP",
"CSTA",
"CXCL8",
"CYP1A2",
"FGFBP1",
"VEGFD",
"FOXA1",
"GSTP1",
"IL1A",
"LRIG2",
"MELTF",
"MMP1",
"MMP3",
"MMP7",
"MMP9",
"NGF",
"NPPA",
"NPY",
"NQO1",
"NTF3",
"NTS",
"OXTR",
"PCK2",
"PDHA1",
"PGR",
"PLAU",
"PLAUR",
"PTGS2",
"SMAD4",
"SPRR3",
"STAR",
"TP53"
)
### Diff exp analysis
cell.types <- c("Nonclassical","Intermediate","Classical","mo-DC")
#myeloids <- JoinLayers(myeloids)
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(myeloids, ident.1 = c(paste0(cell.types[i], "_Nonprogressor"),
paste0(cell.types[i], "_Progressor"))
,ident.2 =
c(
paste0(cell.types[i], "_HC"),
paste0(cell.types[i], "_RRMS")#,
#paste0(cell.types[i], "_Nonprogressor")
)
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.5,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_Prog_NonprogxHC_RRMS_", cell.types[i], "_DEGs.xlsx"))
rm(zk.response0)
}
zk.response0 <- FindMarkers(myeloids, ident.1 = paste0(cell.types[i], "_MS"),
ident.2 = paste0(cell.types[i], "_HC"),
slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.3,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
write.csv(as.data.frame(zk.response0), file=paste0("MAST_MSxHC_", cell.types[i], "_DEGs.csv"))
rm(zk.response0)
# ends here
## Diff exp analysis
bulk <- AggregateExpression(prot.combined, return.seurat = T, slot = "counts", assays = "RNA",
group.by = c("sex", "disease.state"))
prot.combined <- JoinLayers(prot.combined)
prot.combined <- subset(prot.combined, downsample = 10000)
prot.combined[["RNA"]]$counts <- as(object = prot.combined[["RNA"]]$counts, Class = "dgCMatrix")
prot.combined <- RunUMAP(prot.combined, reduction = "pca", dims = 1:30, reduction.name = "umap.pca")
prot.combined <- RunUMAP(prot.combined, reduction = "integrated.cca", dims = 1:30, reduction.name = "umap.cca")
prot.combined <- RunUMAP(prot.combined, assay = "RNA", reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony",
return.model = T)
prot.combined <- FindNeighbors(prot.combined, reduction = "integrated.cca", dims = 1:30)
prot.combined <- FindClusters(prot.combined, resolution = 2, cluster.name = "cca_clusters")
prot.combined <- RunUMAP(prot.combined, reduction = "integrated.cca", dims = 1:30, reduction.name = "umap.cca")
prot.combined <- DimPlot(prot.combined,reduction = "umap.cca", group.by = "GEM", combine = F)
gc()
############################# Nichenet #########################
library(nichenetr)
library(tidyverse)
### subset by highly variable genes
#DefaultAssay(prot.combined) <- "RNA"
#prot.combined <- FindVariableFeatures(prot.combined)
#prot.combined2 <- subset(prot.combined, features = prot.combined@assays[["RNA"]]@var.features)
## Save V5 as V3
saveRDS(object = prot.combined, file = "obj_BP2_harmony_no.unk.Rds")
prot.combined[["RNA"]] <- JoinLayers(prot.combined[["RNA"]])
prot.combined[["RNA"]]$scale.data <- NULL
prot.combined[["RNA3"]] <- as(object = prot.combined[["RNA"]], Class = "Assay")
prot.combined2 <- prot.combined
prot.combined2[["RNA"]] <- prot.combined2[["RNA3"]]
prot.combined2[["RNA3"]] <- NULL
prot.combined2 <- ScaleData(prot.combined2)
## Save and Load data
saveRDS(object = prot.combined, file = "obj_BP2_harmony_no.unk_v3.Rds")
#rm(prot.combined)
#prot.combined <- readRDS("./obj_harmony_patient_reg.out.mito.ncounts.v3.Rds")
gc()
### start analysis
organism = "human"
#if(organism == "human"){
# lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
# ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds"))
# weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))
#} else if(organism == "mouse"){
# lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_mouse_21122021.rds"))
# ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds"))
# weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds"))
#}
## if connection does not work, add references by --->
setwd("/media/patrick/GERVAZIO/Bioinfo/ref_data/nichenet_ref/human/")
lr_network = readRDS("lr_network_human_21122021.rds")
ligand_target_matrix = readRDS("ligand_target_matrix_nsga2r_final.rds")
weighted_networks = readRDS("weighted_networks_nsga2r_final.rds")
setwd("/media/patrick/JAMBERT/Bioinfo/weiner_lab/yota/protollin/scRNAseq_clin_trial/Seurat/BPCells/")
lr_network = lr_network %>% distinct(from, to)
weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network, by = c("from","to"))
## NK ---> Classical Monocytes
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("0","4","8","13"),#"27"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("3","16","20","28","30"#,"31","32"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## NK ---> intermediate Monocytes
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("11"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("3","16","20","28","30",#"31","32"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## NK ---> Nonclassical Monocytes
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("7","22"),#"34"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("3","16","20","28","30",#"31","32"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## NK ---> Monocytes
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("7","22","11","0","4","8","13"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("3","16","20","28","30"#,"31","32"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## Monocytes ---> CD4+ T cells
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("1","2"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("7","22","11","0","4","8","13"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## Monocytes ---> Treg
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("18"),
condition_colname = "treatment", condition_oi = "Protollin", condition_reference = "Untreated",
sender = c("7","22","11","0","4","8","13"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
#### By patient
## Monocytes ---> CD8+ T cells
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("5","6","17"),
condition_colname = "patient_day", condition_oi = "P12_D15", condition_reference = "P12_D0",
sender = c("7","22","11","0","4","8","13"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## Myeloids ---> CD4+ T cells
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("1","2","18"),
condition_colname = "patient_day", condition_oi = "P12_D15", condition_reference = "P12_D0",
sender = c("7","22","11","0","4","8","13","12","21"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## Myeloids ---> CD8+ T cells
## regular pipe
nichenet_output = nichenet_seuratobj_aggregate(
seurat_obj = prot.combined2,
receiver = c("5","6","17"),
condition_colname = "patient_day", condition_oi = "P12_D15", condition_reference = "P12_D0",
sender = c("7","22","11","0","4","8","13","12","21"
),
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks)
## plots
nichenet_output$ligand_activities
write.xlsx(as.data.frame(nichenet_output$ligand_activities), rowNames = T, file="NK_mono.classical_nichenet_ligands.xlsx")
nichenet_output$top_ligands
nichenet_output$ligand_expression_dotplot
nichenet_output$ligand_differential_expression_heatmap
nichenet_output$ligand_target_heatmap
nichenet_output$ligand_target_heatmap + scale_fill_gradient2(low = "whitesmoke",
high = "royalblue", breaks = c(0,0.0045,0.009)) +
xlab("AD response genes in NK") + ylab("Prioritized immmune cell ligands")
DotPlot(prot.combined %>% subset(idents = c("10","23","38","30")), cols = c("royalblue","royalblue"),
features = nichenet_output$top_targets %>% rev(), split.by = "disease") + RotatedAxis()
nichenet_output$ligand_activity_target_heatmap
nichenet_output$ligand_receptor_heatmap
####### Manual pipe
## receiver
receiver = c("80","18","48") ## Monocytes
expressed_genes_receiver = get_expressed_genes(receiver, prot.combined2, pct = 0.10)
background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
## sender
sender_celltypes = c("21","41","1","2","62","20","28","67","57") ## NKs
list_expressed_genes_sender = sender_celltypes %>% unique() %>% lapply(get_expressed_genes, prot.combined2, 0.10) # lapply to get the expressed genes of every sender cell type separately here
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique()
seurat_obj_receiver= subset(prot.combined2, idents = receiver)
seurat_obj_receiver = SetIdent(seurat_obj_receiver, value = seurat_obj_receiver[["disease", drop=TRUE]])
condition_oi = "AD"
condition_reference = "C"
DE_table_receiver = FindMarkers(object = seurat_obj_receiver, ident.1 = condition_oi, ident.2 = condition_reference, min.pct = 0.10) %>% rownames_to_column("gene")
geneset_oi = DE_table_receiver %>% dplyr::filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)
geneset_oi = geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]
ligands = lr_network %>% pull(from) %>% unique()
receptors = lr_network %>% pull(to) %>% unique()
expressed_ligands = intersect(ligands,expressed_genes_sender)
expressed_receptors = intersect(receptors,expressed_genes_receiver)
potential_ligands = lr_network %>% dplyr::filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique()
ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)
ligand_activities = ligand_activities %>% arrange(-aupr_corrected) %>% dplyr::mutate(rank = rank(dplyr::desc(aupr_corrected)))
#write.xlsx(as.data.frame(ligand_activities), rowNames = T, file="NK_mono.classical_nichenet_ligands.xlsx")
best_upstream_ligands = ligand_activities %>% #top_n(-30, aupr_corrected) %>%
arrange(-aupr_corrected) %>% pull(test_ligand) %>% unique()
DotPlot(prot.combined2, features = best_upstream_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis() ### Plot!!!!
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi,
ligand_target_matrix = ligand_target_matrix, n = 200) %>% bind_rows() %>% drop_na()
active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix,
cutoff = 0.33)
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev() %>% make.names()
order_targets = active_ligand_target_links_df$target %>% unique() %>% intersect(rownames(active_ligand_target_links)) %>% make.names()
rownames(active_ligand_target_links) = rownames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23
colnames(active_ligand_target_links) = colnames(active_ligand_target_links) %>% make.names() # make.names() for heatmap visualization of genes like H2-T23
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()
p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Prioritized ligands","Predicted target genes",
color = "purple",legend_position = "top",
x_axis_position = "top",legend_title = "Regulatory potential") +
theme(axis.text.x = element_text(face = "italic")) + scale_fill_gradient2(low = "whitesmoke", high = "purple", breaks = c(0,0.0045,0.0090))
p_ligand_target_network ### Plot!!!!
lr_network_top = lr_network %>% dplyr::filter(from %in% best_upstream_ligands & to %in% expressed_receptors) %>% distinct(from,to)
best_upstream_receptors = lr_network_top %>% pull(to) %>% unique()
lr_network_top_df_large = weighted_networks_lr %>% dplyr::filter(from %in% best_upstream_ligands & to %in% best_upstream_receptors)
lr_network_top_df = lr_network_top_df_large %>% spread("from","weight",fill = 0)
lr_network_top_matrix = lr_network_top_df %>% dplyr::select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df$to)
dist_receptors = dist(lr_network_top_matrix, method = "binary")
hclust_receptors = hclust(dist_receptors, method = "ward.D2")
order_receptors = hclust_receptors$labels[hclust_receptors$order]
dist_ligands = dist(lr_network_top_matrix %>% t(), method = "binary")
hclust_ligands = hclust(dist_ligands, method = "ward.D2")
order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
order_receptors = order_receptors %>% intersect(rownames(lr_network_top_matrix))
order_ligands_receptor = order_ligands_receptor %>% intersect(colnames(lr_network_top_matrix))
vis_ligand_receptor_network = lr_network_top_matrix[order_receptors, order_ligands_receptor]
rownames(vis_ligand_receptor_network) = order_receptors %>% make.names()
colnames(vis_ligand_receptor_network) = order_ligands_receptor %>% make.names()
p_ligand_receptor_network = vis_ligand_receptor_network %>% t() %>% make_heatmap_ggplot("Ligands","Receptors", color = "mediumvioletred",
x_axis_position = "top",legend_title = "Prior interaction potential")
p_ligand_receptor_network ### Plot!!!!
# DE analysis for each sender cell type
# this uses a new nichenetr function - reinstall nichenetr if necessary!
DE_table_all = Idents(prot.combined2) %>% levels() %>% intersect(sender_celltypes) %>%
lapply(get_lfc_celltype, seurat_obj = prot.combined2, condition_colname = "disease",
condition_oi = condition_oi, condition_reference = condition_reference, expression_pct = 0.10, celltype_col = NULL) %>%
reduce(full_join) # use this if cell type labels are the identities of your Seurat object -- if not: indicate the celltype_col properly
DE_table_all[is.na(DE_table_all)] = 0
# Combine ligand activities with DE information
ligand_activities_de = ligand_activities %>% dplyr::select(test_ligand, pearson) %>% dplyr::rename(ligand = test_ligand) %>% left_join(DE_table_all %>% dplyr::rename(ligand = gene))
ligand_activities_de[is.na(ligand_activities_de)] = 0
# make LFC heatmap
lfc_matrix = ligand_activities_de %>% dplyr::select(-ligand, -pearson) %>% as.matrix() %>% magrittr::set_rownames(ligand_activities_de$ligand)
rownames(lfc_matrix) = rownames(lfc_matrix) %>% make.names()
order_ligands = order_ligands[order_ligands %in% rownames(lfc_matrix)]
vis_ligand_lfc = lfc_matrix[order_ligands,]
colnames(vis_ligand_lfc) = vis_ligand_lfc %>% colnames() %>% make.names()
p_ligand_lfc = vis_ligand_lfc %>% make_threecolor_heatmap_ggplot("Prioritized ligands","LFC in Sender",
low_color = "midnightblue",mid_color = "white", mid = median(vis_ligand_lfc),
high_color = "red",legend_position = "top", x_axis_position = "top", legend_title = "LFC") +
theme(axis.text.y = element_text(face = "italic"))
p_ligand_lfc ### Plot!!!!
# change colors a bit to make them more stand out
#p_ligand_lfc = p_ligand_lfc + scale_fill_gradientn(colors = c("midnightblue","blue", "grey95", "grey99","firebrick1","red"),values = c(0,0.1,0.2,0.25, 0.40, 0.7,1), limits = c(vis_ligand_lfc %>% min() - 0.1, vis_ligand_lfc %>% max() + 0.1))
#p_ligand_lfc
# ligand activity heatmap
ligand_aupr_matrix = ligand_activities %>% dplyr::select(aupr_corrected) %>% as.matrix() %>% magrittr::set_rownames(ligand_activities$test_ligand)
rownames(ligand_aupr_matrix) = rownames(ligand_aupr_matrix) %>% make.names()
colnames(ligand_aupr_matrix) = colnames(ligand_aupr_matrix) %>% make.names()
vis_ligand_aupr = ligand_aupr_matrix[order_ligands, ] %>% as.matrix(ncol = 1) %>% magrittr::set_colnames("AUPR")
p_ligand_aupr = vis_ligand_aupr %>%
make_heatmap_ggplot("Prioritized ligands","Ligand activity", color = "darkorange",legend_position = "top",
x_axis_position = "top", legend_title = "AUPR\n(target gene prediction ability)") + theme(legend.text = element_text(size = 9))
#p_ligand_aupr ## Plot!!
# ligand expression Seurat dotplot
order_ligands_adapted <- str_replace_all(order_ligands, "\\.", "-")
rotated_dotplot = DotPlot(prot.combined2 %>% subset(seurat_clusters %in% sender_celltypes), features = order_ligands_adapted,
cols = "RdYlBu") + coord_flip() + theme(legend.text = element_text(size = 10), legend.title = element_text(size = 12)) # flip of coordinates necessary because we want to show ligands in the rows when combining all plots
figures_without_legend = cowplot::plot_grid(
p_ligand_aupr + theme(legend.position = "none", axis.ticks = element_blank()) + theme(axis.title.x = element_text()),
rotated_dotplot + theme(legend.position = "none", axis.ticks = element_blank(), axis.title.x = element_text(size = 12),
axis.text.y = element_text(face = "italic", size = 9),
axis.text.x = element_text(size = 9, angle = 90,hjust = 0)) + ylab("Expression in Sender") + xlab("") + scale_y_discrete(position = "right"),
p_ligand_lfc + theme(legend.position = "none", axis.ticks = element_blank()) + theme(axis.title.x = element_text()) + ylab(""),
p_ligand_target_network + theme(legend.position = "none", axis.ticks = element_blank()) + ylab(""),
align = "hv",
nrow = 1,
rel_widths = c(ncol(vis_ligand_aupr)+6, ncol(vis_ligand_lfc) + 7, ncol(vis_ligand_lfc) + 8, ncol(vis_ligand_target)))
legends = cowplot::plot_grid(
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_aupr)),
ggpubr::as_ggplot(ggpubr::get_legend(rotated_dotplot)),
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_lfc)),
ggpubr::as_ggplot(ggpubr::get_legend(p_ligand_target_network)),
nrow = 1,
align = "h", rel_widths = c(1.5, 1, 1, 1))
combined_plot = cowplot::plot_grid(figures_without_legend, legends, rel_heights = c(10,5), nrow = 2, align = "hv")
combined_plot
############################################ Treated patients analysis
Idents(prot.combined) <- "patient"
prot.combined2 <- subset(prot.combined, idents = c("P12","P13","P15","P16"))
Idents(prot.combined2) <- "seurat_clusters"
prot.combined2 <- JoinLayers(prot.combined2)
## Save and Load data
saveRDS(object = prot.combined2, file = "obj_BP2_unintegrated_P12.13.15.16.Rds")
prot.combined2 <- readRDS("./obj_BP2_unintegrated_P12.13.15.16.Rds")
## Save V5 as V3
prot.combined2[["RNA"]] <- JoinLayers(prot.combined2[["RNA"]])
prot.combined2[["RNA"]]$scale.data <- NULL
prot.combined2[["RNA3"]] <- as(object = prot.combined2[["RNA"]], Class = "Assay")
#prot.combined2 <- prot.combined
prot.combined2[["RNA"]] <- prot.combined2[["RNA3"]]
prot.combined2[["RNA3"]] <- NULL
## Save and Load data
saveRDS(object = prot.combined2, file = "obj_BP2_integrated_P12.13.15.16.v3.Rds")
#rm(prot.combined)
prot.combined <- readRDS("./obj_harmony_patient_reg.out.mito.ncounts.v3.Rds")
prot.combined2 <- ScaleData(prot.combined2)
gc()
Idents(prot.combined) <- "celltypes"
DimPlot(prot.combined2, reduction = "umap", raster = F,
#ncol = 4,
label = T,
repel = T,
#group.by = "seurat_clusters",
#split.by = "patient_day"
)#, combine = F)
FeaturePlot(prot.combined2, features = c("FAM13A","MEGF9","GIMAP7","S100Z"),#,"MX1","IFI44"),
pt.size = 0.1,
ncol = 2,
reduction = "umap", raster = F,
#split.by = "disease_sex"
)
DotPlot(prot.combined2, features = c("TNF","IFNG","KLRC1","NCAM1","IL2RB","IL7R","TBX21","EOMES",# iNKs
"PRF1","GZMM","GZMH","GZMA","GZMB", # mNKs
"LILRB1", "KLRB1", "ZBTB16", # NKT-like
"CD3E","CD3D", # T cells
"CD8A","CD8B","PTPRC","CCL5", #"CD244", # T cells
"CD4","S100A4","SELL", # T cells
"FOXP3", "IL2RA", # Tregs
"TRDC","TRDV1","TRDV2","TRGC2","TRGV9","TRGC1", #gamma delta
"ITGB2","PECAM1","IL3RA","LAMP1",#, #pDC
"CD64", "CD68", "CD71", "CCR5","ITGAM", # Macrophages
"CD1C","ITGAX","FCER1A","CCR7","NRP1", # DC
"CD14","FCGR3A", #monocytes
"PF4", # platelets
"CD19","IGKC","IGHM","CD27","CD1D","CD22","CD86","MS4A1","IGLC2","IGLC3","IGHD","CD79A","CD79B","AIM2", "BANK1","RALGPS2","TNFRSF13B", # B cells
"IL4R","CXCR4", "BTG1", "TCL1A", "YBX3", # Naive B cells
"COCH", "SSPN", "TEX9", "TNFRSF13C", "LINC01781", # Memory B cells
"LINC01857", # mature B cells
"IGHA2","MZB1","TNFRSF17","DERL3","TXNDC5","POU2AF1","CPNE5","NT5DC2"# plasma cells
# "IL10"#,"IL1B","IL15","IL7" #"IL2","IL3","IL6","IL4","IL22"
),
#cols = c("blue","blue"),#"blue"),#"red"),#"green","yellow","gray","pink","brown","lightblue"),
col.max = 20, #idents = #c("Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C"),#"Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"),
#c("Classical Mono_1","Classical Mono_2","Intermediate Mono","Nonclassical Mono",
#"pDC_AD","pDC_C",
#"mo-DC_AD","mo-DC_C"
#),
#idents = "CD8+ TEM",
# c("NK_4_C","NK_4_AD","NK_8_C","NK_8_AD","NK_21_C","NK_21_AD",
# "CD8+ NKT-like_C", "CD8+ NKT-like_AD"#, "NK_AD", "NK_C"
# c("NK_AD","NK_C","Classical Mono_1_AD","Classical Mono_2_AD","Classical Mono_1_C","Classical Mono_2_C","Intermediate Mono_AD","Nonclassical Mono_AD","Intermediate Mono_C","Nonclassical Mono_C"
#c("62","67","49","47"),
dot.scale = 10,
cluster.idents = T, #group.by = "patho",
#scale = F,
#split.by = "disease"
) + RotatedAxis()
prot.combined2 <- RenameIdents(prot.combined2,
`0` = "Classical monocytes_1",
`1` = "Naive CD4+ T cells_1",
`2` = "Memory T CD4+_1",
`3` = "mNK_1",
`4` = "Classical monocytes_2",
`5` = "Effector CD8+ T cells_1",
`6` = "Naive CD8+ T cells_1",
`7` = "Nonclassical monocytes_1",
`8` = "Classical monocytes_3",
`9` = "B cells_1",
`10` = "B cells_2",
`11` = "Intermediate monocytes",
`12` = "mo-DC_1",
`13` = "Classical monocytes_4",
`14` = "Unk_1",
`15` = "B cells_3",
`16` = "NKT-like CD8+_1",
`17` = "Memory T CD8+",
`18` = "Treg",
`19` = "Unk_2",
`20` = "NKT-like CD4+_1",
`21` = "pDC_1",
`22` = "Nonclassical monocytes_3",
`23` = "B cells_4",
`24` = "B cells_5",
`25` = "B cells_6",
`26` = "Unk_3",
`27` = "Classical monocytes_5",
`28` = "iNK_1",
`29` = "Memory T CD4+_2",
`30` = "NKT-like CD8+_2",
`31` = "iNK_2",
`32` = "mNK_4",
`33` = "Unk_4",
`34` = "Nonclassical monocytes_2",
`35` = "Unk_5",
`36` = "Platelets_1",
`38` = "B cells_8",
`40` = "Plasma B cells_1"
)
prot.combined2$celltypes <- sub("(.*)_.*", "\\1", Idents(prot.combined2))
Idents(prot.combined2) <- "celltypes"
cell.types <- c(
"Classical monocytes",
"Memory T CD4+",
"Naive CD4+ T cells",
"mNK",
"Effector CD8+ T cells",
"Memory T CD8+",
"Nonclassical monocytes",
"B cells",
"Intermediate monocytes",
"mo-DC",
"NKT-like CD8+",
"NKT-like CD4+",
"Naive CD8+ T cells",
"Treg",
"pDC",
"iNK",
"Platelets",
"Plasma B cells"
)
prot.combined <- subset(prot.combined2, idents = cell.types)
prot.combined2$celltypes.patient.day <- paste(Idents(prot.combined2), prot.combined2$patient_day, sep = "_")
Idents(prot.combined2) <- "celltypes.patient.day"
prot.combined2 <- JoinLayers(prot.combined2)
cell.types <- c("Nonclassical monocytes","Intermediate monocytes","Classical monocytes","mo-DC")
#myeloids <- JoinLayers(myeloids)
for (i in 1:length(cell.types)) {
zk.response0 <- FindMarkers(prot.combined2, ident.1 = paste0(cell.types[i], "_P16_D15"),
ident.2 = paste0(cell.types[i], "_P16_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.25,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_P16_D15xP16_D0_", cell.types[i], "_DEGs_pct0.25.xlsx"))
rm(zk.response0)
}
prot.combined2$cluster.patient.day <- paste(Idents(prot.combined2), prot.combined2$patient_day, sep = "_")
Idents(prot.combined2) <- "cluster.patient.day"
nums <- as.character(40:40)
#myeloids <- JoinLayers(myeloids)
for (i in 1:length(nums)) {
zk.response0 <- FindMarkers(prot.combined2, ident.1 = paste0(nums[i], "_P15_D15"),
ident.2 = paste0(nums[i], "_P15_D0")
, slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.25,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file=paste0("wilcox_P15_D15xP15_D0_clus.", nums[i], "_DEGs_pct0.25.xlsx"))
rm(zk.response0)
}
# Nonclassical <- subset(prot.combined2, idents = c("Nonclassical monocytes_P12_D0", "Nonclassical monocytes_P12_D15",
# "Nonclassical monocytes_P13_D0", "Nonclassical monocytes_P13_D15",
# "Nonclassical monocytes_P15_D0", "Nonclassical monocytes_P15_D15",
# "Nonclassical monocytes_P16_D0", "Nonclassical monocytes_P16_D15"))
# Idents(Nonclassical) <- factor(Idents(Nonclassical), levels = c("Nonclassical monocytes_P12_D0", "Nonclassical monocytes_P12_D15",
# "Nonclassical monocytes_P13_D0", "Nonclassical monocytes_P13_D15",
# "Nonclassical monocytes_P15_D0", "Nonclassical monocytes_P15_D15",
# "Nonclassical monocytes_P16_D0", "Nonclassical monocytes_P16_D15"
# ))
#
# DotPlot(Nonclassical, features = endocytosis,#"IRF8",
# #cols = c("blue", "red","yellow","green","pink"),
# #col.max = 20,
# #col.min = 5,
# #dot.min = 3,
# dot.scale = 12,
# #cluster.idents = T,
# #scale = F,
# #split.by = "disease.state",
# # idents =
# # c("Nonclassical monocytes_P12_D15", "Nonclassical monocytes_P12_D0", "Nonclassical monocytes_P13_D0", "Nonclassical monocytes_P15_D15",
# # "Nonclassical monocytes_P15_D0", "Nonclassical monocytes_P16_D15", "Nonclassical monocytes_P16_D0", "Nonclassical monocytes_P13_D15"
# # )
# ) + RotatedAxis()
######### celltypes records
library(ggstatsplot)
myeloids <- subset(prot.combined2, idents = c("0","4","8","13","27","11","7","22","12"))
#prot.combined2$disease.patient <- paste(prot.combined2$disease.state, prot.combined2$GEM, sep = "_")
num.cells <- as.data.frame(table(myeloids$patient_day, myeloids$patient_day))
num.cells <- num.cells[num.cells[,3] !=0,][,2:3]
num.cells.celltype <- as.data.frame.matrix(table(myeloids$patient_day, myeloids$seurat_clusters))
freq.num.cells.celltype <- num.cells.celltype / num.cells[,2]*100
tfreq.num.cells.celltype <- t(freq.num.cells.celltype)
tfreq.num.cells.celltype <- tfreq.num.cells.celltype[rowSums(tfreq.num.cells.celltype) > 0,]
write.csv(as.data.frame(tfreq.num.cells.celltype), file="treated_tfreq.num.cells.cluster_myeloids_patient_day.csv")
#tfreq.num.cells.celltype <- read.csv("treated_tfreq.num.cells.cluster_myeloids_patient_day.csv", header = T, row.names = 1)
library(RColorBrewer)
paletteLength <- 40
colors <- colorRampPalette( rev(brewer.pal(11, "Set3")))(paletteLength)
n <- 60
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
#pie(rep(1,n), col=sample(col_vector, n))
barplot(tfreq.num.cells.celltype, col = sample(col_vector, n), legend.text = rownames(tfreq.num.cells.celltype),
xlim = c(0,12), main = "Cell types distribution")
num.cells <- as.data.frame(table(prot.combined2$patient_day, prot.combined2$patient_day))
num.cells <- num.cells[num.cells[,3] !=0,][,2:3]
num.cells.celltype <- as.data.frame.matrix(table(prot.combined2$patient_day, prot.combined2$seurat_clusters))
freq.num.cells.celltype <- num.cells.celltype / num.cells[,2]*100
tfreq.num.cells.celltype <- t(freq.num.cells.celltype)
tfreq.num.cells.celltype <- tfreq.num.cells.celltype[rowSums(tfreq.num.cells.celltype) > 0,]
write.csv(as.data.frame(freq.num.cells.celltype), file="treated_freq.num.cells.cluster_patient_day.csv")
num.cells <- as.data.frame(table(prot.combined$patient_day, prot.combined$patient_day))
num.cells <- num.cells[num.cells[,3] !=0,][,2:3]
num.cells.celltype <- as.data.frame.matrix(table(prot.combined$patient_day, prot.combined$celltypes))
freq.num.cells.celltype <- num.cells.celltype / num.cells[,2]*100
tfreq.num.cells.celltype <- t(freq.num.cells.celltype)
tfreq.num.cells.celltype <- tfreq.num.cells.celltype[rowSums(tfreq.num.cells.celltype) > 0,]
df_freqs <- read.csv("treated_freq.num.cells.cluster_patient_day.csv",header = T)
plt <- ggbetweenstats(data = df_freqs,
x = treatment,
y = X8,
var.equal = F,
#pairwise.comparisons = T,
p.adjust.method = "none",
type = "np"
)
ggsave(filename = "treated_clus8_freq-vlnplot_np_new.pdf",
plot = plt,
width = 4,
height = 8,
device = "pdf")
nums <- c("0","1","2","3","4","5","6","8","9","17")
############# cluster markers
Idents(prot.combined2) <- "seurat_clusters"
myeloids <- subset(prot.combined2, idents = c("0","4","8","13",#"27",
"11","7","22","12"))
## Save V5 as V3
myeloids[["RNA"]] <- JoinLayers(myeloids[["RNA"]])
myeloids[["RNA"]]$scale.data <- NULL
myeloids[["RNA3"]] <- as(object = myeloids[["RNA"]], Class = "Assay")
myeloids[["RNA"]] <- myeloids[["RNA3"]]
myeloids[["RNA3"]] <- NULL
sub_myeloids <- subset(myeloids, #cells = Cells(myeloids[["RNA"]]),
downsample = 5000)
sub_myeloids <- BuildClusterTree(sub_myeloids, reduction = "umap", reorder = T)
markers <- FindAllMarkers(sub_myeloids, only.pos = T) %>% group_by(cluster) %>% dplyr::filter(avg_log2FC >1)
write.xlsx(as.data.frame(markers), rowNames = T,file="wilcox_all.markers_monocyte_clusters.xlsx")
markers %>% group_by(cluster) %>% dplyr::filter(avg_log2FC >1) %>% slice_head(n = 10) %>%
ungroup() -> top5
sub_myeloids <- ScaleData(sub_myeloids, features = top5$gene)
p <- DoHeatmap(sub_myeloids, features = top5$gene, size = 3,
cells = 1:30000) + theme(axis.text = element_text(
size = 7
)) #+ NoLegend()
p
########## Other plots
FeaturePlot(prot.combined, features = "CD48",
pt.size = 0.1,
#ncol = 2,
reduction = "umap", raster = F,
#split.by = "treatment"
)
Idents(prot.combined2) <- "seurat_clusters"
Idents(prot.combined2) <- "celltypes"
p2 <- VlnPlot(prot.combined, features = "CD244",#c("JUN","STAT1", "CCL3", "CCL3L1"),#c("IRF1","IFNG","IFNGR1","IFNGR2"),#c("CD8A","CD4","CD19"),#c("TMEM176A","TMEM176B"),
#split.by = "disease",
group.by = "patient_day",
#group.by = "celltypes",
pt.size = 0.05,
raster = F,
#ncol = 5,
#slot = "counts",
#add.noise = F,
#log = T,
#sort = "increasing",
idents = #c("5"),
"Nonclassical monocytes"
#"8", # moDC
#"2", # Nonclassical
#"6", # Intermediate
# c("0","1","3","4","5","9","17"),#,"13","25"), # Classical monocytes
) + scale_y_continuous(limits = c(0.000, 5.2)) +
stat_summary(fun = mean, geom = "point",size = 30, colour = "black", shape = 95)
p2$layers[[2]]$aes_params$alpha <- 0.1
p2
p2 <- VlnPlot(prot.combined2, features = "TBXAS1",#c("JUN","STAT1", "CCL3", "CCL3L1"),#c("IRF1","IFNG","IFNGR1","IFNGR2"),#c("CD8A","CD4","CD19"),#c("TMEM176A","TMEM176B"),
#split.by = "disease",
group.by = "prog",
pt.size = 0.05,
raster = F,
#ncol = 5,
#slot = "counts",
#add.noise = F,
#log = T,
#sort = "increasing",
idents = #"4",
#"8", # moDC
"2", # Nonclassical
#"6", # Intermediate
#c("0","1","3","4","5","9","13","17","25"), # Classical monocytes
) + scale_y_continuous(limits = c(0.000, 3.6)) +
stat_summary(fun = mean, geom = "point",size = 35, colour = "black", shape = 95)
p2$layers[[2]]$aes_params$alpha <- 0.2
p2
VlnPlot(prot.combined2, features = ccrs,
#"IRF8",
pt.size = 0.01,
ncol = 3,
raster = F,
#idents = "Classical",#c("0","5"),
split.by = "disease.state",
#group.by = "celltypes"
#, combine = F
) + theme(legend.position = 'right')
zk.response0 <- FindMarkers(prot.combined2, ident.1 = c("0","4"),
ident.2 = c("8","13","27"),
slot = "data",
assay = "RNA",
features = NULL,
logfc.threshold = 0,
test.use = "wilcox",
min.pct = 0.25,
min.diff.pct = -Inf,
verbose = TRUE,
only.pos = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
latent.vars = NULL,
min.cells.feature = 3,
min.cells.group = 3,
pseudocount.use = 1,
mean.fxn = NULL,
fc.name = NULL,
base = 2,
densify = FALSE,
recorrect_umi = TRUE
)
zk.response0 <- zk.response0[zk.response0$p_val_adj < 0.1,]
write.xlsx(as.data.frame(zk.response0), rowNames = T,file="wilcox_clus_0.4x8.13.27_DEGs_pct0.25.xlsx")
rm(zk.response0)

V5_prot_clin_trial.R at commit f781e9a, no license · at the source

Overview

Authors: Panayota Kolypetri¹, Patrick da Silva¹, Ronaldo S. Francisco Jr¹, Dan Frenkel², Rachael R. Cecere¹, Pien C. J. Kiliaan¹, Federico Montini¹, Shrishti Saxena¹, William A. Clementi³, Xuejun Liu⁴, Cheng Sun⁵, Regan W. Bergmark⁶, Tarun Singhal¹, Taylor J. Saraceno¹, Joseph Zimmermann⁷, Seth A. Gale¹, Dennis J. Selkoe¹, Tanuja Chitnis¹, Howard L. Weiner¹

ORCID iDs: Panayota Kolypetri, Shrishti Saxena, Seth A. Gale

Institutions: Brigham and Women's Hospital (United States); Tel Aviv University (Israel); I-Mab Biopharma US Limited,Rockville, MD USA; Jiangsu Nhwa Pharmaceutical (China) (China); Harvard University (United States)

Journal: npj aging, volume 12, issue 1, article 126

Dates: received 18 September 2025; accepted 28 April 2026; published online 5 June 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41514-026-00397-3 · PMID 42248874 · PMCID PMC13558635 · OpenAlex W4415973101

Open access: hybrid, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.11.05.25339599

Status: code verified

Categories: human (organism), Alzheimer's / dementia (population), clinical / translational (subfield)

Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Connectivity

Keywords: Diseases, Immunology, Neurology, Neuroscience

Topic: Alzheimer's disease research and treatments (Physiology, Medicine), according to OpenAlex

Funding: Jiangsu Nhwa Pharmaceutical Co., Ltd and IMab Pharmaceuticals; Davis Alzheimer Prevention Program; Ann Romney Center for Neurologic Diseases; NIH (NIA R21AG072486)

Citations: cited by 2 papers (Europe PMC); 76 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 12 matches between paragraphs and lines of code.

ronaldosfjunior/Kolypetri-et-al-2025

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: f781e9af849a60352d89aba78551bbf4e920d484, 12 September 2025

Languages: R (5)

Size: 6 files, 5 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: ggplot2 (5 files), tidyverse (5 files), pheatmap (4 files), car (3 files), DESeq2 (3 files), ggpubr (3 files), limma (3 files), reticulate (3 files), data.table (2 files), Harmony (2 files), patchwork (2 files), reshape2 (2 files), Seurat (2 files), circlize (1 file), cowplot (1 file), Plotly (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

6 files

DESeq2_prot_clin.trial_bulk.R, R, 910 lines, 2 matches
DESeq2_prot_mice_tissue.R, R, 820 lines
DESeq2_prot_pien.R, R, 854 lines, 3 matches
V5_prot_clin_trial.R, R, 2,957 lines, 4 matches
figs_protollin_clin_trial_paper.R, R, 1,851 lines, 3 matches
README.md, Text, 22 lines

At the source: github.com/ronaldosfjunior/Kolypetri-et-al-2025

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: ronaldosfjunior/Kolypetri-et-al-2025

Read it in the paper: doi.org/10.1038/s41514-026-00397-3.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
12 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: ronaldosfjunior/Kolypetri-et-al-2025

Read it in the paper: doi.org/10.1038/s41514-026-00397-3.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 19 authors, 4 keywords, 4 funders, 75 references.

Cite

This paper

Kolypetri, P., da Silva, P., Francisco, R. S., Frenkel, D., Cecere, R. R., Kiliaan, P. C. J., Montini, F., Saxena, S., Clementi, W. A., Liu, X., Sun, C., Bergmark, R. W., Singhal, T., Saraceno, T. J., Zimmermann, J., Gale, S. A., Selkoe, D. J., Chitnis, T., & Weiner, H. L. (2026). Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8+ T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial. npj aging, 12(1), 126. https://doi.org/10.1038/s41514-026-00397-3

BibTeX

@article{kolypetri2026nasal, author = {Kolypetri, Panayota and da Silva, Patrick and Francisco, Ronaldo S. and Frenkel, Dan and Cecere, Rachael R. and Kiliaan, Pien C. J. and Montini, Federico and Saxena, Shrishti and Clementi, William A. and Liu, Xuejun and Sun, Cheng and Bergmark, Regan W. and Singhal, Tarun and Saraceno, Taylor J. and Zimmermann, Joseph and Gale, Seth A. and Selkoe, Dennis J. and Chitnis, Tanuja and Weiner, Howard L.}, title = {{Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8\<sup\>+\</sup\> T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial}}, journal = {npj aging}, year = {2026}, month = jun, volume = {12}, number = {1}, pages = {126}, publisher = {Nature Publishing Group}, issn = {2731-6068}, doi = {10.1038/s41514-026-00397-3}, url = {https://doi.org/10.1038/s41514-026-00397-3}, pmid = {42248874}, pmcid = {PMC13558635} }

RIS

TY - JOUR AU - Kolypetri, Panayota AU - da Silva, Patrick AU - Francisco, Ronaldo S. AU - Frenkel, Dan AU - Cecere, Rachael R. AU - Kiliaan, Pien C. J. AU - Montini, Federico AU - Saxena, Shrishti AU - Clementi, William A. AU - Liu, Xuejun AU - Sun, Cheng AU - Bergmark, Regan W. AU - Singhal, Tarun AU - Saraceno, Taylor J. AU - Zimmermann, Joseph AU - Gale, Seth A. AU - Selkoe, Dennis J. AU - Chitnis, Tanuja AU - Weiner, Howard L. TI - Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8+ T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial T2 - npj aging J2 - NPJ Aging PY - 2026 DA - 2026/06/05 VL - 12 IS - 1 SP - 126 SN - 2731-6068 PB - Nature Publishing Group DO - 10.1038/s41514-026-00397-3 UR - https://doi.org/10.1038/s41514-026-00397-3 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41514-026-00397-3", "type": "article-journal", "title": "Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8+ T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial", "container-title": "npj aging", "author": [ { "family": "Kolypetri", "given": "Panayota" }, { "family": "da Silva", "given": "Patrick" }, { "family": "Francisco", "given": "Ronaldo S." }, { "family": "Frenkel", "given": "Dan" }, { "family": "Cecere", "given": "Rachael R." }, { "family": "Kiliaan", "given": "Pien C. J." }, { "family": "Montini", "given": "Federico" }, { "family": "Saxena", "given": "Shrishti" }, { "family": "Clementi", "given": "William A." }, { "family": "Liu", "given": "Xuejun" }, { "family": "Sun", "given": "Cheng" }, { "family": "Bergmark", "given": "Regan W." }, { "family": "Singhal", "given": "Tarun" }, { "family": "Saraceno", "given": "Taylor J." }, { "family": "Zimmermann", "given": "Joseph" }, { "family": "Gale", "given": "Seth A." }, { "family": "Selkoe", "given": "Dennis J." }, { "family": "Chitnis", "given": "Tanuja" }, { "family": "Weiner", "given": "Howard L." } ], "container-title-short": "NPJ Aging", "volume": "12", "issue": "1", "page": "126", "DOI": "10.1038/s41514-026-00397-3", "PMID": "42248874", "PMCID": "PMC13558635", "ISSN": "2731-6068", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41514-026-00397-3", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 5 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, limma, car, 12 other tools, 3 references
[2] doi:10.1038/s41593-026-02367-0 [code]: A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: Harmony, reticulate, limma, 11 other tools, Alzheimer's / dementia, 2 references
[3] doi:10.1038/s41593-026-02300-5 [code]: Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: Harmony, car, circlize, 11 other tools, 1 reference
[4] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, limma, circlize, 11 other tools, 1 reference
[5] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: reticulate, limma, car, 11 other tools, 1 reference
[6] doi:10.1016/j.isci.2026.115657 [code]: Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: Harmony, limma, car, 11 other tools, clinical / translational
[7] doi:10.1038/s41467-026-70232-6 [code]: Gene expression dynamics of human and mouse craniofacial development at the single-cell level.
Journal: Nature communications
In common: Harmony, reticulate, limma, 10 other tools, 2 references
[8] doi:10.1038/s41586-026-10214-2 [code]: Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, reticulate, car, 11 other tools
[9] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, reticulate, circlize, 10 other tools, 1 reference
[10] doi:10.1038/s41380-026-03629-w [code]: Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: Harmony, limma, circlize, 9 other tools, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8<sup>+</sup> T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repository

ronaldosfjunior/Kolypetri-et-al-2025

Code availability statement

Tracing map

Data

Code and data availability statement

Versions

Version 1, 27 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity

Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8&lt;sup&gt;+&lt;/sup&gt; T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.

Code ↔ Paper

Overview

Abstract

Repository

ronaldosfjunior/Kolypetri-et-al-2025

Code availability statement

Tracing map

Data

Code and data availability statement

Versions

Version 1, 27 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity

Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8<sup>+</sup> T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.