最近,bioconductor团队出版了一本电子书,其整合了其网站上关于单细胞的R包并制定了一套常规的分析流程包括分析,可视化,导入导出。不仅如此,前三章还分别教你如何下载使用R,使用bioconductor网站以及如何设计单细胞实验,对初学者很友好了,哪怕你对R语言一窍不懂,也能跟着走完流程。
能看这本书的都是对单细胞测序有所需求或这有这个意愿去学习相关知识的。这本书主要是整合目前常见的单细胞分析流程并尽可能详细的解释这些流程的每一个步骤,包括原理,所使用的工具以及给出栗子。所以可以根据自己的实际需求取选取合适自己的workflow与其中的步骤。
当然是与单细胞测序分析无关的其他分析,毕竟一口不能吃成大胖子。
因为这本书关于单细胞分析的软件都是用R写的。如果对R完全不了解,你也不知道做出来的到底什么鬼,所以还是简单的介绍一下必备的背景吧。(最有效率的办法还是跟着生信技能树的公众号去学习!)
P.S 最近bioconductor容易出现联网失败的问题,具体的解决方法见Jimmy老师的推文【紧急通知】下载R包却联网失败?初学者的痛
这里有个坑注意一下,这本书的R包要求bioconductor version至少要3.10以上才能用,所以可以一键升级所有R包
#####
rm(list = ls())
options(stringsAsFactors = F)
options(download.file.method = 'libcurl')
options(url.method='libcurl')
chooseBioCmirror()###记得选bioc镜像,清华or中科大
BiocManager::valid() ###测试R包是否过期
BiocManager::install(version = "3.10")###升级bioconductor版本
但是时间比较久,因为要升级全部的R包,建议中午睡个午觉,下午起来就搞定了。
SingleCellExperiment
SingleCellExperiment是一个能够完成常用的70种单细胞测序R包的数据交换的R包,因为它能存储一个sc-seq数据的所有数据在一个class里,可以放单个细胞数据,也可以放不同细胞间基因表达量数据,还可以放不同批次的细胞等等,一个数据就是一个assay。class中分门别类的放好,而且方便调用,一个函数就能去你想要得数据。举个生活中的例子,class就相当于支付宝中的余额宝,数据就相当于你的钱,随时可以从不同银行(不同R包)转入转出,调用的函数就等于你的密码。
本章用的6个功能
这里选用书中最简单的列子来说明,注意的是SingleCellExperiment
输入必须是矩阵,就像你存在支付宝的钱必须是人民币一样
rm(list = ls())
options(stringsAsFactors = F)
options(download.file.method = 'libcurl')
options(url.method='libcurl')
library(scater)
library(SingleCellExperiment)
####
###建立一个矩阵,行名是基因,列名是细胞,注意的是输入必须是矩阵格式
counts_matrix <- data.frame(cell_1 = rpois(10, 10),
cell_2 = rpois(10, 10),
cell_3 = rpois(10, 30))
rownames(counts_matrix) <- paste0("gene_", 1:10)
counts_matrix <- as.matrix(counts_matrix) # must be a matrix object!
head(counts_matrix)
# cell_1 cell_2 cell_3
# gene_1 10 9 29
# gene_2 8 13 36
# gene_3 4 9 36
# gene_4 6 11 42
# gene_5 10 8 26
# gene_6 4 9 35
###之后把矩阵放入class中,记得给list命名,这样方便提取
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))
sce
# class: SingleCellExperiment 是什么class
# dim: 10 3 维度
# metadata(0):
# assays(1): counts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10 行名
# rowData names(0): 基因的元数据的行名
# colnames(3): cell_1 cell_2 cell_3 列名
# colData names(0): 细胞的元数据的列名
# reducedDimNames(0):
# spikeNames(0):
# altExpNames(0):
上面就是存入,对应的就是取出
有2种方法,一种是完整函数,一种是更便捷的办法(类似于密码与指纹支付)
assay(sce, "counts")##取出assay中名为counts的矩阵
counts(sce)##取出assay中名为counts的矩阵
第二种方法名字一定要一致,例如”counts“只能取”counts“,就像别人的指纹不可能转走你的钱,对吧。
SingleCellExperiment
支持直接编辑数据,就像你不需要把余额宝的钱转支付宝,就可以在淘宝购物一样。但有一点不一样的是,不会覆盖你的原始数据,会作为一个新的assay存放在class中
sce <- scater::logNormCounts(sce)##对sce的数据进行log化和正态化
# class: SingleCellExperiment
# dim: 10 3
# metadata(0):
# assays(2): counts logcounts 可以看到这里多了个名字为logcounts的assay
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(0):
# reducedDimNames(0):
# spikeNames(0):
# altExpNames(0):
###提取log数据
logcounts(sce)
# cell_1 cell_2 cell_3
# gene_1 4.473275 4.146048 3.961101
# gene_2 4.167494 4.651269 4.254920
# gene_3 3.245625 4.146048 4.254920
# gene_4 3.778973 4.420662 4.466477
# gene_5 4.473275 3.986273 3.814209
# gene_6 3.245625 4.146048 4.216435
# gene_7 4.328471 3.986273 4.006919
# gene_8 3.778973 4.146048 3.591729
# gene_9 3.986273 3.806574 4.254920
# gene_10 4.836720 3.361811 4.136242
###可以看到数据已经log+正态化
元数据,作为一种注释的数据,注释试验批次,实验条件等等,类似于你的工资条,告诉你你的工资构成。注意元数据是以数据框的形式放入class里
添加细胞的元数据
cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)
cell_metadata
# batch
# cell_1 1
# cell_2 1
# cell_3 2
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
colData = cell_metadata) ##为counts数据附上注释信息
sce
# class: SingleCellExperiment
# dim: 10 3
# metadata(0):
# assays(1): counts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(1): batch ###细胞的元数据的列名为batch
# reducedDimNames(0):
# spikeNames(0):
# altExpNames(0):
##元数据也可以用2种方法取出,一种取出来的是数据框,另一种则是向量
colData(sce)
# DataFrame with 3 rows and 1 column
# batch
# <numeric>
# cell_1 1
# cell_2 1
# cell_3 2
sce$batch
# [1] 1 1 2
当然也可以用函数添加更多的元数据例如质控数据,总数等等。
sce <- scater::addPerCellQC(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
## batch sum detected percent_top_50 percent_top_100
## <numeric> <integer> <integer> <numeric> <numeric>
## cell_1 1 86 10 100 100
## cell_2 1 104 10 100 100
## cell_3 2 280 10 100 100
添加基因的元数据
跟上面大同小异,调用的是不同的函数而已
sce <- scater::addPerFeatureQC(sce)##添加基因质控元数据
rowData(sce)
# DataFrame with 3 rows and 9 columns
# batch sum detected percent_top_50 percent_top_100 percent_top_200
# <numeric> <integer> <integer> <numeric> <numeric> <numeric>
# cell_1 1 77 10 100 100 100
# cell_2 1 88 10 100 100 100
# cell_3 2 325 10 100 100 100
通过PCA t-sne等方法计算counts数据的降维后的数据,一般是二维或者三维。
###PCA
sce <- scater::logNormCounts(sce)
sce <- scater::runPCA(sce) ###通过PCA算法计算降维数据
sce
# class: SingleCellExperiment
# dim: 10 3
# metadata(0):
# assays(2): counts logcounts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(2): mean detected
# colnames(3): cell_1 cell_2 cell_3
# colData names(16): batch sum ... percent_top_500 total
# reducedDimNames(1): PCA 可以看到这里多了通过PCA算法计算的降维数据
# spikeNames(0):
# altExpNames(0):
reducedDim(sce, "PCA")
# PC1 PC2
# cell_1 -1.3487780 -0.1053741
# cell_2 0.8693045 -0.4941895
# cell_3 0.4794735 0.5995637
# attr(,"percentVar")
# [1] 82.02114 17.97886
###tsne
sce <- scater::runTSNE(sce, perplexity = 0.1)
# Perplexity should be lower than K!
sce
# class: SingleCellExperiment
# dim: 10 3
# metadata(0):
# assays(2): counts logcounts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(2): mean detected
# colnames(3): cell_1 cell_2 cell_3
# colData names(16): batch sum ... percent_top_500 total
# reducedDimNames(2): PCA TSNE 多了通过TSNE算法计算的降维数据
# spikeNames(0):
# altExpNames(0):
reducedDim(sce, "TSNE")
# [,1] [,2]
# cell_1 -5395.251 -1806.471
# cell_2 4278.781 -3769.183
# cell_3 1116.470 5575.653
除了自动添加,也可以手动添加降维数据
u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2)
###添加通过Uniform Manifold Approximation and Projection (UMAP) method 计算的降维数据
reducedDim(sce, "UMAP_uwot") <- u
reducedDims(sce) # Now stored in the object.
什么叫可选实验?书中给的定义是 data for a distinct set of features but the same set of samples/cells,就是相同细胞或者样本的不同的features,常见的有spike-in RNA
spike-in RNA
spike-in RNA:人工添加的定量的外源RNA(spike-in RNA),它是不存在第一种因素(真实生物差异)的,所以它们的变化可以直接反映技术噪音,然后将整体的内源基因平均表达量变化与spike-in进行拟合。大致了解一下,后面的步骤会仔细讲解。
###建立一个spike_se的class,只含有feature矩阵
spike_counts <- cbind(cell_1 = rpois(5, 10),
cell_2 = rpois(5, 10),
cell_3 = rpois(5, 30))
rownames(spike_counts) <- paste0("spike_", 1:5)
spike_se <- SummarizedExperiment(list(counts=spike_counts))
spike_se
# class: SummarizedExperiment
# dim: 5 3
# metadata(0):
# assays(1): counts
# rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(0):
###通过altExp函数将spike_se这个class放入sce这个class中去并命名为spike
altExp(sce, "spike") <- spike_se
altExps(sce)
# List of length 1
# names(1): spike 已经成功添加了名为spike的altExps
基本的class部分就含有以上6个部分,还有很多可选项可以根据自己的需求选择(https://osca.bioconductor.org/data-infrastructure.html)
单细胞测序
单细胞测序(英语:Single cell sequencing)采取优化的下一代DNA测序技术(NGS)检测单细胞的序列,可以获得特定微环境下的细胞序列差异以方便研究其功能差异等[1]。对个体细胞的DNA测序可以帮助我们了解例如在癌症中的小范围细胞的变异;对其进行RNA测序可以帮助我们了解和鉴别不同的细胞类型与其表现的基因,对研究发育生物学等有较大裨益。[2] 来源自wiki百科
个人理解:就是把原来的组织更近一步拆分成以细胞为单位,更加精准区别遗传信息的异质性,更有助于人们了解不同细胞间的差异。
更详细介绍可以见单细胞测序扫盲:是什么?为什么?怎么做? 以及公众号单细胞天地
常规的scRNA workflow可以分为一下6个步骤,分别为:实验设计,预处理,导入R,数据处理,统计分析与可视化
开始做单细胞测序之前,需要想好用什么用什么技术去测序,因为现有的测序技术实在是太多了,wiki介绍已经有100多种技术以及被发表(https://en.wikipedia.org/wiki/List_of_single_cell_omics_methods)
1586655684033
所有的scRNA-seq技术可以根据Transcript coverage大致分为以下2种:
表格来源自文章Single-Cell RNA-Seq Technologies and Related Computational Data Analysis
下面介绍一下这2种代表性技术
代表性seq分为scsmart-seq与10xgenome2种,表格来源于[文章](doi: 10.1093/bfgp/elx035 )
1586530477265
两者之间主要的一点不同mRNA逆转录成cDNA,smart用的是全长,而10x只通过3-tag的方法。而选取哪一种主要看你的目的:主要是做isoform or gene fusion,选全长;如果更关注转录组异质性,选10x(因为会有更多细胞)图片来源自文章Single-Cell RNA-Seq Technologies and Related Computational Data Analysis
不好理解的话可以看哔哩哔哩视频【陈巍学基因】视频56:10X分析单细胞表达;以及【陈巍学基因】视频11:单细胞mRNA测序
An external file that holds a picture, illustration, etc. Object name is fgene-10-00317-g001.jpg
根据书的推荐
CellRanger
软件的pipeline(10x公司的配套工具)去获取表达矩阵+ STAR 比对,具体实战例子可以见单细胞天地系列推文[https://mp.weixin.qq.com/s/fP8f4HboMM7m2Nd7AIljlg]后面几部是sc-seq的重头戏
因为不同原因(例如细胞损伤,PCR扩增效率等)造成低质量的库(含有较低的counts,较少的表达基因以及高比例的线粒体或者核RNA)会导致下游分析结果不准确:
所有筛选都是通过阈值来进行筛选
大致步骤:获取QC矩阵-根据自己的需求确定阈值-移除低质量数据
因为技术原因(例如PCR扩增效率等)造成的数据偏差会导致下游分析结果不准确,所以需要进行正态化校正
答案是不,正态化只能考虑技术偏好性,而批次效应还需要考虑不同批次间生物样本的差异
怎么正态化数据?
通过每个细胞的size factor来处理数据
size factor:每个细胞都有自己的偏好性(例如PCR扩增效率不同),通过该细胞的均值比而等比例缩放所有基因。
不同的工具有不同的计算size factor的方法:
最后通过scater包的logNormCounts来log+正态化数据
highly variable genes (HVGs)
还有一种是基于先验信息的方法(如已知细胞的亚型)。比如通过SCDE软件鉴定已知不同细胞亚型间的差异表达基因,然后再基于差异表达基因来聚类分析等。
许多scRNA-seq分析都通过细胞在多个基因中的表达量来比较细胞。例如,聚类通过计算跨基因的欧几里德距离来鉴定具有相似转录组表达丰度的细胞。在这些应用中,每个单独的基因代表数据的一个维度。也就是说,如果一个细胞只有2个基因,我们可以通过二维的图来区分,但是一般细胞中都含有多个基因,就会有多个维度,而人是无法看见3维以上的数据。
所以我们需要将多个特征压缩到一个维中,这样既可以节约工作量,又可以进行可视化。
虽然对于算法来说,推荐10-50个维度,但是不利于人类理解,所以一般都会降到2维或者3维进行可视化。对于降维,通常选用以下3种算法:
1586696993027
1586697021310
1586697058676
小结:t-sne与UMAP是2种主流降维方法,各有利弊。根据自己的需求做选择。
聚类不等同于细胞类型,是人类通过自身的经验来聚集的概念,所以实际上不存在“真正的聚类”这个说法。所以需要通过不同的聚类算法去探索最合适的分类。
在这里插入图片描述
clusterModularity()
函数,设定 as.ratio=TRUE
,会得到配对集群之间观测到的权重与预期权重之比。最后可视化结果得到log后不同cluster之间的观测到的权重与预期权重之比,颜色越深说明聚类效果越好1586697205540
在这里插入图片描述
1586698621562
1586699915066
在这里插入图片描述
1586700066037
1586700537027
1586701030900
实际上我们只需要单边t-test,因为只需要考虑上调的基因做为潜在的marker,通过scran包的findMarkers函数去寻找potential marker。另外,可以调整pval.type参数调整筛选的严格程度。
1586702333552
1586702730628
1586703276441
block
以及design
去除或者作为唯一因素的design矩阵(类似于DE分析的design矩阵)1586744581864
removeBatchEffect()
以及 sva包的 comBat()
,前提是假定细胞群体组成在批次中是已知或相同的。rescaleBatches()
来移除线性的批次效果,因为相当于对每个基因的对数表达值进行线性回归,并进行一些调整以提高性能和效率。对于每个基因,按比例缩小每个批次中的平均表达,直到等于所有批次中的最低均值。我们选择按比例缩小所有表达式的值,因为当批次位于均值方差趋势的不同位置时,这可以缓解方差差异。fastMNN()
来进行MNN校正