一文了解TooManyCells
1. TooManyCells 简介
TooManyCells 是由宾夕法尼亚大学的 Gregory W. Schwartz 等人开发的一种聚类算法,开发的本意是用于 scRNA-seq 数据的分析,当然也可以用于任意的以 observations 为列,以 features 为行的数据
相关文献:TooManyCells identifies and visualizes relationships of single-cell clades(https://www.nature.com/articles/s41592-020-0748-5?tdsourcetag=s_pctim_aiomsg) 发表在 Nature Methods
据我了解,该算法有 2 种包装形式,一种是
too-many-cells
软件,一种是TooManyCellsR
R 包,由于 R 包报错问题尚未解决,所以本教程只涉及too-many-cells
软件参考 github 文档
https://github.com/GregorySchwartz/tooManyCellsR?tdsourcetag=s_pctim_aiomsg
https://github.com/GregorySchwartz/too-many-cells
2. 软件安装
下面仅介绍一种安装方式,更多方式请见 github 文档
可以选择用 conda 安装 curl
用 curl 安装 stack
## 安装 stack
curl -sSL https://get.haskellstack.org/ | sh #这个 curl 可以用 conda 安装
stack setup
用 stack 安装 too-many-cells
## 安装 too-many-cells
git clone https://github.com/GregorySchwartz/too-many-cells.git
cd too-many-cells
stack install
完成后打开
nohup.out
看到最后提示将/blabla/.local/bin
路径加入 PATH 变量,因为这是软件的安装路径
3. 文件架构
个人习惯,project 下建三个文件夹:
scpt,用于存放脚本;
input,用于存放输入信息;
out,用于存放输出信息
4. 关于常规表达矩阵转 cellranger 结果
一开始由于我的疏忽,以为这个软件只能读取 cellranger 结果中的 3 个文件这样格式的输入文件,但我手头只有普通表达矩阵,在向曾老师求助之后他很快发了推文:表达矩阵逆转为10X的标准输出3个文件
从这段非常秀的代码中我学到这么几点:
R 中的“常规写入”
file="matrix.mtx"
sink(file)
cat("%%MatrixMarket matrix coordinate integer general\n")
cat("%\n")
sink()
操作循环的输出结果(如每次循环返回一个小数据框,最后
rbind
到一起)
tmp=do.call(rbind,lapply(1:ncol(ct),function(i){
return(data.frame(row=1:nrow(ct),
col=i,
exp=ct[,i]))
}))
write.csv
中的append
参数,避免覆盖原有内容
5.默认参数
但是,在我又读了一遍 github 文档之后,发现输入既可以是一个文件夹(里面放 cellranger 的 3 个文件),也可以是一个 csv 格式的普通表达矩阵…于是还是采用后者读取数据
除了表达矩阵之外,还需要一个输入文件
labels.csv
,大致长下面这个样子:如果已知细胞有不同的来源,或者数据分析之后对细胞有注释需求都可以通过这个输入文件实现
代码
too-many-cells make-tree \
--matrix-path ../input/expr_count.csv \
--labels-file ../input/labels.csv \
--draw-collection "PieRing" \
--output ../out \
> ../out/clusters.csv
对于我的需求来说最后的输出只有两个文件有用,一个是
clusters.csv
,记录聚类结果,一个是dendrogram.svg
,可视化聚类结果可视化效果:
6. “修剪”树枝
显然,默认参数下的分支太细了,我们可以通过两种方式来调整
直接设置 --min-size 参数为一个值,如 100,以规定最小分支细胞数
设置 --smart-cutoff 参数为一个值,如 4,以规定最小分支细胞数为 4*median absolute deviation,这个我还不太明白,但是我猜测这是根据每个分支的情况决定分支大小,可能会更合理
代码
too-many-cells make-tree \
--prior ../out \
--labels-file ../input/labels.csv \
--smart-cutoff 1 \ #经调试,我的数据最合适的值是1
--min-size 1 \
--draw-collection "PieChart" \
--output ../out_pruned \
> ../out_pruned/clusters_pruned.csv
可视化结果:
只有一种颜色是因为我的 labels 只标了一种,如果 labels 有多种,那么就会呈现这样的效果:
7. 取细胞子集
如果有取出一部分有一定特点的细胞进一步分析,就需要用到
clusters_pruned.csv
了首先查看这个文件的结构
$ head clusters_pruned.csv
cell,cluster,path
AAACGGGAGGTGTTAA.1,9,9/8/7/6/5/4/3/2/1/0
AACACGTTCGGCGGTT.1,9,9/8/7/6/5/4/3/2/1/0
AACCGCGGTATATGAG.1,9,9/8/7/6/5/4/3/2/1/0
ACACCCTTCTGGTTCC.1,9,9/8/7/6/5/4/3/2/1/0
ACCTTTAAGGTGTTAA.1,9,9/8/7/6/5/4/3/2/1/0
ACGAGGACACGTTGGC.1,9,9/8/7/6/5/4/3/2/1/0
AGGGAGTCAGGCTCAC.1,9,9/8/7/6/5/4/3/2/1/0
AGGGATGAGCGATAGC.1,9,9/8/7/6/5/4/3/2/1/0
AGTGGGAAGATGTAAC.1,9,9/8/7/6/5/4/3/2/1/0
看来是通过数字的形式记录 cluster 信息的,并且记录了从小到大的每一个分支,但是这个数字和图怎么对应呢?可以把数字标出来
too-many-cells make-tree \
--prior ../out \
--labels-file ../input/labels.csv \
--smart-cutoff 1 \
--min-size 1 \
--draw-collection "PieChart" \
--draw-node-number \ #只需多加这个参数
--output ../out_pruned \
> ../out_pruned/clusters_pruned.csv
结果:
那么接下来就可以取任意分支的所有细胞的 barcodes 并由此去除该细胞子集的表达信息了
8. 其他可视化选项
由于我没有这些需求,所以这里仅仅是搬运 github 文档
8.1. 分支末端画成饼图
代码
too-many-cells make-tree \
--prior out \
--labels-file labels.csv \
--smart-cutoff 4 \
--min-size 1 \
--draw-collection "PieChart" \
--output out_pruned \
> clusters_pruned.csv
可视化结果:
8.2. 调整树枝宽度
代码
too-many-cells make-tree \
--prior out \
--labels-file labels.csv \
--smart-cutoff 4 \
--min-size 1 \
--draw-collection "PieChart" \
--draw-max-node-size 40 \
--output out_pruned \
> clusters_pruned.csv
可视化结果:
8.3. 不按分支大小缩放
代码
too-many-cells make-tree \
--prior out \
--labels-file labels.csv \
--smart-cutoff 4 \
--min-size 1 \
--draw-collection "PieChart" \
--draw-max-node-size 40 \
--draw-no-scale-nodes \
--output out_pruned \
> clusters_pruned.csv
可视化结果:
9. 体现特定基因表达量
代码
too-many-cells make-tree \
--prior ../out \
--matrix-path ../input/expr_count.csv \
--labels-file ../input/labels.csv \
--smart-cutoff 1 \
--min-size 1 \
--feature-column 2 \
--draw-leaf "DrawItem (DrawThresholdContinuous [(\"gene1\", 0), (\"gene2\", 0)])" \
--draw-colors "[\"#e41a1c\", \"#377eb8\", \"#4daf4a\", \"#eaeaea\"]" \
--draw-scale-saturation 10 \ #如果不加这个参数,很可能表达量普遍较低以至于整张图没有颜色,至于这个值多少比较合适我还没有试过
--output ../out_gene_expression \
> ../out_gene_expression/clusters_pruned.csv
可视化结果:
10. 任意分支之间的差异分析
比如我想比较上图中两根蓝色(也就是 gene1 高表达,gene2 低表达)的分支,一看他们的 cluster 号码分别是 110 和 148(当然可以直接空缺一个数字,表示比较一个分支和其他所有细胞)
代码
too-many-cells differential \
--matrix-path ../input/expr_count.csv \
-n "([110], [148])" \
+RTS -N24
> ../out/differential.csv
11. diversity
比较两个细胞群的多样性(需要先跑过 make-tree 得到结果)
代码
too-many-cells diversity\
--priors ../out1 \
--priors ../out2 \
-o ../out_diversity_stats
12. 伪时序分析
代码
too-many-cells paths\
--prior ../out \
--labels-file ../input/labels.csv \
--bandwidth 3 \
-o ../out_paths
结果:
个人认为还是用 slingshot 吧…