呜呜最近发现我工作效率低的一个原因就是重复性工作没有流程化，一气之下，把seurat分析单套数据的流程封装了起来，步骤包含数据质控、数据标准化、聚类以及初步的细胞类型鉴定。细胞类型鉴定是用每个cluster的top marker来标注的。之后再更新整合多套数据的流程，希望与君分享

1. 用到的所有函数放在了SeuratWrapperFunction.R中了

这个需要用source()函数导入到下面封装好的代码中的

### Time: 20221025 ### Author: zhengyiyi ## load function library(Seurat) library(SingleCellExperiment) library(scater) library(MAST) library(ggplot2) library(patchwork) library(RColorBrewer) library(Seurat) library(openxlsx) library("grid") library("ggplotify") library(magrittr) library(dplyr) library(ggsignif) library("ggplot2") library(cowplot) library(openxlsx) CreateFirstFilterSeuratObject <- function(exp_count, pr_name, Or_ident, min_cell, min_gene, mt_pattern,rb_pattern){ # CreateSeuratObject : filter gene and cell by min.cells(Include features detected in at least this many cells) and min.features(Include cells where at least this many features are detected.) exp_count.seurat <- CreateSeuratObject(exp_count, project = pr_name, min.cells = min_cell, min.features = min_gene) exp_count.seurat@meta.data$orig.ident <- as.factor(Or_ident) exp_count.seurat@active.ident <- exp_count.seurat$orig.ident ## Check the mt gene and rp gene exp_count.seurat[["percent.mt"]] <- PercentageFeatureSet(object = exp_count.seurat, pattern = mt_pattern) exp_count.seurat[["percent.rb"]] <- PercentageFeatureSet(object = exp_count.seurat, pattern = rb_pattern) print(levels(exp_count.seurat@active.ident)) print(table(exp_count.seurat$orig.ident)) print(head(x = exp_count.seurat@meta.data, 5)) print(dim(exp_count.seurat@meta.data)) return(exp_count.seurat) } SeuratQC <- function(exp_count.seurat, name_pre){ folder_name="1.QC" if (file.exists(folder_name)){ print("1.QC existed.") } else{ dir.create(folder_name) } file_name=paste(name_pre,"Violinplot.pdf",sep="_") reads.drop <- isOutlier(as.numeric(exp_count.seurat$nCount_RNA), nmads = 3, type = "lower") feature.drop <- isOutlier(as.numeric(exp_count.seurat$nFeature_RNA), nmads = 3, type = "lower") qc.mito2 <- isOutlier(exp_count.seurat$percent.mt, nmads = 3, type="higher") qc.ribo2 <- isOutlier(exp_count.seurat$percent.rb, nmads = 3, type="higher") ### 保留下那些质量差的细胞，看看这些细胞的情况 low_q_cells <- (reads.drop | feature.drop | qc.mito2 | qc.ribo2) exp_count.seurat$lowquality_cells <- "High_qualitycells" exp_count.seurat$lowquality_cells[low_q_cells] <- "Low_qualitycells" plot1=FeatureScatter(object = exp_count.seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "lowquality_cells") plot2=FeatureScatter(object = exp_count.seurat, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "lowquality_cells") plot3=FeatureScatter(object = exp_count.seurat, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "lowquality_cells") pdf(file = file_name, height = 10, width = 20) print(VlnPlot(object = exp_count.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size=0)) opar <- par(no.readonly = TRUE) par(mfrow=c(2,2)) print(hist(exp_count.seurat$nCount_RNA, breaks = 100, main = "Reads Count Distribution", col = "grey80", xlab = "library size", ylab = "frequency", cex.lab = 1.4, cex.axis = 1.4)) print(abline(v = 10000, col = "red", lwd = 2, lty = 2)) print(hist(exp_count.seurat$nFeature_RNA, breaks = 100, main = "Gene Count Distribution", col = "grey80", xlab = "total mRNA", ylab = "frequency", cex.lab = 1.4, cex.axis = 1.4)) print(abline(v = 1500, col = "red", lwd = 2, lty = 2)) print(hist(exp_count.seurat$percent.mt, breaks = 10, main = "Mitochondrial Percentage Distribution", col = "grey80", xlab = "total mRNA", ylab = "frequency", cex.lab = 1.4, cex.axis = 1.4)) print(abline(v = 20, col = "red", lwd = 2, lty = 2)) print(hist(exp_count.seurat$percent.rb, breaks = 10, main = "Ribosome Percentage Distribution", col = "grey80", xlab = "total mRNA", ylab = "frequency", cex.lab = 1.4, cex.axis = 1.4)) print(abline(v = 20, col = "red", lwd = 2, lty = 2)) par(opar) print(plot1 + plot2 + plot3) dev.off() #copy files to 1.QC result_files <- c(file_name) file.copy(result_files, folder_name ,overwrite = T) file.remove(result_files) print("based on the Violin plot, you can filter the cell by the nFeature_RNA, nCount_RNA and/or percent.mt rb") } SubsetSeuratData <- function(exp_count.seurat, min_gene, max_gene, min_ncountRNA, max_ncountRNA, max_mt_percent, max_rb_percent){ exp_count.seurat <- subset(exp_count.seurat, subset = nCount_RNA < max_ncountRNA & nCount_RNA > min_ncountRNA & nFeature_RNA < max_gene & nFeature_RNA > min_gene & percent.mt < max_mt_percent & percent.rb < max_rb_percent ) # QC qc_pdf_name <- paste("filtered_by_gene", min_gene, "UMI", min_ncountRNA, "percent.mt", max_mt_percent, "percent.rb", max_rb_percent, sep="_") SeuratQC(exp_count.seurat, qc_pdf_name) dim(exp_count.seurat) print(levels(exp_count.seurat$orig.ident)) return(exp_count.seurat) } SelectPCsByLogNormalize <- function(exp_count.seurat, file_prefix, scale_factor, nfeatures, npc_plot){ folder_name="2.RunPCA" if (file.exists(folder_name)){ print("2.RunPCA existed.") } else{ dir.create(folder_name) } DefaultAssay(exp_count.seurat) <- "RNA" # Normalize Find VariableGene and Scale Data exp_count.seurat <- NormalizeData(exp_count.seurat, normalization.method = "LogNormalize", scale.factor = scale_factor) exp_count.seurat <- FindVariableFeatures(exp_count.seurat, selection.method = "vst", nfeatures = nfeatures) exp_count.seurat <- ScaleData(exp_count.seurat, vars.to.regress = c("percent.rb","percent.mt")) # RunPCA exp_count.seurat <- RunPCA(exp_count.seurat, pc.genes = exp_count.seurat@var.genes, npcs = npc_plot, verbose = FALSE,) # Select PCs:1:30pcs selected_pcs_name <- paste0(file_prefix, "_Selected_PCs.pdf") pdf(file = selected_pcs_name, height = 8, width = 8) print(ElbowPlot(object = exp_count.seurat, ndims = npc_plot, reduction = "pca")) dev.off() # Plot pca_plot_name <- paste0(file_prefix, "_PCA_plot",".pdf") pdf(pca_plot_name,8,8) print(DimPlot(exp_count.seurat, reduction="pca", label = TRUE, pt.size=0.5, group.by="ident")) dev.off() #Copy files to 2.RunPCA result_files <- c(selected_pcs_name,pca_plot_name) file.copy(result_files, folder_name ,overwrite = T)#拷贝文件 file.remove(result_files) #移除拷贝完的文件 return(exp_count.seurat) } SelectPCsBySCTransform <- function(exp_count.seurat, file_prefix, npc_plot){ folder_name="2.RunPCA" if (file.exists(folder_name)){ print("2.RunPCA existed.") } else{ dir.create(folder_name) } # Normalize Find VariableGene and Scale Data exp_count.seurat <- SCTransform(exp_count.seurat, vars.to.regress = c("percent.rb","percent.mt")) # RunPCA exp_count.seurat <- RunPCA(exp_count.seurat, verbose = FALSE, npcs = npc_plot) # Select PCs:1:30pcs selected_pcs_name <- paste0(file_prefix, "_Selected_PCs.pdf") pdf(file = selected_pcs_name, height = 8, width = 8) print(ElbowPlot(object = exp_count.seurat, ndims = npc_plot, reduction = "pca")) dev.off() # Plot pca_plot_name <- paste0(file_prefix, "_PCA_plot",".pdf") pdf(pca_plot_name, 8, 8) print(DimPlot(exp_count.seurat,reduction="pca", label = TRUE, pt.size= 0.5,group.by="ident",shape.by="orig.ident")) dev.off() #Copy files to 2.RunPCA result_files <- c(selected_pcs_name,pca_plot_name) file.copy(result_files, folder_name ,overwrite = T)#拷贝文件 file.remove(result_files) #移除拷贝完的文件 return(exp_count.seurat) } # RunUMAP and RunTSNE and Then find cluster PlotCluster <- function(exp_count.seurat, file_prefix, npc_used, k_param, resolution_number){ folder_name="3.PlotCluster" if (file.exists(folder_name)){ print("3.PlotCluster existed.") } else{ dir.create(folder_name) } # FinderNeighbors exp_count.seurat <- FindNeighbors(exp_count.seurat, dims = 1:npc_used, verbose = FALSE, k.param = k_param) exp_count.seurat <- FindClusters(exp_count.seurat, verbose = FALSE, resolution = resolution_number) # RunTSNE and RunUMAP exp_count.seurat <- RunTSNE(exp_count.seurat, dims = 1:npc_used, verbose = FALSE, check_duplicates = FALSE) exp_count.seurat <- RunUMAP(exp_count.seurat, dims = 1:npc_used, verbose = FALSE) # RenameIdents from zero to one levels_define <- as.numeric(levels(exp_count.seurat)) new.cluster.ids <- levels_define + 1 names(new.cluster.ids) <- levels(exp_count.seurat) exp_count.seurat <- RenameIdents(exp_count.seurat, new.cluster.ids) combined_cluster_plotname <- paste0(file_prefix, "_combined_cluster_resolution_",resolution_number,".pdf") pdf(combined_cluster_plotname,7,7) print(DimPlot(exp_count.seurat,reduction="umap",label = TRUE, pt.size = 0.5, label.size = 4.5, repel=TRUE, group.by="ident")) print(DimPlot(exp_count.seurat,reduction="tsne", label = TRUE, pt.size = 0.5, label.size = 4.5, repel=TRUE, group.by="ident")) dev.off() split_clustering_plot_name <- paste0(file_prefix, "_split_clustering_resolution_", resolution_number, ".pdf") pdf(split_clustering_plot_name,14,7) print(DimPlot(exp_count.seurat,reduction="umap",label = TRUE, pt.size = 0.5, label.size = 4.5, group.by="ident", repel=TRUE, split.by="orig.ident")) print(DimPlot(exp_count.seurat,reduction="tsne", label = TRUE, pt.size = 0.5, label.size = 4.5, group.by="ident", repel=TRUE, split.by="orig.ident")) dev.off() # Copy files to 2.Cluster file.copy(combined_cluster_plotname, folder_name ,overwrite = T) file.remove(combined_cluster_plotname) file.copy(split_clustering_plot_name, folder_name ,overwrite = T) file.remove(split_clustering_plot_name) return(exp_count.seurat) } # Find markers for each cluster FindmarkerForCluster <- function(exp_count.seurat, file_prefix, min.pct, logfc.threshold, p_val_adj, mt_rb_pattern){ folder_name="4.MarkersInCluster" if (file.exists(folder_name)){ print("4.MarkersInCluster file existed.") } else{ dir.create(folder_name) } cluster.markers <- FindAllMarkers(object = exp_count.seurat, only.pos = TRUE, min.pct = min.pct, logfc.threshold = logfc.threshold) index <- cluster.markers$p_val_adj < p_val_adj cluster.markers <- cluster.markers[index,] # remove mt and rb gene index <- grep(mt_rb_pattern, cluster.markers$gene) if (length(index) > 0){ cluster.markers <- cluster.markers[-index, ] } save_name <- paste0(file_prefix,"_MarkersInClusters.csv") write.csv(cluster.markers,save_name) file.copy(save_name, folder_name,overwrite = T) file.remove(save_name) return(cluster.markers) } # Top markers for each cluster TopMarkersInCluster <- function(cluster.markers, file_prefix, top_num){ library(dplyr) folder_name="4.MarkersInCluster" if (file.exists(folder_name)){ print("4.MarkersInCluster file existed.") } else{ dir.create(folder_name) } #将readsCountSM.markers传给group_by，group_by按cluster 排序，再将结果传给top_n，top_n按avg_logFC排序，显示每个类中的前两个 top_marker <- cluster.markers %>% group_by(cluster) %>% top_n(n = top_num, wt = avg_log2FC) file_name=paste(file_prefix, "_top_marker", top_num,".csv",sep="") write.csv(top_marker, file =file_name) file.copy(file_name, folder_name,overwrite = T) return(top_marker) file.remove(file_name) } # Rename each cluster with top2 markers MapTop2MarkerEachCluster <- function(exp_count.seurat, cluster.markers, file_prefix){ library(dplyr) library(plyr) exp_count.seurat$seurat_clusters <- exp_count.seurat@active.ident top2 <- cluster.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC) top2 <- top2 %>% group_by(cluster) %>% dplyr::mutate(markers = paste0(gene, collapse = "/")) %>% dplyr::slice(1) marker.names <- top2$markers current.cluster.ids <- as.character(1:(length(unique(Idents(exp_count.seurat))))) new.cluster.ids <- marker.names Idents(exp_count.seurat) <- plyr::mapvalues(Idents(exp_count.seurat), from = current.cluster.ids, to = new.cluster.ids) return(exp_count.seurat) }

2.DataQualityOverviewSingle.R: 数据质控的函数是放在了这个脚本

使用的话需要把source()那行代码改成自己的路径
这边有个特别注意的点就是这个脚本里的第一行
“#!/home/zhengjh/miniconda3/envs/r403/bin/Rscript”
之前我是把脚本放环境变量里了，每次source环境的时候总是就不能直接使用，必须全路径，只要把这个Rscript所在位置换成你用的conda环境中的Rscript就可以了。另外需要chmod +x DataQualityOverviewSingle.R 这个代码，才可以执行

#!/home/zhengjh/miniconda3/envs/r403/bin/Rscript # parameter library(optparse) library(getopt) option_list <- list( make_option(c("-o", "--output"), type = "character", default = FALSE, action = "store", help = "This is the output directory." ), make_option(c("-t", "--Type"), type = "character", default = "10X", action = "store", help = "This is type of exp matrix, default is 10X, the other is .csv which row is gene, column is cell." ), make_option(c("--datapath"), type = "character", default = FALSE, action = "store", help = "This is the path of exp matrix." ), make_option(c("--project_name"), type = "character", default = FALSE, action = "store", help = "This is project name built in seurat object and the pre_fix of files generated by this function." ), make_option(c("--mt_pattern"), type = "character", default = "^mt-", action = "store", help = "This is mitochondria pattern used to calculate its percentage, default is the ^mt-." ), make_option(c("--rb_pattern"), type = "character", default = "^Rpl|^Rps", action = "store", help = "This is ribosome pattern used to calculate its percentage, default is ^Rpl|^Rps" ), make_option(c("--min_cell"), type = "integer", default = 3, action = "store", help = "This is min cells number to filter, default is 3." ), make_option(c("--min_gene"), type = "integer", default = 0, action = "store", help = "This is min gene number to filter, default is 0." ) ) # -help opt = parse_args(OptionParser(option_list = option_list, usage = "This script is general overview data quality of single cell exp matrix!")) print(opt) #### step0.load data #### source("/home/zhengjh/scripts/seurat/function/SeuratWrapperFunction.R") print("Step0: Load Data") ##### step00. load functions #### # read in the functions timestart<-Sys.time() # set the output path cat("Current working dir: ", opt$output) cat("\n") setwd(opt$output) # load ref object #### step01. load exp matrix #### print("Step0: Load exp matrix data") if(opt$Type=="10X"){ print("Read exp matrix from 10X cellranger output") exp_count <- Read10X(data.dir = opt$datapath) print("The rownum and colnum of exp_count is ") print(dim(exp_count)) }else{ print("Read exp matrix from csv file") exp_count <- read.csv(opt$datapath, header = T, row.names = 1) print("The rownum and colnum of normal exp is ") print(dim(exp_count)) } #### step1: Create_seurat_filter_cell ###### print("Step1: Creat Seuart Object with filtering lower genes and lower cells") Or_ident <- rep(opt$project_name, ncol(exp_count)) exp_count.seurat <- CreateFirstFilterSeuratObject(exp_count, opt$project_name, Or_ident, opt$min_cell, opt$min_gene, opt$mt_pattern, opt$rb_pattern) print("The row_num and col_num of the first filtered matrix:") print(dim(exp_count.seurat)) print("The cellnumber of samples:") print(table(exp_count.seurat$orig.ident)) name_pre <- paste0("filter_gene_expressed_", opt$min_cell, "cells", opt$min_gene, "genes") #### step2: Overview of the data quality ###### print("Step2: Overview of the data quality") SeuratQC(exp_count.seurat, name_pre) rds_name <- paste0(opt$project_name, "_dataquality_overview_seurat.rds") saveRDS(exp_count.seurat, file = rds_name) print("Happy~Finished!!!") timeend<-Sys.time() runningtime<-timeend-timestart print(runningtime)

3.LogNormalizeClusteringSingle.R脚本: LogNormalize数据并聚类

#!/home/zhengjh/miniconda3/envs/r403/bin/Rscript # parameter library(optparse) library(getopt) option_list <- list( make_option(c("-o", "--output"), type = "character", default = FALSE, action = "store", help = "This is the output directory." ), make_option(c( "--seuratobject"), type = "character", default = "10X", action = "store", help = "This is absoulte path of seuratobject." ), make_option(c("--min_gene"), type = "integer", default = 50, action = "store", help = "This is min gene number to filter, default is 50." ), make_option(c("--max_gene"), type = "integer", default = 10000, action = "store", help = "This is max gene number to filter, default is 10000." ), make_option(c("--min_ncountRNA"), type = "integer", default = 500, action = "store", help = "This is min ncountRNA to filter, default is 5000." ), make_option(c("--max_ncountRNA"), type = "integer", default = 100000, action = "store", help = "This is min ncountRNA to filter, default is 100000." ), make_option(c("--max_mt_percent"), type = "integer", default = 60, action = "store", help = "This is min mitochondria percentage to filter, default is 60." ), make_option(c("--max_rb_percent"), type = "integer", default = 60, action = "store", help = "This is min ribosome percentage to filter, default is 60." ), make_option(c("--file_prefix"), type = "character", action = "store", help = "This is number of file prefix to used in the following analysis." ), make_option(c("--scale_factor"), type = "integer", default = 10000, action = "store", help = "This is scale factor used in Seurat NormalizeData function, default is 10000." ), make_option(c("--nfeatures"), type = "integer", default = 2000, action = "store", help = "This is number of features used in Seurat FindVariableFeatures function, default is 2000." ), make_option(c("--npc_plot"), type = "integer", default = 50, action = "store", help = "This is number of principles to plot, default is 50." ), make_option(c("--npc_used"), type = "integer", default = 20, action = "store", help = "This is number of principles to use in RunTSNE and RunUMAP function, default is 20." ), make_option(c("--k_parameter"), type = "integer", default = 20, action = "store", help = "This is k parameter used in FindNeighbors function which will influence the clusering number, default is 20." ), make_option(c("--resolution_number"), type = "double", default = 0.5, action = "store", help = "This is resolution number used in FindClusters function , default is 0.5" ), make_option(c("--min_pct"), type = "double", default = 0.25, action = "store", help = "This is min.pct used in FindAllMarkers function , default is 0.25" ), make_option(c("--logfc_threshold"), type = "double", default = 0.25, action = "store", help = "This is logfc.threshold used in FindClusters function , default is 0.25" ), make_option(c("--p_val_adj"), type = "double", default = 0.05, action = "store", help = "This is p_val_adj used to filter marker , default is 0.05" ), make_option(c("--top_num"), type = "integer", default = 10, action = "store", help = "This is top num marker of each cluster, default is 10" ), make_option(c("--mt_rb_pattern"), type = "character", default = "^mt-|^Rpl-|^Rps-", action = "store", help = "This is mt_rb_pattern used to filter mt/rb genes found in cluster markers, default is ^mt-|^Rpl-|^Rps-" ) ) # -help opt = parse_args(OptionParser(option_list = option_list, usage = "This script is general processing of single cell exp matrix!")) print(opt) source("/home/zhengjh/scripts/seurat/function/SeuratWrapperFunction.R") ###### step0. subset seurat object #### timestart<-Sys.time() # set the output path setwd(opt$output) exp_count.seurat <- readRDS(opt$seuratobject) print("Step0: Subset Seurat Object") cat(opt$min_gene, "< gene number <", opt$max_gene, "\n") cat(opt$min_ncountRNA, "< ncountRNA <", opt$max_ncountRNA, "\n") cat("max_mt_percent ", opt$max_mt_percent, "\n") cat("max_rb_percent ", opt$max_rb_percent, "\n") exp_count.seurat <- SubsetSeuratData(exp_count.seurat = exp_count.seurat, min_gene = opt$min_gene, max_gene = opt$max_gene, min_ncountRNA = opt$min_ncountRNA, max_ncountRNA = opt$max_ncountRNA, max_mt_percent = opt$max_mt_percent, max_rb_percent = opt$max_rb_percent) #### step1. normalize data and run pca #### print("Step1: LogNormalize Data and Run PCA") cat("scale factor used in NormalizeData function is ", opt$scale_factor, "\n") cat("nfeatures used in Seurat FindVariableFeatures function is ", opt$nfeatures, "\n") cat("npc_plot used in ElbowPlot is ", opt$npc_plot, "\n") exp_count.seurat <- SelectPCsByLogNormalize(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, scale_factor = opt$scale_factor, nfeatures = opt$nfeatures, npc_plot = opt$npc_plot) #### step2. find clusters #### print("Step2: FindNeighbors and clusters and then run tsne and umap reduction") cat("k_param used in Seurat FindNeighbors function is ", opt$k_param, "\n") cat("resolution_number used in Seurat FindClusters function is ", opt$resolution_number, "\n") cat("npc_used used in RunTSNE and RUNUMAP is ", opt$npc_used, "\n") exp_count.seurat <- PlotCluster(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, npc_used = opt$npc_used, k_param = opt$k_param, resolution_number = opt$resolution_number) #### step3. find top markers ##### print("Step3: find top markers and then map cluster with top2 marke") cat("min.pct used in Seurat FindAllMarkers function is ", opt$min.pct, "\n") cat("logfc.threshold used in Seurat FindAllMarkers function is ", opt$logfc.threshold, "\n") cat("p_val_adj used to filter markers is ", opt$p_val_adj, "\n") cluster.markers <- FindmarkerForCluster(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, min.pct = opt$min_pct, logfc.threshold = opt$logfc_threshold, p_val_adj = opt$p_val_adj, mt_rb_pattern = opt$mt_rb_pattern) TopMarker <- TopMarkersInCluster(cluster.markers = cluster.markers, file_prefix = opt$file_prefix, top_num = opt$top_num) # Top markers heatmap P_heatmap_top_marker <- DoHeatmap(exp_count.seurat, features = TopMarker$gene) filename <- paste0(opt$file_prefix, "_Heatmap_plot_top_markers.pdf") folder_name <- "4.MarkersInCluster" pdf(filename, 30, 15) print(P_heatmap_top_marker) dev.off() file.copy(filename, folder_name, overwrite = T) #copy files file.remove(filename) #### step4. map top2 marker to each cluster ##### print("step4: map top2 marker to each cluster") exp_count.seurat <- MapTop2MarkerEachCluster(exp_count.seurat = exp_count.seurat, cluster.markers = cluster.markers, file_prefix = opt$file_prefix) #### step5. save the clustering plots aftering assign top2 markers ### print("step5. save the clustering plots aftering assign top2 markers") p1 <- DimPlot(exp_count.seurat, reduction = "umap",label = TRUE, pt.size = 0.5, label.size = 4.5, repel = TRUE, group.by = "ident") + NoLegend() p1 <- p1 + labs(title="Clustering UMAP Reduction") p1 <- p1 + theme(plot.title = element_text(hjust = 0.5)) p2 <- DimPlot(exp_count.seurat, reduction = "umap",label = TRUE, pt.size = 0.5, label.size = 4.5, repel = TRUE, group.by = "ident") p2 <- p2 + labs(title="Clustering TSNE Reduction") p2 <- p2 + theme(plot.title = element_text(hjust = 0.5)) plots <- plot_grid(p1, p2, ncol=2,rel_widths = c(1.7, 2.2)) cluster_pdfname <- paste0(opt$file_prefix, "_Clustering_TSNE_UMAP.pdf") folder_name <- "3.PlotCluster" pdf(cluster_pdfname, width = 18, height = 6) print(plots) dev.off() file.copy(cluster_pdfname, folder_name, overwrite = T)#拷贝文件 file.remove(cluster_pdfname) rds_name <- paste0(opt$file_prefix, "_lognormalize_clustering_seurat.rds.rds") saveRDS(exp_count.seurat, file = rds_name) print("Happy~Finished!!") timeend<-Sys.time() runningtime<-timeend-timestart print(runningtime)

4. SCTransformClusteringSingle.R 脚本: SCTtransform数据并聚类

#!/home/zhengjh/miniconda3/envs/r403/bin/Rscript # parameter library(optparse) library(getopt) option_list <- list( make_option(c("-o", "--output"), type = "character", default = FALSE, action = "store", help = "This is the output directory." ), make_option(c( "--seuratobject"), type = "character", default = "10X", action = "store", help = "This is absoulte path of seuratobject." ), make_option(c("--min_gene"), type = "integer", default = 50, action = "store", help = "This is min gene number to filter, default is 50." ), make_option(c("--max_gene"), type = "integer", default = 10000, action = "store", help = "This is max gene number to filter, default is 10000." ), make_option(c("--min_ncountRNA"), type = "integer", default = 500, action = "store", help = "This is min ncountRNA to filter, default is 5000." ), make_option(c("--max_ncountRNA"), type = "integer", default = 100000, action = "store", help = "This is min ncountRNA to filter, default is 100000." ), make_option(c("--max_mt_percent"), type = "integer", default = 60, action = "store", help = "This is min mitochondria percentage to filter, default is 60." ), make_option(c("--max_rb_percent"), type = "integer", default = 60, action = "store", help = "This is min ribosome percentage to filter, default is 60." ), make_option(c("--file_prefix"), type = "character", action = "store", help = "This is number of file prefix to used in the following analysis." ), make_option(c("--npc_plot"), type = "integer", default = 50, action = "store", help = "This is number of principles to plot, default is 50." ), make_option(c("--npc_used"), type = "integer", default = 20, action = "store", help = "This is number of principles to use in RunTSNE and RunUMAP function, default is 20." ), make_option(c("--k_parameter"), type = "integer", default = 20, action = "store", help = "This is k parameter used in FindNeighbors function which will influence the clusering number, default is 20." ), make_option(c("--resolution_number"), type = "double", default = 0.5, action = "store", help = "This is resolution number used in FindClusters function , default is 0.5" ), make_option(c("--min_pct"), type = "double", default = 0.25, action = "store", help = "This is min.pct used in FindAllMarkers function , default is 0.25" ), make_option(c("--logfc_threshold"), type = "double", default = 0.25, action = "store", help = "This is logfc.threshold used in FindClusters function , default is 0.25" ), make_option(c("--p_val_adj"), type = "double", default = 0.05, action = "store", help = "This is p_val_adj used to filter marker , default is 0.05" ), make_option(c("--top_num"), type = "integer", default = 10, action = "store", help = "This is top num marker of each cluster, default is 10" ), make_option(c("--mt_rb_pattern"), type = "character", default = "^mt-|^Rpl-|^Rps-", action = "store", help = "This is mt_rb_pattern used to filter mt/rb genes found in cluster markers, default is ^mt-|^Rpl-|^Rps-" ) ) # -help opt = parse_args(OptionParser(option_list = option_list, usage = "This script is general processing of single cell exp matrix!")) print(opt) source("/home/zhengjh/scripts/seurat/function/SeuratWrapperFunction.R") ###### step0. subset seurat object #### timestart<-Sys.time() # set the output path setwd(opt$output) exp_count.seurat <- readRDS(opt$seuratobject) print("Step0: Subset Seurat Object") cat(opt$min_gene, "< gene number <", opt$max_gene, "\n") cat(opt$min_ncountRNA, "< ncountRNA <", opt$max_ncountRNA, "\n") cat("max_mt_percent ", opt$max_mt_percent, "\n") cat("max_rb_percent ", opt$max_rb_percent, "\n") exp_count.seurat <- SubsetSeuratData(exp_count.seurat = exp_count.seurat, min_gene = opt$min_gene, max_gene = opt$max_gene, min_ncountRNA = opt$min_ncountRNA, max_ncountRNA = opt$max_ncountRNA, max_mt_percent = opt$max_mt_percent, max_rb_percent = opt$max_rb_percent) #### step1. normalize data and run pca #### print("Step1: SCTtransform Data and Run PCA") cat("scale factor used in NormalizeData function is ", opt$scale_factor, "\n") cat("nfeatures used in Seurat FindVariableFeatures function is ", opt$nfeatures, "\n") cat("npc_plot used in ElbowPlot is ", opt$npc_plot, "\n") exp_count.seurat <- SelectPCsBySCTransform(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, npc_plot = opt$npc_plot) #### step2. find clusters #### print("Step2: FindNeighbors and clusters and then run tsne and umap reduction") cat("k_param used in Seurat FindNeighbors function is ", opt$k_param, "\n") cat("resolution_number used in Seurat FindClusters function is ", opt$resolution_number, "\n") cat("npc_used used in RunTSNE and RUNUMAP is ", opt$npc_used, "\n") exp_count.seurat <- PlotCluster(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, npc_used = opt$npc_used, k_param = opt$k_param, resolution_number = opt$resolution_number) #### step3. find top markers ##### print("Step3: find top markers and then map cluster with top2 marke") cat("min.pct used in Seurat FindAllMarkers function is ", opt$min.pct, "\n") cat("logfc.threshold used in Seurat FindAllMarkers function is ", opt$logfc.threshold, "\n") cat("p_val_adj used to filter markers is ", opt$p_val_adj, "\n") cluster.markers <- FindmarkerForCluster(exp_count.seurat = exp_count.seurat, file_prefix = opt$file_prefix, min.pct = opt$min_pct, logfc.threshold = opt$logfc_threshold, p_val_adj = opt$p_val_adj, mt_rb_pattern = opt$mt_rb_pattern) TopMarker <- TopMarkersInCluster(cluster.markers = cluster.markers, file_prefix = opt$file_prefix, top_num = opt$top_num) # Top markers heatmap P_heatmap_top_marker <- DoHeatmap(exp_count.seurat, features = TopMarker$gene) filename <- paste0(opt$file_prefix, "_Heatmap_plot_top_markers.pdf") folder_name <- "4.MarkersInCluster" pdf(filename, 30, 15) print(P_heatmap_top_marker) dev.off() file.copy(filename, folder_name, overwrite = T)# copy files file.remove(filename) #### step4. map top2 marker to each cluster ##### print("step4: map top2 marker to each cluster") exp_count.seurat <- MapTop2MarkerEachCluster(exp_count.seurat = exp_count.seurat, cluster.markers = cluster.markers, file_prefix = opt$file_prefix) #### step5. save the clustering plots aftering assign top2 markers ### print("step5. save the clustering plots aftering assign top2 markers") p1 <- DimPlot(exp_count.seurat, reduction = "umap",label = TRUE, pt.size = 0.5, label.size = 4.5, repel = TRUE, group.by = "ident") + NoLegend() p1 <- p1 + labs(title="Clustering UMAP Reduction") p1 <- p1 + theme(plot.title = element_text(hjust = 0.5)) p2 <- DimPlot(exp_count.seurat, reduction = "umap",label = TRUE, pt.size = 0.5, label.size = 4.5, repel = TRUE, group.by = "ident") p2 <- p2 + labs(title="Clustering TSNE Reduction") p2 <- p2 + theme(plot.title = element_text(hjust = 0.5)) plots <- plot_grid(p1, p2, ncol=2,rel_widths = c(1.7, 2.2)) cluster_pdfname <- paste0(opt$file_prefix, "_Clustering_TSNE_UMAP.pdf") folder_name <- "3.PlotCluster" pdf(cluster_pdfname, width = 18, height = 6) print(plots) dev.off() file.copy(cluster_pdfname, folder_name, overwrite = T)#拷贝文件 file.remove(cluster_pdfname) rds_name <- paste0(opt$file_prefix, "_sctransform_clustering_seurat.rds.rds") saveRDS(exp_count.seurat, file = rds_name) print("Happy~Finished!!") timeend<-Sys.time() runningtime<-timeend-timestart print(runningtime)

5.附上使用说明

1) 把SeuratWrapperFunction.R模块下面的代码复制到文件SeuratWrapperFunction.R，并存储

2) 把DataQualityOverviewSingle.R模块下面的代码复制到文件DataQualityOverviewSingle.R，并存储

3) 把LogNormalizeClusteringSingle.R模块下面的代码复制到文件LogNormalizeClusteringSingle.R，并存储

4) 把SCTransformClusteringSingle.R模块下面的代码复制到文件SCTransformClusteringSingle.R，并存储

5) 可执行权限

chmod +x DataQualityOverviewSingle.R chmod +x LogNormalizeClusteringSingle.R chmod +x SCTransformClusteringSingle.R

6)使用

## 质量控制 Rscript ~/scripts/seurat/pipeline/DataQualityOverviewSingle.R -o /home/zhengjh/other/PinealGland/Results/Zebrafish -t 10X --datapath /home/zhengjh/other/PinealGland/Data/Zebrafish_GSE123778_RAW/GSM3511192_scSeq1 \ --project_name "zebrafish_scseq1" --mt_pattern "^mt-" --rb_pattern "^rpl|^rps" --min_cell 3 --min_gene 100 ## LogNormalize Rscript ~/scripts/seurat/pipeline/LogNormalizeClusteringSingle.R -o ~/other/PinealGland/Results/Zebrafish --seuratobject ~/other/PinealGland/Results/Zebrafish/zebrafish_scseq1_dataquality_overview_seurat.rds \ --min_gene 200 --max_gene 6500 --min_ncountRNA 200 --max_ncountRNA 75000 \ --max_mt_percent 20 --max_rb_percent 20 --file_prefix zebrafish_scseq1_logtransform \ --scale_factor 10000 --nfeatures 2000 --npc_plot 50 --npc_used 30 \ --k_parameter 20 --resolution_number 0.5 --min_pct 0.25 --logfc_threshold 0.25 \ --p_val_adj 0.05 --top_num 10 --mt_rb_pattern "^mt-|rpl|rps" ## SCTransform Rscript ~/scripts/seurat/pipeline/SCTransformClusteringSingle.R -o ~/other/PinealGland/Results/Zebrafish --seuratobject ~/other/PinealGland/Results/Zebrafish/zebrafish_scseq1_dataquality_overview_seurat.rds --min_gene 200 --max_gene 6500 --min_ncountRNA 200 --max_ncountRNA 75000 \ --max_mt_percent 20 --max_rb_percent 20 --file_prefix zebrafish_scseq1_scttransform \ --npc_plot 50 --npc_used 30 \ --k_parameter 20 --resolution_number 0.5 --min_pct 0.25 --logfc_threshold 0.25 \ --p_val_adj 0.05 --top_num 10 --mt_rb_pattern "^mt-|rpl|rps"

7) 结果展示

最后放个图片
祝君开心快乐 健康！
也希望自己一直保持学习 不郁闷呜呜u