如何去学一个R包（上）

序言

FateID是用于定量单细胞转录组数据集中细胞命运偏倚的方法，所述数据集包含从共同祖先（Herman，Sagar和 Grün2018）产生的不同细胞类型。预期祖先群体是数据集的一部分，并且FateID算法被设计用于学习每个祖细胞对一个或多个替代终端命运的预先存在的偏差。该算法的策略是应用迭代随机森林分类（Breiman 2001），以便使用在先前迭代中被分类为训练集的细胞来量化越来越年幼的祖细胞中的命运偏倚。

https://cran.r-project.org/web/packages/FateID/vignettes/FateID.html

软件安装

FateID可以直接从CRAN安装：

1install.packages("FateID")


安装后，加载FateID包：

1library(FateID)


接下来就是应用实例

输入数据

输入数据只需要一个表达矩阵，基因为列，细胞为行。

示例数据可从FateID包中获得。

1data(intestine)


导入的数据集是Grun等人2016发表的文章中的数据：谱系追踪5天后，Lgr5谱系报告基因呈现阳性的小鼠肠上皮细胞的转录物计数矩阵。这些细胞是Lgr5阳性肠干细胞的培养5天后的细胞。

1x <- intestine$x
2head(x[,1:5])
3##                          I5d_3     I5d_4     I5d_6    I5d_8      I5d_9
4## 2200002D01Rik__chr7   3.751327 0.3486954 0.9355945 4.791691  2.4250484
5## 2210407C18Rik__chr11  0.100000 0.5983680 0.3076682 0.100000  1.6469830
6## 2310014L17Rik__chr7   0.100000 0.1000000 0.3076682 0.100000  0.1000000
7## Acsl5__chr19          1.132947 0.5983680 0.5161523 1.651566  0.8719749
8## Actb__chr5            4.813430 2.1173627 3.9683714 2.431937  4.7778317
9## Adh1__chr3           10.256359 7.1820558 0.3076682 7.179934 10.3793249


然后FateID需要对细胞进行分类，这可以通过任何聚类方法生成。例如，本包中的包含的RaceID3算法。注意，这里每个细胞对应的分类要以向量形式作为后续输入。

1y <- intestine$y
2head(y)
3##  I5d_3  I5d_4  I5d_6  I5d_8  I5d_9 I5d_10 
4##      2      3      4      2      5      4


聚类分析可以告知数据集中是否存在成熟细胞类型，其中不同谱系的细胞类型对应于不同的聚类（即不同的分区数）。在示例数据中，cluster6包含肠细胞，其由Alpi基因的高表达标记，而cluster9表示成熟的Paneth细胞（Defa24的高表达）和簇cluster13则是成熟的杯状细胞（Clca3的高表达）。其他稀有细胞类型仅以非常低的数量存在，因此被排除在分析之外。

作为FateID的进一步输入，分化轨迹的端点，即数据集中所有不同细胞谱系的最成熟阶段必须由表示分区中相应簇的整数向量y定义。

1tar <- c(6,9,13)


如果数据中没有具体的细胞分类，则FateID可以基于marker gene 进行分类。运行这步命令，需要提供marker gene ID的列表。该列表的每个组成部分可以包含一个或多个不同谱系的标记基因：

1FMarker <- list(c("Defa20__chr8","Defa24__chr8"), "Clca3__chr3", "Alpi__chr1")
2xf <- getPart(x,FMarker,fthr=NULL,n=5)
3head(xf$part)
4##  I5d_3  I5d_4  I5d_6  I5d_8  I5d_9 I5d_10 
5##      1      1      1      1      1      1
6head(xf$tar)
7## [1] 2 3 4
8tar <- xf$tar
9y <- xf$part


getPart函数提取表达最高谱系之一中表达markers的前n个细胞，并通过这些细胞的平均表达量定义了表达阈值。从而将markers表达高于该阈值的细胞划分为一类。或者可以通过fthr这个参数自定义表达阈值。得到的细胞分类是安照FMarker提供的顺序，从2开始标号。Cluster 1 则是剩下的markers的表达量没有达到阈值的细胞。

确定分类后，FateID可以对剩下的(total - n)细胞进行重新分类。

1rc <- reclassify(x, y, tar, clthr=.75, nbfactor=5, use.dist=FALSE, seed=12345, nbtree=NULL, q=0.9)
2y  <- rc$part


reclassify函数可以重新分类并替代之前结果。此步骤的目的是识别所有具有明显偏向的细胞。此步骤是可选的，但是建议执行，可以在特征选择之前或之后对表达数据执行该函数和其他函数。在样本数据中，x仅包含具有超过由RaceID3推断的变化基因。当然也可以将所有基因用作输入。

1v <- intestine$v
2rc <- reclassify(v, y, tar, clthr=.75, nbfactor=5, use.dist=FALSE, seed=12345, nbtree=NULL, q=0.9)
3y  <- rc$part


重分类函数还基于重要性采样执行特征选择，即保留具有大于给定类的重要性分布的q分位数的所有特征。函数返回简化表达式表，可以替换原始输入表达式数据框：

1x  <- rc$xf


如果输入数据未受任何其他特征选择方法的限制，建议使用此功能。

还可以利用差异基因表达分析来进行特征选择。

1x <- getFeat(v,y,tar,fpv=0.01)


getFeat函数将靶簇的细胞内的基因表达与所有剩余细胞的进行比较，并鉴定出在目标簇中p值低于fpv显着上调的基因。该函数返回一个简化的表达式数据框，可用于后续分析。通常，重新分类功能更适用于特征选择，因为它反映了用于随机森林分类的信息。

计算命运偏见

FateID的核心功能计算数据集中除了目标集群中的细胞外每个细胞的命运偏差。这些细胞被分配给代表相应目标簇的谱系，概率为1，并且在推断所有其他细胞的命运偏差期间该概率不会改变：

1tar <- c(6,9,13)
2x <- intestine$x
3y <- intestine$y
4fb  <- fateBias(x, y, tar, z=NULL, minnr=5, minnrh=10, adapt=TRUE, confidence=0.75, nbfactor=5, use.dist=FALSE, seed=12345, nbtree=NULL)
5## minnr: 5 
6## minnrh: 10 
7## test set size iteration 1 : 5 5 5 
8## randomforest iteration 1 of 15 cells
9………


fateBias是FateID的核心功能，以下就是对该函数的功能进行简单介绍。

fateBias函数的输入对象有特征基因的表达矩阵x，cluster分区向量y，目标clustertar。除此之外还有可选参数z，z用于识别在先前迭代中已被分类为目标簇之一的所有细胞紧邻区域中的非分类细胞的胞间距离矩阵。默认为z=1-cor(x)，但如果有更优的距离度量，则该参数可以用于提供距离矩阵。

FateID中由fateBias函数执行迭代计算。它以表示每个目标集群的一组细胞开始。对于每个目标群集，minnr提取与目标群集中的所有细胞中具有最短中值距离的相邻细胞群。所有目标簇的相邻细胞集会作为下一次迭代的测试集。因此minnr该参数控制算法的步长。在每次迭代中，minnr细胞乘以目标聚类的数量，并且可以在下一次迭代中对训练集做出贡献。

该迭代的训练集包括分配给一个目标簇的所有细胞，并且响应矢量由这些细胞的分区给出。测试集的分类基于随机森林投票完成：如果一个细胞在某个目标簇获得的投票显著比其他簇多，则将其分配给此目标簇并为下一次迭代的训练集做出贡献。没有显着命运偏差的细胞则不会纳入到训练集进行计算。但是，会记录并存储所有细胞的投票比例(可以解释作为命运概率)。另一个重要参数控制哪个cell对给定迭代的训练集有贡献。minnrh则是控制最多对几个目标群集的训练集有贡献。该参数控制考虑测试集分类的分化轨迹上的基因表达范围。如果将minnrh设置为Inf，则之前对某个目标分类具有显著fate bias的细胞都会对训练集有贡献。通常建议将表达范围设置为较小的值，以增加算法的特异性。但是，训练集应该足够大，以保证分类的确定性。minnrh参数应设置为20或更大的值，具体取决于数据集的大小和覆盖范围。如果输入数据集很大并且数据集中有大量细胞可用于覆盖整个分化轨迹的所有谱系，则可以增大该参数。如果一个细胞集中包含数百个细胞的大多数数据集，我们使用minnr=5和minnrh=20。

FateID还允许动态的测试集大小，其中minnr参数针对每个目标集群单独调整，基于前一次迭代中的分类成功。如果adapt=TRUE，在每次迭代中成功分类的细胞的数量是确定的，由confidence参数为目标簇给出的最低票数的细胞的数量会影响测试集中的细胞。然后通过在添加伪数量1之后将该数量除以所有目标聚类的最大值来导出权重。在下一次迭代中通过相应权重对每个聚类重新调整测试集大小。

如果本地邻域的先前分类成功率低，则这导致测试集大小减小，因此分类较慢。如果在不同的分支上存在很大差异，这种自适应方案很重要。如果在短的分支并触及naïve compartment，其分类成功性可以预期减少，则在未来迭代中该分支上的进展减慢，而在更成熟的分布密集的分支上计算仍然全速进行。

作为替代方法，FateID算法还可以基于到距离来提供分类。当use.dist设置为时TRUE，则距离矩阵z（或1-cor(x)）被解释为特征矩阵。其余参数是随机森林算法的控制参数，通常不必进行调整。

命运偏差的可视化

各种降维方法通常用于单细胞转录组分析，以便可视化细胞分群。FateID包通过主曲线计算维度减少数目，以实现命运偏差和伪时间排序的可视化。

1dr  <- compdr(x, z=NULL, m=c("tsne","cmd","dm","lle","umap"), k=c(2,3), lle.n=30, dm.sigma="local", dm.distance="euclidean", tsne.perplexity=30, seed=12345)
2## finding neighbours
3## calculating weights
4## computing coordinates
5## finding neighbours
6## calculating weights
7## computing coordinates


前两个输入参数与fateBias函数的输入参数相同。该参数k表示减少到的维度(通常二维或者三维)，但是也可以将维度设置在三维以上，然后可视化投影到维度子集之后的数据。其余参数是用于降维的各种算法的主要控制参数。该函数使用Rtsne包进行t-SNE映射计算（Maaten和Hinton 2008），使用stats的cmdscale函数进行经典多维缩放，使用lle包中的lle函数进行局部线性嵌入计算，用destiny进行diffusion maps（安格勒等人。2016），使用umap包计算Umaps 。为了加速计算，可以仅选择维度的子集作为输入参数m。

所有结果都可以通过plotFateMap函数绘制。

1plotFateMap(y,dr,k=2,m="tsne")


以三维方式绘制可打开交互式RGL设备，以允许旋转绘图并放大和缩小：

1plotFateMap(y,dr,k=3,m="tsne")


通过增加参数(t+簇号)提供目标簇的名称和fateBias的输出在降维中突出显示命运偏差

1plotFateMap(y,dr,k=2,m="tsne",fb=fb,g="t6")

在降维的空间中计算principal curves以近似分化轨迹 。

对于给定目标聚类的principal curves的计算，只用那些在目标聚类具有大于trthr的随机森林投票的一小部分的细胞。如果没有给出这个参数，那么只用对这个目标聚类有显着偏差(p<0.05)的细胞。如果参数prc为TRUE则principal curves绘制将在图中。

1pr <- plotFateMap(y,dr,k=2,m="tsne",trthr=.33,fb=fb,prc=TRUE)


如果已知所有细胞类型的共同上游祖先的最小节点簇，则该簇的数量也可以作为输入参数start给出，以便初始化principal curve 的计算，其中principal curve 从起始簇开始并以一个目标集群结束。为了强制主要曲线遍历表示相应目标聚类的细胞群，初始目标聚类内的细胞权重增加10倍，用于主曲线计算。

如果principal curves 已经算过，则plotFateMap的返回对象是两个对象的列表。第一个参数pr，是由principal.curve函数计算的每个目标集群的 principal curve 对象列表。第二个参数trc是对于每个包含所有具有显著命运偏差或比trthr大的细胞IDs的目标簇的 principal curve的cell IDs 的 list。cell IDs 按照它们从曲线的最短距离点获得的沿主曲线的位置排序。

还可以在降维图中突出显示指定的基因或者基因集。

1v <- intestine$v
2pr <-plotFateMap(y, dr, k=2, m="tsne", g=c("Defa20__chr8", "Defa24__chr8"), n="Defa", x=v)


最后，如果g="E"和函数fateBias返回的列表fb作为输入，则函数将基于不同目标聚类的概率绘制命运偏差的熵。命运偏差的熵水平将指示对应于更多个多能细胞状态的更大的命运偏倚熵的多能性水平。如果执行g="E"，函数将返回一个具有每个cell的命运偏差熵的向量：

1E <- plotFateMap(y,dr,k=2,m="tsne",g="E",fb=fb)


1head(E)
2##     I5d_3     I5d_4     I5d_6     I5d_8     I5d_9    I5d_10 
3## 0.3947993 0.5427533 0.4772478 0.6879617 0.6941222 0.6986183


也可以直接计算主曲线而无需调用plotFateMap函数：

1pr  <- prcurve(y,fb,dr,k=2,m="tsne",trthr=0.4,start=3)


此函数具有与函数相同的输入参数，plotFateMap并由plotFateMap函数调用。