#' @description: proProcess, cluster and remove the batch effect

# load package

library(Seurat)

library(harmony)

library(clustree)

library(ggpubr)

library(dplyr)

library(patchwork)

library(ComplexHeatmap)

library(circlize)

library(vegan)

set.seed(101)

library(future)

plan("multiprocess", workers = 10) 

options(future.globals.maxSize = 100000 * 1024^2) # set 50G RAM

setwd(dir = "/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA")

source(file = "/home/longzhilin/Analysis_Code/Visualization/colorPalettes.R")

#### four samples

T1.data <- Read10X(data.dir = "/data/raw_data/lzl/RenalTumor-20200713/cellRanger_result/scRNA-V5.0/T1/outs/filtered_feature_bc_matrix")

colnames(T1.data) <- paste0("T1_", colnames(T1.data))

#36601*8849

T2.data <- Read10X(data.dir = "/data/raw_data/lzl/RenalTumor-20200713/cellRanger_result/scRNA-V5.0/T2/outs/filtered_feature_bc_matrix")

colnames(T2.data) <- paste0("T2_", colnames(T2.data))

#36601*15440

T3.data <- Read10X(data.dir = "/data/raw_data/lzl/RenalTumor-20200713/cellRanger_result/scRNA-V5.0/T3/outs/filtered_feature_bc_matrix")

colnames(T3.data) <- paste0("T3_", colnames(T3.data))

#36601*14644

T4.data <- Read10X(data.dir = "/data/raw_data/lzl/RenalTumor-20200713/cellRanger_result/scRNA-V5.0/T4/outs/filtered_feature_bc_matrix")

colnames(T4.data) <- paste0("T4_", colnames(T4.data))

#36601*11016

#### Preliminary filtration

#min.cell >3 & min.features >200                        

T1 <- CreateSeuratObject(counts = T1.data,

                         project = "T1",

                         min.cells = 3,

                         min.features = 200)

#23137*8823                         

T2 <- CreateSeuratObject(counts = T2.data,

                         project = "T2",

                         min.cells = 3,

                         min.features = 200) 

#23487*15437                         

T3 <- CreateSeuratObject(counts = T3.data,

                         project = "T3",

                         min.cells = 3,

                         min.features = 200)

#24106*14602                         

T4 <- CreateSeuratObject(counts = T4.data,

                         project = "T4",

                         min.cells = 3,

                         min.features = 200) 

#23065*10792

#### Mitochondrial gene ratio

T1[["percent.mt"]] <- PercentageFeatureSet(T1, pattern = "^MT-")

T2[["percent.mt"]] <- PercentageFeatureSet(T2, pattern = "^MT-")

T3[["percent.mt"]] <- PercentageFeatureSet(T3, pattern = "^MT-")

T4[["percent.mt"]] <- PercentageFeatureSet(T4, pattern = "^MT-")

#### Draw a statistical graph of the number of genes/count number/proportion of mitochondrial genes

pdf(file = "1.QualityControl/count.feature.mt.pdf")

VlnPlot(T1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(T1, feature1 = "nCount_RNA", feature2 = "percent.mt")

ploT1 <- FeatureScatter(T1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1 + ploT1

VlnPlot(T2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(T2, feature1 = "nCount_RNA", feature2 = "percent.mt")

ploT1 <- FeatureScatter(T2, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1 + ploT1

VlnPlot(T3, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(T3, feature1 = "nCount_RNA", feature2 = "percent.mt")

ploT1 <- FeatureScatter(T3, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1 + ploT1

VlnPlot(T4, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(T4, feature1 = "nCount_RNA", feature2 = "percent.mt")

ploT1 <- FeatureScatter(T4, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1 + ploT1

ggdensity(T1@meta.data, x = "nCount_RNA", title = "T1")

ggdensity(T1@meta.data, x = "nFeature_RNA", title = "T1")

ggdensity(T1@meta.data, x = "percent.mt", title = "T1")

ggdensity(T2@meta.data, x = "nCount_RNA", title = "T2")

ggdensity(T2@meta.data, x = "nFeature_RNA", title = "T2")

ggdensity(T2@meta.data, x = "percent.mt", title = "T2")

ggdensity(T3@meta.data, x = "nCount_RNA", title = "T3")

ggdensity(T3@meta.data, x = "nFeature_RNA", title = "T3")

ggdensity(T3@meta.data, x = "percent.mt", title = "T3")

ggdensity(T4@meta.data, x = "nCount_RNA", title = "T4")

ggdensity(T4@meta.data, x = "nFeature_RNA", title = "T4")

ggdensity(T4@meta.data, x = "percent.mt", title = "T4")

dev.off()

######################### Detect the resolution parameters of each sample cluster. After the parameters are determined, you can block them without executing [test]

set.resolutions <- seq(0.5, 2, by = 0.1)

pdf(file = "1.QualityControl/PCA-test.pdf")

#### T1

T1.pro <- subset(T1, subset = nFeature_RNA > 200 & percent.mt < 10 & nCount_RNA > 1000 & nFeature_RNA < 6000) 

T1.pro <- SCTransform(T1.pro, vars.to.regress = c("nCount_RNA", "percent.mt"), verbose = FALSE)

T1.pro <- RunPCA(T1.pro, npcs = 100, verbose = FALSE)

ElbowPlot(object = T1.pro, ndims = 100)

T1.pro  <- FindNeighbors(object = T1.pro , dims = 1:50, verbose = FALSE)

T1.pro  <- FindClusters(object = T1.pro , resolution = set.resolutions, verbose = FALSE) 

clustree(T1.pro)

T1.pro  <- RunUMAP(T1.pro , dims = 1:50)

T1.res <- sapply(set.resolutions, function(x){

    p <- DimPlot(object = T1.pro, reduction = 'umap',label = TRUE, group.by = paste0("SCT_snn_res.", x))

    print(p)

})

#### T2

T2.pro <- subset(T2, subset = nFeature_RNA > 200 & percent.mt < 10 & nCount_RNA > 1000 & nFeature_RNA < 6000)

T2.pro <- SCTransform(T2.pro, vars.to.regress = c("nCount_RNA", "percent.mt"), verbose = FALSE)

T2.pro <- RunPCA(T2.pro, npcs = 100, verbose = FALSE)

ElbowPlot(object = T2.pro, ndims = 100)

T2.pro  <- FindNeighbors(object = T2.pro , dims = 1:50, verbose = FALSE)

T2.pro  <- FindClusters(object = T2.pro , resolution = set.resolutions, verbose = FALSE) 

clustree(T2.pro)

T2.pro  <- RunUMAP(T2.pro , dims = 1:50)

T2.res <- sapply(set.resolutions, function(x){

    p <- DimPlot(object = T2.pro, reduction = 'umap',label = TRUE, group.by = paste0("SCT_snn_res.", x))

    print(p)

})

#### T3

T3.pro <- subset(T3, subset = nFeature_RNA > 200 & percent.mt < 10 & nCount_RNA > 1000 & nFeature_RNA < 6000) 

T3.pro <- SCTransform(T3.pro, vars.to.regress = c("nCount_RNA", "percent.mt"), verbose = FALSE)

T3.pro <- RunPCA(T3.pro, npcs = 100, verbose = FALSE)

ElbowPlot(object = T3.pro, ndims = 100)

T3.pro  <- FindNeighbors(object = T3.pro , dims = 1:50, verbose = FALSE)

T3.pro  <- FindClusters(object = T3.pro , resolution = set.resolutions, verbose = FALSE) 

clustree(T3.pro)

T3.pro  <- RunUMAP(T3.pro , dims = 1:50)

T3.res <- sapply(set.resolutions, function(x){

    p <- DimPlot(object = T3.pro, reduction = 'umap',label = TRUE, group.by = paste0("SCT_snn_res.", x))

    print(p)

})

#### T4

T4.pro <- subset(T4, subset = nFeature_RNA > 200 & percent.mt < 10 & nCount_RNA > 1000 & nFeature_RNA < 6000) 

T4.pro <- SCTransform(T4.pro, vars.to.regress = c("nCount_RNA", "percent.mt"), verbose = FALSE)

T4.pro <- RunPCA(T4.pro, npcs = 100, verbose = FALSE)

ElbowPlot(object = T4.pro, ndims = 100)

T4.pro  <- FindNeighbors(object = T4.pro , dims = 1:50, verbose = FALSE)

T4.pro  <- FindClusters(object = T4.pro , resolution = set.resolutions, verbose = FALSE) 

clustree(T4.pro)

T4.pro  <- RunUMAP(T4.pro , dims = 1:50)

T4.res <- sapply(set.resolutions, function(x){

    p <- DimPlot(object = T4.pro, reduction = 'umap',label = TRUE, group.by = paste0("SCT_snn_res.", x))

    print(p)

})

dev.off()

#### remove doublet

library(DoubletFinder) # Require cleanup of low-quality cells in advance

source(file = "/home/longzhilin/Analysis_Code/SingleCell/doubletDetect.R")

pdf("1.QualityControl/doublet.pdf")

T1.pro1 <- doubletDetect(Seurat.object = T1.pro, PCs = 1:50, doublet.rate = 0.061, annotation = "SCT_snn_res.0.7", sct = T) #7893 ~8000

T2.pro1 <- doubletDetect(Seurat.object = T2.pro, PCs = 1:50, doublet.rate = 0.106, annotation = "SCT_snn_res.1.7", sct = T) #13992 ~14000

T3.pro1 <- doubletDetect(Seurat.object = T3.pro, PCs = 1:50, doublet.rate = 0.091, annotation = "SCT_snn_res.0.7", sct = T) #11973 ~12000

T4.pro1 <- doubletDetect(Seurat.object = T4.pro, PCs = 1:50, doublet.rate = 0.061, annotation = "SCT_snn_res.0.5", sct = T) #8054 ~8000

dev.off()

saveRDS(T1.pro1, file = "T1.pro1.rds")

saveRDS(T2.pro1, file = "T2.pro1.rds")

saveRDS(T3.pro1, file = "T3.pro1.rds")

saveRDS(T4.pro1, file = "T4.pro1.rds")

pdf("1.QualityControl/doublet.cell.pdf")

DimPlot(object = T1.pro1, reduction = 'umap', group.by = "Doublet")

DimPlot(object = T2.pro1, reduction = 'umap', group.by = "Doublet")

DimPlot(object = T3.pro1, reduction = 'umap', group.by = "Doublet")

DimPlot(object = T4.pro1, reduction = 'umap', group.by = "Doublet")

dev.off()

T1.pro2 <- subset(T1.pro1, subset = Doublet == "Singlet") #21247*7449

T2.pro2 <- subset(T2.pro1, subset = Doublet == "Singlet") #21605*12574

T3.pro2 <- subset(T3.pro1, subset = Doublet == "Singlet") #21947*10937

T4.pro2 <- subset(T4.pro1, subset = Doublet == "Singlet") #20459*7605

saveRDS(T1.pro2, file = "T1.pro2.rds")

saveRDS(T2.pro2, file = "T2.pro2.rds")

saveRDS(T3.pro2, file = "T3.pro2.rds")

saveRDS(T4.pro2, file = "T4.pro2.rds")

############################################## merge data and correct the batch effect

DefaultAssay(T1.pro2) <- "RNA"

DefaultAssay(T2.pro2) <- "RNA"

DefaultAssay(T3.pro2) <- "RNA"

DefaultAssay(T4.pro2) <- "RNA"

source(file = "/home/longzhilin/Analysis_Code/SingleCell/variableFeatureSelection.R")

Renal.list <- list(T1 = T1.pro2, T2 = T2.pro2, T3 = T3.pro2, T4 = T4.pro2)

Renal.list.Stardard <- variableFeatureSelection(seurat.lists = Renal.list, method = "Stardard", nfeatures = 3000)

saveRDS(Renal.list.Stardard, file = "Renal.list.Stardard.3000.rds")

Renal.list.SCT <- variableFeatureSelection(seurat.lists = Renal.list, method = "SCT", nfeatures = 3000)

saveRDS(Renal.list.SCT, file = "Renal.list.SCT.3000.rds")

#### 

# assay=SCT

data.merge <- merge(Renal.list.SCT[[1]], y = Renal.list.SCT[2:length(Renal.list.SCT)], project = "Renal")

DefaultAssay(data.merge) <- "SCT"

seurat.features.SCT <- SelectIntegrationFeatures(object.list = Renal.list.SCT, nfeatures = 3000)

VariableFeatures(data.merge) <- seurat.features.SCT

# Remove previous clustering results

index <- match(paste0("SCT_snn_res.", seq(0.5, 2, by=0.1)), colnames(data.merge@meta.data))

data.merge@meta.data <- data.merge@meta.data[,-index]

# assay=RNA

seurat.features.RNA <- SelectIntegrationFeatures(object.list = Renal.list.Stardard, nfeatures = 3000)

DefaultAssay(data.merge) <- "RNA"

VariableFeatures(data.merge) <- seurat.features.RNA

data.merge <- NormalizeData(data.merge, verbose = FALSE)

data.merge <- ScaleData(data.merge, verbose = FALSE, vars.to.regress = c("nCount_RNA", "percent.mt"), features = rownames(data.merge@assays$RNA@data))

DefaultAssay(data.merge) <- "SCT"

saveRDS(data.merge, file = "data.merge.rds")

##############################################2.Evaluation of cellcycle and patient bias

DefaultAssay(data.merge) <- "SCT"

pdf("1.QualityControl/filtered.statistics.pdf")

VlnPlot(object = data.merge, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "orig.ident", cols = Palettes$group_pal[1:length(unique(data.merge@meta.data$orig.ident))], pt.size = 0)

FeatureScatter(object = data.merge, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "orig.ident")

FeatureScatter(object = data.merge, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "orig.ident")

dev.off()

# Draw the distribution of the number of samples

cell.number <- as.data.frame(table(data.merge$orig.ident))

pdf("1.QualityControl/highQuality.pdf")

ggbarplot(cell.number, x="Var1", y="Freq", fill = "Var1", color = "Var1", palette = Palettes$group_pal[1:length(unique(data.merge@meta.data$orig.ident))],

          sort.by.groups=FALSE, #不按组排序

          label = T, xlab = "", ylab = "Cell Number") + theme(legend.position="none") 

dev.off()

#### Assess the cell cycle effect

s.genes <- cc.genes$s.genes

g2m.genes <- cc.genes$g2m.genes

data.merge <- CellCycleScoring(data.merge, s.features = s.genes, g2m.features = g2m.genes, set.ident = F)

data.merge <- RunPCA(data.merge, features = c(s.genes, g2m.genes))

pdf("1.QualityControl/cellCycle.afterMerge.pdf")

DimPlot(data.merge, dims = c(1, 2), reduction = "pca", group.by = "Phase")

DimPlot(data.merge, dims = c(1, 3), reduction = "pca", group.by = "Phase")

DimPlot(data.merge, dims = c(2, 3), reduction = "pca", group.by = "Phase")

dev.off()

set.resolutions <- seq(0.2, 1.2, by = 0.1)

#### First observe whether the clustering effect will depend on the sample

a <- data.merge

a <- RunPCA(a, npcs = 100, verbose = T)

pdf("1.QualityControl/merge.observe.batch.pdf")

ElbowPlot(object = a, ndims = 100)

a <- FindNeighbors(a, dims = 1:50, verbose = T)

a <- FindClusters(object = a, resolution = set.resolutions, verbose = T) 

clustree(a)

a <- RunUMAP(a, dims = 1:50)

DimPlot(object = a, reduction = 'umap',label = TRUE, group.by = "orig.ident")

DimPlot(object = a, reduction = 'umap',label = TRUE, group.by = "Phase")

merge.res <- sapply(set.resolutions, function(x){

    p <- DimPlot(object = a, reduction = 'umap',label = TRUE, group.by = paste0("SCT_snn_res.", x))

    print(p)

})

dev.off()

##############################################3.correct batch effect

source(file = "/home/longzhilin/Analysis_Code/SingleCell/scRNA.Integrate.multipleSample.R")

pdf("2.Cluster/SCT.Harmony.Integration.PC40.feature3000.pdf")

data.merge.harmony.PC40.SCT <- Harmony.integration.reduceDimension(seurat.object = data.merge, assay = "SCT", set.resolutions = seq(0.2, 1.2, by = 0.1), PC = 40, nfeatures = 3000, npcs = 50)

dev.off()

saveRDS(data.merge.harmony.PC40.SCT, file = "data.merge.harmony.PC40.SCT.feature3000.rds")