单细胞02.数据质控，单细胞数据过滤
关于生信学习，我们的建议是纯生信 → R 语言 → 单细胞、空间组及多组学三步走策略。费曼学习法可能是生信入门的好方法，但是真正提高还是要靠在课题进展中不断地实践！本次单细胞 21天教程以《癌生物学》《肿瘤免疫12讲2.0》为基础，并结合经典前沿综述和单细胞文章为基础，手把手展示单细胞数据分析流程。
第一周内容安排
单细胞02.数据质控，浅谈单细胞数据过滤！
单细胞03.数据构建，Seurat对象那些事儿！
单细胞04.数据降维，PCA+UMAP或tSNE！
单细胞05.细胞聚类，单细胞分群那些事儿！
单细胞06.数据整合，单细胞数据批次处理！
单细胞07.分群鉴定，细胞亚群注释那些事儿！
单细胞02.数据质控，浅谈单细胞数据过滤
学习内容和流程
1.阅读推文。了解单细胞分析的思路，通过Seurat官网了解单细胞数据分析流程，完成推文中的代码操作，理解代码的含义。
2.观看视频。根据视频解读，重点掌握单细胞数据分析的基本流程和相关概念。
3.发朋友圈。继续打卡，将推文转发至朋友圈，一句话总结学习感受或收获。
如果阅读了比较多10x单细胞转录组的文献，我们可能会发现单细胞数据分析的思路。归纳下来，单细胞数据大致分三层次：
①基本分析，包括测序统计信息（测序的饱和度、捕获细胞数、平均测序量和平均基因检测数），细胞过滤（检测基因数、检测的UMI数和线粒体风度百分比等）和降维聚类，有时还要涉及数据的合并与批次矫正。
②常规分析，包括细胞注释、差异分析、富集分析和细胞通讯等。常规分析更注重目标亚群的鉴定和分析。其中，有些内容属于高级分析，更注重与样本其他维度的信息进行衔接，侧重于细胞亚群的功能、通路和表型。
③高级分析，这部分往往是个性化的，会更注重与临床、分子机制及其功能的联系，也更注重与湿实验的相互验证。比如结合FISH或多重免疫荧光做验证，Marker基因与临床或表型的相关分析，预后分析等等。
单细胞测序可以获得每个样本中成千上万个细胞的信息，但不是所有的细胞数据都能满足分析要求。因此，在进行降维聚类等操作前，需要根据一定标准，过滤掉不合格的细胞。具体到10x Genomics测序，单细胞数据往往来自一个油包水磁珠为单元的单细胞扩增。但是，在10x的芯片管道中，有些管道捕获了裂解的细胞，这时只有少量的RNA被扩增，基因数会比较少；有些管道捕获了2个甚至2个以上的细胞，这时会有较多的RNA被扩增，基因数比较多，而且可能整合不同类型的细胞，造成干扰；有些管道捕获了衰老或凋亡的细胞，这些细胞往往不具有研究价值，除非研究对象就是衰老或凋亡细胞。
因此，在进行单细胞测序分析时，样本中细胞的活率、状态等非常重要！单细胞数据分析过程中，我们要过滤掉表达基因过多、过少或线粒体基因表达量过高的细胞（衰老或凋亡细胞的线粒体基因超表达）。在实操层面上，Seurat官网已经提供非常全面的单细胞数据质控的代码，内容包括数据前处理、数据标准化、鉴定高变基因和数据scaling。
代码（根据Seurat官网整理，最好收藏下来）
网址：https://satijalab.org/seurat/articles/pbmc3k_tutorial
################单细胞02.数据质控，单细胞数据过滤###################
rm(list = ls())
gc()
getwd()
#### 读取数据，Seurat数据对象构建 ####
library(dplyr)
library(Seurat)
library(patchwork)
# 注意更改个人路径
pbmc.data <- Read10X(data.dir = "~/pbmc3k_filtered_gene_bc_matrices/hg19/")
# 构建Seurat数据对象 (未标准化的Raw data)
pbmc <- CreateSeuratObject
(counts = pbmc.data, #单细胞表达矩阵数据的来源
project = "pbmc3k", #自定义的项目名称
min.cells = 3, #三个以上细胞表达的基因保留
min.features = 200) #保留表达200个以上基因的细胞
pbmc
An object of class Seurat
# 13714 features across 2700 samples within 1 assay，13174个基因，2700个细胞
# Active assay: RNA (13714 features, 0 variable features)
Seurat对象创建之后，我们就需要进行数据质控操作了。正如前面所说，我们要过滤掉表达基因过多、过少或线粒体基因表达量过高的细胞，官网也是这么解释的！
在实操中，Seurat官网教程首先处理线粒体基因，并通过双中括号[[ ]]操作符号，函数PercentageFeatureSet()来增加列名。操作后，metadata中增加percent.mt了。
# 通过[[ ]]操作增加metadata的列名pbmc[["percent.mt"]] = PercentageFeatureSet(pbmc, pattern = "^MT-")
通过小提琴图查看基因表达情况、细胞分布和线粒体基因。
# 用小提琴图可视化质控参数VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
通过散点图探究三个参数之间的关系。
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 = FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 = FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot3 = FeatureScatter(pbmc, feature1 = "nFeature_RNA", feature2 = "percent.mt")plot1 + plot2 + plot3
# 过滤掉基因数多于2,500或少于 200的细胞
# 过滤掉线粒体基因百分比大于5%的细胞pbmc = subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
鉴定高变基因。为什么进行这一步操作呢？可以理解为经验带来的启示。We and others have found that focusing on these genes in downstream analysis helps to highlight biological signal in single-cell datasets.（参考文献：Accounting for technical noise in single-cell RNA-seq experiments）。
###鉴定高变基因，FindVariableFeatures()函数
pbmc = FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 = head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 = VariableFeaturePlot(pbmc, pt.size = 2, raster = T)
plot2 = LabelPoints(plot = plot1, points = top10,
xnudge = 0,
ynudge = 0)
plot1 + plot2
数据Scaling。其中用到ScaleData()函数。为什么进行这一步操作呢？Scaling is an essential step in the Seurat workflow, but only on genes that will be used as input to PCA.
###数据Scaling
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- NormalizeData(pbmc)
all.genes = rownames(pbmc)
pbmc = ScaleData(pbmc, features = all.genes)
上述所有操作都是为了得到高质量的细胞表达信息，以用于后续分析。有些参数无需修改即可用于各种单细胞数据的处理，但是有些参数要适当调整，比如线粒体基因百分比等。这就需要我们在实际分析中，多做摸索和尝试，具体情况具体分析了。
最后编辑于 04-01 · 浏览 758
