单细胞Seurat包升级之2,700 PBMCs分析（上）
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆写于19.7.18+28
官方总共有7个数据集，涵盖了新版本的不同亮点
这一次跟着教程https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html走一遍最简单的教程，记录自己的学习轨迹
隔了这么久，中间肯定发生了什么故事嘿嘿😜
飞机+高铁确实很适合学习，能稳下心来回顾自己的收获(还要一会才能到站，趁电脑还剩一点电量，车上完成下车轻松)
过去的10天，学到了很多东西，技能树单细胞培训、第一届生信人才发展论坛、今天来到北京，准备明天为期5天的北大龙星计划助教，不知道有没有认识的小伙伴参加。
前言
这个数据集是10X放出的2700个外周血单核细胞(PBMC，Peripheral Blood Mononuclear Cells )公共数据，使用 Illumina NextSeq 500测序，原始数据在：https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
这个教程更新于2019-6-24，它会做这么几件事：QC and data filtration, calculation of high-variance genes, dimensional reduction, graph-based clustering, and the identification of cluster markers.
从读取数据开始
10X的数据可以使用Read10X这个函数，会返回一个UMI count矩阵，其中的每个值表示每个基因(行表示)在每个细胞(列表示)的分子数量
加载PBMC数据
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
这个矩阵长什么样？和我们平常见到的bulk转录组矩阵一样吗？
找几个基因，看看前30个细胞的情况
> pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
[[ suppressing 30 column names ‘AAACATACAACCAC’, ‘AAACATTGAGCTAC’, ‘AAACATTGATCAGC’ ... ]]
CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5
TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . .
MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .
从这个结果可以得出两个结论：它是sparse Matrix；它很杂乱
但其实仔细看，这个.就表示空着，数值为0，也就是没有检测到该分子。因为单细胞数据和bulk转录组数据还不太一样，scRNA的大部分数值都是0，原因一是由于建库时细胞的起始反应模板浓度就很低；二是测序存在dropout现象(本来基因在这个细胞有表达量，但没有测到)。
这里Seurat采用了一种聪明的再现方式，比原来用0表示的矩阵大大减小了空间占用，可以对比一下：
dense.size <- object.size(as.matrix(pbmc.data))
dense.size # 709548272 bytes
sparse.size <- object.size(pbmc.data)
sparse.size # 29861992 bytes
dense.size/sparse.size #空间缩小23.8倍
然后利用这个矩阵构建`Seurat`对象
这个对象是一个容器，里面装了数据(比如表达矩阵)和分析结果(比如PCA、聚类)。另外关于这个对象的详细解释，见：https://github.com/satijalab/seurat/wiki/Seurat#slots
# 用原始数据（非标准化）先构建一个对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
> pbmc
An object of class Seurat
13714 features across 2700 samples within 1 assay
Active assay: RNA (13714 features)
发现这里有个active assay，它其中存的就是表达矩阵，用pbmc[["RNA"]]@counts就可以访问
使用两个中括号是对assay的快速访问，例如访问metadata：
pbmc[[]]就如同object@meta.data
另外看上面的图，seurat对象还有很多信息：
# 可以利用str来查看
> str(pbmc)
Formal class 'Seurat' [package "Seurat"] with 12 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay' [package "Seurat"] with 7 slots
.. .. .. ..@ counts :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..
# 后面还有很多
# 如果要查看其他的信息，可以用@，例如
> pbmc@version
[1] ‘3.0.0.9000’
预处理阶段
标准流程包括了根据QC 结果进行的细胞选择和过滤，标准化过程，检测显著差异基因
质控（QC）及筛选细胞
这篇文章：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/介绍了一些QC的标准
例如：
检测每个细胞中检测到的unique genes数量(这种情况，一般低质量的细胞或空的液滴"droplet"中基因数量也会较少；如果一个液滴中有两个细胞"doublets"或者存在多个细胞"multiplets"，这样会导致检测到的基因数量出奇的高)
利用细胞中检测到的molecules总量也可以(它的结果和unique genes方法高度相关)
比对到线粒体基因组的reads数量，因为一般低质量或死亡的细胞中会广泛存在线粒体基因组污染；可以利用
PercentageFeatureSet函数计算，顾名思义这个函数会计算来自某一个feature的reads count数量，需要做的就是挑出来以MT-开头的feature对应的所有基因。
# 质控筛选, [[符号可以在pbmc这个对象中添加metadata，利用这个方法还可以挑出其他的feature
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# 检查下metadata
> head(pbmc@meta.data, 3)
orig.ident nCount_RNA nFeature_RNA percent.mito percent.mt
AAACATACAACCAC 10X_PBMC 2419 779 0.030177759 3.0177759
AAACATTGAGCTAC 10X_PBMC 4903 1352 0.037935958 3.7935958
AAACATTGATCAGC 10X_PBMC 3147 1129 0.008897363 0.8897363
# 然后进行一个可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
这张图可以给出的提示就是过滤信息：过滤的时候可以过滤掉feature大于2500和小于200的(就是避免前面所说基因数过多/过少的情况)；然后线粒体占比图显示大部分都集中在5%以下，因此超过5%的可以过滤掉
# 另外可以用一个散点图（FeatureScatter）来绘制两组feature信息的相关性，最后组合在一起（CombinePlots）
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
这个图给出的意思是：(左边👈)有表达量的reads和在线粒体的reads数量一般成反比，线粒体reads占比很高的情况一般有表达量的很少；相反，如果是真正表达的基因reads，很少会来自线粒体基因组；(右边👉)就是刚才所说的比对结果中的基因数量(feature)一般会和测序得到的reads count值成正比
一切就绪，最后进行一个过滤操作
# 最后进行一个过滤操作
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
数据标准化
关于Seurat标准化，可以看这一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf
当过滤掉一部分细胞以后，就要会数据进行标准化。默认是进行一个全局的LogNormalize操作，具体方法是：
normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result
翻译成公式就是：log1p(value/colSums[cell-idx] *scale_factor)，其中log1p指的是log(x + 1)这里也有解释：https://github.com/satijalab/seurat/issues/833
当然，除了这一种默认的算法，还有：
CLR: Applies a centered log ratio transformation
RC: Relative counts. Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. No log-transformation is applied. For counts per million (CPM) set
scale.factor = 1e6
得到的结果存储在：pbmc[["RNA"]]@data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
# 当然其中都是默认参数，于是直接写这种也是可以的
pbmc <- NormalizeData(pbmc)
这一步要因人而异，可以看：https://bioinformatics.stackexchange.com/questions/5115/seurat-with-normalized-count-matrix 假入有一个TPM的count矩阵，那么就可以不需要使用Seurat::NormalizeData()操作了，因为TPM已经根据测序深度进行了标准化，只需要进行log降一下维度即可。如果后续进行ScaleData操作，程序会检测是否使用了Seurat提供的标准化方法，如果我们使用自己的的标准化数据，那么就可能出现一个warning提醒，不过到时候不想被提醒，可以设置check.for.norm =F
鉴定差异基因HVGs(highly variable features)
意思就是：在细胞与细胞间进行比较，选择表达量差别最大的(也即是同一个基因在有的细胞中表达量很高，同时在部分细胞中表达很少)，利用FindVariableFeatures函数，会计算一个mean-variance结果，也就是给出表达量均值和方差的关系并且得到top variable features。
那么关于如何计算这些top variable features，这个函数是有三种算法的，分别是：
(默认) vst：对方差(variance)和均值(mean)都取log值+loess线性拟合
mean.var.plot方法：average expression & dispersion (for each feature) + z-scoresdispersion方法：the highest dispersion values
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(pbmc), 10)
下一次就是Scale、Cluster等一系列操作
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~