Skip to content

GitLab

Commit dba01155 authored by Stefano Beretta's avatar Stefano Beretta
Browse files

Added analysis scripts and notebooks.

parent de05f37d
Hide whitespace changes
Inline Side-by-side
SampleSheet.csv 0 → 100644
View file @ dba01155
Lane,Sample,Index,Group
*,FB-10-A4,SI-GA-G1,AB_Ctrl
*,FB-11-C0,SI-GA-G2,CD_High
*,FB-12-F1,SI-GA-G3,EF_Low
*,FB-16-B1,SI-GA-G4,AB_Ctrl
*,FB-17-E3,SI-GA-G5,EF_Low
*,FB-18-F2,SI-GA-G6,EF_Low
*,FB-4-B3,SI-GA-G7,AB_Ctrl
*,FB-5-C1,SI-GA-G8,CD_High
*,FB-6-D1,SI-GA-G9,CD_High
scVelo-Notebooks/Macro_Controls.ipynb 0 → 100644
View file @ dba01155
This source diff could not be displayed because it is too large. You can view the blob instead.
scVelo-Notebooks/Macro_IFNa.ipynb 0 → 100644
View file @ dba01155
This source diff could not be displayed because it is too large. You can view the blob instead.
scVelo-Notebooks/Macro_IFNaDHFR.ipynb 0 → 100644
View file @ dba01155
This source diff could not be displayed because it is too large. You can view the blob instead.
scripts/1_fastMNN_analysis.R 0 → 100644
View file @ dba01155
library(Seurat)
library(SeuratWrappers)
library(openxlsx)
wdir <- "Birocchi_SciTraMed2022/results"
full_min <- readRDS(paste(wdir, "Full", "Full_minimal.rds", sep = "/"))
# Compute mito %
full_min[["percent.mt"]] <- PercentageFeatureSet(full_min, pattern = "^mt-")
png(filename = paste(wdir, "Full", "01-fastMNN_2", paste("Full", "VlnQCplots.png", sep = "_"), sep = "/"),
units = "in", width = 12, height = 9, res = 100)
VlnPlot(full_min, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.2)
dev.off()
# Filter cells
full_min <- subset(x = full_min, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 10)
# Normalize
full_min <- NormalizeData(full_min)
full_min <- FindVariableFeatures(full_min, selection.method = "vst",
nfeatures = ceiling(dim(full_min)[1] * 0.2),
verbose = T)
# Run FastMNN
full_min <- RunFastMNN(object.list = SplitObject(full_min, split.by = "RNA_Group"))
full_min <- RunUMAP(full_min, reduction = "mnn", dims = 1:50)
full_min <- FindNeighbors(full_min, reduction = "mnn", dims = 1:50)
for (res in c(1.0)) {
full_min <- FindClusters(full_min, resolution = res)
}
png(filename = paste(wdir, "Full", "01-fastMNN_2", "Full_RNA_Group_res1.2_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min,
group.by = c("RNA_Group", "RNA_snn_res.1.2"),
ncol = 2, label = T)
dev.off()
png(filename = paste(wdir, "Full", "01-fastMNN_2", "Full_RNA_Group_res1.8_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min,
group.by = c("RNA_Group", "RNA_snn_res.1.8"),
ncol = 2, label = T)
dev.off()
# Compute Markers
df_markers <- list()
for (md in c("RNA_snn_res.1")) { # "RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"
Idents(full_min) <- md
png(paste(wdir, "Full", "01-fastMNN_2", paste0("Full_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = full_min, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = full_min, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "Full", "01-fastMNN_2", paste("Full_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
full_min@misc[["markers"]] <- df_markers
# Save final object
saveRDS(object = full_min, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
##########################
### Filter Macrophages ###
##########################
Idents(full_min) <- "RNA_snn_res.1"
full_min_macro <- subset(x = full_min, idents = c(0, 1, 4, 5, 12, 14, 18, 19, 22))
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "Full_Macrophages_subset_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap", group.by = "RNA_snn_res.1", label = T)
dev.off()
full_min_macro <- FindVariableFeatures(full_min_macro, selection.method = "vst",
nfeatures = ceiling(dim(full_min_macro)[1] * 0.2),
verbose = T)
full_min_macro <- RunFastMNN(object.list = SplitObject(full_min_macro, split.by = "RNA_Group"))
full_min_macro <- RunUMAP(full_min_macro, reduction = "mnn", dims = 1:30)
full_min_macro <- FindNeighbors(full_min_macro, reduction = "mnn", dims = 1:30)
for (res in c(0.6, 1.2, 1.8)) {
full_min_macro <- FindClusters(full_min_macro, resolution = res)
}
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
width = 14, height = 12, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
ncol = 2, label = T)
dev.off()
df_markers <- list()
for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(full_min_macro) <- md
png(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = full_min_macro, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = full_min_macro, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2",
paste("MacrophagesFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
full_min_macro@misc[["markers_full"]] <- df_markers
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap", group.by = c("RNA_Group", "orig.ident"), ncol = 2, label = T)
dev.off()
# cell-cycle scoring
cc.genes <- readRDS(file = "regev_lab_cell_cycle_mouse.rds")
s.genes <- cc.genes[["s.genes"]]
g2m.genes <- cc.genes[["g2m.genes"]]
full_min_macro@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
full_min_macro <- CellCycleScoring(object = full_min_macro, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
full_min_macro@meta.data$CC.Difference <- full_min_macro@meta.data$S.Score - full_min_macro@meta.data$G2M.Score
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
group.by = c("Phase"),
ncol = 2)
dev.off()
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
group.by = c("Phase"), split.by = "RNA_Group",
ncol = 3)
dev.off()
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
group.by = c("Phase"), split.by = "Phase",
ncol = 3)
dev.off()
# Save final object
saveRDS(object = full_min_macro,
file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
## Filter Monocytes
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
Idents(macro_obj) <- "RNA_snn_res.0.6"
monoc <- subset(macro_obj, idents = c(7, 10))
monoc <- FindVariableFeatures(monoc, selection.method = "vst", nfeatures = ceiling(dim(monoc)[1] * 0.2), verbose = T)
monoc <- RunFastMNN(object.list = SplitObject(monoc, split.by = "RNA_Group"))
monoc <- RunUMAP(monoc, reduction = "mnn", dims = 1:20)
monoc <- FindNeighbors(monoc, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2, 1.8)) {
monoc <- FindClusters(monoc, resolution = res)
}
png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
width = 14, height = 12, units = "in", res = 200)
DimPlot(monoc, reduction = "umap",
group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
ncol = 2, label = T)
dev.off()
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(monoc) <- md
png(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", paste0("MonocytesFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = monoc, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = monoc, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2",
paste("MonocytesFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
monoc@misc[["markers_full"]] <- df_markers
png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(monoc, reduction = "umap",
group.by = c("RNA_Group", "orig.ident"),
ncol = 2, label = T)
dev.off()
# Save object
saveRDS(object = monoc, file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/"))
# Filter out Clu 3
Idents(monoc) <- "RNA_snn_res.0.6"
monoc_filt <- subset(monoc, idents = c(0, 1, 2))
monoc_filt <- FindVariableFeatures(monoc_filt,
selection.method = "vst", nfeatures = ceiling(dim(monoc_filt)[1] * 0.2), verbose = T)
monoc_filt <- RunFastMNN(object.list = SplitObject(monoc_filt, split.by = "RNA_Group"))
monoc_filt <- RunUMAP(monoc_filt, reduction = "mnn", dims = 1:20)
monoc_filt <- FindNeighbors(monoc_filt, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2, 1.8)) {
monoc_filt <- FindClusters(monoc_filt, resolution = res)
}
png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
width = 14, height = 12, units = "in", res = 200)
DimPlot(monoc_filt, reduction = "umap",
group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
ncol = 2, label = T)
dev.off()
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(monoc_filt) <- md
png(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", paste0("MonocytesFiltered_No3_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = monoc_filt, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = monoc_filt, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2",
paste("MonocytesFiltered_No3_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
monoc_filt@misc[["markers_full"]] <- df_markers
png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(monoc_filt, reduction = "umap",
group.by = c("RNA_Group", "orig.ident"),
ncol = 2, label = T)
dev.off()
saveRDS(object = monoc_filt,
file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/"))
######################
### Filter T cells ###
######################
Idents(full_min) <- "RNA_snn_res.1"
full_min_tcell <- subset(x = full_min, idents = c(2, 3, 6, 7, 11, 17))
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Full_Tcell_subset_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_tcell, reduction = "umap", group.by = "RNA_snn_res.1", label = T)
dev.off()
full_min_tcell <- FindVariableFeatures(full_min_tcell, selection.method = "vst",
nfeatures = ceiling(dim(full_min_tcell)[1] * 0.2),
verbose = T)
full_min_tcell <- RunFastMNN(object.list = SplitObject(full_min_tcell, split.by = "RNA_Group"))
full_min_tcell <- RunUMAP(full_min_tcell, reduction = "mnn", dims = 1:30)
full_min_tcell <- FindNeighbors(full_min_tcell, reduction = "mnn", dims = 1:30)
for (res in c(0.6, 1.2, 1.8)) {
full_min_tcell <- FindClusters(full_min_tcell, resolution = res)
}
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
width = 12, height = 8, units = "in", res = 300)
FeaturePlot(full_min_tcell, reduction = "umap", features = c("Cd4", "Cd8a"))
dev.off()
for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8", "RNA_snn_res.2.2", "RNA_snn_res.2.4")) {
Idents(full_min_tcell) <- md
png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = full_min_tcell, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_vln_plot.png"), sep = "/"),
width = 12, height = 8, units = "in", res = 300)
print(VlnPlot(object = full_min_tcell, features = c("Cd4", "Cd8a"), pt.size = 0.2))
dev.off()
}
df_markers <- list()
for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(full_min_tcell) <- md
png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = full_min_macro, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = full_min_tcell, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2",
paste("TcellFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
full_min_tcell@misc[["markers_full"]] <- df_markers
saveRDS(object = full_min_tcell,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
### Subset of Cd4 and Cd8
full_min_tcell <- FindClusters(full_min_tcell, resolution = 8)
## Cd4
Idents(full_min_tcell) <- "RNA_snn_res.8"
cd4_sub <- subset(x = full_min_tcell, idents = c(0, 1, 3, 4, 5, 7, 9, 10, 13, 14, 16, 17, 18, 21, 22, 24, 27, 30, 31, 33, 34, 36, 38, 43, 44, 45, 46, 47, 48, 50, 51, 52, 53, 55, 56, 57))
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
width = 12, height = 8, units = "in", res = 300)
FeaturePlot(cd4_sub, reduction = "umap", features = c("Cd4", "Cd8a"))
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub_fastMNN_umap_plot.png", sep = "/"),
width = 12, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap", pt.size = 1.5, label = T)
dev.off()
cd4_sub <- NormalizeData(cd4_sub)
cd4_sub <- FindVariableFeatures(cd4_sub, selection.method = "vst",
nfeatures = ceiling(dim(cd4_sub)[1] * 0.2),
verbose = T)
cd4_sub <- RunFastMNN(object.list = SplitObject(cd4_sub, split.by = "RNA_Group"))
cd4_sub <- RunUMAP(cd4_sub, reduction = "mnn", dims = 1:15)
cd4_sub <- FindNeighbors(cd4_sub, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
cd4_sub <- FindClusters(cd4_sub, resolution = res)
}
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(cd4_sub) <- md
png(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", paste0("Cd4sub_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = cd4_sub, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = cd4_sub, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub",
paste("Cd4sub_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
cd4_sub@misc[["markers"]] <- df_markers
# cell-cycle scoring
cd4_sub@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
cd4_sub <- CellCycleScoring(object = cd4_sub, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
cd4_sub@meta.data$CC.Difference <- cd4_sub@meta.data$S.Score - cd4_sub@meta.data$G2M.Score
png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
group.by = c("Phase"),
ncol = 2)
dev.off()
png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
group.by = c("Phase"), split.by = "RNA_Group",
ncol = 3)
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
group.by = c("Phase"), split.by = "Phase",
ncol = 3)
dev.off()
# Save object
saveRDS(object = cd4_sub,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/"))
## Cd8
cd8_sub <- subset(x = full_min_tcell, idents = c(2, 6, 8, 11, 12, 15, 19, 20, 23, 25, 26, 28, 29, 32, 35, 37, 39, 40, 41, 42, 49, 54, 58, 59, 60, 61))
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
width = 12, height = 8, units = "in", res = 300)
FeaturePlot(cd8_sub, reduction = "umap", features = c("Cd4", "Cd8a"), order = T)
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub_fastMNN_umap_plot.png", sep = "/"),
width = 12, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap", pt.size = 1.5, label = T, order= T)
dev.off()
cd8_sub <- NormalizeData(cd8_sub)
cd8_sub <- FindVariableFeatures(cd8_sub, selection.method = "vst",
nfeatures = ceiling(dim(cd8_sub)[1] * 0.2),
verbose = T)
cd8_sub <- RunFastMNN(object.list = SplitObject(cd8_sub, split.by = "RNA_Group"))
cd8_sub <- RunUMAP(cd8_sub, reduction = "mnn", dims = 1:15)
cd8_sub <- FindNeighbors(cd8_sub, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
cd8_sub <- FindClusters(cd8_sub, resolution = res)
}
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(cd8_sub) <- md
png(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", paste0("Cd8sub_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = cd8_sub, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = cd8_sub, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub",
paste("Cd8sub_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
cd8_sub@misc[["markers"]] <- df_markers
# cell-cycle scoring
cd8_sub@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
cd8_sub <- CellCycleScoring(object = cd8_sub, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
cd8_sub@meta.data$CC.Difference <- cd8_sub@meta.data$S.Score - cd8_sub@meta.data$G2M.Score
png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
group.by = c("Phase"),
ncol = 2)
dev.off()
png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
group.by = c("Phase"), split.by = "RNA_Group",
ncol = 3)
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
group.by = c("Phase"), split.by = "Phase",
ncol = 3)
dev.off()
# Save object
saveRDS(object = cd8_sub,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/"))
##############################
### Filter Denditric cells ###
##############################
full_obj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(full_obj) <- "RNA_snn_res.1.2"
DimPlot(object = full_obj, reduction = "umap", label = T)
full_min_dc <- subset(x = full_obj, idents = c(10, 11, 16, 22))
png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "Full_Dendritic_subset_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap", group.by = "RNA_snn_res.1.2", label = T)
dev.off()
full_min_dc <- NormalizeData(full_min_dc)
full_min_dc <- FindVariableFeatures(full_min_dc, selection.method = "vst",
nfeatures = ceiling(dim(full_min_dc)[1] * 0.2), verbose = T)
full_min_dc <- RunFastMNN(object.list = SplitObject(full_min_dc, split.by = "RNA_Group"))
full_min_dc <- RunUMAP(full_min_dc, reduction = "mnn", dims = 1:15)
full_min_dc <- FindNeighbors(full_min_dc, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
full_min_dc <- FindClusters(full_min_dc, resolution = res)
}
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
Idents(full_min_dc) <- md
png(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", paste0("DendriticFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = full_min_dc, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = full_min_dc, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2",
paste("DendriticFilteredNEW_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
full_min_dc@misc[["markers"]] <- df_markers
# cell-cycle scoring
full_min_dc@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
full_min_dc <- CellCycleScoring(object = full_min_dc, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
full_min_dc@meta.data$CC.Difference <- full_min_dc@meta.data$S.Score - full_min_dc@meta.data$G2M.Score
png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
group.by = c("Phase"),
ncol = 2)
dev.off()
png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
group.by = c("Phase"), split.by = "RNA_Group",
ncol = 3)
dev.off()
png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
group.by = c("Phase"), split.by = "Phase",
ncol = 3)
dev.off()
# Save object
saveRDS(object = full_min_dc,
file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
######################
### Filter B cells ###
######################
fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(fmnn_full) <- "RNA_snn_res.0.6"
fmnn_bcells <- subset(fmnn_full, idents = c(12))
fmnn_bcells <- FindVariableFeatures(fmnn_bcells, selection.method = "vst", nfeatures = ceiling(dim(fmnn_bcells)[1] * 0.2), verbose = T)
fmnn_bcells <- RunFastMNN(object.list = SplitObject(fmnn_bcells, split.by = "RNA_Group"))
fmnn_bcells <- RunUMAP(fmnn_bcells, reduction = "mnn", dims = 1:20)
fmnn_bcells <- FindNeighbors(fmnn_bcells, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2)) {
fmnn_bcells <- FindClusters(fmnn_bcells, resolution = res)
}
png(filename = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
width = 14, height = 12, units = "in", res = 200)
DimPlot(fmnn_bcells, reduction = "umap",
group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2"),
ncol = 2, label = T)
dev.off()
df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2")) {
Idents(fmnn_bcells) <- md
png(paste(wdir, "BcellFiltered", "01-fastMNN_2", paste0("BcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
width = 10, height = 8, units = "in", res = 300)
print(DimPlot(object = fmnn_bcells, reduction = "umap", pt.size = 1.5, label = T))
dev.off()
df_markers[[md]] <- FindAllMarkers(object = fmnn_bcells, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
wb_markers <- createWorkbook()
for (clu in unique(df_markers[[md]]$cluster)) {
df_markers_clu <- subset(df_markers[[md]], cluster == clu)
addWorksheet(wb_markers, clu)
writeData(wb_markers, sheet = clu, x = df_markers_clu)
}
saveWorkbook(wb = wb_markers, overwrite = T,
file = paste(wdir, "BcellFiltered", "01-fastMNN_2",
paste("BcellFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
fmnn_bcells@misc[["markers_full"]] <- df_markers
png(filename = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
width = 14, height = 8, units = "in", res = 300)
DimPlot(fmnn_bcells, reduction = "umap",
group.by = c("RNA_Group", "orig.ident"),
ncol = 2, label = T)
dev.off()
# Save object
saveRDS(object = fmnn_bcells, file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/"))
scripts/2_fastMNN_DGE.R 0 → 100644
View file @ dba01155
library(Seurat)
library(SeuratWrappers)
library(ggplot2)
library(dplyr)
library(ggrepel)
library(reshape)
library(openxlsx)
library(stringr)
library(RColorBrewer)
library(SingleR)
library(ComplexHeatmap)
library(circlize)
plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr) {
data <- data.frame(gene = row.names(dge_res), pval = -log10(dge_res$p_val_adj), lfc = dge_res$avg_logFC)
data <- na.omit(data)
data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed",
data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed",
abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant"))
data <- data[order(data$pval, decreasing = TRUE),]
highl <- head(subset(data, color != "NonSignificant"), 20)
vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
theme_bw(base_size = 12) +
theme(legend.position = "right") +
xlab("log2 Fold Change") +
ylab(expression(-log[10]("adjusted p-value"))) +
geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
geom_point(size = 2, alpha = 0.8, na.rm = T) +
scale_color_manual(name = "Expression",
values = c(Overexpressed = "#CD4F39",
Underexpressed = "#008B00",
NonSignificant = "darkgray")) +
geom_label_repel(data = highl, aes(label = gene), show.legend = FALSE)
return(vol)
}
wdir <- "Birocchi_SciTraMed2022/results"
adj.pval.thr <- 1e-06
logfc.thr <- 0
# Full
fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(fmnn_full) <- "RNA_Group"
fmnn_full_mark_low <- FindMarkers(object = fmnn_full, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_full_mark_high <- FindMarkers(object = fmnn_full, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_full_mark_highlow <- FindMarkers(object = fmnn_full, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_full_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_full_mark_low,
"CD_High-AB_Ctrl" = fmnn_full_mark_high,
"CD_High-EF_Low" = fmnn_full_mark_highlow)
fmnn_full@misc[["dgemarkers"]] <- fmnn_full_dge_mark
write.xlsx(x = fmnn_full_dge_mark,
file = paste(wdir, "Full", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_full,
file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_full_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "Full", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Macrophages
fmnn_macro <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
Idents(fmnn_macro) <- "RNA_Group"
fmnn_macro_mark_low <- FindMarkers(object = fmnn_macro, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_macro_mark_high <- FindMarkers(object = fmnn_macro, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_macro_mark_highlow <- FindMarkers(object = fmnn_macro, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_macro_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_macro_mark_low,
"CD_High-AB_Ctrl" = fmnn_macro_mark_high,
"CD_High-EF_Low" = fmnn_macro_mark_highlow)
fmnn_macro@misc[["dgemarkers"]] <- fmnn_macro_dge_mark
write.xlsx(x = fmnn_macro_dge_mark,
file = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_macro,
file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_macro_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# T cell
fmnn_tcell <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
Idents(fmnn_tcell) <- "RNA_Group"
fmnn_tcell_mark_low <- FindMarkers(object = fmnn_tcell, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_tcell_mark_high <- FindMarkers(object = fmnn_tcell, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_tcell_mark_highlow <- FindMarkers(object = fmnn_tcell, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_tcell_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_tcell_mark_low,
"CD_High-AB_Ctrl" = fmnn_tcell_mark_high,
"CD_High-EF_Low" = fmnn_tcell_mark_highlow)
fmnn_tcell@misc[["dgemarkers"]] <- fmnn_tcell_dge_mark
write.xlsx(x = fmnn_tcell_dge_mark,
file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_tcell,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_tcell_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Cd4 - T cell
fmnn_cd4 <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/"))
Idents(fmnn_cd4) <- "RNA_Group"
fmnn_cd4_mark_low <- FindMarkers(object = fmnn_cd4, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_cd4_mark_high <- FindMarkers(object = fmnn_cd4, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_cd4_mark_highlow <- FindMarkers(object = fmnn_cd4, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_cd4_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_cd4_mark_low,
"CD_High-AB_Ctrl" = fmnn_cd4_mark_high,
"CD_High-EF_Low" = fmnn_cd4_mark_highlow)
fmnn_cd4@misc[["dgemarkers"]] <- fmnn_cd4_dge_mark
write.xlsx(x = fmnn_cd4_dge_mark,
file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd4sub", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_cd4,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_cd4_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd4sub",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Cd8 - T cell
fmnn_cd8 <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/"))
Idents(fmnn_cd8) <- "RNA_Group"
fmnn_cd8_mark_low <- FindMarkers(object = fmnn_cd8, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_cd8_mark_high <- FindMarkers(object = fmnn_cd8, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_cd8_mark_highlow <- FindMarkers(object = fmnn_cd8, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_cd8_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_cd8_mark_low,
"CD_High-AB_Ctrl" = fmnn_cd8_mark_high,
"CD_High-EF_Low" = fmnn_cd8_mark_highlow)
fmnn_cd8@misc[["dgemarkers"]] <- fmnn_cd8_dge_mark
write.xlsx(x = fmnn_cd8_dge_mark,
file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd8sub", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_cd8,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_cd8_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd8sub",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Dendritic NEW
fmnn_dc <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
Idents(fmnn_dc) <- "RNA_Group"
fmnn_dc_mark_low <- FindMarkers(object = fmnn_dc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_dc_mark_high <- FindMarkers(object = fmnn_dc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_dc_mark_highlow <- FindMarkers(object = fmnn_dc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_dc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_dc_mark_low,
"CD_High-AB_Ctrl" = fmnn_dc_mark_high,
"CD_High-EF_Low" = fmnn_dc_mark_highlow)
fmnn_dc@misc[["dgemarkers"]] <- fmnn_dc_dge_mark
write.xlsx(x = fmnn_dc_dge_mark,
file = paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_dc,
file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_dc_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# T cell
fmnn_tc <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
Idents(fmnn_tc) <- "RNA_Group"
fmnn_tc_mark_low <- FindMarkers(object = fmnn_tc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_tc_mark_high <- FindMarkers(object = fmnn_tc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_tc_mark_highlow <- FindMarkers(object = fmnn_tc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_tc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_tc_mark_low,
"CD_High-AB_Ctrl" = fmnn_tc_mark_high,
"CD_High-EF_Low" = fmnn_tc_mark_highlow)
fmnn_tc@misc[["dgemarkers"]] <- fmnn_tc_dge_mark
write.xlsx(x = fmnn_tc_dge_mark,
file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_tc,
file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_tc_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# T cell
fmnn_tcell <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
dim(fmnn_tcell)
cd4_ssub <- subset(x = fmnn_tcell, Cd4 > Cd8a & Cd4 > 0)
dim(cd4_ssub)
FeaturePlot(object = cd4_ssub, features = c("Cd4", "Cd8a"), reduction = "umap", split.by = "RNA_Group")
Idents(cd4_ssub) <- "RNA_Group"
cd4_ssub_mark_low <- FindMarkers(object = cd4_ssub, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
cd4_ssub_mark_high <- FindMarkers(object = cd4_ssub, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
cd4_ssub_mark_highlow <- FindMarkers(object = cd4_ssub, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
cd4_ssub_dge_mark <- list("EF_Low-AB_Ctrl" = cd4_ssub_mark_low,
"CD_High-AB_Ctrl" = cd4_ssub_mark_high,
"CD_High-EF_Low" = cd4_ssub_mark_highlow)
cd4_ssub@misc[["dgemarkers"]] <- cd4_ssub_dge_mark
for (contr in names(cd4_ssub@misc$dgemarkers)) {
deg_res <- cd4_ssub@misc$dgemarkers[[contr]]
deg_res.logfc <- as.vector(deg_res$avg_logFC)
names(deg_res.logfc) <- row.names(deg_res)
deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
print(deg_res.logfc[1:5])
# GSEA
contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = exah.gmt, nPerm = 10000)
if (nrow(contr.enricher.gsea) > 0) {
# Write res
contr.enricher.gsea.res <- list(contr.enricher.gsea)
names(contr.enricher.gsea.res) <- c(contr)
write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, "Cd4_Exhaustion_GSEA.xlsx", sep = "/"))
# Plot
for (t in seq(1, nrow(contr.enricher.gsea))) {
d <- contr.enricher.gsea$Description[t]
p <- gseaplot2(contr.enricher.gsea, geneSetID = t,
#title = gsub("_", " " , d),
color = c("forestgreen"),
pvalue_table = T)
fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, paste0("Cd4_", d, "_GSEA.pdf"), sep = "/")
ggsave(filename = fn, plot = p, width = 10, height = 7)
}
}
}
# for (contr in names(cd4_ssub@misc$dgemarkers)) {
# deg_res <- cd4_ssub@misc$dgemarkers[[contr]]
# print(nrow(deg_res))
# deg_res_top <- subset(deg_res, p_val_adj < 1e-06 & avg_logFC > 0)
# cat(cdsub, contr, nrow(deg_res_top), "\n")
# # Enrichment
# contr.enrichement <- enricher(rownames(deg_res_top), TERM2GENE = exah.gmt, pvalueCutoff = 0.05, universe = rownames(cd4_ssub))
# #print(contr.enrichement)
# if (nrow(contr.enrichement) > 0) {
# print(contr.enrichement@result[,1:7])
# # # Write res
# # contr.enricher.gsea.res <- list(contr.enricher.gsea)
# # names(contr.enricher.gsea.res) <- c(contr)
# # write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, "Exhaustion_GSEA.xlsx", sep = "/"))
# # # Plot
# # for (t in seq(1, nrow(contr.enricher.gsea))) {
# # d <- contr.enricher.gsea$Description[t]
# # p <- gseaplot2(contr.enricher.gsea, geneSetID = t,
# # #title = gsub("_", " " , d),
# # color = c("forestgreen"),
# # pvalue_table = T)
# # fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, paste0(d, "_GSEA.pdf"), sep = "/")
# # ggsave(filename = fn, plot = p, width = 10, height = 7)
# # }
# }
# }
saveRDS(object = cd4_ssub, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "TcellFiltered_fastMNN_Cd4subTest_final.rds", sep = "/"))
cd8_ssub <- subset(x = fmnn_tcell, Cd8a > Cd4 & Cd8a > 0)
dim(cd8_ssub)
FeaturePlot(object = cd8_ssub, features = c("Cd4", "Cd8a"), reduction = "umap", split.by = "RNA_Group")
Idents(cd8_ssub) <- "RNA_Group"
cd8_ssub_mark_low <- FindMarkers(object = cd8_ssub, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
cd8_ssub_mark_high <- FindMarkers(object = cd8_ssub, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
cd8_ssub_mark_highlow <- FindMarkers(object = cd8_ssub, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
cd8_ssub_dge_mark <- list("EF_Low-AB_Ctrl" = cd8_ssub_mark_low,
"CD_High-AB_Ctrl" = cd8_ssub_mark_high,
"CD_High-EF_Low" = cd8_ssub_mark_highlow)
cd8_ssub@misc[["dgemarkers"]] <- cd8_ssub_dge_mark
for (contr in names(cd8_ssub@misc$dgemarkers)) {
deg_res <- cd8_ssub@misc$dgemarkers[[contr]]
deg_res.logfc <- as.vector(deg_res$avg_logFC)
names(deg_res.logfc) <- row.names(deg_res)
deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
print(deg_res.logfc[1:5])
# GSEA
contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = exah.gmt, nPerm = 10000)
if (nrow(contr.enricher.gsea) > 0) {
# Write res
contr.enricher.gsea.res <- list(contr.enricher.gsea)
names(contr.enricher.gsea.res) <- c(contr)
write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, "Cd8_Exhaustion_GSEA.xlsx", sep = "/"))
# Plot
for (t in seq(1, nrow(contr.enricher.gsea))) {
d <- contr.enricher.gsea$Description[t]
p <- gseaplot2(contr.enricher.gsea, geneSetID = t,
#title = gsub("_", " " , d),
color = c("forestgreen"),
pvalue_table = T)
fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, paste0("Cd8_", d, "_GSEA.pdf"), sep = "/")
ggsave(filename = fn, plot = p, width = 10, height = 7)
}
}
}
# for (contr in names(cd8_ssub@misc$dgemarkers)) {
# deg_res <- cd8_ssub@misc$dgemarkers[[contr]]
# print(nrow(deg_res))
# deg_res_top <- subset(deg_res, p_val_adj < 1e-06 & avg_logFC > 0)
# cat(cdsub, contr, nrow(deg_res_top), "\n")
# # Enrichment
# contr.enrichement <- enricher(rownames(deg_res_top), TERM2GENE = exah.gmt, pvalueCutoff = 0.05, universe = rownames(cd8_ssub))
# #print(contr.enrichement)
# if (nrow(contr.enrichement) > 0) {
# print(contr.enrichement@result[,1:7])
# # # Write res
# # contr.enricher.gsea.res <- list(contr.enricher.gsea)
# # names(contr.enricher.gsea.res) <- c(contr)
# # write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, "Exhaustion_GSEA.xlsx", sep = "/"))
# # # Plot
# # for (t in seq(1, nrow(contr.enricher.gsea))) {
# # d <- contr.enricher.gsea$Description[t]
# # p <- gseaplot2(contr.enricher.gsea, geneSetID = t,
# # #title = gsub("_", " " , d),
# # color = c("forestgreen"),
# # pvalue_table = T)
# # fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, paste0(d, "_GSEA.pdf"), sep = "/")
# # ggsave(filename = fn, plot = p, width = 10, height = 7)
# # }
# }
# }
saveRDS(object = cd8_ssub, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "TcellFiltered_fastMNN_Cd8subTest_final.rds", sep = "/"))
# Monocytes
fmnn_monoc <- readRDS(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/"))
Idents(fmnn_monoc) <- "RNA_Group"
fmnn_monoc_mark_low <- FindMarkers(object = fmnn_monoc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_monoc_mark_high <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_monoc_mark_highlow <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_monoc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_monoc_mark_low,
"CD_High-AB_Ctrl" = fmnn_monoc_mark_high,
"CD_High-EF_Low" = fmnn_monoc_mark_highlow)
fmnn_monoc@misc[["dgemarkers"]] <- fmnn_monoc_dge_mark
write.xlsx(x = fmnn_monoc_dge_mark,
file = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_monoc,
file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_monoc_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Monocytes No3
fmnn_monoc <- readRDS(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/"))
Idents(fmnn_monoc) <- "RNA_Group"
fmnn_monoc_mark_low <- FindMarkers(object = fmnn_monoc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_monoc_mark_high <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_monoc_mark_highlow <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_monoc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_monoc_mark_low,
"CD_High-AB_Ctrl" = fmnn_monoc_mark_high,
"CD_High-EF_Low" = fmnn_monoc_mark_highlow)
fmnn_monoc@misc[["dgemarkers"]] <- fmnn_monoc_dge_mark
write.xlsx(x = fmnn_monoc_dge_mark,
file = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers_No3.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_monoc,
file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_monoc_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "_No3_volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# Full cluster markers
fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(fmnn_full) <- "RNA_snn_res.0.6"
clu_deg <- list()
for (clu in unique(fmnn_full$RNA_snn_res.0.6)) {
print(clu)
clu_mark <- list()
fmnn_full_clu <- subset(fmnn_full, ident = clu)
Idents(fmnn_full_clu) <- "RNA_Group"
clu_mark[["EF_Low-AB_Ctrl"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
clu_mark[["CD_High-AB_Ctrl"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
clu_mark[["CD_High-EF_Low"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
clu_deg[[clu]] <- clu_mark
}
fmnn_full@misc[["deg_clumarkers"]] <- clu_deg
saveRDS(object = fmnn_full, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
wb <- createWorkbook()
for (clu in names(clu_deg)) {
addWorksheet(wb, clu)
writeData(wb, clu, "EF_Low-AB_Ctrl", startCol = 1, startRow = 1)
writeData(wb, clu, clu_deg[[clu]][["EF_Low-AB_Ctrl"]], startCol = 1, startRow = 2, rowNames = TRUE)
writeData(wb, clu, "CD_High-AB_Ctrl", startCol = 10, startRow = 1)
writeData(wb, clu, clu_deg[[clu]][["CD_High-AB_Ctrl"]], startCol = 10, startRow = 2, rowNames = TRUE)
}
## Save workbook to working directory
saveWorkbook(wb, paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_deg_clumarkers.xlsx", sep = "/"), overwrite = TRUE)
# SingleR Annotation
fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
organism.db <- list(ImmGenData = ImmGenData(),
MouseRNAseqData = MouseRNAseqData())
for(el in names(organism.db)){
# Main
lname <- paste("SingleR", el, sep = "_")
fmnn_full@misc[[lname]] <- SingleR(test = fmnn_full@assays$RNA@data,
ref = organism.db[[el]],
labels = organism.db[[el]]$label.main)
lname_lab <- paste(lname, "labels", sep = "_")
fmnn_full@meta.data[[lname_lab]] <- fmnn_full@misc[[lname]]$labels
# png(filename = paste(plot_dir, paste(sample, lname, "umap", "plot.png", sep = "_"), sep = "/"),
# units = "in", width = 12, height = 9, res = 200)
# print(DimPlot(object = fmnn_full, reduction = "umap", group.by = lname_lab, label = T))
# dev.off()
# Refined
lname_ref <- paste("SingleRrefined", el, sep = "_")
fmnn_full@misc[[lname_ref]] <- SingleR(test = fmnn_full@assays$RNA@data,
ref = organism.db[[el]],
labels = organism.db[[el]]$label.fine)
lname_ref_lab <- paste(lname_ref, "labels", sep = "_")
fmnn_full@meta.data[[lname_ref_lab]] <- fmnn_full@misc[[lname_ref]]$labels
# png(filename = paste(plot_dir, paste(sample, lname_ref, "umap", "plot.png", sep = "_"), sep = "/"),
# units = "in", width = 12, height = 9, res = 200)
# print(DimPlot(object = fmnn_full, reduction = "umap", group.by = lname_ref_lab, label = T))
# dev.off()
}
saveRDS(object = fmnn_full, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
# B cells
fmnn_bcells <- readRDS(file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/"))
Idents(fmnn_bcells) <- "RNA_Group"
fmnn_bcells_mark_low <- FindMarkers(object = fmnn_bcells, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_bcells_mark_high <- FindMarkers(object = fmnn_bcells, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
fmnn_bcells_mark_highlow <- FindMarkers(object = fmnn_bcells, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
fmnn_bcells_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_bcells_mark_low,
"CD_High-AB_Ctrl" = fmnn_bcells_mark_high,
"CD_High-EF_Low" = fmnn_bcells_mark_highlow)
fmnn_bcells@misc[["dgemarkers"]] <- fmnn_bcells_dge_mark
write.xlsx(x = fmnn_bcells_dge_mark,
file = paste(wdir, "BcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers_No3.xlsx", sep = "/"), row.names = TRUE)
saveRDS(object = fmnn_bcells,
file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/"))
for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) {
res_dge <- fmnn_bcells_dge_mark[[contr]]
print(head(res_dge))
vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "BcellFiltered", "02-fastMNN_2-post-analysis",
paste(contr, "_volcano.png", sep = "-"), sep = "/"),
plot = vol, width = 10, height = 9, dpi = 600)
}
# # Macrophages cluster markers
# fmnn_macro <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
# Idents(fmnn_macro) <- "RNA_snn_res.0.6"
# clu_deg <- list()
# for (clu in unique(fmnn_macro$RNA_snn_res.0.6)) {
# print(clu)
# clu_mark <- list()
# fmnn_macro_clu <- subset(fmnn_macro, ident = clu)
# Idents(fmnn_macro_clu) <- "RNA_Group"
# clu_mark[["EF_Low-AB_Ctrl"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0)
# clu_mark[["CD_High-AB_Ctrl"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0)
# clu_mark[["CD_High-EF_Low"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0)
# clu_deg[[clu]] <- clu_mark
# }
# fmnn_macro@misc[["deg_clumarkers"]] <- clu_deg
# saveRDS(object = fmnn_macro,
# file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
#
# full_res <- data.frame()
# for (clu in names(clu_deg)) {
# for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) {
# clu_name <- paste0(clu, "-", contr)
# res <- subset(clu_deg[[clu]][[contr]], p_val_adj < 0.05 & abs(avg_logFC) >= 1)
# res_filt <- res[, "avg_logFC", drop = FALSE]
# colnames(res_filt) <- clu_name
# full_res <- merge(full_res, res_filt, by = 0, all = T)
# rownames(full_res) <- full_res$Row.names
# full_res$Row.names <- NULL
# }
# }
# nrow(full_res)
# #full_res[is.na(full_res)] <- 0
#
# wb <- createWorkbook()
# for (clu in names(clu_deg)) {
# addWorksheet(wb, clu)
# writeData(wb, clu, "EF_Low-AB_Ctrl", startCol = 1, startRow = 1)
# writeData(wb, clu, clu_deg[[clu]][["EF_Low-AB_Ctrl"]], startCol = 1, startRow = 2, rowNames = TRUE)
# writeData(wb, clu, "CD_High-AB_Ctrl", startCol = 10, startRow = 1)
# writeData(wb, clu, clu_deg[[clu]][["CD_High-AB_Ctrl"]], startCol = 10, startRow = 2, rowNames = TRUE)
# }
# ## Save workbook to working directory
# saveWorkbook(wb, paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_deg_clumarkers.xlsx", sep = "/"), overwrite = TRUE)
#
\ No newline at end of file
scripts/3_macroGTvsCtrlPost.R 0 → 100644
View file @ dba01155
library(Seurat)
library(dplyr)
library(openxlsx)
library(clusterProfiler)
library(ReactomePA)
library(org.Mm.eg.db)
library(dplyr)
library(stringr)
library(ComplexHeatmap)
library(RColorBrewer)
library(grid)
library(scales)
library(ggrepel)
library(cowplot)
library(gridExtra)
library(circlize)
plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr) {
data <- data.frame(gene = row.names(dge_res), pval = -log10(dge_res$p_val_adj), lfc = dge_res$avg_logFC)
data <- na.omit(data)
data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed",
data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed",
abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant"))
data <- data[order(data$pval, decreasing = TRUE),]
highl <- head(subset(data, color != "NonSignificant"), 20)
vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
theme_bw(base_size = 12) +
theme(legend.position = "right") +
xlab("log2 Fold Change") +
ylab(expression(-log[10]("adjusted p-value"))) +
geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
geom_point(size = 2, alpha = 0.8, na.rm = T) +
scale_color_manual(name = "Expression",
values = c(Overexpressed = "#CD4F39",
Underexpressed = "#008B00",
NonSignificant = "darkgray")) +
geom_label_repel(data = highl, aes(label = gene), show.legend = FALSE)
return(vol)
}
wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
md <- macro_obj@meta.data
md <- mutate(md, RNA_Treatment = ifelse(md$RNA_Group == "AB_Ctrl", "Control", "GeneTherapy"))
rownames(md) <- colnames(macro_obj)
head(md)
treat <- md$RNA_Treatment
names(treat) <- rownames(md)
macro_obj <- AddMetaData(object = macro_obj, metadata = treat, col.name = "RNA_Treatment")
Idents(macro_obj) <- "RNA_Treatment"
treat_dge <- FindMarkers(object = macro_obj, ident.1 = "GeneTherapy", ident.2 = "Control", logfc.threshold = 0, min.pct = 0)
head(treat_dge)
write.xlsx(x = list("GeneTherapy-Control" = treat_dge),
file = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control-DGE.xlsx", sep = "/"), row.names = T)
adj.pval.thr <- 1e-06
logfc.thr <- 0
vol <- plotVolcano(treat_dge, adj.pval.thr, logfc.thr)
ggsave(filename = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control-volcano.png", sep = "/"),
plot = vol, width = 10, height = 9, dpi = 200)
org_db <- org.Mm.eg.db
gmt_obj <- read.gmt("h.all.v7.1.symbols_mouse.gmt")
c7_gmt <- read.gmt("c7.all.v7.1.symbols_mouse.gmt")
pval_thr <- 0.05
out_dir <- paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", sep = "/")
contr <- "GeneTherapy-Control"
dge_res_con <- treat_dge
dge_res_con$symbol <- rownames(dge_res_con)
dge_res_con$entrez <- mapIds(x = org_db,
keys = dge_res_con$symbol,
column = "ENTREZID",
keytype = "SYMBOL",
multiVals = "first")
head(dge_res_con)
# Sort by pval
gene.res.logfc <- as.vector(dge_res_con$avg_logFC)
names(gene.res.logfc) <- dge_res_con$symbol
gene.res.logfc <- sort(gene.res.logfc, decreasing = TRUE)
gene.res.logfc <- gene.res.logfc[!duplicated(names(gene.res.logfc))]
print(gene.res.logfc[1:10])
gene.res.logfc.entrez <- as.vector(dge_res_con$avg_logFC)
names(gene.res.logfc.entrez) <- dge_res_con$entrez
gene.res.logfc.entrez <- sort(gene.res.logfc.entrez, decreasing = TRUE)
gene.res.logfc.entrez <- gene.res.logfc.entrez[!duplicated(names(gene.res.logfc.entrez))]
print(gene.res.logfc.entrez[1:10])
# Top genes
top.genes <- subset(dge_res_con, p_val_adj < 0.05)
top.genes.names <- as.character(top.genes$symbol)
print(top.genes.names[1:10])
top.genes.entrez <- as.character(top.genes$entrez)
print(top.genes.entrez[1:10])
# Top genes pos
top.genes.pos <- subset(dge_res_con, p_val_adj < 0.05 & avg_logFC > 0)
top.genes.pos.names <- as.character(top.genes.pos$symbol)
print(top.genes.pos.names[1:10])
top.genes.pos.entrez <- as.character(top.genes.pos$entrez)
print(top.genes.pos.entrez[1:10])
# Top genes neg
top.genes.neg <- subset(dge_res_con, p_val_adj < 0.05 & avg_logFC < 0)
top.genes.neg.names <- as.character(top.genes.neg$symbol)
print(top.genes.neg.names[1:10])
top.genes.neg.entrez <- as.character(top.genes.neg$entrez)
print(top.genes.neg.entrez[1:10])
##################
### Enrichment ###
##################
enr_res <- list()
## enrich GO
contr.enrich_go <- enrichGO(gene = top.genes.names,
universe = as.character(dge_res_con$symbol),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res[["enr_GO_BP"]] <- contr.enrich_go
## enrich KEGG
contr.enrich_kegg <- enrichKEGG(gene = top.genes.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg.symb <- setReadable(contr.enrich_kegg, org_db, keyType = "ENTREZID")
enr_res[["enr_KEGG"]] <- contr.enrich_kegg.symb
## enrich Pathway
contr.enrich_path <- enrichPathway(gene = top.genes.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path.symb <- setReadable(contr.enrich_path, org_db, keyType = "ENTREZID")
enr_res[["enr_Pathway"]] <- contr.enrich_path.symb
## Enrich
contr.hallmark <- enricher(top.genes.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res[["enr_HALLMARK"]] <- contr.hallmark
contr.c7 <- enricher(top.genes.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = c7_gmt)
enr_res[["enr_C7"]] <- contr.hallmark
write.xlsx(x = enr_res, file = paste(out_dir, paste0(contr, "_enrichment.xlsx"), sep = "/"))
######################
### Enrichment Pos ###
######################
enr_res_pos <- list()
## enrich GO
contr.enrich_go.pos <- enrichGO(gene = top.genes.pos.names,
universe = as.character(dge_res_con$symbol),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res_pos[["enr_pos_GO_BP"]] <- contr.enrich_go.pos
## enrich KEGG
contr.enrich_kegg_pos <- enrichKEGG(gene = top.genes.pos.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg_pos.symb <- setReadable(contr.enrich_kegg_pos, org_db, keyType = "ENTREZID")
enr_res_pos[["enr_pos_KEGG"]] <- contr.enrich_kegg_pos.symb
## enrich Pathway
contr.enrich_path_pos <- enrichPathway(gene = top.genes.pos.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path_pos.symb <- setReadable(contr.enrich_path_pos, org_db, keyType = "ENTREZID")
enr_res_pos[["enr_pos_Pathway"]] <- contr.enrich_path_pos.symb
## Enrich
contr.hallmark_pos <- enricher(top.genes.pos.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res_pos[["enr_pos_HALLMARK"]] <- contr.hallmark_pos
contr.c7_pos <- enricher(top.genes.pos.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = c7_gmt)
enr_res_pos[["enr_pos_C7"]] <- contr.c7_pos
write.xlsx(x = enr_res_pos, file = paste(out_dir, paste0(contr, "_enrichment_pos.xlsx"), sep = "/"))
######################
### Enrichment Neg ###
######################
enr_res_neg <- list()
## enrich GO
contr.enrich_go.neg <- enrichGO(gene = top.genes.neg.names,
universe = as.character(dge_res_con$symbol),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res_neg[["enr_neg_GO_BP"]] <- contr.enrich_go.neg
## enrich KEGG
contr.enrich_kegg_neg <- enrichKEGG(gene = top.genes.neg.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg_neg.symb <- setReadable(contr.enrich_kegg_neg, org_db, keyType = "ENTREZID")
enr_res_neg[["enr_neg_KEGG"]] <- contr.enrich_kegg_neg.symb
## enrich Pathway
contr.enrich_path_neg <- enrichPathway(gene = top.genes.neg.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path_neg.symb <- setReadable(contr.enrich_path_neg, org_db, keyType = "ENTREZID")
enr_res_neg[["enr_neg_Pathway"]] <- contr.enrich_path_neg.symb
## Enrich
contr.hallmark_neg <- enricher(top.genes.neg.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res_neg[["enr_neg_HALLMARK"]] <- contr.hallmark_neg
contr.c7_neg <- enricher(top.genes.neg.names,
universe = as.character(dge_res_con$symbol),
pvalueCutoff = pval_thr,
TERM2GENE = c7_gmt)
enr_res_neg[["enr_neg_C7"]] <- contr.c7_neg
write.xlsx(x = enr_res_neg, file = paste(out_dir, paste0(contr, "_enrichment_neg.xlsx"), sep = "/"))
############
### GSEA ###
############
gsea_res <- list()
## gsea GO
cat("# gsea GO", "\n")
contr.gsea_go <- gseGO(geneList = gene.res.logfc.entrez,
OrgDb = org_db,
ont = "BP",
nPerm = 10000,
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_go.symb <- setReadable(contr.gsea_go, org_db, keyType = "auto")
gsea_res[["gsea_GO_BP"]] <- contr.gsea_go.symb
## GSEA KEGG
cat("# gsea KEGG", "\n")
contr.gsea_kegg <- gseKEGG(geneList = gene.res.logfc.entrez,
organism = "mmu",
nPerm = 10000,
minGSSize = 120,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_kegg.symb <- setReadable(contr.gsea_kegg, org_db, keyType = "ENTREZID")
gsea_res[["gsea_KEGG"]] <- contr.gsea_kegg.symb
## GSEA Pathway
cat("# gsea Pathway", "\n")
contr.gsea_path <- gsePathway(geneList = gene.res.logfc.entrez,
organism = "mouse",
nPerm = 10000,
minGSSize = 120,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_path.symb <- setReadable(contr.gsea_path, org_db, keyType = "ENTREZID")
gsea_res[["gsea_Pathway"]] <- contr.gsea_path.symb
## GSEA HALLMARK
contr.enricher.gsea <- GSEA(gene.res.logfc, TERM2GENE = gmt_obj, nPerm = 10000, pvalueCutoff = pval_thr)
gsea_res[["gsea_HALLMARK"]] <- contr.enricher.gsea
contr.enricher.gsea.c7 <- GSEA(gene.res.logfc, TERM2GENE = c7_gmt, nPerm = 10000, pvalueCutoff = pval_thr)
gsea_res[["gsea_C7"]] <- contr.enricher.gsea.c7
write.xlsx(x = gsea_res, file = paste(out_dir, paste0(contr, "_GSEA.xlsx"), sep = "/"))
#######################################################################################
# No cluster 6
wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
clu6_mark <- subset(macro_obj@misc$markers_full$RNA_snn_res.0.6, cluster == "6")
nrow(clu6_mark)
head(clu6_mark)
org_db <- org.Mm.eg.db
gmt_obj <- read.gmt("h.all.v7.1.symbols_mouse.gmt")
pval_thr <- 0.05
out_dir <- paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", sep = "/")
contr <- "Clu6Markers"
dge_res_con <- clu6_mark
dge_res_con$entrez <- mapIds(x = org_db,
keys = dge_res_con$gene,
column = "ENTREZID",
keytype = "SYMBOL",
multiVals = "first")
head(dge_res_con)
# Sort by pval
gene.res.logfc <- as.vector(dge_res_con$avg_logFC)
names(gene.res.logfc) <- dge_res_con$gene
gene.res.logfc <- sort(gene.res.logfc, decreasing = TRUE)
gene.res.logfc <- gene.res.logfc[!duplicated(names(gene.res.logfc))]
print(gene.res.logfc[1:10])
gene.res.logfc.entrez <- as.vector(dge_res_con$avg_logFC)
names(gene.res.logfc.entrez) <- dge_res_con$entrez
gene.res.logfc.entrez <- sort(gene.res.logfc.entrez, decreasing = TRUE)
gene.res.logfc.entrez <- gene.res.logfc.entrez[!duplicated(names(gene.res.logfc.entrez))]
print(gene.res.logfc.entrez[1:10])
# Top genes
top.genes <- subset(dge_res_con, p_val_adj < 1e-06)
top.genes.names <- as.character(top.genes$gene)
print(top.genes.names[1:10])
top.genes.entrez <- as.character(top.genes$entrez)
print(top.genes.entrez[1:10])
# Top genes pos
top.genes.pos <- subset(dge_res_con, p_val_adj < 1e-06 & avg_logFC > 0)
top.genes.pos.names <- as.character(top.genes.pos$gene)
print(top.genes.pos.names[1:10])
top.genes.pos.entrez <- as.character(top.genes.pos$entrez)
print(top.genes.pos.entrez[1:10])
# Top genes neg
top.genes.neg <- subset(dge_res_con, p_val_adj < 1e-06 & avg_logFC < 0)
top.genes.neg.names <- as.character(top.genes.neg$gene)
print(top.genes.neg.names[1:10])
top.genes.neg.entrez <- as.character(top.genes.neg$entrez)
print(top.genes.neg.entrez[1:10])
##################
### Enrichment ###
##################
enr_res <- list()
## enrich GO
contr.enrich_go <- enrichGO(gene = top.genes.names,
universe = as.character(dge_res_con$gene),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res[["enr_GO_BP"]] <- contr.enrich_go
## enrich KEGG
contr.enrich_kegg <- enrichKEGG(gene = top.genes.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg.symb <- setReadable(contr.enrich_kegg, org_db, keyType = "ENTREZID")
enr_res[["enr_KEGG"]] <- contr.enrich_kegg.symb
## enrich Pathway
contr.enrich_path <- enrichPathway(gene = top.genes.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path.symb <- setReadable(contr.enrich_path, org_db, keyType = "ENTREZID")
enr_res[["enr_Pathway"]] <- contr.enrich_path.symb
## Enrich
contr.hallmark <- enricher(top.genes.names,
universe = as.character(dge_res_con$gene),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res[["enr_HALLMARK"]] <- contr.hallmark
write.xlsx(x = enr_res, file = paste(out_dir, paste0(contr, "_enrichment.xlsx"), sep = "/"))
######################
### Enrichment Pos ###
######################
enr_res_pos <- list()
## enrich GO
contr.enrich_go.pos <- enrichGO(gene = top.genes.pos.names,
universe = as.character(dge_res_con$gene),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res_pos[["enr_pos_GO_BP"]] <- contr.enrich_go.pos
## enrich KEGG
contr.enrich_kegg_pos <- enrichKEGG(gene = top.genes.pos.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg_pos.symb <- setReadable(contr.enrich_kegg_pos, org_db, keyType = "ENTREZID")
enr_res_pos[["enr_pos_KEGG"]] <- contr.enrich_kegg_pos.symb
## enrich Pathway
contr.enrich_path_pos <- enrichPathway(gene = top.genes.pos.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path_pos.symb <- setReadable(contr.enrich_path_pos, org_db, keyType = "ENTREZID")
enr_res_pos[["enr_pos_Pathway"]] <- contr.enrich_path_pos.symb
## Enrich
contr.hallmark_pos <- enricher(top.genes.pos.names,
universe = as.character(dge_res_con$gene),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res_pos[["enr_pos_HALLMARK"]] <- contr.hallmark_pos
write.xlsx(x = enr_res_pos, file = paste(out_dir, paste0(contr, "_enrichment_pos.xlsx"), sep = "/"))
######################
### Enrichment Neg ###
######################
enr_res_neg <- list()
## enrich GO
contr.enrich_go.neg <- enrichGO(gene = top.genes.neg.names,
universe = as.character(dge_res_con$gene),
OrgDb = org_db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = pval_thr,
keyType = 'SYMBOL')
enr_res_neg[["enr_neg_GO_BP"]] <- contr.enrich_go.neg
## enrich KEGG
contr.enrich_kegg_neg <- enrichKEGG(gene = top.genes.neg.entrez,
organism = "mmu",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_kegg_neg.symb <- setReadable(contr.enrich_kegg_neg, org_db, keyType = "ENTREZID")
enr_res_neg[["enr_neg_KEGG"]] <- contr.enrich_kegg_neg.symb
## enrich Pathway
contr.enrich_path_neg <- enrichPathway(gene = top.genes.neg.entrez,
organism = "mouse",
universe = as.character(dge_res_con$entrez),
pvalueCutoff = pval_thr)
contr.enrich_path_neg.symb <- setReadable(contr.enrich_path_neg, org_db, keyType = "ENTREZID")
enr_res_neg[["enr_neg_Pathway"]] <- contr.enrich_path_neg.symb
## Enrich
contr.hallmark_neg <- enricher(top.genes.neg.names,
universe = as.character(dge_res_con$gene),
pvalueCutoff = pval_thr,
TERM2GENE = gmt_obj)
enr_res_neg[["enr_neg_HALLMARK"]] <- contr.hallmark_neg
write.xlsx(x = enr_res_neg, file = paste(out_dir, paste0(contr, "_enrichment_neg.xlsx"), sep = "/"))
############
### GSEA ###
############
gsea_res <- list()
## gsea GO
cat("# gsea GO", "\n")
contr.gsea_go <- gseGO(geneList = gene.res.logfc.entrez,
OrgDb = org_db,
ont = "BP",
nPerm = 10000,
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_go.symb <- setReadable(contr.gsea_go, org_db, keyType = "auto")
gsea_res[["gsea_GO_BP"]] <- contr.gsea_go.symb
## GSEA KEGG
cat("# gsea KEGG", "\n")
contr.gsea_kegg <- gseKEGG(geneList = gene.res.logfc.entrez,
organism = "mmu",
nPerm = 10000,
minGSSize = 120,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_kegg.symb <- setReadable(contr.gsea_kegg, org_db, keyType = "ENTREZID")
gsea_res[["gsea_KEGG"]] <- contr.gsea_kegg.symb
## GSEA Pathway
cat("# gsea Pathway", "\n")
contr.gsea_path <- gsePathway(geneList = gene.res.logfc.entrez,
organism = "mouse",
nPerm = 10000,
minGSSize = 120,
pvalueCutoff = pval_thr,
verbose = FALSE)
contr.gsea_path.symb <- setReadable(contr.gsea_path, org_db, keyType = "ENTREZID")
gsea_res[["gsea_Pathway"]] <- contr.gsea_path.symb
## GSEA HALLMARK
contr.enricher.gsea <- GSEA(gene.res.logfc, TERM2GENE = gmt_obj, nPerm = 10000, pvalueCutoff = pval_thr)
gsea_res[["gsea_HALLMARK"]] <- contr.enricher.gsea
write.xlsx(x = gsea_res, file = paste(out_dir, paste0(contr, "_GSEA.xlsx"), sep = "/"))
scripts/4_monocle.R 0 → 100644
View file @ dba01155
library(monocle3)
library(Seurat)
library(SeuratData)
library(SeuratWrappers)
library(ggplot2)
library(dplyr)
wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
out_dir <- paste(wdir, "MacrophagesFiltered", "03-fastMNN_2-pseudotime", sep = "/")
# Load the data
expression_matrix <- macro_obj@assays$RNA@counts
cell_metadata <- macro_obj@meta.data
cell_metadata$orig.ident <- factor(cell_metadata$orig.ident, levels = unique(cell_metadata$orig.ident))
gene_annotation <- data.frame("gene_short_name" = rownames(macro_obj))
rownames(gene_annotation) <- gene_annotation$gene_short_name
# Make the CDS object
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 100)
cds <- align_cds(cds, num_dim = 50, alignment_group = "RNA_Group")
cds <- reduce_dimension(cds)
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_Group_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds, color_cells_by = "RNA_Group", label_cell_groups = FALSE, group_label_size = 8, cell_size = 0.5)
dev.off()
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds, color_cells_by = "RNA_snn_res.0.6", label_cell_groups = T, group_label_size = 8, cell_size = 0.5)
dev.off()
cds <- cluster_cells(cds = cds,resolution = 0.001)
plot_cells(cds, cell_size = 1, color_cells_by = "partition")
plot_cells(cds, genes=c("Pglyrp1", "Ace"))
marker_test_res <- top_markers(cds, group_cells_by="partition", reference_cells=1000, cores=8)
head(marker_test_res)
marker_test_res$cell_group
top_specific_markers <- marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>% filter(cell_group == "6") %>%
top_n(20, pseudo_R2)
marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>% filter(cell_group == "3") %>%
top_n(20, pseudo_R2) %>% select(gene_short_name) %>% list()
tt3 <- as.data.frame(top_specific_markers)
str(tt3)
tt3$cell_group
subset(tt3, tt3$cell_group %in% c("3", "4", "5", "6"))
top_specific_marker_ids <- unique(top_specific_markers %>% pull(gene_id))
cds <- learn_graph(cds)
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
color_cells_by = "RNA_snn_res.0.6",
label_groups_by_cluster=FALSE,
label_leaves=F,
label_branch_points=F,
graph_label_size=1.5,
group_label_size = 6, cell_size = 0.5)
dev.off()
cds <- order_cells(cds)
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_pseudotime_umap_trajectory_plot3.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
color_cells_by = "pseudotime",
label_groups_by_cluster=FALSE,
label_leaves=T,
label_branch_points=T,
graph_label_size=2.5,
group_label_size = 8, cell_size = 0.5)
dev.off()
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="7"){
cell_ids <- which(colData(cds)[, "RNA_snn_res.0.6"] == time_bin)
closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_pseudotime_umap_trajectory_plot2.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
color_cells_by = "pseudotime",
label_groups_by_cluster=FALSE,
label_leaves=T,
label_branch_points=T,
graph_label_size=2.5,
group_label_size = 8, cell_size = 0.5)
dev.off()
cds_pr_test_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
head(cds_pr_test_res, n = 10)
nrow(cds_pr_test_res)
pr_deg_ids <- subset(cds_pr_test_res, q_value < 0.05)
head(pr_deg_ids, n = 10)
nrow(pr_deg_ids)
plot_cells(cds, genes = rownames(head(pr_deg_ids[order(pr_deg_ids$q_value, decreasing = F),], n = 9)),
show_trajectory_graph = T,
label_cell_groups = FALSE,
label_leaves = FALSE)
# plot_cells(cds, genes = rownames(pr_deg_ids)[1:6],
# show_trajectory_graph = T,
# label_cell_groups = FALSE,
# label_leaves = FALSE)
gene_module_df <- find_gene_modules(cds[rownames(pr_deg_ids),], resolution=c(10^-10, 10^seq(-6,-1)))
cell_group_df <- tibble::tibble(cell=row.names(colData(cds)),
cell_group=colData(cds)$RNA_snn_res.0.6)
agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
pheatmap::pheatmap(agg_mat,angle_col = 0, scale="column", clustering_method="ward.D2")
# mod <- 5 # Cluster 6
# mod <- 33 # Clsuter 10
# pr_deg_ids_mod <- pr_deg_ids[gene_module_df[gene_module_df$module == mod,][["id"]],]
# plot_cells(cds, genes = rownames(head(pr_deg_ids_mod[order(pr_deg_ids_mod$q_value, decreasing = F),], n = 9)),
# show_trajectory_graph = T,
# label_cell_groups = FALSE,
# label_leaves = FALSE)
plot_cells(cds,
genes = gene_module_df %>% filter(module %in% c(7, 25, 54, 14, 32, 12, 19, 38, 24, 58, 15, 5, 16, 43, 18, 29)),
label_cell_groups = FALSE,
show_trajectory_graph = FALSE)
# Analyzing branches in single-cell trajectories
cds_subset <- choose_cells(cds)
subset_pr_test_res <- graph_test(cds_subset, neighbor_graph="principal_graph", cores=4)
pr_deg_ids <- row.names(subset(subset_pr_test_res, q_value < 0.05))
gene_module_df <- find_gene_modules(cds_subset[pr_deg_ids,], resolution=0.001)
agg_mat <- aggregate_gene_expression(cds_subset, gene_module_df)
module_dendro <- hclust(dist(agg_mat))
gene_module_df$module <- factor(gene_module_df$module, levels = row.names(agg_mat)[module_dendro$order])
plot_cells(cds_subset,
genes=gene_module_df,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
clu6_genes <- c("Mrc1", "Ccl8")
clu6_lineage_cds <- cds[rowData(cds)$gene_short_name %in% clu6_genes,
colData(cds)$RNA_snn_res.0.6 %in% c("6", "7", "10")]
plot_genes_in_pseudotime(clu6_lineage_cds, color_cells_by="RNA_snn_res.0.6", min_expr=0.5)
#####################################################
#####################################################
Idents(macro_obj) <- "RNA_Group"
# Ctrl
macro_obj_ctrl <- subset(macro_obj, idents = "AB_Ctrl")
expression_matrix_ctrl <- macro_obj_ctrl@assays$RNA@counts
cell_metadata_ctrl <- macro_obj_ctrl@meta.data
cell_metadata_ctrl$orig.ident <- factor(cell_metadata_ctrl$orig.ident, levels = unique(cell_metadata_ctrl$orig.ident))
gene_annotation_ctrl <- data.frame("gene_short_name" = rownames(macro_obj_ctrl))
rownames(gene_annotation_ctrl) <- gene_annotation_ctrl$gene_short_name
cds_ctrl <- new_cell_data_set(expression_data = expression_matrix_ctrl,
cell_metadata = cell_metadata_ctrl,
gene_metadata = gene_annotation_ctrl)
cds_ctrl <- preprocess_cds(cds_ctrl, num_dim = 50)
cds_ctrl <- reduce_dimension(cds_ctrl)
### Construct and assign the made up partition
recreate.partition_ctrl <- c(rep(1, length(cds_ctrl@colData@rownames)))
names(recreate.partition_ctrl) <- cds_ctrl@colData@rownames
recreate.partition_ctrl <- as.factor(recreate.partition_ctrl)
cds_ctrl@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_ctrl
### Assign the cluster info
list_cluster_ctrl <- macro_obj_ctrl@meta.data$RNA_snn_res.0.6
names(list_cluster_ctrl) <- macro_obj_ctrl@assays[["RNA"]]@data@Dimnames[[2]]
cds_ctrl@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_ctrl
### Could be a space-holder, but essentially fills out louvain parameters
cds_ctrl@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_ctrl)[["UMAP"]] <- macro_obj_ctrl@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
##########cds_ctrl@preprocess_aux$gene_loadings <- macro_obj_ctrl@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_ctrl@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#########colnames(reducedDims(cds_ctrl)[["UMAP"]]) <- NULL
plot_cells(cds = cds_ctrl, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_ctrl <- learn_graph(cds_ctrl, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl,
color_cells_by = "RNA_snn_res.0.6",
label_groups_by_cluster=FALSE,
label_leaves=F,
label_branch_points=F,
graph_label_size=1.5,
group_label_size = 6, cell_size = 1)
dev.off()
cds_ctrl <- order_cells(cds_ctrl)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl,
color_cells_by = "pseudotime",
label_groups_by_cluster=FALSE,
label_leaves=T,
label_branch_points=T,
graph_label_size=2.5,
group_label_size = 8, cell_size = 1)
dev.off()
####
cds_ctrl_ps_res = graph_test(cds_ctrl, neighbor_graph = "principal_graph", cores = 4)
cds_ctrl_ps_res_sig = row.names(subset(cds_ctrl_ps_res, q_value == 0))
length(cds_ctrl_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl, genes = cds_ctrl_ps_res_sig[0:25],
show_trajectory_graph = T,
label_cell_groups = FALSE,
label_leaves = FALSE)
dev.off()
gene_module_df_ctrl <- find_gene_modules(cds_ctrl[cds_ctrl_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_ctrl <- tibble::tibble(cell=row.names(colData(cds_ctrl)),
cell_group=colData(cds_ctrl)$RNA_snn_res.0.6)
agg_mat_ctrl <- aggregate_gene_expression(cds_ctrl, gene_module_df_ctrl, cell_group_df_ctrl)
row.names(agg_mat_ctrl) <- stringr::str_c("Module ", row.names(agg_mat_ctrl))
pheatmap::pheatmap(agg_mat_ctrl, scale = "column", clustering_method = "ward.D2",
filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_Modules.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_ctrl, cell_size = 1,
genes = gene_module_df_ctrl,
label_cell_groups = FALSE,
show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_ctrl)
nrow(gene_module_df_ctrl)
for (mod in levels(gene_module_df_ctrl$module)) {
tt <- as.data.frame(subset(gene_module_df_ctrl, module == mod))
print(nrow(tt))
lineage_cds_ctrl <- cds_ctrl[rowData(cds_ctrl)$gene_short_name %in% tt$id,]
png(filename = paste(out_dir, paste0("MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
width = 10, height = 12, units = "in", res = 300)
print(plot_genes_in_pseudotime(cds_subset = lineage_cds_ctrl, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
dev.off()
}
# High
macro_obj_high <- subset(macro_obj, idents = "CD_High")
expression_matrix_high <- macro_obj_high@assays$RNA@counts
cell_metadata_high <- macro_obj_high@meta.data
cell_metadata_high$orig.ident <- factor(cell_metadata_high$orig.ident, levels = unique(cell_metadata_high$orig.ident))
gene_annotation_high <- data.frame("gene_short_name" = rownames(macro_obj_high))
rownames(gene_annotation_high) <- gene_annotation_high$gene_short_name
cds_high <- new_cell_data_set(expression_data = expression_matrix_high,
cell_metadata = cell_metadata_high,
gene_metadata = gene_annotation_high)
cds_high <- preprocess_cds(cds_high, num_dim = 50)
cds_high <- reduce_dimension(cds_high)
### Construct and assign the made up partition
recreate.partition_high <- c(rep(1, length(cds_high@colData@rownames)))
names(recreate.partition_high) <- cds_high@colData@rownames
recreate.partition_high <- as.factor(recreate.partition_high)
cds_high@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_high
### Assign the cluster info
list_cluster_high <- macro_obj_high@meta.data$RNA_snn_res.0.6
names(list_cluster_high) <- macro_obj_high@assays[["RNA"]]@data@Dimnames[[2]]
cds_high@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_high
### Could be a space-holder, but essentially fills out louvain parameters
cds_high@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_high)[["UMAP"]] <- macro_obj_high@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
#########cds_high@preprocess_aux$gene_loadings <- macro_obj_high@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_high@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#########colnames(cds_high@reducedDims$UMAP) <- NULL
plot_cells(cds = cds_high, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_high <- learn_graph(cds_high, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high,
color_cells_by = "RNA_snn_res.0.6",
label_groups_by_cluster=FALSE,
label_leaves=F,
label_branch_points=F,
graph_label_size=1.5,
group_label_size = 6, cell_size = 1)
dev.off()
cds_high <- order_cells(cds_high)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high,
color_cells_by = "pseudotime",
label_groups_by_cluster=FALSE,
label_leaves=T,
label_branch_points=T,
graph_label_size=2.5,
group_label_size = 8, cell_size = 1)
dev.off()
####
cds_high_ps_res = graph_test(cds_high, neighbor_graph = "principal_graph", cores = 4)
cds_high_ps_res_sig = row.names(subset(cds_high_ps_res, q_value == 0))
length(cds_high_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high, genes = cds_high_ps_res_sig[0:25],
show_trajectory_graph = T,
label_cell_groups = FALSE,
label_leaves = FALSE)
dev.off()
gene_module_df_high <- find_gene_modules(cds_high[cds_high_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_high <- tibble::tibble(cell=row.names(colData(cds_high)),
cell_group=colData(cds_high)$RNA_snn_res.0.6)
agg_mat_high <- aggregate_gene_expression(cds_high, gene_module_df_high, cell_group_df_high)
row.names(agg_mat_high) <- stringr::str_c("Module ", row.names(agg_mat_high))
pheatmap::pheatmap(agg_mat_high, scale = "column", clustering_method = "ward.D2",
filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_Modules.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_high, cell_size = 1,
genes = gene_module_df_high,
label_cell_groups = FALSE,
show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_high)
nrow(gene_module_df_high)
for (mod in levels(gene_module_df_high$module)) {
tt <- as.data.frame(subset(gene_module_df_high, module == mod))
print(nrow(tt))
lineage_cds_high <- cds_high[rowData(cds_high)$gene_short_name %in% tt$id,]
png(filename = paste(out_dir, paste0("MacrophagesFiltered_High_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
width = 10, height = 12, units = "in", res = 300)
print(plot_genes_in_pseudotime(cds_subset = lineage_cds_high, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
dev.off()
}
# Low
macro_obj_low <- subset(macro_obj, idents = "EF_Low")
expression_matrix_low <- macro_obj_low@assays$RNA@counts
cell_metadata_low <- macro_obj_low@meta.data
cell_metadata_low$orig.ident <- factor(cell_metadata_low$orig.ident, levels = unique(cell_metadata_low$orig.ident))
gene_annotation_low <- data.frame("gene_short_name" = rownames(macro_obj_low))
rownames(gene_annotation_low) <- gene_annotation_low$gene_short_name
cds_low <- new_cell_data_set(expression_data = expression_matrix_low,
cell_metadata = cell_metadata_low,
gene_metadata = gene_annotation_low)
cds_low <- preprocess_cds(cds_low, num_dim = 50)
cds_low <- reduce_dimension(cds_low)
### Construct and assign the made up partition
recreate.partition_low <- c(rep(1, length(cds_low@colData@rownames)))
names(recreate.partition_low) <- cds_low@colData@rownames
recreate.partition_low <- as.factor(recreate.partition_low)
cds_low@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_low
### Assign the cluster info
list_cluster_low <- macro_obj_low@meta.data$RNA_snn_res.0.6
names(list_cluster_low) <- macro_obj_low@assays[["RNA"]]@data@Dimnames[[2]]
cds_low@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_low
### Could be a space-holder, but essentially fills out louvain parameters
cds_low@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_low)[["UMAP"]] <- macro_obj_low@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
############cds_low@preprocess_aux$gene_loadings <- macro_obj_low@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_low@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#colnames(cds_low@reducedDims$UMAP) <- NULL
#plot_cells(cds = cds_low, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_low <- learn_graph(cds_low, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low,
color_cells_by = "RNA_snn_res.0.6",
label_groups_by_cluster=FALSE,
label_leaves=F,
label_branch_points=F,
graph_label_size=1.5,
group_label_size = 6, cell_size = 1)
dev.off()
cds_low <- order_cells(cds_low)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low,
color_cells_by = "pseudotime",
label_groups_by_cluster=FALSE,
label_leaves=T,
label_branch_points=T,
graph_label_size=2.5,
group_label_size = 8, cell_size = 1)
dev.off()
####
cds_low_ps_res = graph_test(cds_low, neighbor_graph = "principal_graph", cores = 4)
cds_low_ps_res_sig = row.names(subset(cds_low_ps_res, q_value == 0))
length(cds_low_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low, genes = cds_low_ps_res_sig[0:25],
show_trajectory_graph = T,
label_cell_groups = FALSE,
label_leaves = FALSE)
dev.off()
gene_module_df_low <- find_gene_modules(cds_low[cds_low_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_low <- tibble::tibble(cell=row.names(colData(cds_low)),
cell_group=colData(cds_low)$RNA_snn_res.0.6)
agg_mat_low <- aggregate_gene_expression(cds_low, gene_module_df_low, cell_group_df_low)
row.names(agg_mat_low) <- stringr::str_c("Module ", row.names(agg_mat_low))
pheatmap::pheatmap(agg_mat_low, scale = "column", clustering_method = "ward.D2",
filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_Modules.png", sep = "/"),
width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_low, cell_size = 1,
genes = gene_module_df_low,
label_cell_groups = FALSE,
show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_low)
nrow(gene_module_df_low)
for (mod in levels(gene_module_df_low$module)) {
tt <- as.data.frame(subset(gene_module_df_low, module == mod))
print(nrow(tt))
lineage_cds_low <- cds_low[rowData(cds_low)$gene_short_name %in% tt$id,]
png(filename = paste(out_dir, paste0("MacrophagesFiltered_Low_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
width = 10, height = 12, units = "in", res = 300)
print(plot_genes_in_pseudotime(cds_subset = lineage_cds_low, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
dev.off()
}
cds_obj <- list("Full" = cds, "AB_Ctrl" = cds_ctrl, "CD_High" = cds_high, "EF_Low" = cds_low)
saveRDS(object = cds_obj, file = paste(out_dir, "Monocle_cds_objects.rds", sep = "/"))
scripts/5_prepareSCVelo.R 0 → 100644
View file @ dba01155
library(SeuratWrappers)
library(Seurat)
library(SeuratDisk)
ts <- c("FB_10-A4_T", "FB_12-F1_T", "FB_16-B1_T", "FB_18-F2_T", "FB_4-B3_T", "FB_6-D1_T", "FB_11-C0_T", "FB_17-E3_T", "FB_5-C1_T")
veldata <- list()
velobj <- list()
for (i in ts) {
veldata[[i]] <- ReadVelocity(file = paste0("Birocchi_SciTraMed2022/velocyto_loom/", i, ".loom"))
velobj[[i]] <- as.Seurat(x = veldata[[i]])
velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "merged", replacement = "", x = colnames(velobj[[i]])))
velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "x$", perl = T, replacement = "", x = colnames(velobj[[i]])))
velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = ":", perl = T, replacement = "_", x = colnames(velobj[[i]])))
velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "_", perl = T, replacement = "-", x = colnames(velobj[[i]])))
velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "-T-", perl = T, replacement = "_", x = colnames(velobj[[i]])))
velobj[[i]]@meta.data$orig.ident <- i
}
velobj_ctrl <- merge(merge(velobj[["FB_10-A4_T"]], velobj[["FB_16-B1_T"]]), velobj[["FB_4-B3_T"]])
velobj_high <- merge(merge(velobj[["FB_11-C0_T"]], velobj[["FB_5-C1_T"]]), velobj[["FB_6-D1_T"]])
velobj_low <- merge(merge(velobj[["FB_12-F1_T"]], velobj[["FB_18-F2_T"]]), velobj[["FB_17-E3_T"]])
velo_full <- merge(merge(velobj_ctrl, velobj_high), velobj_low)
wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
colnames(macro_obj)[1:5]
# write.csv(Cells(macro_obj), file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_cellID_obs.csv", sep = "/"), row.names = FALSE)
# write.csv(Embeddings(macro_obj, reduction = "umap"), file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_cell_embeddings.csv", sep = "/"))
# colnames(macro_obj@meta.data)
# write.csv(macro_obj@meta.data[,c("orig.ident", "RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1", "RNA_snn_res.1.2")],
# file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_mdata.csv", sep = "/"))
#
# #write.csv(Cells(obj), file = "cellID_obs.csv")
# macro_obj_embedding <- Embeddings(macro_obj, reduction = "umap")
# head(macro_obj_embedding)
#
# #write.csv(Embeddings(obj, reduction = "umap"), file = "cell_embeddings.csv")
# #my_clusters <- obj@meta.data$RNA_snn_res.0.6
# macro_obj_clu <- FetchData(object = macro_obj, vars = "RNA_snn_res.0.6")
# colnames(macro_obj_clu) <- c("Seurat_RNA_snn_res.0.6")
# head(macro_obj_clu)
# #write.csv(obj@meta.data$seurat_clusters, file = "clusters.csv")
# macro_obj_groups <- FetchData(object = macro_obj, "RNA_Group")
# colnames(macro_obj_groups) <- c("RNA_Group")
# head(macro_obj_groups)
#
# cells <- intersect(Cells(velo_full), Cells(macro_obj))
# length(Cells(velo_full))
# length(Cells(macro_obj))
# length(cells)
# velo_full <- subset(velo_full, cells = cells)
#
# velo_full[["RNA"]] <- velo_full[["spliced"]]
# #########velo_full[["RNA"]] <- macro_obj@assays$RNA
# # adding clusters and groups
# velo_full <- AddMetaData(object = velo_full, metadata = macro_obj_clu)
# velo_full <- AddMetaData(object = velo_full, metadata = macro_obj_groups)
#
# colnames(velo_full@meta.data)
#
# names(velo_full@assays)
# str(velo_full)
# # creating fake umap slots in order to create the reduction slot associated
# velo_full <- SCTransform(velo_full)
# velo_full <- RunPCA(velo_full)
# velo_full <- RunUMAP(velo_full, dims = 1:30, n.components = 2)
# velo_full <- FindNeighbors(velo_full, dims = 1:30)
# velo_full <- FindClusters(velo_full)
# DefaultAssay(velo_full) <- "RNA"
# # restore original UMAPs
# head(velo_full@reductions$umap@cell.embeddings)
# head(macro_obj_embedding[rownames(velo_full@reductions$umap@cell.embeddings),])
# velo_full@reductions$umap@cell.embeddings <- macro_obj_embedding[rownames(velo_full@reductions$umap@cell.embeddings),]
# SaveH5Seurat(velo_full, filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2","MacrophagesFiltered_velo.h5Seurat", sep = "/"), overwrite = T)
# Convert(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2","MacrophagesFiltered_velo.h5Seurat", sep = "/"), dest = "h5ad", overwrite = T)
# Groups
Idents(macro_obj) <- "RNA_Group"
velobj_groups <- list("AB_Ctrl" = velobj_ctrl, "CD_High" = velobj_high, "EF_Low" = velobj_low)
for (group in names(velobj_groups)) {
print(group)
velo_sub <- velobj_groups[[group]]
macro_obj_sub <- subset(macro_obj, idents = group)
macro_obj_sub_embedding <- Embeddings(macro_obj_sub, reduction = "umap")
head(macro_obj_sub_embedding)
macro_obj_sub_clu <- FetchData(object = macro_obj_sub, vars = "RNA_snn_res.0.6")
head(macro_obj_sub_clu)
cells_sub <- intersect(Cells(velo_sub), Cells(macro_obj_sub))
length(Cells(velo_sub))
length(Cells(macro_obj_sub))
length(cells_sub)
velo_sub <- subset(velo_sub, cells = cells_sub)
velo_sub[["RNA"]] <- velo_sub[["spliced"]]
# adding clusters
velo_sub <- AddMetaData(object = velo_sub, metadata = macro_obj_sub_clu)
# creating fake umap slots in order to create the reduction slot associated
velo_sub <- SCTransform(velo_sub)
velo_sub <- RunPCA(velo_sub)
velo_sub <- RunUMAP(velo_sub, dims = 1:30)
velo_sub <- FindNeighbors(velo_sub, dims = 1:30)
velo_sub <- FindClusters(velo_sub)
DefaultAssay(velo_sub) <- "RNA"
# restore original UMAPs
head(velo_sub@reductions$umap@cell.embeddings)
head(macro_obj_sub_embedding[rownames(velo_sub@reductions$umap@cell.embeddings),])
velo_sub@reductions$umap@cell.embeddings <- macro_obj_sub_embedding[rownames(velo_sub@reductions$umap@cell.embeddings),]
SaveH5Seurat(velo_sub, filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered-", group, "_velo.h5Seurat"), sep = "/"), overwrite = T)
Convert(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered-", group, "_velo.h5Seurat"), sep = "/"), dest = "h5ad", overwrite = T)
}
scripts/6_prepare_plots.R 0 → 100644
View file @ dba01155
library(Seurat)
library(ggplot2)
library(ComplexHeatmap)
library(dplyr)
library(circlize)
library(RColorBrewer)
library(gridExtra)
library(cowplot)
library(stringr)
library(scales)
library(clusterProfiler)
library(ggrepel)
library(openxlsx)
# Functions
plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr, hgenes = NULL) {
data <- data.frame(gene = row.names(dge_res), pval = -log10(dge_res$p_val_adj), lfc = dge_res$avg_logFC)
data <- na.omit(data)
data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed",
data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed",
abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant"))
data <- data[order(data$pval, decreasing = TRUE),]
if (is.null(hgenes)) {
highl <- head(subset(data, color != "NonSignificant"), 20)
} else {
highl <- subset(data, gene %in% hgenes)
}
vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
theme_bw(base_size = 20) +
theme(legend.position = "none",
panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
xlab("log2 Fold Change") +
ylab(expression(-log[10]("adjusted p-value"))) +
geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
geom_point(size = 2, alpha = 0.8, na.rm = T) +
scale_color_manual(name = "Expression",
values = c(Overexpressed = "#CD4F39",
Underexpressed = "#008B00",
NonSignificant = "darkgray")) +
geom_label_repel(data = subset(highl, lfc > 0), aes(label = gene), size = 8,
show.legend = FALSE, direction = "y", nudge_x = 1, nudge_y = 40, hjust = 0) +
geom_label_repel(data = subset(highl, lfc <= 0), aes(label = gene), size = 8,
show.legend = FALSE, direction = "y", nudge_x = -.5)
return(vol)
}
doSCHeatmapMean <- function(fname, clustering, ntop, num_color, high_genes, clu_lab, duplicated_labs = F) {
col_scales <- list(colorRampPalette(rev(brewer.pal(n = 9, name ="RdBu")))(100),
colorRamp2(c(-4, -3, -2, -1, 0, 1, 2, 3, 4), rev(brewer.pal(n = 9, name ="RdBu"))),
colorRamp2(c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2), rev(brewer.pal(n = 9, name ="RdBu"))),
colorRamp2(c(-1.5, -1, -0.5, 0, 0.5, 1, 1.5), rev(brewer.pal(n = 7, name ="RdBu"))),
colorRamp2(c(-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1), rev(brewer.pal(n = 9, name ="RdBu"))))
if (num_color > length(col_scales)) {
col_sc <- col_scales[[length(col_scales)]]
} else {
col_sc <- col_scales[[num_color]]
}
obj <- readRDS(fname)
top_genes <- obj@misc$markers[[clustering]] %>% group_by(cluster) %>% top_n(ntop, avg_logFC)
if (nrow(obj@assays$RNA@scale.data) == 0) {
hm_data <- as.matrix(obj@assays$RNA@data[top_genes$gene,])
#print(hm_data[1:3, 1:5])
hm_data <- t(scale(t(hm_data)))
#print(hm_data[1:3, 1:5])
} else {
hm_data <- obj@assays$RNA@scale.data[top_genes$gene,]
}
mm <- merge(obj@meta.data, t(hm_data), by = 0)
rownames(mm) <- mm$Row.names
mm$Row.names <- NULL
mm_filt <- mm[, c(clustering, top_genes$gene)]
mm_matrix <- mm_filt %>% group_by_at(clustering) %>% summarise_all(list(mean)) %>% data.frame()
#print(mm_matrix)
c19 <- hue_pal()(nrow(mm_matrix))
if (length(clu_lab) > 0) {
rownames(mm_matrix) <- clu_lab
names(c19) <- clu_lab
} else {
rownames(mm_matrix) <- mm_matrix[[clustering]]
names(c19) <- mm_matrix[[clustering]]
}
mm_matrix[[clustering]] <- NULL
hm_matrix <- as.matrix(mm_matrix)
if (duplicated_labs) {
colnames(hm_matrix) <- str_split_fixed(colnames(hm_matrix), pattern = "[.]", n = 2)[,1] # !!!
}
gg <- intersect(colnames(hm_matrix), high_genes)
gg_index <- which(colnames(hm_matrix) %in% gg)
labels <- colnames(hm_matrix)[gg_index]
ha_tre <- columnAnnotation(link = anno_mark(at = gg_index, labels = labels,
labels_rot = 45, padding = 0.75,
labels_gp = gpar(fontsize = 14),
link_width = unit(15, "mm"),
side = "bottom"))
if (length(clu_lab) > 0) {
rsplit <- factor(clu_lab, levels = clu_lab)
} else {
rsplit <- as.numeric(levels(obj@meta.data[[clustering]]))
}
ha_col <- rowAnnotation(df = data.frame(Cluster = rsplit), show_annotation_name = FALSE,
show_legend = FALSE, col = list(Cluster = c19),
annotation_legend_param = list(Cluster = list(grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"),
title_gp = gpar(col = "black", fontsize = 14),
labels_gp = gpar(col = "black", fontsize = 12))))
hm = Heatmap(matrix = hm_matrix,
name = "Expr.",
col = col_sc,
cluster_columns = FALSE,
cluster_rows = FALSE,
rect_gp = gpar(col = "grey", lwd = 0),
border = FALSE,
show_row_names = FALSE,
show_column_names = FALSE,
column_title_gp = gpar(fontsize = 0),
row_split = rsplit,
#cluster_row_slices = TRUE,
#column_split = obj@meta.data[[clustering]],
column_split = rep(0:(nrow(hm_matrix) - 1), each = ntop),
column_gap = unit(4, "mm"),
row_title_gp = gpar(fontsize = 18),
row_title_rot = 0,
left_annotation = ha_col,
bottom_annotation = ha_tre)
return(hm)
}
ggtable <- function(d, p = NULL) {
as.ggplot(tableGrob2(d, p))
}
tableGrob2 <- function(d, p = NULL, labels = NULL) {
d <- d[order(rownames(d)),]
if (is.null(labels)) {
tp <- tableGrob(d, theme=ttheme_default(base_size = 24))
} else {
tp <- tableGrob(d, rows = labels, parse = TRUE, theme=ttheme_default(base_size = 24))
}
if (is.null(p)) {
return(tp)
}
pcol <- unique(ggplot_build(p)$data[[1]][["colour"]])
j <- which(tp$layout$name == "rowhead-fg")
for (i in seq_along(pcol)) {
tp$grobs[j][[i+1]][["gp"]] = gpar(col = pcol[i], fontsize = 24)
}
return(tp)
}
gseaScores <- getFromNamespace("gseaScores", "DOSE")
gsInfo <- function(object, geneSetID) {
geneList <- object@geneList
if (is.numeric(geneSetID))
geneSetID <- object@result[geneSetID, "ID"]
geneSet <- object@geneSets[[geneSetID]]
exponent <- object@params[["exponent"]]
df <- gseaScores(geneList, geneSet, exponent, fortify=TRUE)
df$ymin=0
df$ymax=0
pos <- df$position == 1
h <- diff(range(df$runningScore))/20
df$ymin[pos] <- -h
df$ymax[pos] <- h
df$geneList <- geneList
df$Description <- object@result[geneSetID, "Description"]
return(df)
}
gseaplot2 <- function(x, geneSetID, title = "", color="green", base_size = 14,
rel_heights=c(1.5, .5, 1), subplots = 1:3, pvalue_table = FALSE, ES_geom="line", nopadj = FALSE, labels = NULL) {
ES_geom <- match.arg(ES_geom, c("line", "dot"))
geneList <- position <- NULL ## to satisfy codetool
if (length(geneSetID) == 1) {
gsdata <- gsInfo(x, geneSetID)
} else {
gsdata <- do.call(rbind, lapply(geneSetID, gsInfo, object = x))
}
p <- ggplot(gsdata, aes_(x = ~x)) + xlab(NULL) +
#theme_classic(base_size) +
theme_classic(base_size) +
theme(panel.grid.major = element_line(colour = "grey92"),
panel.grid.minor = element_line(colour = "grey92"),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()) +
scale_x_continuous(expand=c(0,0))
if (ES_geom == "line") {
es_layer <- geom_line(aes_(y = ~runningScore, color= ~Description), size=1)
} else {
es_layer <- geom_point(aes_(y = ~runningScore, color= ~Description), size=1, data = subset(gsdata, position == 1))
}
p.res <- p + es_layer +
theme(legend.position = c(.8, .8), legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"))
p.res <- p.res + ylab("Running Enrichment Score") +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x=element_blank(),
plot.margin=margin(t=.2, r = .2, b=0, l=.2, unit="cm"))
i <- 0
for (term in unique(gsdata$Description)) {
idx <- which(gsdata$ymin != 0 & gsdata$Description == term)
gsdata[idx, "ymin"] <- i
gsdata[idx, "ymax"] <- i + 1
i <- i + 1
}
p2 <- ggplot(gsdata, aes_(x = ~x)) +
geom_linerange(aes_(ymin=~ymin, ymax=~ymax, color=~Description)) +
xlab(NULL) + ylab(NULL) + theme_classic(base_size) +
theme(legend.position = "none",
plot.margin = margin(t=-.1, b=0,unit="cm"),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.line.x = element_blank()) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0))
if (length(geneSetID) == 1) {
## geneList <- gsdata$geneList
## j <- which.min(abs(geneList))
## v1 <- quantile(geneList[1:j], seq(0,1, length.out=6))[1:5]
## v2 <- quantile(geneList[j:length(geneList)], seq(0,1, length.out=6))[1:5]
## v <- sort(c(v1, v2))
## inv <- findInterval(geneList, v)
v <- seq(1, sum(gsdata$position), length.out=9)
inv <- findInterval(rev(cumsum(gsdata$position)), v)
if (min(inv) == 0) inv <- inv + 1
col=c(rev(brewer.pal(5, "Blues")), brewer.pal(5, "Reds"))
ymin <- min(p2$data$ymin)
yy <- max(p2$data$ymax - p2$data$ymin) * .3
xmin <- which(!duplicated(inv))
xmax <- xmin + as.numeric(table(inv)[unique(inv)])
d <- data.frame(ymin = ymin, ymax = yy,
xmin = xmin,
xmax = xmax,
col = col[unique(inv)])
p2 <- p2 + geom_rect(
aes_(xmin=~xmin,
xmax=~xmax,
ymin=~ymin,
ymax=~ymax,
fill=~I(col)),
data=d,
alpha=.9,
inherit.aes=FALSE)
}
## p2 <- p2 +
## geom_rect(aes(xmin=x-.5, xmax=x+.5, fill=geneList),
## ymin=ymin, ymax = ymin + yy, alpha=.5) +
## theme(legend.position="none") +
## scale_fill_gradientn(colors=color_palette(c("blue", "red")))
df2 <- p$data #data.frame(x = which(p$data$position == 1))
df2$y <- p$data$geneList[df2$x]
p.pos <- p + geom_segment(data=df2, aes_(x=~x, xend=~x, y=~y, yend=0), color="grey")
p.pos <- p.pos + ylab("Ranked list metric") +
xlab("Rank in Ordered Dataset") +
theme(plot.margin=margin(t = -.1, r = .2, b=.2, l=.2, unit="cm"))
#if (!is.null(title) && !is.na(title) && title != "") {
p.res <- p.res + ggtitle(title) + theme(plot.title = element_text(size = 40, face = "bold", hjust = 0.5))
#}
if (length(color) == length(geneSetID)) {
p.res <- p.res + scale_color_manual(values=color)
if (length(color) == 1) {
p.res <- p.res + theme(legend.position = "none")
p2 <- p2 + scale_color_manual(values = "black")
} else {
p2 <- p2 + scale_color_manual(values = color)
}
}
if (pvalue_table) {
if (nopadj) {
pd <- x[geneSetID, c("Description", "pvalue", "NES")]
} else {
pd <- x[geneSetID, c("Description", "pvalue", "p.adjust", "NES")]
}
pd <- pd[order(pd[,1], decreasing=FALSE),]
#pd$Description <- gsub("HALLMARK_", "", pd$Description)
#rownames(pd) <- pd$Description
pd <- pd[,-1]
pd <- round(pd, 4)
#tp <- tableGrob2(pd, p.res)
if (is.null(labels)) {
tp <- tableGrob(pd, rows = NULL, theme=ttheme_default(base_size = 26))
p.res <- p.res + theme(legend.position = "none") +
annotation_custom(tp,
xmin = quantile(p.res$data$x, .65),
xmax = quantile(p.res$data$x, .95),
ymin = quantile(p.res$data$runningScore, .5),
ymax = quantile(p.res$data$runningScore, .9))
} else {
tp <- tableGrob2(pd, p.res, labels)
p.res <- p.res + theme(legend.position = "none") +
annotation_custom(tp,
xmin = quantile(p.res$data$x, .6),
xmax = quantile(p.res$data$x, .9),
ymin = quantile(p.res$data$runningScore, .8),
ymax = quantile(p.res$data$runningScore, .95))
}
}
plotlist <- list(p.res, p2, p.pos)[subplots]
n <- length(plotlist)
plotlist[[n]] <- plotlist[[n]] +
theme(axis.line.x = element_line(),
axis.ticks.x=element_line(),
axis.text.x = element_text())
if (length(subplots) == 1) {
return(plotlist[[1]] + theme(plot.margin=margin(t=.2, r = .2, b=.2, l=.2, unit="cm")))
}
if (length(rel_heights) > length(subplots)) {
rel_heights <- rel_heights[subplots]
}
plot_grid(plotlist = plotlist, ncol=1, align="v", rel_heights=rel_heights)
}
#################
wdir <- "Birocchi_SciTraMed2022/results"
out_dir <- "Birocchi_SciTraMed2022/Paper"
##################
#### Figure 1 ####
##################
# 1 - UMAP cluster res 0.6
fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
full_umap_plot <- cbind(fobj@reductions$umap@cell.embeddings,
fobj@meta.data[rownames(fobj@reductions$umap@cell.embeddings),])
label <- data.frame(
UMAP_1 = c(-9.2, -7.5, -8, -6, -5.4, -4.8, -4, -2.5, -2.5, -0.9, -2, 1, 1.5, 4, 4, 7, 4, 5.2, 9),
UMAP_2 = c(1, -1, 6, 1, -1.5, -3.3, 6.5, -1.7, -5.5, -3.3, -8, -6.5, -10, -3.5, -1, 2.5, 3, 6.5, 3),
RNA_snn_res.0.6 = c(9, 3, 10, 4, 14, 5, 12, 18, 6, 17, 13, 8, 7, 11, 2, 1, 0, 15, 16),
label = c("Naive~T~cells", "CD4", "NK", "CD8", "T~regs", "Prolif.~T~cells", "B~cells", "Mast~cells",
"Unknown", "pDC", "DC1", "DC2", "mregDC", "Monocytes", "Macro~3", "Macro~2", "Macro~1", "Resid.like~Macro", "Prolif.~Macro")
)
p <- ggplot(full_umap_plot, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.6), size = 1, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none") +
#geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), parse = T, size = 5) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("")
ggsave(filename = paste(out_dir, "Fig1_1-Full_fastMNN_Res06.png", sep = "/"),
plot = p, width = 10, height = 10, dpi = 100)
# 2 - Dotplot of enrichment pos HALLMARK
ds_high <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/"))
ds_high_hall_enrich <- ds_high$results$enrichment$enrichment_h.all.v6.1.symbols_mouse_pos@result[,c("ID", "GeneRatio", "p.adjust", "Count")]
ds_high_hall_enrich$Contr <- "IFN-alpha"
ds_low <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/"))
ds_low_hall_enrich <- ds_low$results$enrichment$enrichment_h.all.v6.1.symbols_mouse_pos@result[,c("ID", "GeneRatio", "p.adjust", "Count")]
ds_low_hall_enrich$Contr <- "IFN-alpha-DHFR"
ds_hall_enrich <- rbind(ds_high_hall_enrich, ds_low_hall_enrich)
gsize <- as.numeric(sub("/\\d+$", "", as.character(ds_hall_enrich$GeneRatio)))
gcsize <- as.numeric(sub("^\\d+/", "", as.character(ds_hall_enrich$GeneRatio)))
ds_hall_enrich$GeneRatio = gsize/gcsize
ds_hall_enrich <- subset(ds_hall_enrich, p.adjust < 0.05)
ds_hall_enrich$ID <- gsub("HALLMARK_", "", ds_hall_enrich$ID)
#ds_hall_enrich$ID <- tolower(gsub("_", " ", ds_hall_enrich$ID))
print(head(ds_hall_enrich))
id_order <- ds_hall_enrich %>% select(ID, GeneRatio) %>% group_by(ID) %>% summarize_all(sum) %>% data.frame()
id_order <- id_order[order(id_order$GeneRatio, decreasing = F),]
ds_hall_enrich$ID <- factor(ds_hall_enrich$ID, levels = id_order$ID)
p <- ggplot(ds_hall_enrich, aes(x = GeneRatio, y = ID, size = Count, color = p.adjust)) +
theme_bw(base_size = 20) +
scale_color_gradient(low = "red", high = "blue", guide = guide_colorbar(reverse = T)) +
geom_point() +
ylab("") +
facet_grid(.~Contr, labeller = label_parsed)
ggsave(plot = p,
filename = paste(out_dir, "Fig1_2-Full_Hallmark_enrichment_pos.png", sep = "/"),
width = 12, height = 8, dpi = 100)
# 3/4 - GSEA hallmark ifna response High vs CTRL o Low vs CTRL
ds_high <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/"))
ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[1,]
p <- gseaplot2(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = 1, base_size = 18,
title = expression(IFN-alpha),
color = c("red"),
pvalue_table = T)
ggsave(filename = paste(out_dir, "Fig1_3-CD_High-AB_Ctrl_HALLMARK_INT_ALPHA_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
ds_low <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/"))
p <- gseaplot2(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = 1, base_size = 18,
title = expression(IFN-alpha-DHFR),
color = c("red"),
pvalue_table = T)
ggsave(filename = paste(out_dir, "Fig1_4-EF_Low-AB_Ctrl_HALLMARK_INT_ALPHA_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
# 5 - Heatmap DEGs HIGH or LOW vs CTRL for each cluster
fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
clu_deg <- fobj@misc[["deg_clumarkers"]]
full_res_genes <- c()
for (clu in names(clu_deg)) {
for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) {
clu_name <- paste0(clu, "-", contr)
res <- subset(clu_deg[[clu]][[contr]], p_val_adj < 0.05 & abs(avg_logFC) >= 1)
full_res_genes <- union(full_res_genes, rownames(res))
}
}
full_res <- data.frame()
for (clu in as.character(seq(0, 18, 1))) {
for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) {
clu_name <- paste0(clu, "-", contr)
res <- clu_deg[[clu]][[contr]]
res_filt <- res[rownames(res) %in% full_res_genes, "avg_logFC", drop = FALSE]
colnames(res_filt) <- clu_name
full_res <- merge(full_res, res_filt, by = 0, all = T)
rownames(full_res) <- full_res$Row.names
full_res$Row.names <- NULL
}
}
nrow(full_res)
full_res[is.na(full_res)] <- 0
sinfo <- data.frame(Dataset = colnames(full_res))
sinfo <- cbind(sinfo, str_split_fixed(sinfo$Dataset, "-", 3))
colnames(sinfo) <- c("Dataset", "Cluster", "Dose", "Ctrl")
rownames(sinfo) <- sinfo$Dataset
sinfo$Cluster <- factor(sinfo$Cluster, levels = as.character(seq(0, 18, 1)))
sinfo$Ctrl <- NULL
sinfo$Dataset <- NULL
sinfo$Treatment <- ifelse(sinfo$Dose == "EF_Low", "IFN-a-DHFR", "IFN-a")
sinfo$Dose <- NULL
ctypes <- c("Macro1","Macro2", "Macro3","CD4","CD8","ProlifT","Unk",
"mregDC","DC2","NaiveT","NK","Mono","Bcells","DC1","Tregs",
"ResMacro","ProlifMacro","pDC","MastCells")
levels(sinfo$Cluster) <- ctypes
inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1",
"Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr",
"Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7",
"Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b",
"Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9",
"Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1",
"Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars")
ifn_font <- ifelse(rownames(full_res) %in% inf_alpha, 2, 1)
ifn_col <- ifelse(rownames(full_res) %in% inf_alpha, "red", "black")
c19 <- hue_pal()(19)
names(c19) <- levels(sinfo$Cluster)
ha_col <- columnAnnotation(df = sinfo,
show_annotation_name = FALSE,
col = list(Treatment = c("IFN-a-DHFR" = "#6699FF", "IFN-a" = "#CC0066"), Cluster = c19),
annotation_name_side = "right",
annotation_legend_param = list(Treatment = list(at = c("IFN-a-DHFR", "IFN-a"),
grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"),
title_gp = gpar(col = "black", fontsize = 14),
labels = c(expression(paste("IFN-", alpha, "-DHFR ")),
expression(paste("IFN-", alpha, " "))),
# labels = gt_render(c("IFN-&alpha-DHFR ", "IFN-alpha- ")),
labels_gp = gpar(col = "black", fontsize = 12)),
Cluster = list(grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"),
title_gp = gpar(col = "black", fontsize = 14),
labels_gp = gpar(col = "black", fontsize = 12))))
hm = Heatmap(matrix = as.matrix(full_res),
col = colorRampPalette(rev(brewer.pal(n = 9, name ="RdYlBu")))(100),
cell_fun = function(j, i, x, y, width, height, fill) {
if(abs(as.matrix(full_res)[i, j]) > 0.5)
grid.text(sprintf("%.1f", as.matrix(full_res)[i, j]), x, y, gp = gpar(fontsize = 9))
},
cluster_columns = FALSE,
cluster_rows = TRUE,
rect_gp = gpar(col = "grey", lwd = 0),
border = FALSE,
name = "LogFC",
show_row_names = TRUE,
row_names_gp = gpar(fontsize = 11, col = ifn_col),
show_column_names = F,
column_title_gp = gpar(col = "black", fontsize = 11),
heatmap_legend_param = list(legend_height = unit(3, "cm"),
title_gp = gpar(col = "black", fontsize = 14),
labels_gp = gpar(col = "black", fontsize = 12)),
column_gap = unit(1, "mm"),
column_split = sinfo$Cluster,
top_annotation = ha_col,
#split = 2,
row_dend_reorder = TRUE,
width = 18, height = 10)
png(paste(out_dir, "Fig1_5-Full_fastMNN_deg_clumarkers_abslogFC1_fdr05_05Lab_noColNames.png", sep = "/"), width = 1800, height = 1000, res = 100)
draw(hm, merge_legends = T)
dev.off()
#########################
#### Figure 1 Suppl. ####
#########################
# 1 - Heatmap of top marker mean expression
hm <- doSCHeatmapMean(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"),
"RNA_snn_res.0.6", 10, 3, c("Saa3", "Chil3", "Ly6i", "C1qc", "Selenop", "Ms4a7", "Il1b", "Dusp1", "Socs3", "Furin", "Penk", "Trbc2",
"Cd8b1", "Cd8a", "Klrd1", "Ube2c", "Dut", "Stmn1", "Fscn1", "Il2ra", "Pkib", "Ifi27l2a", "Tmem123",
"Socs2", "Ccl17", "Il1r2", "Plet1", "S1pr1", "Klf2", "Dusp10", "Prf1", "Car2", "Klrb1c", "Cxcl10",
"Ifitm6", "Clec4e", "Igkc", "Igha", "Cd79a", "Cd24a", "Xcr1", "Naaa", "Izumo1r", "Cd2", "Ctla4", "Hexb",
"Cd81", "Psap", "Cks1b", "Top2a", "Siglech", "Ccr9", "Rnase6", "Mcpt2", "Cma1", "Cpa3"),
c("Macro1","Macro2", "Macro3","CD4","CD8","Prolif. T","Unknown", "mregDC","DC2","Naive T","NK","Monocytes",
"B cells","DC1","T regs","Res. Macro","Prolif. Macro","pDC","Mast Cells"))
png(paste(out_dir, "Fig1Suppl_1-Full_fastMNN_HM_top10_c3_clumean_topnames_2.png", sep = "/"),
width = 1500, height = 800, res = 100)
draw(hm)
dev.off()
# 2 - HALLMARK_IFNA_RESPONSE modulescore expression
fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1",
"Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr",
"Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7",
"Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b",
"Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9",
"Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1",
"Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars")
inf_alpha_filt <- inf_alpha[inf_alpha %in% rownames(fobj)]
fobj <- AddModuleScore(object = fobj, features = list(IFN_ALPHA = inf_alpha_filt), ctrl = length(inf_alpha_filt), name = "IFN_ALPHA")
umap_plot <- cbind(fobj@reductions$umap@cell.embeddings,
fobj@meta.data[rownames(fobj@reductions$umap@cell.embeddings),])
umap_plot$RNA_Group2 <- factor(umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
str(umap_plot)
summary(umap_plot$IFN_ALPHA1)
p <- ggplot(umap_plot, aes(x = UMAP_1, y = UMAP_2)) +
geom_point(data = subset(umap_plot, IFN_ALPHA1 > 0.2 & IFN_ALPHA1 < 0.5), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) +
geom_point(data = subset(umap_plot, IFN_ALPHA1 <= 0.2 & IFN_ALPHA1 > 0.1), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) +
geom_point(data = subset(umap_plot, IFN_ALPHA1 >= 0.5 & IFN_ALPHA1 < 0.75), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) +
geom_point(data = subset(umap_plot, IFN_ALPHA1 <= 0.1), aes(colour = IFN_ALPHA1), size = 1, alpha = 1) +
geom_point(data = subset(umap_plot, IFN_ALPHA1 >= 0.75), aes(colour = IFN_ALPHA1), size = 1, alpha = 1) +
scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"),
name = "Expr.",
#values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)),
values = rescale(c(-.1,.05,.2,.4,.7)),
guide = "colorbar") +
#theme_bw(base_size = 20) +
#theme_nothing(font_size = 20) +
theme_void(base_size = 20) +
#ggtitle("INTERFERON ALPHA RESPONSE") +
theme(legend.position = "right",
strip.text.x = element_text(size = 40)) +
#xlab("UMAP 1") +
#ylab("UMAP 2") +
facet_grid(.~RNA_Group2, labeller = label_parsed)
ggsave(filename = paste(out_dir, "Fig1Suppl_2-Full_fastMNN_IFNa_split.png", sep = "/"), plot = p, width = 20, height = 8, dpi = 100)
##################
#### Figure 2 ####
##################
# 1 - UMAP plot T cells res 0.6 with labels
tobj <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
t_umap_plot <- cbind(tobj@reductions$umap@cell.embeddings,
tobj@meta.data[rownames(tobj@reductions$umap@cell.embeddings),])
label <- data.frame(
UMAP_1 = c(0.5, 0, -5, 0, 3.5, 4, 2.2, -2.5),
UMAP_2 = c(1.5, -3.5, -0.3, 4.3, -3, 2.5, 0.5, 0.5),
RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7),
label = c("CD4 1", "CD8 1", "Proliferating", "CD4 2", "CD8 2", "Naive", "T regs", "CD4 3")
)
p <- ggplot(t_umap_plot, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none") +
#geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), size = 8) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("")
ggsave(filename = paste(out_dir, "Fig2_1-TcellFiltered_fastMNN_Res06.png", sep = "/"),
plot = p, width = 10, height = 8, dpi = 100)
# 2 - GSEA DEG High vs. CTRL e Low vs. CTRL Tcells, di signature Exhaustion_GSEA_new (Wherry et al., 2007)
tobj <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))
t_gmt <- read.gmt(paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "PaperSig.gmt", sep = "/"))
for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) {
deg_res <- tobj@misc$dgemarkers[[contr]]
deg_res.logfc <- as.vector(deg_res$avg_logFC)
names(deg_res.logfc) <- row.names(deg_res)
deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
print(deg_res.logfc[1:5])
# GSEA
contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = t_gmt, nPerm = 10000, pvalueCutoff = 0.2)
num_id <- 0
for (p in seq(1, nrow(contr.enricher.gsea@result))) {
if (as.character(contr.enricher.gsea@result[p, "ID"]) == "Exhaustion") {
num_id <- p
}
}
if (contr == "CD_High-AB_Ctrl") {
p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha),
color = c("red"),
pvalue_table = T, nopadj = T)
} else {
p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha-DHFR),
color = c("red"),
pvalue_table = T, nopadj = T)
}
ggsave(filename = paste(out_dir, paste0("Fig2_2-TcellFiltered_", contr, "_Exhaustion_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
}
# 3 - UMAP plot DC with labels res 0.2
dcobj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
dc_umap_plot <- cbind(dcobj@reductions$umap@cell.embeddings,
dcobj@meta.data[rownames(dcobj@reductions$umap@cell.embeddings),])
dc_umap_plot$RNA_snn_res.0.1 <- factor(dc_umap_plot$RNA_snn_res.0.1, levels = c(0, 1, 2, 3))
#levels(dc_umap_plot$RNA_snn_res.0.1) <- c("0", "1", "2", "3")
dc_umap_plot$RNA_snn_res.0.1
label <- data.frame(
UMAP_1 = c(-5, 5, -1, 1),
UMAP_2 = c(2.5, 0, -7.5, 10),
RNA_snn_res.0.1 = c(0, 1, 2, 3),
label = c("mregDC", "DC2", "DC1", "pDC")
)
p <- ggplot(dc_umap_plot, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.1), size = 1.5, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none") +
#geom_text(aes(label = RNA_snn_res.0.1), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), size = 8) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("")
ggsave(filename = paste(out_dir, "Fig2_3-DendriticFilteredNEW_fastMNN_Res01.png", sep = "/"),
plot = p, width = 10, height = 8, dpi = 100)
# 4/5 - GSEA DEG High vs. CTRL e Low vs. CTRL DC di signatures “Unstimulated DC vs. 4h LPS DN” e “Hallmark Oxidative phosphorylation”
dcobj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
c7_gmt <- read.gmt(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "c7.all.v6.1.symbols_mouse.gmt", sep = "/"))
for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) {
deg_res <- dcobj@misc$dgemarkers[[contr]]
deg_res.logfc <- as.vector(deg_res$avg_logFC)
names(deg_res.logfc) <- row.names(deg_res)
deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
print(deg_res.logfc[1:5])
# GSEA
contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = c7_gmt, nPerm = 10000, pvalueCutoff = 0.2)
num_id <- 0
for (p in seq(1, nrow(contr.enricher.gsea@result))) {
if (as.character(contr.enricher.gsea@result[p, "ID"]) == "GSE14000_UNSTIM_VS_4H_LPS_DC_DN") {
num_id <- p
}
}
if (contr == "CD_High-AB_Ctrl") {
p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha),
color = c("red"),
pvalue_table = T, nopadj = T)
} else {
p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha-DHFR),
color = c("red"),
pvalue_table = T, nopadj = T)
}
ggsave(filename = paste(out_dir, paste0("Fig2_4-DendriticFilteredNEW_", contr, "_GSE14000_UNSTIM_VS_4H_LPS_DC_DN_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
}
# HALLMARK_OXIDATIVE_PHOSPHORYLATION
ds_high <- readRDS(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/"))
num_id <- 0
for (p in seq(1, nrow(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result))) {
if (as.character(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[p, "ID"]) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION") {
num_id <- p
}
}
ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[num_id,3:8]
p <- gseaplot2(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha),
color = c("red"),
pvalue_table = T, nopadj = T)
ggsave(filename = paste(out_dir, "Fig2_5-DendriticFilteredNEW_CD_High-AB_Ctrl_HALLMARK_OXIDATIVE_PHOSPHORYLATION_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
# HALLMARK_OXIDATIVE_PHOSPHORYLATION
ds_low <- readRDS(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/"))
num_id <- 0
for (p in seq(1, nrow(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result))) {
if (as.character(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[p, "ID"]) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION") {
num_id <- p
}
}
ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[num_id,3:8]
p <- gseaplot2(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = num_id, base_size = 18,
title = expression(IFN-alpha-DHFR),
color = c("red"),
pvalue_table = T, nopadj = T)
ggsave(filename = paste(out_dir, "Fig2_5-DendriticFilteredNEW_EF_Low-AB_Ctrl_HALLMARK_OXIDATIVE_PHOSPHORYLATION_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
#########################
#### Figure 2 Suppl. ####
#########################
# 1 - Heatmap marker genes Tcells with labels
hm <- doSCHeatmapMean(fname = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"),
clustering = "RNA_snn_res.0.6", ntop = 15, num_color = 3,
high_genes = c("Furin","Cd4","Cxcr6","Bhlhe40","Gzma","Cd8a","Irf8","Prf1","Stmn1","Pclaf","Hist1h2ap",
"Tubb5","Ccl1","Penk","Xcl1","Ifng","Cd7","Gzmk","Irf7","Eomes","S1pr1","Dusp10","Ramp3",
"Tcf7","Ikzf2","Izumo1r","Foxp3","Ctla4","Nme1","Lgals1","Pgk1","Lgals7"),
clu_lab = c("CD4 1", "CD8 1", "Proliferating", "CD4 2", "CD8 2", "Naive", "T regs", "CD4 3"))
png(paste(out_dir, "Fig2Suppl_1-TcellFiltered_fastMNN_RNA_snn_res.0.6_top15-color3-wider2.png", sep = "/"),
width = 1500, height = 700, res = 150)
draw(hm)
dev.off()
# 2 - Heatmap marker genes DC with labels
# obj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
# Idents(obj) <- "RNA_snn_res.0.1"
# marks <- FindAllMarkers(object = obj, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
# min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
# obj@misc$markers[["RNA_snn_res.0.1"]] <- marks
# saveRDS(object = obj, file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))
hm <- doSCHeatmapMean(fname = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"),
clustering = "RNA_snn_res.0.1", ntop = 15, num_color = 3,
high_genes = c("Fscn1", "Tmem123", "Ccr7", "Ms4a6c", "S100a6", "Fcer1g", "Naaa", "Xcr1", "Cd24a", "Ly6d", "Ccr9", "Siglech"),
clu_lab = c("mregDC", "DC2", "DC1", "pDC"))
png(paste(out_dir, "Fig2Suppl_2-DendriticFilteredNEW_fastMNN_RNA_snn_res.0.1_top15-color3-wider2.png", sep = "/"),
width = 1500, height = 700, res = 150)
draw(hm)
dev.off()
##################
#### Figure 3 ####
##################
# 1 - UMAP plot Mono/macrophages split by treatment, res 0.6
mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
m_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings,
mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),])
label <- data.frame(
UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5),
UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0),
RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")
)
p <- ggplot(m_umap_plot, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none") +
#geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), size = 8) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("")
ggsave(filename = paste(out_dir, "Fig3_1-MacrophagesFiltered_fastMNN_Res06.png", sep = "/"),
plot = p, width = 10, height = 8, dpi = 100)
# 2 - Vulcano plot DEG Gene Therapy vs Control
hgenes <- c("Ccl8", "Pf4", "Ccl7", "Selenop", "Clec4a1", "Cd72", "Cd81", "Eno1", "Mrc1", "Msr1",
"Ly6c2", "Isg15", "Ier3", "Iigp1", "Dusp1", "Ier5", "Jun", "Gm4951")
res_dge <- read.xlsx(xlsxFile = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control-DGE.xlsx", sep = "/"),
sheet = "GeneTherapy-Control", rowNames = T)
vol <- plotVolcano(dge_res = res_dge, adj.pval.thr = 0.01, logfc.thr = 0, hgenes = hgenes)
ggsave(filename = paste(out_dir, "Fig3_2-MacrophagesFiltered_GeneTherapy_Ctrl_Volcano.png", sep = "/"), plot = vol, width = 10, height = 9, dpi = 100)
# 3 - UMAP plot Mono/macrophages split by treatment, res 0.6
mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
m_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings,
mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),])
m_umap_plot$RNA_Group2 <- factor(m_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(m_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
label <- data.frame(
UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5),
UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0),
RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")
)
p <- ggplot(m_umap_plot, aes(UMAP_1, UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none",
strip.text.x = element_text(size = 40)) +
#geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), size = 8) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("") +
facet_grid(.~RNA_Group2, labeller = label_parsed)
ggsave(filename = paste(out_dir, "Fig3_3-MacrophagesFiltered_fastMNN_Res06_split.png", sep = "/"),
plot = p, width = 28, height = 8, dpi = 100)
# 4 - Featureplot genes
# genes <- c("Pglyrp1", "Ace", "Plac8", "Ifitm6", "Stmn1", "Top2a", "Chil3", "Arg1", "Il1b", "Dusp1", "Tnf", "Ccl8", "Pf4", "Mrc1")
# png(filename = paste(out_dir, "Fig3_4-MacrophagesFiltered_fastMNN_FeatureplotGenes.png", sep = "/"), width = 1200, height = 1200, res = 100)
# FeaturePlot(object = mobj, features = genes, reduction = "umap", pt.size = 1, ncol = 4, cols = c("blue", "white", "red"), order = T) + NoLegend()
# dev.off()
#########################
#### Figure 3 Suppl. ####
#########################
# 1 - GeneTherapy Dotplot of enrichment pos HALLMARK
gt_enr_pos <- read.xlsx(xlsxFile = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control_enrichment_pos.xlsx", sep = "/"),
sheet = "enr_pos_HALLMARK")
ds_hall_enrich <- gt_enr_pos[,c("ID", "GeneRatio", "p.adjust", "Count")]
gsize <- as.numeric(sub("/\\d+$", "", as.character(ds_hall_enrich$GeneRatio)))
gcsize <- as.numeric(sub("^\\d+/", "", as.character(ds_hall_enrich$GeneRatio)))
ds_hall_enrich$GeneRatio = gsize/gcsize
ds_hall_enrich <- subset(ds_hall_enrich, p.adjust < 0.05)
ds_hall_enrich$ID <- gsub("HALLMARK_", "", ds_hall_enrich$ID)
#ds_hall_enrich$ID <- tolower(gsub("_", " ", ds_hall_enrich$ID))
print(head(ds_hall_enrich))
id_order <- ds_hall_enrich %>% select(ID, GeneRatio) %>% group_by(ID) %>% summarize_all(sum) %>% data.frame()
id_order <- id_order[order(id_order$GeneRatio, decreasing = F),]
ds_hall_enrich$ID <- factor(ds_hall_enrich$ID, levels = id_order$ID)
p <- ggplot(ds_hall_enrich, aes(x = GeneRatio, y = ID, size = Count, color = p.adjust)) +
theme_bw(base_size = 20) +
scale_color_gradient(low = "red", high = "blue", guide = guide_colorbar(reverse = T)) +
geom_point() +
ylab("")
ggsave(plot = p,
filename = paste(out_dir, "Fig3Suppl_1-MacrophagesFiltered_GT_Ctrl_Hallmark_enrichment_pos.png", sep = "/"),
width = 10, height = 8, dpi = 150)
# 2 - GSEA MacroPolarization
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
gmt.obj <- read.gmt("/TIGET/coltella_sc_new/reference/MacroPolarizationSig.gmt")
# DGE Marker
for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) {
deg_res <- macro_obj@misc$dgemarkers[[contr]]
deg_res.logfc <- as.vector(deg_res$avg_logFC)
names(deg_res.logfc) <- row.names(deg_res)
deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
print(deg_res.logfc[1:5])
# GSEA
contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = gmt.obj, nPerm = 10000)
contr.enricher.gsea@result[,1:10]
ptitle <- expression(IFN-alpha)
if (contr == "EF_Low-AB_Ctrl") { ptitle <- expression(IFN-alpha-DHFR) }
p <- gseaplot2(contr.enricher.gsea, geneSetID = c(1,2), base_size = 18,
title = ptitle,
color = c("red", "forestgreen"),
pvalue_table = T, nopadj = T, labels = c("M1 Macro. IL-4", bquote("M2 Macro. IFN"~gamma~"+ LPS")))
ggsave(filename = paste(out_dir, paste0("Fig3Suppl_2-MacrophagesFiltered_", contr, "_MacroPolarization_GSEA.png"), sep = "/"),
plot = p, width = 16, height = 8, dpi = 100)
}
# 3 - Heatmap marker genes Macrophages with labels
hm <- doSCHeatmapMean(fname = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"),
clustering = "RNA_snn_res.0.6", ntop = 15, num_color = 3,
high_genes = c("Saa3", "Ccl5", "Ly6i", "Tmem176b", "Cx3cr1", "Hexb", "Dusp1", "Socs3", "Il1b", "Plac8", "Il1r2", "Mgl2",
"Ccl8", "Cd63", "Il18bp", "Chil3", "Arg1", "Ccl24", "Pf4", "Mrc1", "Ms4a4c", "Ly6c2", "Cd74", "C1qb",
"Stmn1", "Pclaf", "Birc5", "Pglyrp1", "Gsr", "Clec4e"),
clu_lab = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), duplicated_labs = T)
png(paste(out_dir, "Fig3Suppl_3-MacrophagesFiltered_fastMNN_RNA_snn_res.0.6_top15-color3-wider2.png", sep = "/"),
width = 1800, height = 700, res = 150)
draw(hm)
dev.off()
# 4 - Feature plot hypoxia (modulescore) from https://www.gsea-msigdb.org/gsea/msigdb/cards/SEMENZA_HIF1_TARGETS
mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
hypoxia_sig <- c("Adm","Adra1b","Ak3","Aldoart2","Aldoc","Bnip3","Car9","Cdkn1a","Cited2","Cp","Edn1","Eno1","Epo","Flt1","Gm20425",
"Hk1","Hk2","Hmox1","Igf2","Igfbp1","Igfbp2","Igfbp3","Ldha","Nos2","P4ha2","Pfkl","Pgk1","Serpine1","Slc2a1","Slc2a3",
"Tfrc","Tgfb3","Tpi1","Vegfa")
hypoxia_sig_filt <- hypoxia_sig[hypoxia_sig %in% rownames(mobj)]
mobj <- AddModuleScore(object = mobj, features = list(HYPOXIA = hypoxia_sig_filt), ctrl = length(hypoxia_sig_filt), name = "HYPOXIA")
macro_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings,
mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),])
macro_umap_plot$RNA_Group2 <- factor(macro_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(macro_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
str(umap_plot)
p <- ggplot(macro_umap_plot, aes(x = UMAP_1, y = UMAP_2)) +
#geom_point(aes(colour = HYPOXIA1), size = 1, alpha = 0.8) +
geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= -0.25 & HYPOXIA1 < 0.25), aes(colour = HYPOXIA1), size = 1.5, alpha = 0.8) +
geom_point(data = subset(macro_umap_plot, HYPOXIA1 < -0.25), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) +
geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= 0.25 & HYPOXIA1 < 0.5), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) +
geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= 0.5), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) +
# scale_colour_gradient2(low="#2166ac",
# mid="gray95",
# high="#b2182b",
# name = "Expr.",
# midpoint = -0.05) +
scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"),
name = "Expr.",
#values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)),
values = rescale(c(-.5,-.4,-.1,.2,.7)),
guide = "colorbar") +
theme_void(base_size = 20) +
#ggtitle("INTERFERON ALPHA RESPONSE") +
theme(legend.position = "right",
strip.text.x = element_text(size = 40))
#xlab("UMAP 1") +
#ylab("UMAP 2") +
#facet_grid(.~RNA_Group2, labeller = label_parsed)
ggsave(filename = paste(out_dir, "Fig3Suppl_4-MacrophagesFiltered_Hypoxia.png", sep = "/"), plot = p, width = 10, height = 8, dpi = 100)
##################
#### Figure 4 ####
##################
# 1 - Monocle pseudotime
cds_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "03-fastMNN_2-pseudotime", "Monocle_cds_objects.rds", sep = "/"))
label <- data.frame(
UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5),
UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0),
RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")
)
# # Single plots
# prepare_plot <- function(cds) {
# S_matrix <- reducedDims(cds)$UMAP
# data_df <- data.frame(S_matrix[,c(1,2)])
# data_df$sample_name <- row.names(data_df)
# data_df <- as.data.frame(cbind(data_df, colData(cds)))
# data_df$RNA_Group2 <- factor(data_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
# levels(data_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
# #str(data_df)
# ica_space_df <- t(cds@principal_graph_aux$UMAP$dp_mst) %>%
# as.data.frame() %>%
# dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
# dplyr::mutate(sample_name = rownames(.),
# sample_state = rownames(.))
# dp_mst <- cds@principal_graph$UMAP
# edge_df <- dp_mst %>%
# igraph::as_data_frame() %>%
# dplyr::select_(source = "from", target = "to") %>%
# dplyr::left_join(ica_space_df %>%
# dplyr::select_(
# source="sample_name",
# source_prin_graph_dim_1="prin_graph_dim_1",
# source_prin_graph_dim_2="prin_graph_dim_2"),
# by = "source") %>%
# dplyr::left_join(ica_space_df %>%
# dplyr::select_(
# target="sample_name",
# target_prin_graph_dim_1="prin_graph_dim_1",
# target_prin_graph_dim_2="prin_graph_dim_2"),
# by = "target")
# p <- ggplot(data_df, aes(x = UMAP_1, y = UMAP_2)) +
# geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) +
# theme_nothing(font_size = 20) +
# #theme_bw(base_size = 20) +
# theme(legend.position = "none",
# strip.text.x = element_text(size = 40)) +
# geom_segment(data = edge_df,
# aes_string(x="source_prin_graph_dim_1",
# y="source_prin_graph_dim_2",
# xend="target_prin_graph_dim_1",
# yend="target_prin_graph_dim_2"),
# size = 1,
# color = "black",
# linetype = "solid",
# na.rm = TRUE) +
# #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
# geom_label(data = label, aes(label = label), size = 8) +
# #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
# xlab("") +
# ylab("") +
# facet_grid(.~RNA_Group2, labeller = label_parsed)
# return(p)
# }
# p <- prepare_plot(cds_obj$AB_Ctrl)
# ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Ctrl_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"),
# plot = p, width = 10, height = 8, dpi = 100)
# p <- prepare_plot(cds_obj$CD_High)
# ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_High_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"),
# plot = p, width = 10, height = 8, dpi = 100)
# p <- prepare_plot(cds_obj$EF_Low)
# ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Low_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"),
# plot = p, width = 10, height = 8, dpi = 100)
# UMAP Data
data_df_ctrl <- data.frame(reducedDims(cds_obj$AB_Ctrl)$UMAP[,c(1,2)])
data_df_ctrl$sample_name <- row.names(data_df_ctrl)
data_df_ctrl <- as.data.frame(cbind(data_df_ctrl, colData(cds_obj$AB_Ctrl)))
data_df_high <- data.frame(reducedDims(cds_obj$CD_High)$UMAP[,c(1,2)])
data_df_high$sample_name <- row.names(data_df_high)
data_df_high <- as.data.frame(cbind(data_df_high, colData(cds_obj$CD_High)))
data_df_low <- data.frame(reducedDims(cds_obj$EF_Low)$UMAP[,c(1,2)])
data_df_low$sample_name <- row.names(data_df_low)
data_df_low <- as.data.frame(cbind(data_df_low, colData(cds_obj$EF_Low)))
data_df <- rbind(data_df_ctrl, rbind(data_df_high, data_df_low))
data_df$RNA_Group2 <- factor(data_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(data_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
# Graph Data
ica_space_df_ctrl <- t(cds_obj$AB_Ctrl@principal_graph_aux$UMAP$dp_mst) %>%
as.data.frame() %>%
dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
dplyr::mutate(sample_name = rownames(.),
sample_state = rownames(.))
dp_mst_ctrl <- cds_obj$AB_Ctrl@principal_graph$UMAP
edge_df_ctrl <- dp_mst_ctrl %>%
igraph::as_data_frame() %>%
dplyr::select_(source = "from", target = "to") %>%
dplyr::left_join(ica_space_df_ctrl %>%
dplyr::select_(
source="sample_name",
source_prin_graph_dim_1="prin_graph_dim_1",
source_prin_graph_dim_2="prin_graph_dim_2"),
by = "source") %>%
dplyr::left_join(ica_space_df_ctrl %>%
dplyr::select_(
target="sample_name",
target_prin_graph_dim_1="prin_graph_dim_1",
target_prin_graph_dim_2="prin_graph_dim_2"),
by = "target")
ica_space_df_high <- t(cds_obj$CD_High@principal_graph_aux$UMAP$dp_mst) %>%
as.data.frame() %>%
dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
dplyr::mutate(sample_name = rownames(.),
sample_state = rownames(.))
dp_mst_high <- cds_obj$CD_High@principal_graph$UMAP
edge_df_high <- dp_mst_high %>%
igraph::as_data_frame() %>%
dplyr::select_(source = "from", target = "to") %>%
dplyr::left_join(ica_space_df_high %>%
dplyr::select_(
source="sample_name",
source_prin_graph_dim_1="prin_graph_dim_1",
source_prin_graph_dim_2="prin_graph_dim_2"),
by = "source") %>%
dplyr::left_join(ica_space_df_high %>%
dplyr::select_(
target="sample_name",
target_prin_graph_dim_1="prin_graph_dim_1",
target_prin_graph_dim_2="prin_graph_dim_2"),
by = "target")
ica_space_df_low <- t(cds_obj$EF_Low@principal_graph_aux$UMAP$dp_mst) %>%
as.data.frame() %>%
dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
dplyr::mutate(sample_name = rownames(.),
sample_state = rownames(.))
dp_mst_low <- cds_obj$EF_Low@principal_graph$UMAP
edge_df_low <- dp_mst_low %>%
igraph::as_data_frame() %>%
dplyr::select_(source = "from", target = "to") %>%
dplyr::left_join(ica_space_df_low %>%
dplyr::select_(
source="sample_name",
source_prin_graph_dim_1="prin_graph_dim_1",
source_prin_graph_dim_2="prin_graph_dim_2"),
by = "source") %>%
dplyr::left_join(ica_space_df_low %>%
dplyr::select_(
target="sample_name",
target_prin_graph_dim_1="prin_graph_dim_1",
target_prin_graph_dim_2="prin_graph_dim_2"),
by = "target")
edge_df_ctrl$RNA_Group <- "AB_Ctrl"
edge_df_high$RNA_Group <- "CD_High"
edge_df_low$RNA_Group <- "EF_Low"
edge_df <- rbind(edge_df_ctrl, rbind(edge_df_high, edge_df_low))
edge_df$RNA_Group2 <- factor(edge_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(edge_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
p <- ggplot(data_df, aes(x = UMAP_1, y = UMAP_2)) +
geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) +
theme_nothing(font_size = 20) +
#theme_bw(base_size = 20) +
theme(legend.position = "none",
strip.text.x = element_text(size = 40)) +
geom_segment(data = edge_df,
aes_string(x="source_prin_graph_dim_1",
y="source_prin_graph_dim_2",
xend="target_prin_graph_dim_1",
yend="target_prin_graph_dim_2"),
size = 1,
color = "black",
linetype = "solid",
na.rm = TRUE, inherit.aes=FALSE) +
#geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) +
geom_label(data = label, aes(label = label), size = 8) +
#guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) +
xlab("") +
ylab("") +
facet_grid(.~RNA_Group2, labeller = label_parsed)
ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"),
plot = p, width = 28, height = 8, dpi = 100)
# 3 - Feature plot IFN_ALPHA (modulescore).
mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1",
"Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr",
"Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7",
"Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b",
"Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9",
"Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1",
"Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars")
inf_alpha_filt <- inf_alpha[inf_alpha %in% rownames(mobj)]
mobj <- AddModuleScore(object = mobj, features = list(IFNa = inf_alpha_filt), ctrl = length(inf_alpha_filt), name = "IFNa")
macro_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings,
mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),])
macro_umap_plot$RNA_Group2 <- factor(macro_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low"))
levels(macro_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR")
summary(macro_umap_plot$IFNa1)
p <- ggplot(macro_umap_plot, aes(x = UMAP_1, y = UMAP_2)) +
#geom_point(aes(colour = IFNa1), size = 1.5, alpha = 1) +
geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.1 & IFNa1 < 0.3), aes(colour = IFNa1), size = 2, alpha = 0.8) +
geom_point(data = subset(macro_umap_plot, IFNa1 < 0.1), aes(colour = IFNa1), size = 2, alpha = 1) +
#geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.25 & IFNa1 < 0.5), aes(colour = IFNa1), size = 1.5, alpha = 1) +
geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.3), aes(colour = IFNa1), size = 2, alpha = 1) +
# scale_colour_gradient2(low="#2166ac",
# mid="gray95",
# high="#b2182b",
# name = "Expr.",
# midpoint = 0.25) +
# scale_color_gradientn(colours = c("blue","blue","grey95","red","red"),
# #values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)),
# values = rescale(c(-.1,0,.2,.5,.65)),
# guide = "colorbar") +
scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"),
name = "Expr.",
#values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)),
values = rescale(c(-.1,0,.2,.35,.65)),
guide = "colorbar") +
theme_void(base_size = 20) +
#ggtitle("INTERFERON ALPHA RESPONSE") +
theme(legend.position = "right",
strip.text.x = element_text(size = 40)) +
#xlab("UMAP 1") +
#ylab("UMAP 2") +
facet_grid(.~RNA_Group2, labeller = label_parsed)
ggsave(filename = paste(out_dir, "Fig4_3-MacrophagesFiltered_IFNa.png", sep = "/"), plot = p, width = 28, height = 8, dpi = 100)
# # 2 - GSEA signatures M1 e M2 DEG High vs. CTRL e Low vs. CTRL (macro)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/CD_High-AB_Ctrl/Alternatively_Activated_IL-4_Signature_GSEA)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/CD_High-AB_Ctrl/Classically_Activated_LPS-IFN_Signature_GSEA)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/EF_Low-AB_Ctrl/Alternatively_Activated_IL-4_Signature_GSEA)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/EF_Low-AB_Ctrl/Classically_Activated_LPS-IFN_Signature_GSEA)
# mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
# pol_gmt <- read.gmt(paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "MacroPolarizationSig.gmt", sep = "/"))
# for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) {
# deg_res <- mobj@misc$dgemarkers[[contr]]
# deg_res.logfc <- as.vector(deg_res$avg_logFC)
# names(deg_res.logfc) <- row.names(deg_res)
# deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
# deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
# print(deg_res.logfc[1:5])
# # GSEA
# contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = pol_gmt, nPerm = 10000, pvalueCutoff = 0.2)
# for (p in seq(1, nrow(contr.enricher.gsea@result))) {
# d <- as.character(contr.enricher.gsea@result[p, "ID"])
# if (contr == "CD_High-AB_Ctrl") {
# p <- gseaplot2(contr.enricher.gsea, geneSetID = p, base_size = 18,
# title = bquote(.(d)~": " ~ IFN-alpha ~ " vs. Control"),
# color = c("red"),
# pvalue_table = T)
# } else {
# p <- gseaplot2(contr.enricher.gsea, geneSetID = p, base_size = 18,
# title = bquote(.(d)~": " ~ IFN-alpha-DHFR ~ " vs. Control"),
# color = c("red"),
# pvalue_table = T)
# }
# ggsave(filename = paste(out_dir, paste0("Fig4_2_MacrophagesFiltered_", contr, "_", d, "_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
# }
# }
# # 4 - GSEA signatures M1 e M2 on pos marker genes of cluster 6 (res 0.6)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/ RNA_snn_res.0.6_Macro_polarization_GSEA/ clu6_Alternatively_Activated_IL-4_Signature_GSEA)
# #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/ RNA_snn_res.0.6_Macro_polarization_GSEA/clu6_Classically_Activated_LPS-IFN_Signature_GSEA)
# mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
# pol_gmt <- read.gmt(paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "MacroPolarizationSig.gmt", sep = "/"))
# deg_res <- subset(mobj@misc$markers_full$RNA_snn_res.0.6, cluster == "6")
# deg_res.logfc <- as.vector(deg_res$avg_logFC)
# names(deg_res.logfc) <- deg_res$gene
# deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE)
# deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))]
# print(deg_res.logfc[1:5])
# # GSEA
# contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = pol_gmt, nPerm = 10000, pvalueCutoff = 0.2)
# for (cc in c(1, 2)) {
# d <- as.character(contr.enricher.gsea@result[cc, "ID"])
# p <- gseaplot2(contr.enricher.gsea, geneSetID = cc, base_size = 18,
# title = bquote(.(d)),
# color = c("red"),
# pvalue_table = T)
# ggsave(filename = paste(out_dir, paste0("Fig4_4_MacrophagesFiltered_Clu6_", d, "_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100)
# }
# # 3 - Featureplot Ccl8, Mrc1
# mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
# ccl8_expr_plot <- cbind(as.data.frame(mobj@reductions$umap@cell.embeddings),
# "Expr" = as.numeric(mobj@assays$RNA@data["Ccl8", rownames(mobj@reductions$umap@cell.embeddings)]))
# p <- ggplot(ccl8_expr_plot, aes(x = UMAP_1, y = UMAP_2)) +
# geom_point(aes(colour = Expr), size = 1, alpha = 0.8) +
# scale_colour_gradient2(low="#2166ac",
# mid="grey",
# high="#b2182b",
# name = "Expr.",
# midpoint = (((max(ccl8_expr_plot$Expr)-min(ccl8_expr_plot$Expr))/2)-0.7)) +
# theme_bw(base_size = 20) +
# #ggtitle("Ccl8 Expression") +
# theme(legend.position = "right") +
# xlab("UMAP 1") +
# ylab("UMAP 2")
# ggsave(filename = paste(out_dir, "Fig5_3_MacrophagesFiltered_Ccl8_Expr.png", sep = "/"),
# plot = p, width = 10, height = 8, dpi = 100)
# mrc1_expr_plot <- cbind(as.data.frame(mobj@reductions$umap@cell.embeddings),
# "Expr" = mobj@assays$RNA@data["Mrc1", rownames(mobj@reductions$umap@cell.embeddings)])
# p <- ggplot(mrc1_expr_plot, aes(x = UMAP_1, y = UMAP_2)) +
# geom_point(data = subset(mrc1_expr_plot, Expr < 1), aes(colour = Expr), size = 1, alpha = 0.8) +
# geom_point(data = subset(mrc1_expr_plot, Expr >= 1), aes(colour = Expr), size = 1.1, alpha = 1) +
# scale_colour_gradient2(low="#2166ac",
# mid="grey",
# high="#b2182b",
# name = "Expr.",
# midpoint = (((max(mrc1_expr_plot$Expr)-min(mrc1_expr_plot$Expr))/2)-0.4)) +
# theme_bw(base_size = 20) +
# #ggtitle("Mrc1 Expression") +
# theme(legend.position = "right") +
# xlab("UMAP 1") +
# ylab("UMAP 2")
# ggsave(filename = paste(out_dir, "Fig5_3_MacrophagesFiltered_Mrc1_Expr.png", sep = "/"),
# plot = p, width = 10, height = 8, dpi = 100)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment