企业官网建设流程全解析

##系统报错改为英文 Sys.setenv(LANGUAGE = "en") ##禁止转化为因子 options(stringsAsFactors = FALSE) ##清空环境 rm(list=ls()) library(Seurat) library(tidyverse) library(dplyr) library(ggplot2) library(patchwork) library(cowplot) library(SingleR) library(celldex) library(clustree) library(scDblFinder) loaded_pkgs <- sessionInfo()$otherPkgs data.frame( Package = names(loaded_pkgs), Version = sapply(loaded_pkgs, function(x) x$Version), row.names = NULL ) .libPaths() setwd("/mnt/DATA/home/jiangxingyu109/scRNAseq/GSExxxxx-rawdata/") ##储存GEO上下载的数据到目标目录下，10x需要的3个文件 ####数据读取#### dir <- c("BC2/" , "BC3/" , "BC5/" , "BC6/", "BC10/", "BC11/", "BC16/", "BC17/", "BC20/", "BC21/" ,"BC22/") names(dir) <- c("BC2", "BC3", "BC5" ,"BC6" , "BC10", "BC11" ,"BC16", "BC17", "BC20" ,"BC21" ,"BC22") scRNAlist <- list() for (i in seq_along(dir)) { counts <- Read10X(data.dir = dir[i]) sample_name <- basename(dir[i]) scRNAlist[[i]] <- CreateSeuratObject( counts = counts, project = sample_name, min.cells = 3, min.features = 300 ) scRNAlist[[i]] <- RenameCells( scRNAlist[[i]], add.cell.id = sample_name ) } setwd("/mnt/DATA/home/jiangxingyu109/scRNAseq/") save(scRNAlist, file = "scRNAlist260522new.Rdata") ##26年5月更新版本，原始数据汇总，没有进行质量控制，这一版是用新版的seurat做的 ##首先加载scRNAlist2605.Rdata ####质量控制#### # 2) 逐样本 QC 指标（含 HB） for (i in 1:length(scRNAlist)) { sc <- scRNAlist[[i]] sc[["percent.mt"]] <- PercentageFeatureSet(sc, pattern = "^MT-") HB_genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ") HB_m <- match(HB_genes, rownames(sc)) HB_genes <- rownames(sc)[HB_m] HB_genes <- HB_genes[!is.na(HB_genes)] sc[["percent.HB"]] <- PercentageFeatureSet(sc, features = HB_genes) scRNAlist[[i]] <- sc rm(sc) } ##先NormalizeData for (i in seq_along(scRNAlist)) { scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]], verbose = FALSE) } ##查看的小提琴图 violin_before <- list() for(i in 1:length(scRNAlist)){ violin_before[[i]] <- VlnPlot(scRNAlist[[i]], features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.HB"), pt.size = 0.01, ncol = 4) + plot_annotation(title = paste("Sample", i, "- Before QC"), theme = theme(plot.title = element_text(hjust = 0.5, face = "bold"))) } violin_before violin_before[[1]] violin_before[[2]] violin_before[[3]] violin_before[[4]] violin_before[[5]] violin_before[[6]] violin_before[[7]] violin_before[[8]] violin_before[[9]] violin_before[[10]] violin_before[[11]] #### 保存 质量控制 前小提琴图为 SCI TIFF #### dir.create("QC_violin_before_TIFF", showWarnings = FALSE) names(scRNAlist) <- c("BC2", "BC3", "BC5", "BC6", "BC10", "BC11", "BC16", "BC17", "BC20", "BC21", "BC22") violin_before <- list() for(i in 1:length(scRNAlist)){ sample_name <- names(scRNAlist)[i] p <- VlnPlot( scRNAlist[[i]], features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.HB"), pt.size = 0.01, ncol = 4 ) + plot_annotation( title = paste0(sample_name, " - Before QC"), theme = theme( plot.title = element_text( hjust = 0.5, face = "bold", size = 16 ) ) ) & theme_classic(base_size = 14) & theme( axis.title = element_text(size = 13, face = "bold"), axis.text = element_text(size = 11, color = "black"), plot.title = element_text(size = 13, face = "bold", hjust = 0.5), legend.position = "none" ) violin_before[[sample_name]] <- p tiff( filename = paste0("QC_violin_before_TIFF/", sample_name, "_QC_before.tiff"), width = 12, height = 4, units = "in", res = 600, compression = "lzw" ) print(p) dev.off() } #### 进行样本过滤#### ##逐样本过滤（保持你的阈值） scRNAlist <- lapply(scRNAlist, function(x) { subset(x, subset = nFeature_RNA > 300 & percent.mt < 10) }) ####去除双细胞#### set.seed(123) for(i in 1:length(scRNAlist)){ # 转换为SingleCellExperiment sce <- as.SingleCellExperiment(scRNAlist[[i]]) # 运行scDblFinder sce <- scDblFinder(sce) # 查看该样本的双细胞统计 cat("样本", i, "双细胞统计:\n") print(table(sce$scDblFinder.class)) # 添加回Seurat对象 scRNAlist[[i]]$doublet_class <- sce$scDblFinder.class scRNAlist[[i]]$doublet_score <- sce$scDblFinder.score # 记录过滤信息 before <- ncol(scRNAlist[[i]]) scRNAlist[[i]] <- subset(scRNAlist[[i]], doublet_class == "singlet") after <- ncol(scRNAlist[[i]]) cat(paste("样本", i, "过滤:", before, "→", after, "细胞, 去除", before - after, "个双细胞\n\n")) } ##这一步时间较长，会迭代很多次 ##合并样本（保留 orig.ident） dir <- c("BC2/" , "BC3/" , "BC5/" , "BC6/", "BC10/", "BC11/", "BC16/", "BC17/", "BC20/", "BC21/" ,"BC22/") names(dir) <- c("BC2", "BC3", "BC5" ,"BC6" , "BC10", "BC11" ,"BC16", "BC17", "BC20" ,"BC21" ,"BC22") scRNAlist_merge <- merge(x = scRNAlist[[1]], y = scRNAlist[-1], add.cell.ids = names(dir)) #### 合并后 QC 小提琴图 #### p_qc_merge <- VlnPlot( object = scRNAlist_merge, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.HB"), group.by = "orig.ident", pt.size = 0.01, ncol = 4 ) + theme_classic(base_size = 14) + theme( axis.title = element_text(size = 13, face = "bold"), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(size = 11, color = "black"), plot.title = element_text(size = 13, face = "bold", hjust = 0.5), legend.position = "none" ) p_qc_merge tiff( filename = "Merged_samples_QC_violin.tiff", width = 14, height = 6, units = "in", res = 600, compression = "lzw" ) print(p_qc_merge) dev.off() save(scRNAlist_merge, file = "scRNAlist_merge2605new.Rdata") ##处理过程详细记录，是新版 ##这里为止只保存两个处理数据，数据经过NormalizeData、nFeature_RNA > 300、percent.mt 过滤、scDblFinder 去双细胞 ##新版的S4对象处理跟以前不相同，以后的包肯定越来越新，下方用的scRNAlist_merge2605new.Rdata读取并命名的scRNA1 ##系统报错改为英文 Sys.setenv(LANGUAGE = "en") ##禁止转化为因子 options(stringsAsFactors = FALSE) ##清空环境 rm(list=ls()) ##防呆设计 scRNA1 = scRNAlist_merge meta <- scRNA1@meta.data qc_ncount_summary <- do.call(rbind, lapply(split(meta, meta$orig.ident), function(df) { x <- df$nCount_RNA data.frame( sample = unique(df$orig.ident), min = min(x, na.rm = TRUE), q25 = as.numeric(quantile(x, 0.25, na.rm = TRUE)), median = median(x, na.rm = TRUE), q75 = as.numeric(quantile(x, 0.75, na.rm = TRUE)), q95 = as.numeric(quantile(x, 0.95, na.rm = TRUE)), q99 = as.numeric(quantile(x, 0.99, na.rm = TRUE)), q995 = as.numeric(quantile(x, 0.995, na.rm = TRUE)), max = max(x, na.rm = TRUE) ) })) qc_ncount_summary ##这里注意有细胞RNA数拖尾，后续差异分析要注意，就是ncount 99.5%和最大值差距很大，这一部分细胞需要去除 ##harmony整合 scRNA_harmony = scRNA1 scRNA_harmony <- NormalizeData(scRNA_harmony) scRNA_harmony <- FindVariableFeatures(scRNA_harmony) scRNA_harmony <- ScaleData(scRNA_harmony) scRNA_harmony <- RunPCA(scRNA_harmony) ##检查PCA，选择dims值 ElbowPlot(scRNA_harmony, ndims = 50) ##还是选择30 dims=30 set.seed(123) scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident", dims = 1:30) scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:30) scRNA_harmony <- FindClusters(scRNA_harmony, resolution = 0.15) ##我记得是0.15刚好是原文中的组数，后续进一步区分 scRNA_harmony <- RunUMAP(scRNA_harmony, dims = 1:30, reduction = "harmony") ##应该是13或14个 plot1 = DimPlot(scRNA_harmony, reduction = "umap", group.by = c("orig.ident")) plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by = c( "seurat_clusters")) plotc = plot1 + plot2 plotc ##tsne scRNA_harmony <- RunTSNE(scRNA_harmony, reduction = "harmony", dims = 1:30) plot3 = DimPlot(scRNA_harmony, reduction = "tsne", group.by = c("orig.ident")) plot4 = DimPlot(scRNA_harmony, reduction = "tsne", group.by = c( "seurat_clusters")) plotd = plot3 + plot4 plotd

这是第一部分，截至到实现文章中的分组数