太猛了！万字长文单细胞分析全流程讲解，看完就能发文章！建议收藏！（附代码）
Hello~大家好，这里是晨曦的单细胞笔记，上次的联合多数据集不知道大家学习的如何呢，考虑到各位小伙伴们希望扎扎实实的学习单细胞分析技术，这里特别准备了单细胞基础流程讲解，来为大家打下坚实的基础，同时为大家充分的打开单细胞分析的大门，只要大家扎扎实实的走下来，都会成为自己课题组单细胞分析的得力干将哦~
首次揭秘！不做实验也能发10 SCI，CNS级别空间转录组套路全解析（附超详细代码！）
过关神助！99%审稿人必问，多数据集联合分析，你注意到这点了吗？
本笔记的示例数据可以从Seurat官网获得
library(dplyr)library(Seurat)library(patchwork)# Load the PBMC datasetpbmc.data <- Read10X(data.dir = '../data/pbmc3k/filtered_gene_bc_matrices/hg19/')# Initialize the Seurat object with the raw (non-normalized data).pbmc <- CreateSeuratObject(counts = pbmc.data, project = 'pbmc3k', min.cells = 3, min.features = 200)pbmc晨曦解读
晨曦解读
然后我们来看一下我们构建的Seurat对象是什么样子的（其中assays存储的为表达矩阵；meta.data存储的为细胞信息；active.ident存储的为表达矩阵细胞名；active.assay存储的为表达矩阵的名）
#使用pbmc.data函数可以读取表达矩阵，这里我们简单看一下这三个基因在表达矩阵中的样子pbmc.data[c('CD3D', 'TCL1A', 'MS4A1'), 1:30]## 3 x 30 sparse Matrix of class 'dgCMatrix'#### 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 . . . .#可以看到上面的是稀疏矩阵的形式#然后我们简单探索一下如果不是稀疏矩阵，运算速度会是怎样的？dense.size <- object.size(as.matrix(pbmc.data))dense.size## 709591472 bytessparse.size <- object.size(pbmc.data)sparse.size## 29905192 bytesdense.size/sparse.size## 23.7 bytes#这里可以很清楚的看到稀疏矩阵确实可以明显的帮助我们减少内存的消耗
晨曦解读
稀疏矩阵指的是，我们认为细胞与细胞之间的表达大部分基因都应该是相似的，滑稽一下，如果大部分都不相似，可能真是变异了~然后这里我们只突出显示表达不同的基因表达，这样可以减少内存的消耗，这就是稀疏矩阵~（这个表达矩阵和bulk-seq表达矩阵类似，都是横坐标为Gene symbol，纵坐标为细胞名）
#上面构建完Seurat对象后，我们接下来第一步就是对我们的数据进行质控#首先我们计算一下线粒体Gene的比例（这里可以思考一下，为什么我们要计算线粒体基因的比例）pbmc[['percent.mt']] <- PercentageFeatureSet(pbmc, pattern = '^MT-')#计算完线粒体基因的比例后，我们可视化一下VlnPlot(pbmc, features = c('nFeature_RNA', 'nCount_RNA', 'percent.mt'), ncol = 3)晨曦解读
第一张图为不同样本上基因的表达数量，其中一个点代表一个细胞；第二张图为不同细胞RNA表达的数量；第三个图则代表线粒体的百分比，如果线粒体基因占比过高则说明有问题，细胞要么是坏死要么就是特殊细胞（心肌细胞或者神经细胞）根据需求往往是10%一下，但是根据特殊情况25%以下也有可能，所以需要我们认清楚我们的样本来源~（总而言之，这个可视化的目的就是为了评估我们的样本，也就是这些细胞的质量如何）
#上面这么看可能不够直观，下面再进行相关性的可视化plot1 <- FeatureScatter(pbmc, feature1 = 'nCount_RNA', feature2 = 'percent.mt')plot2 <- FeatureScatter(pbmc, feature1 = 'nCount_RNA', feature2 = 'nFeature_RNA')plot1 plot2
晨曦解读
第一张图，是样本中RNA的表达数量和线粒体基因的比例，两者之间没什么关系，这个才是正常的样本，因为线粒体基因的数量大约是300多个，而且十分稳定，基本上细胞坏死，RNA降解掉，线粒体基因的数量都不会大幅度改变，所以如果RNA数量降低，而线粒体基因升高，则提示细胞可能存在坏死；第二章图，基因的数量和RNA的数量呈正相关，这个也是很好理解的，所以根据这个可视化评估，我们的细胞是正常且合格的（其实这个图就是看趋势，上面的那三个小提琴图则是看数量，一般来说线粒体基因比例在10%一下，然后基因的数量和RNA的数量只要大部分比较密集即可，因为这个根据样本的变化，这两个数值也会有不同的变化，比如说有的要求基因的数量是4000左右，RNA的数量为20000-25000左右，和这个不一样则只能说明样本的不同）
#然后这里我们设置如下标准，作为我们的质控指标pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)晨曦解读
既然我们对数据进行了改动，那么我们下一步自然是对数据进行标准化，来平衡数据之间的改动，这里其实涉及到更多数据处理的原因，但是咱们毕竟是生信分析，这样理解更好方便记忆，感兴趣的同学可以去查阅相关资料（那么扩展一下，后面也有一步对数据进行挑选，随后当然也得接上数据标准化的处理）
pbmc <- NormalizeData(pbmc, normalization.method = 'LogNormalize', scale.factor = 10000)#当然这里我们也可以简化代码，变成以下的样子pbmc <- NormalizeData(pbmc)#然后接下来我们就要寻找高可变基因了，关于高可变基因的解读我们后面会进行解答pbmc <- FindVariableFeatures(pbmc, selection.method = 'vst', nfeatures = 2000)#这里默认寻找2000个高可变基因，当然自己也可以人为设置，一般是1000-5000之间# 这里我们展示一下前10个改变最大的基因top10 <- head(VariableFeatures(pbmc), 10)# 对高可变基因进行可视化plot1 <- VariableFeaturePlot(pbmc)plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)plot1 plot2#下面可视化的解读则是，左边和右边的图表达的是一个东西，只不过右边的图把表达差异在各个样本中最显著的基因表示了出来，然后高可变基因用红色表示
晨曦解读
高可变基因即，我们认为细胞与细胞之间的基因表达其实趋于相似，就是说表达你高一点，我少一点其实并不能表现功能，只有那种表达有显著差异的才会被我们认为是后续研究的主要基因（毕竟scRNA-seq的数据量太大，这样的目的其实是为了减少我们分析的数据量，减少我们“无用”的数据运算）
#前面我们提供如果对数据进行了挑选和改变就要对数据进行标准化来平衡数据之间的改动~all.genes <- rownames(pbmc)pbmc <- ScaleData(pbmc, features = all.genes)#这里默认对所有基因进行线性归一化，这一步如果没有服务器往往运行比较慢，我们也可以使用其默认形式，只对高可变基因进行数据标准化的操作pbmc <- ScaleData(pbmc)晨曦解读
这个时候我们回过头来，再来观察一下我们一开始构建的Seurat对象，为什么这里我们要从头观察，不知道大家有没有这样一个疑惑，就是我这里每一步都在对同样一个变量名进行赋值，那么按照我们以前学习的经验，变量的赋值是不断的叠加的，也就是说新的赋值会代替旧的赋值，那么对于Seurat对象是这样吗？
大家从下面这个图可以很清楚的看到，我这里从新展示了经过一系列操作的Seurat对象，在assays这里的表达矩阵被规整的整理成了data即原始表达矩阵，以及scale.data经过标准化的表达矩阵，所以说并没有被变量的重复赋值给替换掉，这里可能会有同学问，这是因为什么？
这是因为Seurat的数据格式是HDF5类型，当然知道这些对我们的分析没有太多的用处，大家只需要知道，我们每一步得到的数据并没有消失，至于数据的具体格式，大家感兴趣可以自行搜索学习~
#接下来我们就要进行PCA降维，具体解释我们会在后续解答~pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
晨曦解读
所谓的PCA降维即，我们可以把一个样本看成是一个维度，我们都知道我们所处的世界是三维世界，对于拥有上万个样本，也就是维度的矩阵，我们压根就无法理解这么高维度的信息，所以我们要把其进行降维，可以简单理解为一个总结
举个例子：我们对于又高又瘦皮肤白的男生，我这么形容消耗了我很多口舌，但是我可以用一个字来概括，就是帅~这就是降维，这就是PCA
#然后针对降维的结果表示为PC的形式，针对展示PC的可视化这里也提供了两种，但是为了教程的连贯性这里就不展现了，感兴趣的同学可以照着运行即可得到可视化结果VizDimLoadings(pbmc, dims = 1:2, reduction = 'pca')DimPlot(pbmc, reduction = 'pca')#当然一个PC内包含的基因，我们可以用热图的形式来展现出来DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)#接下来我们既然已经对庞大的数据进行了降维（也就是聚堆）的形式，那么我们究竟要选择几个PC来代表这么庞大的数据呢，肯定不能都选，否则我的数据量还是这些，就没有我们前面一直渗透的缩小缩小再缩小的含义了pbmc <- JackStraw(pbmc, num.replicate = 100)pbmc <- ScoreJackStraw(pbmc, dims = 1:20)#上面两个代码是通过不同的方式来帮助我们选择PC的数目，并且分别都对应不同的可视化JackStrawPlot(pbmc, dims = 1:15)
晨曦解读
上面这个图我们主要是看右边的图例，也就是P值，P值在PC11之后有一个剧烈的降低，所以我们选择PC为11作为我们后续分析的一个参数
ElbowPlot(pbmc)晨曦解读
这个图的我们主要看拐点，我们可以很清楚的看到，这个图在PC为10的时候图形趋于平稳，所以这里我们选择PC为10作为我们后续分析的一个参数（这里其实我们也可以知道不同的判断方法选择的PC数是不一样的，而且理论上来说保存更多的PC可以保存更多的生物学差异，所以这里我们灵活选择即可，因为都不算错~）
当然这里Seurat官网鼓励我们使用不同数量的PC（10、15甚至50）来进行后续的分析，而且往往后续的分析不会有太大的差异~
#既然我们确定了我们的PC数，也就是我们用这些PC来指代我们的整体#为什么这里我们要确定降维的维数→可以简单理解为只有确定了这个维数，我们才可以知道后续要把哪些主成分进进行聚类#下面我们开始进行细胞聚类，这里是把我们看似一个整体的大量细胞分出亚群pbmc <- FindNeighbors(pbmc, dims = 1:10)#这里的参数10就是我们前面确定的维数pbmc <- FindClusters(pbmc, resolution = 0.5)#这里的0.5决定了细胞亚群的多少，Seurat官网推荐0.4-1.2都可以进行尝试
晨曦解读
这里的resolution参数的选择有一个小技巧，比如我们一开始设置0.5，然后我们进行后续的操作，当我们tSNE可视化细胞亚群结果后，发现有的细胞亚群太大，也就是细胞太多，这个时候我们可以回过头来把这个参数调大点，这样来回修改操作，最后得到我们满意的tSNE可视化结果~
#让我们看一下细胞亚群分类的结果head(Idents(pbmc), 5)## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 ## 2 3 2 1 ## AAACCGTGTATGCG-1 ## 6 ## Levels: 0 1 2 3 4 5 6 7 8晨曦解读
可以看到我们的细胞已经被编上了亚群的编号，并且一共分了8个亚群
#前面既然已经把亚群分好了，但是光秃秃的数据肯定是没有感召力的，我们就要把我们的亚群进行可视化#这里介绍两种方法， 一种是UMAP，一种是tSNE，方法没有好坏之分，文章中也可以都放上去，总的来说一句话，挑你看的顺眼而且结果好的就可以（Ps：一般来说UMAP运算速度快，tSNE比较常用）pbmc <- RunUMAP(pbmc, dims = 1:10)DimPlot(pbmc, reduction = 'umap')
#另一种可视化方法-tSNEpbmc <- RunTSNE(pbmc, dims = 1:10)DimPlot(pbmc, reduction = 'tsne', label = TRUE)晨曦解读
个人感觉，我会选择tSNE~
接下来我们既然已经把细胞亚群给分类了出来，下一步，我们自然要开始定义我们细胞亚群的属性也就是功能，也就是众多工具频出的一个步骤——细胞注释
#细胞注释我粗略的整理出来三种方法，可以供大家选择#1.首先是依赖标记基因的细胞注释即该标记基因可以代表细胞亚群的功能和属性#2.通过R包（SingleR包进行自动细胞注释）#3.人工注释（当然前提你需要知道你研究的组织总各种细胞亚群的marker基因是什么，然后通过数据库查找判断细胞类型，这个方法其实包含在第一种之中）cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)#这个比较是基于cluster2和其它所有亚群进行对比head(cluster2.markers, n = 5)## p_val avg_log2FC pct.1 pct.2 p_val_adj## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61
晨曦解读
运行上面的代码可以得到cluster2亚群的所有差异基因，这个差异基因的选择是基于cluster2亚群最低可以检测到的（min.pct参数）百分率（0.25）来设定的，运行后得到下面的表格
晨曦解读
我们来举个例子大家就知道这个表格怎么解读了，拿IL32基因，pvalue代表该细胞亚群的IL32这个基因同除该细胞亚群外的这个基因去做差异分析得到的P值，也就是说跟bulk-seq类似，小于0.05具有统计学意义，然后FC就是越高特异性越强，pct.1则是代表该基因在对应的cluster所有细胞中表达的比例，pct.2代表该基因在其它cluster所有细胞中表达的比例，这两个值的比例表达基因在细胞中表达的广度差异，就是说该亚群表达比例高，在别的亚群表达比例低，那么就是我的专属marker这个意思
然后在这个表格我们一般关心两个数值，一个是FC，一个就是adjp~
#接下来找一下cluster5和cluster0、1、2、3的markercluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)head(cluster5.markers, n = 5)## p_val avg_log2FC pct.1 pct.2 p_val_adj## FCGR3A 2.150929e-209 4.267579 0.975 0.039 2.949784e-205## IFITM3 6.103366e-199 3.877105 0.975 0.048 8.370156e-195## CFD 8.891428e-198 3.411039 0.938 0.037 1.219370e-193## CD68 2.374425e-194 3.014535 0.926 0.035 3.256286e-190## RP11-290F20.3 9.308287e-191 2.722684 0.840 0.016 1.276538e-186#寻找每一个cluster里与其它所有细胞相比之后的差异marker，并且每个亚群只展示FC最大的两位pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)## # A tibble: 18 x 7## # Groups: cluster [9]## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene ## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> ## 1 1.74e-109 1.07 0.897 0.593 2.39e-105 0 LDHB ## 2 1.17e- 83 1.33 0.435 0.108 1.60e- 79 0 CCR7 ## 3 0 5.57 0.996 0.215 0 1 S100A9 ## 4 0 5.48 0.975 0.121 0 1 S100A8 ## 5 7.99e- 87 1.28 0.981 0.644 1.10e- 82 2 LTB ## 6 2.61e- 59 1.24 0.424 0.111 3.58e- 55 2 AQP3 ## 7 0 4.31 0.936 0.041 0 3 CD79A ## 8 9.48e-271 3.59 0.622 0.022 1.30e-266 3 TCL1A ## 9 1.17e-178 2.97 0.957 0.241 1.60e-174 4 CCL5 ## 10 4.93e-169 3.01 0.595 0.056 6.76e-165 4 GZMK ## 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A ## 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1 ## 13 1.05e-265 4.89 0.986 0.071 1.44e-261 6 GZMB ## 14 6.82e-175 4.92 0.958 0.135 9.36e-171 6 GNLY ## 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A ## 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1## 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4 ## 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP晨曦解读
接下来我们既然得到了细胞亚群中的marker基因，那么下面我们就要把其进行可视化，毕竟可视化能够更加直观一点
VlnPlot(pbmc, features = c('MS4A1', 'CD79A'))FeaturePlot(pbmc, features = c('MS4A1', 'GNLY', 'CD3E', 'CD14', 'FCER1A', 'FCGR3A', 'LYZ', 'PPBP', 'CD8A'))top10 <- pbmc.markers %>%group_by(cluster) %>%top_n(n = 10, wt = avg_log2FC)DoHeatmap(pbmc, features = top10$gene) NoLegend()
晨曦解读
其实到了这一步我们的细胞亚群也找到了，marker基因也有了，我们接下来就要进行细胞注释了，其实大家可以多找找自己领域的scRNA-seq的文章往往可能会找到作者整理好的，marker基因对应细胞亚群的表格，比如下面这个表格
然后找到上面这个表格后，我们就可以根据我们细胞亚群的marker基因进行选择，然后把我们之前的细胞亚群的数字命名进行重命名（这里和大家说一下并不是说你的细胞marker是这个基因它的亚群就一定是这个，需要进行表达量、组织、大概身份猜测、多途径、文献等等才可以确定一个细胞亚群的身份，而且我们并不需要知道所有细胞亚群的身份，只需要知道我们后续感兴趣，能用到的细胞亚群足以）
new.cluster.ids <- c('Naive CD4 T', 'CD14 Mono', 'Memory CD4 T', 'B', 'CD8 T', 'FCGR3A Mono', 'NK', 'DC', 'Platelet')names(new.cluster.ids) <- levels(pbmc)pbmc <- RenameIdents(pbmc, new.cluster.ids)DimPlot(pbmc, reduction = 'umap', label = TRUE, pt.size = 0.5) NoLegend()晨曦解读
这样就构成了我们scRNA-seq文章中常看到的Figure1，然后保存一下我们的数据，以用来后续的分析
saveRDS(pbmc, file = '../output/pbmc3k_final.rds')那么晨曦的单细胞笔记（3）到这里就结束啦，截止到这一步，我们对我们的细胞进行了分群且同时进行了注释并且保存了文件以便后续的分析~
我是晨曦，我们下次再见~
回复晨曦03，即可获得今天的范例数据和代码哦~
欢迎大家关注解螺旋生信频道-挑圈联靠公号~