单细胞数据挖掘 P2.2 构建Seurat对象，质控、绘图

「生信技能树」单细胞数据挖掘_哔哩哔哩_bilibili生信技能树又一单细胞系列课程上线啦~小伙伴们快来学习呀https://www.bilibili.com/video/BV1pa4y1s76J?p=5该章节承接P2.1单细胞课程

一、筛选思路及方法

# 2.1 构建seurat对象，质控
#In total, 2,343 cells from tumor cores were included in this analysis.
#quality controlstandards:
#1) genes detected in < 3 cells were excluded; 筛选基因
#2) cells with < 50 total detected genes were excluded; 筛选细胞
#3) cells with ≥ 5% of mitochondria-expressed genes were excluded. 筛选细胞

先设定筛选标准，需确定以下内容：

1、获取可用的细胞数

2、筛选基因：将表达太少的基因去掉（信息不够）

3、筛选细胞：将表达基因信息数量太少的细胞去掉（该基因信息不全）

4、筛选细胞：将线粒体基因表达过多的细胞去掉，认为检测到的有效基因类型不够，认为都是线粒体基因。

筛选细胞及可分析的基因。

sessionInfo()library("Seurat")?CreateSeuratObjectsce.meta <- data.frame(Patient_ID=group$Patient_ID,row.names = group$sample)head(sce.meta)table(sce.meta$Patient_ID)# 这个函数 CreateSeuratObject 有多种多样的执行方式scRNA = CreateSeuratObject(counts=a.filt,meta.data = sce.meta,min.cells = 3, min.features = 50)#counts:a matrix-like object with unnormalized data with cells as columns and features as rows #meta.data:Additional cell-level metadata to add to the Seurat object#min.cells: features detected in at least this many cells. #min.features:cells where at least this many features are detected.head(scRNA@meta.data)#nCount_RNA：the number of cell total counts#nFeature_RNA：the number of cell's detected genesummary(scRNA@meta.data)scRNA@assays$RNA@counts[1:4,1:4]# 可以看到，之前的counts矩阵存储格式发生了变化：4 x 4 sparse Matrix of class "dgCMatrix"dim(scRNA)#  20047  2342 仅过滤掉一个细胞

？可以查看seurat这个构建对象的语法的用法。这个可以自行查看。

构建creatseurateobject之前，需要先构建一个表达矩阵注释的metadata信息，这里主要是分组信息。本语句中命名为sce.meta。

creatseurateobject()语句本质上是需要一个counts矩阵，以及meta的信息

help中的文件笔记：

1、names.field以及names.delim可以在初期即设置细胞类型，需要选择相应的选择，具体需要查看详细的文档。

2、min.cells以及min.features是质控筛选的指标。

3、min.cells是筛选基因的，即该基因至少在3个细胞中表达，才认为个基因有意义。

4、min.features是筛选细胞的，即要求该种细胞必须要表达50个基因才认为是质量满意的细胞。纳入后续分析。

构建出的seurat对象是个多嵌套层级结构的对象，里面包含了许多内容。

可以通过语句查看各项信息，具体是个稀疏矩阵，包括 scRNA@assays$RNA@counts[1:4,1:4] 查看前4个样本的信息。

dim()可以看到底有多少细胞最终被筛选出来，可以看到其实只是筛选了一个细胞出去。n由2343变成了2342。

筛选线粒体基因（MT/ERCC为标记的两种基因）

 #接下来根据线粒体基因表达筛选低质量细胞#Calculate the proportion of transcripts mapping to mitochondrial genestable(grepl("^MT-",rownames(scRNA)))#FALSE #20050 没有染色体基因scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")head(scRNA@meta.data)summary(scRNA@meta.data)#结果显示没有线粒体基因，因此这里过滤也就没有意义，但是代码留在这里# 万一大家的数据里面有线粒体基因，就可以如此这般进行过滤啦。pctMT=5 #≥ 5% of mitochondria-expressed genesscRNA <- subset(scRNA, subset = percent.mt < pctMT)dim(scRNA)

MT开头往往是线粒体基因的标志，因此先看看有多少是MT开头的。这里使用grepl()函数来进行判断是否有MT开头的基因名。

对MT基因所占的比例进行赋值。本例中没有，但如果有就可以用其中的语言进行过滤了。

除了MT为标记的线粒体基因，还可能是ERCC为标记的基因为线粒体基因。

因此类似的方法可以再使用一遍。

table(grepl("^ERCC-",rownames(scRNA)))#FALSE  TRUE #19961    86  发现是有ERCC基因#External RNA Control Consortium，是常见的已知浓度的外源RNA分子spike-in的一种#指标含义类似线粒体含量，ERCC含量大，则说明total sum变小scRNA[["percent.ERCC"]] <- PercentageFeatureSet(scRNA, pattern = "^ERCC-")head(scRNA@meta.data)summary(scRNA@meta.data)rownames(scRNA)[grep("^ERCC-",rownames(scRNA))]#可以看到有不少ERCC基因sum(scRNA$percent.ERCC< 40)#较接近原文过滤数量2149，但感觉条件有点宽松了，先做下去看看#网上看了相关教程，一般ERCC占比不高于10%sum(scRNA$percent.ERCC< 10)   #就只剩下460个cell，明显低于文献中的数量pctERCC=40scRNA <- subset(scRNA, subset = percent.ERCC < pctERCC)dim(scRNA)# 20047  2142   原文为19752  2149dim(a.filt)#23460  2343 未过滤前

给scRNA ERCC标记基因给予一个percentage的参数，观察一下ERCC基因的含量%有多少。

用rownames考察有多少个ERCC相关的染色体基因。

分别用10和40来做cutoff看看，如果ERCC基因的数量<10和40分别有多少个基因，发现太严格之后呢就很少，因此选择40作为cutoff，尽量保持和原文一致的基因数量。

Tip：

grep()与grepl()函数，这两个函数最大的区别在于grep返回找到的位置，grepl返回是否包含要查找的内容。

具体见下：

>     #创建一个字符串向量
>     x <- c("d", "a", "c", "abba")
>     #查找包含a的元素所在的位置
>     grep("a", x)
[1] 2 4
>     #判断每个元素是否包含a，返回的是逻辑向量
>     grepl("a", x)
[1] FALSE  TRUE FALSE  TRUE
>
>     #同时匹配多个内容，查找包含a或者c的元素所在的位置
>     grep("a|c", x)
[1] 2 3 4
>
>     #同时匹配多个内容，判断每个元素是否包含a或者c，返回的是逻辑向量
>     grepl("a|c", x)
[1] FALSE  TRUE  TRUE  TRUE

三、简单可视化

  # 2.2 可视化#图A：观察不同组cell的counts、feature分布col.num <- length(unique(scRNA@meta.data$Patient_ID))library(ggplot2)p1_1.1 <- VlnPlot(scRNA,features = c("nFeature_RNA"),group.by = "Patient_ID",cols =rainbow(col.num)) +theme(legend.position = "none") +labs(tag = "A")p1_1.1p1_1.2 <- VlnPlot(scRNA,features = c("nCount_RNA"),group.by = "Patient_ID",cols =rainbow(col.num)) +theme(legend.position = "none") p1_1.2p1_1 <- p1_1.1 | p1_1.2p1_1VlnPlot(scRNA,features = c("nFeature_RNA","nCount_RNA","percent.ERCC"))#图B：nCount_RNA与对应的nFeature_RNA关系p1_2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",group.by = "Patient_ID",pt.size = 1.3) +labs(tag = "B")p1_2FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.ERCC")sessionInfo()

A图：小提琴图，拼接一下|

B图：散点图。