（一）单细胞数据分析——单细胞数据预处理

admin 管理员组

文章数量: 1086943

（一）单细胞数据分析——单细胞数据预处理

由于毕业设计是单细胞数据的处理，所以把整个过程所用到的方法进行一个整理，这是第一个部分，对得到的单细胞数据进行质控、降维、聚类等预处理。下面开始：

第一步：导入R包（部分R包可能用不到，因为做课题的时候需要就全部导入了，无伤大雅！）

library(scibet)
library(Seurat)
library(scater)
library(scran)
library(dplyr)
library(Matrix)
library(cowplot)
library(ggplot2)
library(harmony)
rm(list = ls())

第二步：读入所需要的数据

这里是我使用的一个乳腺癌单细胞测序数据

file.dir <- './GSE161529_RAW'
samplefile <- list.files('./GSE161529_RAW')
samplefile.mtx <- samplefile[grepl(pattern = 'matrix.mtx.gz',x = samplefile)]
samplefile.barcodes <- samplefile[grepl(pattern = 'barcodes.tsv.gz',x = samplefile)]
names(samplefile.mtx) <- substr(x = samplefile.mtx,start = 1,stop = 10)
names(samplefile.barcodes) <- substr(x = samplefile.barcodes,start = 1,stop = 10)
feature.file <- 'GSE161529_features.tsv.gz'
feature.names <- read.table(file = feature.file, header = FALSE, sep = "\t", row.names = NULL)
select_data <- c("GSM******") #这里的******指的是样本编号
GSE161529List <- lapply(select_data, function(sampleName){print(sampleName)tmp.barcode <- paste(file.dir, samplefile.barcodes[sampleName], sep = "/")tmp.mtx <- paste(file.dir, samplefile.mtx[sampleName], sep = "/")cell.barcodes <- read.table(file = tmp.barcode, header = FALSE, sep = "\t", row.names = NULL)data <- readMM(file = tmp.mtx)colnames(data) <- cell.barcodes[,1]rownames(data) <- feature.names[,2]CreateSeuratObject(counts = data, project = sampleName)
})
names(GSE161529List) <- select_data
BRCA_GSE161529 <- merge(GSE161529List[[1]], GSE161529List[-1], add.cell.ids = names(GSE161529List))
BRCA_GSE161529$sampleNames <- BRCA_GSE161529$orig.ident
head(BRCA_GSE161529)

当出现以下弹出消息则表示正在读入：

 最终head()的结果如下：

第三步：质量控制、细胞周期异质性减弱与细胞的降维聚类

因为这些步骤都是普遍通用的，这里封装为SeuratPreTreatment()方法，通过四个质控指标nCount_RNA(每个细胞的read counts数目)≥500、nFeature_RNA(每个细胞所检测到的基因数目) ≥200、mitoRatio（线粒体基因比例）<0.2、GenesPerUMI（单位read counts检测到的基因数目的比例）>0.8，进行质量控制；使用CellCycleScoring()减弱细胞周期异质性对下游分析的影响；之后就是标化、降维PCA/Harmony（Harmony是针对需要消除批次效应时使用，不一定都要使用）与聚类。

SeuratPreTreatment <- function(obj=NULL,QC=FALSE,nCount_cut=500,nFeature_cut=200,log10GenesPerUMI_cut=0.8,mitoRatio_cut=0.2){if(is.null(obj)){print("obj is null!")return(NULL)}obj$mitoRatio <- PercentageFeatureSet(object = obj, pattern = "^MT-")obj$mitoRatio <- obj@meta.data$mitoRatio / 100obj$log10GenesPerUMI <- log10(obj$nFeature_RNA)/log10(obj$nCount_RNA)obj$riboRatio <- PercentageFeatureSet(object = obj, pattern = "^RP[SL]")obj$riboRatio <- obj@meta.data$mitoRatio / 100obj <- CellCycleScoring(obj, g2m.features = g2m_genes, s.features = s_genes)if(QC){filtered_obj <- subset(x = obj, subset= (nCount_RNA >= nCount_cut) & (nFeature_RNA >= nFeature_cut) & (log10GenesPerUMI > log10GenesPerUMI_cut) & (mitoRatio < mitoRatio_cut))print(dim(filtered_obj))filtered_obj <- NormalizeData(filtered_obj, normalization.method = "LogNormalize")filtered_obj <- FindVariableFeatures(filtered_obj, selection.method = "vst", nfeatures = 3000)filtered_obj <- ScaleData(filtered_obj, vars.to.regress = c('nCount_RNA'))filtered_obj <- RunPCA(filtered_obj)filtered_obj <- RunUMAP(filtered_obj,reduction = "pca",dims = 1:30,seed.use = 12345)filtered_obj <- FindNeighbors(filtered_obj,reduction = 'pca', dims = 1:30, verbose = FALSE)filtered_obj <- FindClusters(filtered_obj,resolution = 0.5, verbose = FALSE,random.seed=20220727)return(filtered_obj)}else{obj <- NormalizeData(obj, normalization.method = "LogNormalize")obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 3000)obj <- ScaleData(obj, vars.to.regress = c('nCount_RNA'))obj <- RunPCA(obj)obj <- RunUMAP(obj,reduction = "pca",dims = 1:30,seed.use = 12345)obj <- FindNeighbors(obj,reduction = 'pca', dims = 1:30, verbose = FALSE)obj <- FindClusters(obj,resolution = 0.5, verbose = FALSE,random.seed=20220727)return(obj)}
}

方法封装完成后，质控参数QC设置为TRUE，使用该方法预处理数据

BRCA_GSE161529_QC <- SeuratPreTreatment(BRCA_GSE161529, QC = TRUE)

第四步：细胞注释

使用经典的细胞代表的基因集对数据进行注释，基因集详见代码；同时将细胞注释的代码封装为方法cellAnnotation()

cellAnnotation <- function(obj,markerList,assay='SCT',slot='data'){markergene <- unique(do.call(c,markerList))Idents(obj) <- 'seurat_clusters'cluster.averages <- AverageExpression(obj,assays=assay, slot = slot,features=markergene, return.seurat = TRUE)if(assay=='SCT'){scale.data <- cluster.averages@assays$SCT@data}else{scale.data <- cluster.averages@assays$RNA@data}print(scale.data)cell_score <- sapply(names(markergeneList),function(x){tmp <- scale.data[rownames(scale.data)%in%markergeneList[[x]],]if(is.matrix(tmp)){if(nrow(tmp)>=2){res <- apply(tmp,2,max)return(res)}else{return(rep(-2,ncol(tmp)))}}else{return(tmp)}})print(cell_score)celltypeMap <- apply(cell_score,1,function(x){colnames(cell_score)[which(x==max(x))]},simplify = T)obj@meta.data$cellType_auto <- plyr::mapvalues(x = obj@active.ident, from = names(celltypeMap), to = celltypeMap)return(obj)
}
lymphocyte <- c('CD3D','CD3E','CD79A','MS4A1','MZB1')
myeloid <- c('CD68','CD14','TPSAB1' , 'TPSB2','CD1E','CD1C','LAMP3', 'IDO1')
EOC <- c('EPCAM','KRT19','CD24')
fibo_gene <- c('DCN','FAP','COL1A2')
endo_gene <- c('PECAM1','VWF')
markergeneList <- list(lymphocyte=lymphocyte,myeloid=myeloid,Epi=EOC,fibo=fibo_gene,endo=endo_gene)
BRCA_GSE161529_QC <- cellAnnotation(obj=BRCA_GSE161529_QC,assay = 'RNA',markerList=markergeneList)

可以通过DimPlot()查看预处理后的单细胞UMAP图

DimPlot(object = BRCA_GSE161529_QC,reduction = 'umap',group.by = c('seurat_clusters','cellType_auto'), label = TRUE)

第五步：存入预处理后的单细胞数据

saveRDS(object = BRCA_GSE161529_QC,file = './BRCA_GSE161529_obj.RDS')

第（一）部分单细胞数据预处理就完成了，之后还有对细胞亚群的分群，伪时序分析等教程哦！关注不迷路！

本文标签： （一）单细胞数据分析单细胞数据预处理