1. TooManyCells 简介
- TooManyCells 是由宾夕法尼亚大学的 Gregory W. Schwartz 等人开发的一种聚类算法,开发的本意是用于 scRNA-seq 数据的分析,当然也可以用于任意的以 observations 为列,以 features 为行的数据
- 相关文献:TooManyCells identifies and visualizes relationships of single-cell clades 发表在 Nature Methods
- 据我了解,该算法有 2 种包装形式,一种是
too-many-cells
软件,一种是TooManyCellsR
R 包,由于 R 包报错问题尚未解决,所以本教程只涉及too-many-cells
软件 - 参考 github 文档
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 个文件这样格式的输入文件,但我手头只有普通表达矩阵,在向曾老师求助之后他很快发了推文:https://mp.weixin.qq.com/s/NaZ5kz3ew2O01cFEnK8sXg
-
从这段非常秀的代码中我学到这么几点
- 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 吧...
友情宣传
- 全国巡讲全球听(买一得五),第二期 ,你的生物信息学入门课
- 生信技能树的2019年终总结 ,你的生物信息学成长宝藏
- 2020学习主旋律,B站74小时免费教学视频为你领路
网友评论