Skip to content

GitLab

Commit 918c0314 authored by Stefano Beretta's avatar Stefano Beretta
Browse files

Added initial scripts and data.

parent b32186fe
Hide whitespace changes
Inline Side-by-side
0_PreProcessData.R 0 → 100644
View file @ 918c0314
library(SeuratObject)
library(Seurat)
library(harmony)
library(DoubletFinder)
library(grid)
library(cowplot)
library(gridExtra)
library(RColorBrewer)
library(patchwork)
library(scales)
#########################
### Utility Functions ###
#########################
define_region <- function(row, col){
viewport(layout.pos.row = row, layout.pos.col = col)
}
########################
### General Settings ###
########################
# Working dir
wdir <- "squadrito_livertumor2022_scrnaseq"
# Raw data dir (download from GEO)
raw_dir <- "GEO_data"
# Results dir
out_dir <- paste(wdir, "results", sep = "/")
dir.create(path = out_dir, showWarnings = F)
# Samples
sample_info <- data.frame("samples" = c("Sample3", "Sample4", "Sample7", "Sample10", "Sample11", "Sample14", "Sample19", "Sample22"),
"groups" = c("CTRL", "IFN_R", "CTRL", "IFN_PR", "IFN_PR", "IFN_R", "CTRL", "IFN_PR"))
# cell-cycle scoring - NicheNet convertion of Seurat cc.genes.updated.2019 lists
s.genes <- c("Mcm5","Pcna","Tyms","Fen1","Mcm7","Mcm4","Rrm1","Ung","Gins2","Mcm6",
"Cdca7","Dtl","Prim1","Uhrf1","Cenpu","Hells","Rfc2","Polr1b","Nasp",
"Rad51ap1","Gmnn","Wdr76","Slbp","Ccne2","Ubr7","Pold3","Msh2","Atad2",
"Rad51","Rrm2","Cdc45","Cdc6","Exo1","Tipin","Dscc1","Blm","Casp8ap2",
"Usp1","Clspn","Pola1","Chaf1b","Mrpl36","E2f8")
g2m.genes <- c("Hmgb2","Cdk1","Nusap1","Ube2c","Birc5","Tpx2","Top2a","Ndc80",
"Cks2","Nuf2","Cks1b","Mki67","Tmpo","Cenpf","Tacc3","Pimreg",
"Smc4","Ccnb2","Ckap2l","Ckap2","Aurkb","Bub1","Kif11","Anp32e",
"Tubb4b","Gtse1","Kif20b","Hjurp","Cdca3","Jpt1","Cdc20","Ttk",
"Cdc25c","Kif2c","Rangap1","Ncapd2","Dlgap5","Cdca2","Cdca8","Ect2",
"Kif23","Hmmr","Aurka","Psrc1","Anln","Lbr","Ckap5","Cenpe","Ctcf",
"Nek2","G2e3","Gas2l3","Cbx5","Cenpa")
reg_vars <- c("CC.Difference", "percent.mt", "nCount_RNA")
#############################
### Create Seurat Objects ###
#############################
for (sample in sample_info$samples) {
print(sample)
group <- as.character(sample_info[sample_info$samples == sample, "groups"])
mtx_file <- paste(raw_dir, paste0(sample, "_matrix.mtx.gz"), sep = "/")
if (!file.exists(mtx_file)) {
stop(paste("Cannot find", sample, "MTX file:", mtx_file))
}
bcode_file <- paste(raw_dir, paste0(sample, "_barcodes.tsv.gz"), sep = "/")
if (!file.exists(bcode_file)) {
stop(paste("Cannot find", sample, "barcode file:", bcode_file))
}
feature_file <- paste(raw_dir, paste0(sample, "_features.tsv.gz"), sep = "/")
if (!file.exists(feature_file)) {
stop(paste("Cannot find", sample, "features file:", feature_file))
}
sample_mtx <- ReadMtx(mtx = mtx_file, cells = bcode_file, features = feature_file)
colnames(x = sample_mtx) <- paste(sample, colnames(x = sample_mtx), sep = '_') # add id tag to cellnames
sample_obj <- CreateSeuratObject(counts = sample_mtx, project = sample)
sample_obj[["RNA_Group"]] <- group
sample_dir <- paste(out_dir, sample, sep = "/")
dir.create(path = sample_dir, showWarnings = F)
saveRDS(object = sample_obj, file = paste(sample_dir, paste0(sample, "_final.rds"), sep = "/"))
}
####################################################
### Exploratory analysis for parameter selection ###
####################################################
# Perform DoubletFinder on each samples with different
# combinations of parameters
for (sample in sample_info$samples) {
print(sample)
sample_dir <- paste(out_dir, sample, sep = "/")
plot_dir <- paste(sample_dir, "plots", sep = "/")
dir.create(path = plot_dir, showWarnings = F)
obj <- readRDS(paste(sample_dir, paste0(sample, "_final.rds"), sep = "/"))
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^mt-")
obj <- CellCycleScoring(object = obj, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
obj@meta.data$CC.Difference <- obj@meta.data$S.Score - obj@meta.data$G2M.Score
Idents(obj) <- "orig.ident"
png(paste(plot_dir, paste0(sample, "_prefilter.png"), sep = "/"), width = 1200, height = 1000)
print(VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.1))
dev.off()
obj@misc$analysis_params$min.feature <- 1000
obj@misc$analysis_params$max.feature <- 6000
obj@misc$analysis_params$max.pc.mito <- 10
obj <- subset(obj, subset = nFeature_RNA > obj@misc$analysis_params$min.feature &
nFeature_RNA < obj@misc$analysis_params$max.feature &
percent.mt < obj@misc$analysis_params$max.pc.mito)
png(paste(plot_dir, paste0(sample, "_postfilter.png"), sep = "/"), width = 1200, height = 1000)
print(VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.1))
dev.off()
obj <- SCTransform(obj, vars.to.regress = c("CC.Difference", "percent.mt", "nCount_RNA"), verbose = FALSE)
DefaultAssay(object = obj) <- "SCT"
obj@misc$analysis_params$max_pca <- 50
obj <- RunPCA(object = obj, npcs = obj@misc$analysis_params$max_pca)
png(paste(plot_dir, paste0(sample, "_PCA_elbowplot_SCT.png"), sep = "/"), width = 1200, height = 1000)
print(ElbowPlot(object = obj, ndims = obj@misc$analysis_params$max_pca, reduction = "pca"))
dev.off()
obj@misc$analysis_params$num_pc <- 30
obj <- RunUMAP(object = obj, dims = 1:obj@misc$analysis_params$num_pc)
## pK Identification (no ground-truth) --------------------------------------------------
sweep.res.obj <- paramSweep_v3(obj, PCs = 1:obj@misc$analysis_params$num_pc, sct = TRUE)
sweep.stats_obj <- summarizeSweep(sweep.res.obj, GT = FALSE)
bcmvn_obj <- find.pK(sweep.stats_obj)
dev.off()
pK <- as.numeric(as.character(bcmvn_obj$pK))
BCmetric <- bcmvn_obj$BCmetric
pK_choose <- pK[which(BCmetric %in% max(BCmetric))]
png(paste(plot_dir, paste0(sample, "_DF_SCT_pK_selection.png"), sep = "/"), width = 1200, height = 1000)
par(mar = c(5,4,4,8)+1, cex.main = 1.2, font.main = 2)
plot(x = pK, y = BCmetric, pch = 16, type = "b", col = "blue", lty = 1)
abline(v = pK_choose, lwd = 2, col = 'red', lty = 2)
title("The BCmvn distributions")
text(pK_choose, max(BCmetric), as.character(pK_choose), pos = 4, col = "red")
dev.off()
for (pK_select in pK) {
for (doub_pct in c(0.05, 0.07, 0.09)) {
nExp <- round(ncol(obj) * doub_pct) #expected % doublet - approx
obj <- doubletFinder_v3(obj, pN = 0.25, pK = pK_select, nExp = nExp, PCs = 1:obj@misc$analysis_params$num_pc, sct = TRUE)
df_ident <- grep("DF.classification", colnames(obj@meta.data), value = T)[1]
Idents(obj) <- df_ident
new_cn <- paste0("DF_SCT_nExp", doub_pct, "_pK", pK_select)
colnames(obj@meta.data)[colnames(obj@meta.data) == df_ident] <- new_cn
str(obj@meta.data)
dp_df <- DimPlot(object = obj, reduction = "umap", pt.size = 1, group.by = new_cn, split.by = new_cn, ncol= 2) + NoLegend()
dp_phase <- DimPlot(object = obj, reduction = "umap", pt.size = 1, group.by = "Phase", split.by = new_cn, ncol= 3)
fp1 <- FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Cd19"), split.by = new_cn, order = T, ncol= 2)
fp2 <- FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Cd3g"), split.by = new_cn, order = T, ncol= 2)
fp3 <- FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Itgam"), split.by = new_cn, order = T, ncol = 2)
fp4 <- FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Alas2"), split.by = new_cn, order = T, ncol = 2)
png(paste(plot_dir, paste0("DF_SCT_nExp", doub_pct, "_pK", pK_select, ".png"), sep = "/"), height = 1500, width = 1700)
plot.new()
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 2)))
print(dp_phase, vp = define_region(row = 1, col = 1))
print(dp_df, vp = define_region(row = 1, col = 2))
print(fp1, vp = define_region(row = 2, col = 1))
print(fp2, vp = define_region(row = 2, col = 2))
print(fp3, vp = define_region(row = 3, col = 1))
print(fp4, vp = define_region(row = 3, col = 2))
dev.off()
}
}
saveRDS(object = obj, paste(sample_dir, paste0(sample, "_afterDoubletFinder.rds"), sep = "/"))
}
# Prepare some plots on DoubletFinder results
for (sample in sample_info$samples) {
print(sample)
sample_dir <- paste(out_dir, sample, sep = "/")
plot_dir <- paste(sample_dir, "plots", sep = "/")
obj = readRDS(paste(sample_dir, paste0(sample, "_afterDoubletFinder.rds"), sep = "/"))
for (pK_select in c(0.005, seq(0.01, 0.3, 0.01))) {
for (doub_pct in c(0.05, 0.07, 0.09)) {
nExp <- round(ncol(obj) * doub_pct)
new_cn <- paste0("DF_SCT_nExp", doub_pct, "_pK", pK_select)
subtitle <- list("Full" = "", "Cd19" = "", "Cd3g" = "", "Itgam" = "", "Alas2" = "")
num_title <- 0
for (md_val in names(table(obj[[new_cn]]))) {
if (num_title == 0) {
num_title <- 1
subtitle[["Full"]] <- paste(subtitle[["Full"]], paste(md_val, table(obj[[new_cn]])[[md_val]], sep = " : "), sep = "\n")
for (gg in c("Cd19", "Cd3g", "Itgam", "Alas2")) {
subtitle[[gg]] <- paste0(gg, ">0 :\t", subtitle[[gg]], paste(md_val, table(subset(obj, subset = !!sym(gg) > 0)[[new_cn]])[[md_val]], sep = "= "))
}
} else {
subtitle[["Full"]] <- paste(subtitle[["Full"]], paste(md_val, table(obj[[new_cn]])[[md_val]], sep = " : "), sep = "\t|\t")
for (gg in c("Cd19", "Cd3g", "Itgam", "Alas2")) {
subtitle[[gg]] <- paste(subtitle[[gg]], paste(md_val, table(subset(obj, subset = !!sym(gg) > 0)[[new_cn]])[[md_val]], sep = "= "), sep = "\t| ")
}
}
num_title <- num_title + 1
}
Idents(obj) <- new_cn
dp_df <- DimPlot(object = obj, reduction = "umap", pt.size = 1, group.by = new_cn, split.by = new_cn, ncol= 2) + NoLegend() + ggtitle(subtitle[["Full"]])
dp_phase <- DimPlot(object = obj, reduction = "umap", pt.size = 1, group.by = "Phase", split.by = new_cn, ncol= 3)
fp1 <- plot_grid(ggdraw() + draw_label(subtitle[["Cd19"]]),
FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Cd19"), split.by = new_cn, order = T, ncol= 2),
ncol = 1, rel_heights = c(0.1, 1))
fp2 <- plot_grid(ggdraw() + draw_label(subtitle[["Cd3g"]]),
FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Cd3g"), split.by = new_cn, order = T, ncol= 2),
ncol = 1, rel_heights = c(0.1, 1))
fp3 <- plot_grid(ggdraw() + draw_label(subtitle[["Itgam"]]),
FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Itgam"), split.by = new_cn, order = T, ncol = 2),
ncol = 1, rel_heights = c(0.1, 1))
fp4 <- plot_grid(ggdraw() + draw_label(subtitle[["Alas2"]]),
FeaturePlot(object = obj, reduction = "umap", pt.size = 1, features = c("Alas2"), split.by = new_cn, order = T, ncol = 2),
ncol = 1, rel_heights = c(0.1, 1))
png(paste(plot_dir, paste0("DF_SCT_nExp", doub_pct, "_pK", pK_select, ".png"), sep = "/"), height = 1500, width = 1700)
plot.new()
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 2)))
print(dp_phase, vp = define_region(row = 1, col = 1))
print(dp_df, vp = define_region(row = 1, col = 2))
print(fp1, vp = define_region(row = 2, col = 1))
print(fp2, vp = define_region(row = 2, col = 2))
print(fp3, vp = define_region(row = 3, col = 1))
print(fp4, vp = define_region(row = 3, col = 2))
dev.off()
}
}
}
# Select the best combinations of parameters of DoubletFinder
# for each sample
samples <- list("Sample3" = c("nExp" = 0.07, "pK" = 0.005),
"Sample4" = c("nExp" = 0.09, "pK" = 0.005),
"Sample7" = c("nExp" = 0.07, "pK" = 0.005),
"Sample10" = c("nExp" = 0.09, "pK" = 0.01),
"Sample11" = c("nExp" = 0.07, "pK" = 0.005),
"Sample14" = c("nExp" = 0.09, "pK" = 0.005),
"Sample19" = c("nExp" = 0.05, "pK" = 0.005),
"Sample22" = c("nExp" = 0.09, "pK" = 0.02)
)
# Add metadata information on doublets in each sample
full_DoubFinder <- data.frame()
for (sample in names(samples)) {
print(sample)
sample_dir <- paste(out_dir, sample, sep = "/")
obj <- readRDS(paste(sample_dir, paste0(sample, "_afterDoubletFinder.rds"), sep = "/"))
nExp <- round(ncol(obj) * samples[[sample]][["nExp"]]) #expected % doublet - approx
obj@misc$analysis_params$doubletFinder <- list(pN = 0.25,
pK = samples[[sample]][["pK"]],
nExp = nExp,
nExpPerc = samples[[sample]][["nExp"]],
PCs = obj@misc$analysis_params$num_pc)
new_cn <- paste0("DF_SCT_nExp", obj@misc$analysis_params$doubletFinder$nExpPerc, "_pK", obj@misc$analysis_params$doubletFinder$pK)
df_sel_classif <- obj@meta.data[, new_cn, drop = F]
print(head(df_sel_classif))
colnames(df_sel_classif) <- "DoubFinderClassif"
full_DoubFinder <- rbind(full_DoubFinder, df_sel_classif)
saveRDS(object = obj, paste(sample_dir, paste0(sample, "_afterDoubletFinder.rds"), sep = "/"))
}
############################################
### Create Full (Combined) Seurat Object ###
############################################
samples_obj <- list()
for (sample in names(samples)) {
print(sample)
sample_dir <- paste(out_dir, sample, sep = "/")
samples_obj[[sample]] <- readRDS(file = paste(sample_dir, paste0(sample, "_final.rds"), sep = "/"))
}
Full_obj <- merge(x = samples_obj[[1]], y = samples_obj[2:length(samples_obj)])
dir.create(path = paste(out_dir, "Full", sep = "/"))
saveRDS(object = Full_obj, paste(out_dir, "Full", "Full_final.rds", sep = "/"))
plot_dir <- paste(out_dir, "Full", "plots", sep = "/")
dir.create(path = plot_dir, showWarnings = F)
# Assign MT genes
Full_obj[["percent.mt"]] <- PercentageFeatureSet(Full_obj, pattern = "^mt-")
# Identify tresholds visualy and apply tresholds
Idents(Full_obj) <- "RNA_Group"
png(filename = paste(plot_dir, "Full_Prefilter_Vlnplot.png", sep = "/"), width = 1500, height = 1000)
VlnPlot(Full_obj, features = reg_vars, ncol = 3, pt.size = 0)
dev.off()
png(filename = paste(plot_dir, "Full_Prefilter_ScatterMT.png", sep = "/"), width = 1200, height = 800)
FeatureScatter(Full_obj, feature1 = "nCount_RNA", feature2 = "percent.mt")
dev.off()
png(filename = paste(plot_dir, "Full_Prefilter_ScatterFeatures.png", sep = "/"), width = 1200, height = 800)
FeatureScatter(Full_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
dev.off()
Full_obj@misc$analysis_params$min.feature <- 1000
Full_obj@misc$analysis_params$max.feature <- 6000
Full_obj@misc$analysis_params$max.pc.mito <- 10
table(Full_obj$orig.ident)
dim(Full_obj)
Full_obj <- subset(Full_obj, subset = nFeature_RNA > Full_obj@misc$analysis_params$min.feature &
nFeature_RNA < Full_obj@misc$analysis_params$max.feature &
percent.mt < Full_obj@misc$analysis_params$max.pc.mito)
table(Full_obj$orig.ident)
dim(Full_obj)
png(filename = paste(plot_dir, "Full_Postfilter_Vlnplot.png", sep = "/"), width = 1500, height = 1000)
VlnPlot(Full_obj, features = reg_vars, ncol = 3, pt.size = 0)
dev.off()
png(filename = paste(plot_dir, "Full_Postfilter_ScatterMT.png", sep = "/"), width = 1200, height = 800)
FeatureScatter(Full_obj, feature1 = "nCount_RNA", feature2 = "percent.mt")
dev.off()
png(filename = paste(plot_dir, "Full_Postfilter_ScatterFeatures.png", sep = "/"), width = 1200, height = 800)
FeatureScatter(Full_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
dev.off()
## Import filtered object and process (SCTransform, runPCA)
# normalize, Scaledata are not computed, since we are not doing diff expression at this point
Full_obj <- SCTransform(Full_obj, vars.to.regress = reg_vars, verbose = TRUE)
DefaultAssay(object = Full_obj) <- "SCT"
Full_obj@misc$analysis_params$max_pca <- 50
Full_obj <- RunPCA(object = Full_obj, npcs = Full_obj@misc$analysis_params$max_pca)
png(filename = paste(plot_dir, "Full_ElbowPlot.png", sep = "/"), width = 1200, height = 800)
ElbowPlot(Full_obj, ndims = Full_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
Full_obj@misc$analysis_params$num_pc <- 40
## Import Doublets
nrow(full_DoubFinder)
dim(Full_obj)
full_DoubFinder_md <- full_DoubFinder$DoubFinderClassif
names(full_DoubFinder_md) <- rownames(full_DoubFinder)
Full_obj <- AddMetaData(object = Full_obj, metadata = full_DoubFinder_md, col.name = "DoubletFinder_classification")
colnames(Full_obj@meta.data)
subtitle <- ""
num_title <- 0
for (md_val in names(table(Full_obj[["DoubletFinder_classification"]]))) {
if (num_title == 0) {
num_title <- 1
subtitle <- paste(subtitle, paste(md_val, table(Full_obj[["DoubletFinder_classification"]])[[md_val]], sep = " : "), sep = "\n")
} else {
subtitle <- paste(subtitle, paste(md_val, table(Full_obj[["DoubletFinder_classification"]])[[md_val]], sep = " : "), sep = "\t|\t")
}
num_title <- num_title + 1
}
png(paste(plot_dir, "Full_DoubletFinder_classification.png", sep = "/"), height = 1200, width = 1400, res = 100)
DimPlot(object = Full_obj, reduction = "umap", pt.size = 1, group.by = "DoubletFinder_classification") + ggtitle(subtitle) +
annotation_custom(tableGrob(table(Full_obj$orig.ident, Full_obj$DoubletFinder_classification)), xmin=10, ymin=10)
dev.off()
png(paste(plot_dir, "Full_DoubletFinder_classification_split.png", sep = "/"), height = 1200, width = 1800)
DimPlot(object = Full_obj, reduction = "umap", pt.size = 1, group.by = "DoubletFinder_classification", split.by = "DoubletFinder_classification") + ggtitle(subtitle)
dev.off()
#####################
## Remove Doublets ##
#####################
Idents(Full_obj) <- "DoubletFinder_classification"
table(Full_obj$orig.ident, Full_obj$DoubletFinder_classification)
dim(Full_obj)
Full_obj <- subset(Full_obj, idents = "Singlet")
table(Full_obj$orig.ident)
dim(Full_obj)
## Process (normalize, Scaledata, SCTransform, runPCA)
# Scale data will be used when commputing top marker (we think)
Full_obj <- NormalizeData(object = Full_obj)
Full_obj <- ScaleData(object = Full_obj, vars.to.regress = reg_vars,
display.progress = T, features = rownames(Full_obj))
Full_obj <- SCTransform(Full_obj, vars.to.regress = reg_vars, verbose = TRUE)
DefaultAssay(object = Full_obj) <- "SCT"
Full_obj <- RunPCA(object = Full_obj, npcs = Full_obj@misc$analysis_params$max_pca)
ElbowPlot(Full_obj, ndims = Full_obj@misc$analysis_params$max_pca, reduction = "pca")
Full_obj <- RunUMAP(Full_obj, dims = 1:Full_obj@misc$analysis_params$num_pc, seed.use = 1,reduction = "pca")
Full_obj <- FindNeighbors(Full_obj ,reduction = "pca", dims = 1:Full_obj@misc$analysis_params$num_pc,
force.recalc = T, n.trees = 5, annoy.metric = "euclidean")
# Perform clustering on SCT
for (res in c(0.4, 0.6, 1.2)) {
Full_obj <- FindClusters(Full_obj , resolution = res)
}
# We use the SCT_snn_res.1.2. that was already in the Full object as calculated by you.
# Annotate new labels on FULL OBJ after doublet removal
Full_obj@meta.data$CellID <- Full_obj@assays$RNA@counts@Dimnames[[2]]
Full_obj@meta.data$Newlabels <- rep("Undefined", length(Full_obj@meta.data$orig.ident))
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "0"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "1"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "2"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "3"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "4"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "5"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "6"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "7"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "8"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "9"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "10"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "11"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "12"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "13"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "14"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "15"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "16"] <- "B cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "17"] <- "Neutrophils"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "18"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "19"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "20"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "21"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "22"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "23"] <- "Neutrophils"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "24"] <- "Erythroblasts"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "25"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "26"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "27"] <- "LSECs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "28"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "29"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "30"] <- "T and NK cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "31"] <- "B1a-like"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "32"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "33"] <- "Mast cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "34"] <- "Neutrophils"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "35"] <- "Cancer cells"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "36"] <- "APCs"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "37"] <- "Hepatocytes"
Full_obj@meta.data$Newlabels[Full_obj@meta.data$SCT_snn_res.1.2 == "38"] <- "Undefined"
##########################
## Export Table results ##
##########################
data_dir <- paste(wdir, "data", sep = "/")
dir.create(path = data_dir, showWarnings = F)
# Counts
full_counts <- Full_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "Full_final_DBR_labeled_counts.csv.gz", sep = "/"), "w")
write.csv(full_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
full_md <- Full_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_metadata.csv.gz", sep = "/"), "w")
write.csv(x = full_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
full_umap <- Full_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "Full_final_DBR_labeled_umap.csv.gz", sep = "/"), "w")
write.csv(x = full_umap, gz_out_umap)
close(gz_out_umap)
# Save final object
saveRDS(object = Full_obj, paste(out_dir, "Full", "Full_final_DBR_labeled.rds", sep = "/"))
\ No newline at end of file
1_Analysis_TNK.R 0 → 100644
View file @ 918c0314
library(SeuratObject)
library(Seurat)
library(harmony)
library(DoubletFinder)
library(grid)
library(cowplot)
library(gridExtra)
library(RColorBrewer)
library(patchwork)
library(scales)
# Working dir
wdir <- "squadrito_livertumor2022_scrnaseq"
# Data dir
data_dir <- paste(wdir, "data", sep = "/")
# Results dir
out_dir <- paste(wdir, "results", sep = "/")
dir.create(path = out_dir, showWarnings = F)
# Counts
gz_in_counts <- gzfile(paste(data_dir, "Full_final_DBR_labeled_counts.csv.gz", sep = "/"), "rt")
full_counts <- read.csv(gz_in_counts, row.names = 1)
close(gz_in_counts)
# Meta Data
gz_in_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_metadata.csv.gz", sep = "/"), "rt")
full_md <- read.csv(gz_in_md, row.names = 1)
close(gz_in_md)
# UMAP Coordinates
gz_in_umap <- gzfile(paste(data_dir, "Full_final_DBR_labeled_umap.csv.gz", sep = "/"), "rt")
full_umap <- read.csv(gz_in_umap, row.names = 1)
close(gz_in_umap)
Full_obj <- CreateSeuratObject(counts = full_counts, meta.data = full_md)
##########################
## T cells and NK cells ##
##########################
Full_obj <- SetIdent(object = Full_obj, value = "Newlabels")
TNK_obj <- subset(Full_obj, idents = c("T and NK cells"))
out_dir_TNK <- paste(out_dir, "T_NK_cells", sep = "/")
dir.create(path = out_dir_TNK, showWarnings = F)
saveRDS(TNK_obj, file = paste(out_dir_TNK, "T_NK_cells_DBR.rds", sep = "/"))
# Processing of the subset T cells, (Normalization, SCTransform, PCA and UMAP)
## We are not computing diff expression here so scale data is not done, we will do diff expression on the full obj
TNK_obj@meta.data$CellID <- TNK_obj@assays$RNA@counts@Dimnames[[2]] #CellID -> Will be used latter to be sure IDs are matching with full object
DefaultAssay(object = TNK_obj) <- "RNA"
TNK_obj <- NormalizeData(object = TNK_obj)
TNK_obj <- SCTransform(TNK_obj, vars.to.regress = c("nFeature_RNA", "nCount_RNA", "percent.mt"), verbose = TRUE)
DefaultAssay(object = TNK_obj) <- "SCT"
TNK_obj@misc$analysis_params$max_pca <- 50
TNK_obj <- RunPCA(object = TNK_obj, npcs = TNK_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_TNK, "T_and_NK_Elbowplot.png", sep = "/"))
ElbowPlot(TNK_obj, ndims = TNK_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
TNK_obj@misc$analysis_params$num_pc <- 35
TNK_obj <- RunUMAP(TNK_obj, dims = 1:TNK_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
TNK_obj <- FindNeighbors(TNK_obj, reduction = "pca", dims = 1:TNK_obj@misc$analysis_params$num_pc, force.recalc = T)
TNK_obj <- FindClusters(TNK_obj, resolution = 0.8)
png(filename = paste(out_dir_TNK, "T_and_NK_UMAP_res0.8.png", sep = "/"), width = 1200, height = 1000)
DimPlot(object = TNK_obj, reduction = "umap", group.by = "SCT_snn_res.0.8", pt.size = 1, label = T, label.size = 10)
dev.off()
#Annotate new labels rough on T cells
TNK_obj@meta.data$Newlabels_rough <- rep("Undefined",length(TNK_obj@meta.data$orig.ident))
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "0"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "1"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "2"] <- "NK"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "3"] <- "NKT"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "4"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "5"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "6"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "7"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "8"] <- "ILCs"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "9"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "10"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "11"] <- "CD4 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "12"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "13"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "14"] <- "ILCs"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "15"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "16"] <- "Undefined"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "17"] <- "CD8 T cells"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "18"] <- "NKT"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "19"] <- "Undefined"
TNK_obj@meta.data$Newlabels_rough[TNK_obj@meta.data$SCT_snn_res.0.8 == "20"] <- "gd T cells"
png(filename = paste(out_dir_TNK, "T_and_NK_UMAP_res0.8_Labels.png", sep = "/"), width = 1800, height = 1000)
DimPlot(object = TNK_obj, reduction = "umap", group.by = "SCT_snn_res.0.8", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = TNK_obj, reduction = "umap", group.by = "Newlabels_rough", pt.size = 1)
dev.off()
# Save final object
saveRDS(TNK_obj, file = paste(out_dir_TNK, "T_NK_cells_DBR.rds", sep = "/"))
## NKT - detailed annotation NKT cells ##
# extract NKT cells
TNK_obj <- SetIdent(object = TNK_obj, value = "Newlabels_rough")
TNKonly_obj <- subset(TNK_obj, idents = c("NKT"))
# There are very few cells left in the sample, after trying it many times we realised that it works much better.
# keeping the SCTransform from the T cell and NK object.
# Also the NKT belong to the T cells, so we need them to be somehow taken into consideration when spliting NKT subpopulaitons.
DefaultAssay(object = TNKonly_obj) <- "SCT"
TNKonly_obj@misc$analysis_params$max_pca <- 50
TNKonly_obj <- RunPCA(object = TNKonly_obj, npcs = TNKonly_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_TNK, "NKT_Elbowplot.png", sep = "/"))
ElbowPlot(TNKonly_obj, ndims = TNKonly_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
TNKonly_obj@misc$analysis_params$num_pc <- 30
TNKonly_obj <- RunUMAP(TNKonly_obj, dims = 1:TNKonly_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
TNKonly_obj <- FindNeighbors(TNKonly_obj ,reduction = "pca", dims = 1:TNKonly_obj@misc$analysis_params$num_pc, force.recalc = T)
TNKonly_obj <- FindClusters(TNKonly_obj , resolution = 0.6)
png(filename = paste(out_dir_TNK, "NKT_UMAP_res0.6.png", sep = "/"), width = 1200, height = 1000)
DimPlot(object = TNKonly_obj, reduction = "umap", group.by = "SCT_snn_res.0.6", pt.size = 2, label = T, label.size = 10)
dev.off()
# Cluster Annotation
TNKonly_obj@meta.data$Newlabels_detailed <- rep("Undefined",length(TNKonly_obj@meta.data$orig.ident))
TNKonly_obj@meta.data$Newlabels_detailed[TNKonly_obj@meta.data$SCT_snn_res.0.6 == "0"] <- "MAIT" #Trdc negative, Tcrg-C4 --> clearly alpha beta+ + signs of cytotoxic effect (Gzma/b+ and Prf1+)
TNKonly_obj@meta.data$Newlabels_detailed[TNKonly_obj@meta.data$SCT_snn_res.0.6 == "1"] <- "iNKT" #Sell negative CD8a high CD160 high
TNKonly_obj@meta.data$Newlabels_detailed[TNKonly_obj@meta.data$SCT_snn_res.0.6 == "2"] <- "MAIT" #alpha-beta+ and Ccr2+ Cxcr6+ but Cd7 neg
TNKonly_obj@meta.data$Newlabels_detailed[TNKonly_obj@meta.data$SCT_snn_res.0.6 == "3"] <- "gd NKT"
TNKonly_obj@meta.data$Newlabels_detailed[TNKonly_obj@meta.data$SCT_snn_res.0.6 == "4"] <- "Gzma NKT" #Trac neg, Tcrg-C4/Trdc pos
saveRDS(TNKonly_obj, file = paste(out_dir_TNK, "NKT_DBR.rds", sep = "/"))
# Reintegrate new NKT cell labels into T_NK_cell object
# extract labels from NKT cells
labels_NKT <- data.frame(ToSelect = TNKonly_obj@meta.data$CellID,
Classification = TNKonly_obj@meta.data$Newlabels_detailed)
# integrate NK T cell annotation into T_NK_:cell object
TNK_obj@meta.data$Newlabels_detailed <- TNK_obj@meta.data$Newlabels_rough
alllabels <- labels_NKT
TNK_obj@meta.data$Newlabels_detailed <- TNK_obj@meta.data$Newlabels_rough
orderDF <- data.frame(allIDs = TNK_obj@meta.data$CellID,
order=c(1:length(TNK_obj@meta.data$CellID)))
head(alllabels)
head(orderDF)
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
TNK_obj@meta.data$Newlabels_detailed[TNK_obj@meta.data$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Save final object (after new labels integration)
saveRDS(TNK_obj, file = paste(out_dir_TNK, "T_NK_cells_DBR.rds", sep = "/"))
## CD8 - detailed annotation CD8 T cells ##
# generate CD8 T cell subset
TNK_obj <- SetIdent(object = TNK_obj, value = "Newlabels_rough")
CD8_obj <- subset(TNK_obj, idents = c("CD8 T cells"))
# process CD8 T cell subset
DefaultAssay(object = CD8_obj) <- "SCT"
CD8_obj@misc$analysis_params$max_pca <- 50
CD8_obj <- RunPCA(object = CD8_obj, npcs = CD8_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_TNK, "CD8_Elbowplot.png", sep = "/"))
ElbowPlot(CD8_obj, ndims = CD8_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
CD8_obj@misc$analysis_params$num_pc <- 35
CD8_obj <- RunUMAP(CD8_obj, dims = 1:CD8_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
CD8_obj <- FindNeighbors(CD8_obj, reduction = "pca", dims = 1:CD8_obj@misc$analysis_params$num_pc, force.recalc = T)
CD8_obj <- FindClusters(CD8_obj, resolution = 0.3)
png(filename = paste(out_dir_TNK, "CD8_UMAP_res0.3.png", sep = "/"), width = 1200, height = 1000)
DimPlot(object = CD8_obj, reduction = "umap", group.by = "SCT_snn_res.0.3", pt.size = 2, label = T, label.size = 10)
dev.off()
# CD8 T cell annotation
CD8_obj@meta.data$Newlabels_detailed <- rep("Undefined",length(CD8_obj@meta.data$orig.ident))
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "0"] <- "CD8 Tsm-like cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "1"] <- "CD8 Teff1 cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "2"] <- "CD8 Tex cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "3"] <- "CD8 Tsm-like cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "4"] <- "T prol cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "5"] <- "CD8 Tm cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "6"] <- "CD8 Teff2 cells"
CD8_obj@meta.data$Newlabels_detailed[CD8_obj@meta.data$SCT_snn_res.0.3 == "7"] <- "CD8 Teff1 cells"
png(filename = paste(out_dir_TNK, "CD8_UMAP_res0.3_Labels.png", sep = "/"), width = 1800, height = 1000)
DimPlot(object = CD8_obj, reduction = "umap", group.by = "SCT_snn_res.0.3", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = CD8_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1)
dev.off()
# save CD8 T cell object
saveRDS(CD8_obj, file = paste(out_dir_TNK, "CD8Tcells_DBR.rds", sep = "/"))
# extract labels from CD8 T cells
labels_CD8 <- data.frame(ToSelect = CD8_obj@meta.data$CellID,
Classification = CD8_obj@meta.data$Newlabels_detailed)
# reintegrate CD8 T cell labels into T_NK_cell object
alllabels <- labels_CD8
orderDF <- data.frame(allIDs = TNK_obj@meta.data$CellID, order = c(1:length(TNK_obj@meta.data$CellID)))
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
TNK_obj@meta.data$Newlabels_detailed[TNK_obj@meta.data$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Save final object (after new labels integration)
saveRDS(TNK_obj, file = paste(out_dir_TNK, "T_NK_cells_DBR.rds", sep = "/"))
## CD4 - detailed annotation CD4 T cells ##
# generation of CD4 T cell object
TNK_obj <- SetIdent(object = TNK_obj, value = "Newlabels_rough")
CD4_obj <- subset(TNK_obj, idents = c("CD4 T cells"))
#process CD4 T cell subset
DefaultAssay(object = CD4_obj) <- "SCT"
CD4_obj@misc$analysis_params$max_pca <- 50
CD4_obj <- RunPCA(object = CD4_obj, npcs = CD4_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_TNK, "CD4_Elbowplot.png", sep = "/"))
ElbowPlot(CD4_obj, ndims = CD4_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
CD4_obj@misc$analysis_params$num_pc <- 35
CD4_obj <- RunUMAP(CD4_obj, dims = 1:CD4_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
CD4_obj <- FindNeighbors(CD4_obj, reduction = "pca", dims = 1:CD4_obj@misc$analysis_params$num_pc, force.recalc = T)
CD4_obj <- FindClusters(CD4_obj, resolution = 0.3)
png(filename = paste(out_dir_TNK, "CD4_UMAP_res0.3.png", sep = "/"), width = 1200, height = 1000)
DimPlot(object = CD4_obj, reduction = "umap", group.by = "SCT_snn_res.0.3", pt.size = 2, label = T, label.size = 10)
dev.off()
#CD4 T cell annotation
CD4_obj@meta.data$Newlabels_detailed <- rep("Undefined",length(CD4_obj@meta.data$orig.ident))
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "0"] <- "CD4 Tr1-like cells"
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "1"] <- "CD4 Tsm-like cells"
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "2"] <- "CD4 Treg cells"
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "3"] <- "CD4 Tr1-like cells"
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "4"] <- "CD4 Tex cells"
CD4_obj@meta.data$Newlabels_detailed[CD4_obj@meta.data$SCT_snn_res.0.3 == "5"] <- "CD4 Teff-like cells"
png(filename = paste(out_dir_TNK, "CD4_UMAP_res0.3_Labels.png", sep = "/"), width = 1800, height = 1000)
DimPlot(object = CD4_obj, reduction = "umap", group.by = "SCT_snn_res.0.3", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = CD4_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1)
dev.off()
# save CD4 T cell object
saveRDS(CD4_obj, file = paste(out_dir_TNK, "CD4Tcells_DBR.rds", sep = "/"))
#extract labels from CD4 T cells
labels_CD4 <- data.frame(ToSelect = CD4_obj@meta.data$CellID,
Classification = CD4_obj@meta.data$Newlabels_detailed)
#reintegrate CD4 T cell labels into T_NK_cell object
alllabels <- labels_CD4
orderDF <- data.frame(allIDs = TNK_obj@meta.data$CellID, order = c(1:length(TNK_obj@meta.data$CellID)))
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
TNK_obj@meta.data$Newlabels_detailed[TNK_obj@meta.data$CellID%in%alllabels$ToSelect] <- alllabels$Classification
png(filename = paste(out_dir_TNK, "T_and_NK_UMAP_res0.8_Labels_Detailed.png", sep = "/"), width = 2700, height = 1000)
DimPlot(object = TNK_obj, reduction = "umap", group.by = "SCT_snn_res.0.8", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = TNK_obj, reduction = "umap", group.by = "Newlabels_rough", pt.size = 1) |
DimPlot(object = TNK_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1, label = T, label.size = 6)
dev.off()
# Create files to reintegrate T_NK_cell labels into full object
# extract labels from T_NK_cell label object
T_NK_cell_labels_rough <- data.frame(ToSelect = TNK_obj@meta.data$CellID,
Classification = TNK_obj@meta.data$Newlabels_rough)
T_NK_cell_labels_detailed <- data.frame(ToSelect = TNK_obj@meta.data$CellID,
Classification = TNK_obj@meta.data$Newlabels_detailed)
saveRDS(T_NK_cell_labels_rough, file = paste(out_dir_TNK, "T_NK_cell_labels_rough.rds", sep = "/"))
saveRDS(T_NK_cell_labels_detailed, file = paste(out_dir_TNK, "T_NK_cell_labels_detailed.rds", sep = "/"))
##########################
## Export Table results ##
##########################
# Counts
TNK_counts <- TNK_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_counts.csv.gz", sep = "/"), "w")
write.csv(TNK_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
TNK_md <- TNK_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_metadata.csv.gz", sep = "/"), "w")
write.csv(x = TNK_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
TNK_umap <- TNK_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_umap.csv.gz", sep = "/"), "w")
write.csv(x = TNK_umap, gz_out_umap)
close(gz_out_umap)
# Save final object (after new labels integration)
saveRDS(TNK_obj, file = paste(out_dir_TNK, "T_NK_cells_DBR.rds", sep = "/"))
2_Analysis_APC.R 0 → 100644
View file @ 918c0314
library(SeuratObject)
library(Seurat)
library(harmony)
library(DoubletFinder)
library(grid)
library(cowplot)
library(gridExtra)
library(RColorBrewer)
library(patchwork)
library(scales)
# Working dir
wdir <- "squadrito_livertumor2022_scrnaseq"
# Data dir
data_dir <- paste(wdir, "data", sep = "/")
# Results dir
out_dir <- paste(wdir, "results", sep = "/")
dir.create(path = out_dir, showWarnings = F)
# Counts
gz_in_counts <- gzfile(paste(data_dir, "Full_final_DBR_labeled_counts.csv.gz", sep = "/"), "rt")
full_counts <- read.csv(gz_in_counts, row.names = 1)
close(gz_in_counts)
# Meta Data
gz_in_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_metadata.csv.gz", sep = "/"), "rt")
full_md <- read.csv(gz_in_md, row.names = 1)
close(gz_in_md)
# UMAP Coordinates
gz_in_umap <- gzfile(paste(data_dir, "Full_final_DBR_labeled_umap.csv.gz", sep = "/"), "rt")
full_umap <- read.csv(gz_in_umap, row.names = 1)
close(gz_in_umap)
Full_obj <- CreateSeuratObject(counts = full_counts, meta.data = full_md)
#######################################
## Macrophages and DCs (termed APCs) ##
#######################################
Full_obj <- SetIdent(object = Full_obj, value = "Newlabels")
APC_obj <- subset(Full_obj, idents = c("APCs"))
out_dir_APC <- paste(out_dir, "APCs", sep = "/")
dir.create(path = out_dir_APC, showWarnings = F)
saveRDS(APC_obj, file = paste(out_dir_APC, "APCs_DBR.rds", sep = "/"))
# Processing of the subset Macrophages, (Normalization, SCTransform, PCA and UMAP)
# We are not computing diff expression here so scale data is not done, we will do diff expression on the full obj
APC_obj@meta.data$CellID <- APC_obj@assays$RNA@counts@Dimnames[[2]] #CellID -> Will be used latter to be sure IDs are matching with full object
DefaultAssay(object = APC_obj) <- "RNA"
APC_obj <- NormalizeData(object = APC_obj)
APC_obj <- SCTransform(APC_obj, vars.to.regress = c("nFeature_RNA", "nCount_RNA", "percent.mt"), verbose = TRUE)
DefaultAssay(object = APC_obj) <- "SCT"
APC_obj@misc$analysis_params$max_pca <- 50
APC_obj <- RunPCA(object = APC_obj, npcs = APC_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_APC, "APCs_DBR_Elbowplot.png", sep = "/"))
ElbowPlot(APC_obj, ndims = APC_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
APC_obj@misc$analysis_params$num_pc <- 35
APC_obj <- RunUMAP(APC_obj, dims = 1:APC_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
APC_obj <- FindNeighbors(APC_obj, reduction = "pca", dims = 1:APC_obj@misc$analysis_params$num_pc, force.recalc = T)
APC_obj <- FindClusters(APC_obj, resolution = 1.2)
png(filename = paste(out_dir_APC, "APCs_DBR_UMAP_res1.2.png", sep = "/"), width = 1200, height = 1000)
DimPlot(object = APC_obj, reduction = "umap", group.by = "SCT_snn_res.1.2", pt.size = 2, label = T, label.size = 10)
dev.off()
#Annotate new labels rough on macrophages
APC_obj@meta.data$Newlabels_rough <- rep("Undefined",length(APC_obj@meta.data$orig.ident))
APC_obj@meta.data$CellID <- APC_obj@assays$RNA@counts@Dimnames[[2]]
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "0"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "1"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "2"] <- "Mo DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "3"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "4"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "5"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "6"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "7"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "8"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "9"] <- "Mo DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "10"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "11"] <- "Monocyte & Macrophages"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "12"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "13"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "14"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "15"] <- "Undefined"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "16"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "17"] <- "Undefined"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "18"] <- "Undefined"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "19"] <- "DCs"
APC_obj@meta.data$Newlabels_rough[APC_obj@meta.data$SCT_snn_res.1.2 == "20"] <- "Undefined"
png(filename = paste(out_dir_APC, "APCs_DBR_UMAP_res1.2_Labels.png", sep = "/"), width = 1800, height = 1000)
DimPlot(object = APC_obj, reduction = "umap", group.by = "SCT_snn_res.1.2", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = APC_obj, reduction = "umap", group.by = "Newlabels_rough", pt.size = 1)
dev.off()
#Annotate new labels detailed on Macrophages
APC_obj@meta.data$Newlabels_detailed <- rep("Undefined",length(APC_obj@meta.data$orig.ident))
APC_obj@meta.data$CellID <- APC_obj@assays$RNA@counts@Dimnames[[2]]
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "0"] <- "IFNa TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "1"] <- "TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "2"] <- "Mo DCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "3"] <- "IFNa TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "4"] <- "IFNa TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "5"] <- "Monocytes"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "6"] <- "IFNa TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "7"] <- "IFNa TAMs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "8"] <- "CD8 cDC1"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "9"] <- "Mo DCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "10"] <- "pDCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "11"] <- "KCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "12"] <- "CDPs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "13"] <- "pre DCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "14"] <- "Ccr7 DCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "15"] <- "Undefined"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "16"] <- "pre DCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "17"] <- "Undefined"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "18"] <- "Undefined"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "19"] <- "pDCs"
APC_obj@meta.data$Newlabels_detailed[APC_obj@meta.data$SCT_snn_res.1.2 == "20"] <- "Undefined"
png(filename = paste(out_dir_APC, "APCs_DBR_UMAP_res1.2_Labels_Detailed.png", sep = "/"), width = 1800, height = 1000)
DimPlot(object = APC_obj, reduction = "umap", group.by = "SCT_snn_res.1.2", pt.size = 1, label = T, label.size = 10) |
DimPlot(object = APC_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1, label = T, label.size = 8)
dev.off()
###Create files to reintegrate T_NK_cell labels into full object
#extract labels from T_NK_cell label object
labels_APC_rough <- data.frame(ToSelect = APC_obj@meta.data$CellID,
Classification = APC_obj@meta.data$Newlabels_rough)
labels_APC_detailed <- data.frame(ToSelect = APC_obj@meta.data$CellID,
Classification = APC_obj@meta.data$Newlabels_detailed)
saveRDS(labels_APC_rough, file = paste(out_dir_APC, "labels_APC_rough.rds", sep = "/"))
saveRDS(labels_APC_detailed, file = paste(out_dir_APC, "labels_APC_detailed.rds", sep = "/"))
##########################
## Export Table results ##
##########################
# Counts
APC_counts <- APC_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_counts.csv.gz", sep = "/"), "w")
write.csv(APC_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
APC_md <- APC_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_metadata.csv.gz", sep = "/"), "w")
write.csv(x = APC_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
APC_umap <- APC_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_umap.csv.gz", sep = "/"), "w")
write.csv(x = APC_umap, gz_out_umap)
close(gz_out_umap)
# Save final object
saveRDS(APC_obj, file = paste(out_dir_APC, "APCs_DBR.rds", sep = "/"))
3_Analysis_UndRem.R 0 → 100644
View file @ 918c0314
library(ggplot2)
library(Seurat)
library(openxlsx)
########################
### General Settings ###
########################
# Working dir
wdir <- "~/Downloads/TIGET/Squadrito/scSquadrito/squadrito_livertumor2022_scrnaseq"
# Data dir
data_dir <- paste(wdir, "data", sep = "/")
# Results dir
out_dir <- paste(wdir, "results", sep = "/")
dir.create(path = out_dir, showWarnings = F)
# Full dir
full_dir <- paste(out_dir, "Full", sep = "/")
dir.create(path = full_dir, showWarnings = F)
# Plot dir
plot_dir <- paste(full_dir, "plots", sep = "/")
dir.create(path = plot_dir, showWarnings = F)
# Full Counts
gz_in_counts <- gzfile(paste(data_dir, "Full_final_DBR_labeled_counts.csv.gz", sep = "/"), "rt")
full_counts <- read.csv(gz_in_counts, row.names = 1)
close(gz_in_counts)
# Full Meta Data
gz_in_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_metadata.csv.gz", sep = "/"), "rt")
full_md <- read.csv(gz_in_md, row.names = 1)
close(gz_in_md)
###################################
## Reintegrate APCS and T-NK-NKT ##
## cells labels into full object ##
###################################
# TNK Meta Data
gz_in_md <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_metadata.csv.gz", sep = "/"), "rt")
TNK_md <- read.csv(gz_in_md, row.names = 1)
close(gz_in_md)
T_NK_cell_labels_rough <- data.frame(ToSelect = TNK_md$CellID, Classification = TNK_md$Newlabels_rough)
T_NK_cell_labels_detailed <- data.frame(ToSelect = TNK_md$CellID, Classification = TNK_md$Newlabels_detailed)
# APC Meta Data
gz_in_md <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_metadata.csv.gz", sep = "/"), "rt")
APC_md <- read.csv(gz_in_md, row.names = 1)
close(gz_in_md)
labels_APC_rough <- data.frame(ToSelect = APC_md$CellID, Classification = APC_md$Newlabels_rough)
labels_APC_detailed <- data.frame(ToSelect = APC_md$CellID, Classification = APC_md$Newlabels_detailed)
# Reintegrate T_NK_cell_rough labels into full object
alllabels <- T_NK_cell_labels_rough
full_md$Newlabels_rough <- full_md$Newlabels
orderDF <- data.frame(allIDs = full_md$CellID, order = c(1:length(full_md$CellID)))
head(alllabels)
head(orderDF)
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
full_md$Newlabels_rough[full_md$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Reintegrate T_NK_cell_detailed labels into full object
alllabels <- T_NK_cell_labels_detailed
full_md$Newlabels_detailed <- full_md$Newlabels
orderDF <- data.frame(allIDs = full_md$CellID, order = c(1:length(full_md$CellID)))
head(alllabels)
head(orderDF)
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
full_md$Newlabels_detailed[full_md$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Reintegrate rough APC cell labels into full object
alllabels <- labels_APC_rough
orderDF <- data.frame(allIDs = full_md$CellID, order = c(1:length(full_md$CellID)))
head(alllabels)
head(orderDF)
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
full_md$Newlabels_rough[full_md$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Reintegrate detailed APC cell labels into full object
alllabels <- labels_APC_detailed
orderDF <- data.frame(allIDs = full_md$CellID, order = c(1:length(full_md$CellID)))
head(alllabels)
head(orderDF)
alllabels <- merge(alllabels, orderDF, by = 1)
alllabels <- alllabels[order(alllabels$order),]
full_md$Newlabels_detailed[full_md$CellID%in%alllabels$ToSelect] <- alllabels$Classification
# Save Full Meta data
gz_out_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_metadata.csv.gz", sep = "/"), "w")
write.csv(x = full_md, gz_out_md)
close(gz_out_md)
############################
## Remove Undefined Cells ##
############################
Full_obj <- CreateSeuratObject(counts = full_counts, meta.data = full_md)
Full_obj <- SetIdent(object = Full_obj, value = "Newlabels_detailed")
Full_obj <- subset(Full_obj, idents = "Undefined", invert = T)
# Reperform clustering after removal of undefined cells
Full_obj@meta.data$CellID <- Full_obj@assays$RNA@counts@Dimnames[[2]] #CellID -> Will be used latter to be sure IDs are matching with full object
DefaultAssay(object = Full_obj) <- "RNA"
Full_obj <- NormalizeData(object = Full_obj)
Full_obj <- SCTransform(Full_obj, vars.to.regress = c("CC.Difference", "percent.mt", "nCount_RNA"), verbose = TRUE)
DefaultAssay(object = Full_obj) <- "SCT"
Full_obj@misc$analysis_params$max_pca <- 50
Full_obj <- RunPCA(object = Full_obj, npcs = Full_obj@misc$analysis_params$max_pca)
png(filename = paste(plot_dir, "Full_und_rem_Elbowplot.png", sep = "/"))
ElbowPlot(Full_obj, ndims = Full_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
Full_obj@misc$analysis_params$num_pc <- 35
Full_obj <- RunUMAP(Full_obj, dims = 1:Full_obj@misc$analysis_params$num_pc, seed.use = 20, reduction = "pca")
Full_obj <- FindNeighbors(Full_obj, reduction = "pca", dims = 1:Full_obj@misc$analysis_params$num_pc, force.recalc = T)
### Table Markers Fullobj ###
DefaultAssay(object = Full_obj) <- "RNA"
Full_obj <- SetIdent(Full_obj, value = "Newlabels")
tableFullobj <- FindAllMarkers(object = Full_obj, logfc.threshold = 0.25, min.pct = 0.1, test.use = "wilcox")
write.xlsx(x = list("Allcells_Newlabels_AllMarkersF" = tableFullobj), file = paste(data_dir, "scRNAseq_DiFgenes_Fullobj_FC0.25.xlsx", sep = "/"))
# Save full object
saveRDS(object = Full_obj, paste(full_dir, "Full_final_DBR_labeled_Und_rem.rds", sep = "/"))
# Recluster T and NK cells
Full_obj <- SetIdent(object = Full_obj, value = "Newlabels")
TNK_obj <- subset(Full_obj, idents = "T and NK cells")
out_dir_TNK <- paste(out_dir, "T_NK_cells", sep = "/")
dir.create(out_dir_TNK, showWarnings = F)
TNK_obj@meta.data$CellID <- TNK_obj@assays$RNA@counts@Dimnames[[2]] #CellID -> Will be used latter to be sure IDs are matching with full object
DefaultAssay(object = TNK_obj) <- "RNA"
TNK_obj <- NormalizeData(object = TNK_obj)
TNK_obj <- SCTransform(TNK_obj, vars.to.regress = c("CC.Difference", "percent.mt", "nCount_RNA"), verbose = TRUE)
DefaultAssay(object = TNK_obj) <- "SCT"
TNK_obj@misc$analysis_params$max_pca <- 50
TNK_obj <- RunPCA(object = TNK_obj, npcs = TNK_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_TNK, "T_and_NK_Und_rem_Elbowplot.png", sep = "/"))
ElbowPlot(TNK_obj, ndims = TNK_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
TNK_obj@misc$analysis_params$num_pc <- 35
TNK_obj <- RunUMAP(TNK_obj, dims = 1:TNK_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
TNK_obj <- FindNeighbors(TNK_obj, reduction = "pca", dims = 1:TNK_obj@misc$analysis_params$num_pc, force.recalc = T)
### Table Markers T&NKcells ###
DefaultAssay(object = TNK_obj) <- "RNA"
TNK_obj <- SetIdent(TNK_obj, value = "Newlabels_detailed")
tableTNKcells <- FindAllMarkers(TNK_obj, logfc.threshold = 0.25, min.pct = 0.1, test.use = "wilcox")
write.xlsx(x = list("T_NK_cells_Newlabels_detailed_A" = tableTNKcells), file = paste(data_dir, "scRNAseq_DIFgenes_T&NKcells_FC0.25.xlsx", sep = "/"))
# Save T_NK object
saveRDS(object = TNK_obj, paste(out_dir_TNK, "T_NK_final_DBR_labeled_Und_rem.rds", sep = "/"))
# Recluster APCs
Full_obj <- SetIdent(object = Full_obj, value = "Newlabels")
APC_obj <- subset(Full_obj, idents = "APCs")
out_dir_APC <- paste(out_dir,"APCs", sep = "/")
dir.create(out_dir_APC, showWarnings = F)
APC_obj@meta.data$CellID <- APC_obj@assays$RNA@counts@Dimnames[[2]] #CellID -> Will be used latter to be sure IDs are matching with full object
DefaultAssay(object = APC_obj) <- "RNA"
APC_obj <- NormalizeData(object = APC_obj)
APC_obj <- SCTransform(APC_obj, vars.to.regress = c("CC.Difference", "percent.mt", "nCount_RNA"), verbose = TRUE)
DefaultAssay(object = APC_obj) <- "SCT"
APC_obj@misc$analysis_params$max_pca <- 50
APC_obj <- RunPCA(object = APC_obj, npcs = APC_obj@misc$analysis_params$max_pca)
png(filename = paste(out_dir_APC, "APCs_DBR_Elbowplot.png", sep = "/"))
ElbowPlot(APC_obj, ndims = APC_obj@misc$analysis_params$max_pca, reduction = "pca")
dev.off()
APC_obj@misc$analysis_params$num_pc <- 35
APC_obj <- RunUMAP(APC_obj, dims = 1:APC_obj@misc$analysis_params$num_pc, seed.use = 1, reduction = "pca")
APC_obj <- FindNeighbors(APC_obj, reduction = "pca", dims = 1:APC_obj@misc$analysis_params$num_pc, force.recalc = T)
### Table Markers APCs ###
DefaultAssay(object = APC_obj) <- "RNA"
APC_obj <- SetIdent(APC_obj, value = "Newlabels_detailed")
tableAPCs <- FindAllMarkers(object = APC_obj, logfc.threshold = 0.25, min.pct = 0.1, test.use = "wilcox")
write.xlsx(x = list("APCs_Newlabels_detailed_AllMark" = tableAPCs), file = paste(data_dir, "scRNAseq_DIFgenes_APCs_FC0.25.xlsx", sep = "/"))
# Save APC
saveRDS(object = APC_obj, paste(out_dir_APC, "APC_final_DBR_labeled_Und_rem.rds", sep = "/"))
##########################
## Export Table results ##
##########################
data_dir <- paste(wdir, "data", sep = "/")
dir.create(path = data_dir, showWarnings = F)
## Full ##
# Counts
full_counts <- Full_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "Full_final_DBR_labeled_Und_rem_counts.csv.gz", sep = "/"), "w")
write.csv(full_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
full_md <- Full_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "Full_final_DBR_labeled_Und_rem_metadata.csv.gz", sep = "/"), "w")
write.csv(x = full_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
full_umap <- Full_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "Full_final_DBR_labeled_Und_rem_umap.csv.gz", sep = "/"), "w")
write.csv(x = full_umap, gz_out_umap)
close(gz_out_umap)
## TNK ##
# Counts
TNK_counts <- TNK_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_Und_rem_counts.csv.gz", sep = "/"), "w")
write.csv(TNK_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
TNK_md <- TNK_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_Und_rem_metadata.csv.gz", sep = "/"), "w")
write.csv(x = TNK_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
TNK_umap <- TNK_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "T_NK_cells_DBR_labeled_Und_rem_umap.csv.gz", sep = "/"), "w")
write.csv(x = TNK_umap, gz_out_umap)
close(gz_out_umap)
## APCs ##
# Counts
APC_counts <- APC_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_Und_rem_counts.csv.gz", sep = "/"), "w")
write.csv(APC_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
APC_md <- APC_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_Und_rem_metadata.csv.gz", sep = "/"), "w")
write.csv(x = APC_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
APC_umap <- APC_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "APCs_cells_DBR_labeled_Und_rem_umap.csv.gz", sep = "/"), "w")
write.csv(x = APC_umap, gz_out_umap)
close(gz_out_umap)
data/.DS_Store 0 → 100644
View file @ 918c0314
File added
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