专栏首页生信菜鸟团利用 Timescape 做肿瘤进化鱼图

利用 Timescape 做肿瘤进化鱼图

写在前面

前面我们使用 pyclone 分析了肿瘤样本的 clusters 结构,接下来我们进一步分析肿瘤进化,画一个鱼图,需要用到的工具是 citup 和 Timescape

参考:

  • https://github.com/sfu-compbio/citup
  • https://shahlab.ca/projects/citup/
  • https://www.cnblogs.com/xlij1205/p/10848917.html
  • https://bioconductor.org/packages/release/bioc/vignettes/timescape/inst/doc/timescape_vignette.html

citup 软件的安装和使用

CITUP 全称 Conality Inference in Tumors Using Phylogeny,是一种使用系统发育论进行肿瘤的克隆推断生物信息学工具,可用于从单个患者获得的多个样本来推断肿瘤异质性。

给定每个样本的突变频率,CITUP 使用基于优化的算法来找到最能解释数据的进化树。但是,前提条件是,CITUP 也是针对于深度测序数据的,这次勉强使用,为的是学习一下这个工具。

官网推荐我们使用 conda 来安装

# 添加channels
conda config --add channels http://conda.anaconda.org/dranew
conda config --add channels https://conda.anaconda.org/IBMDecisionOptimization/linux-64

# 创建小环境,然后安装,其实可以安装在前面的 pyclone 小环境里
conda create -n citup
conda activate citup
conda install -y citup h5py

# 安装完成后可以调用软件的帮助文档
run_citup_iter.py -h
run_citup_qip.py  -h

实际运行需要准备两个文件,一个是突变位点在不同样本的突变频率,行是突变位点,列是样本。另一个是突变的 cluster,只有单列,记录每个突变位点所在的 cluster 。两个文件的突变位点顺序要一致。具体可以参考:https://shahlab.ca/projects/citup/

在这里,我用前面 pyclone 结果进行处理:

for case in case{1..6}
do
cat ./9.pyclone/${case}_pyclone_analysis/tables/loci.tsv | cut -f 6 | sed '1d' | paste - - - -  >./9.pyclone/${case}_pyclone_analysis/freq.txt
cat ./9.pyclone/${case}_pyclone_analysis/tables/loci.tsv | cut -f 3 | sed '1d' | paste - - - - |cut -f 1 >./9.pyclone/${case}_pyclone_analysis/cluster.txt
## 获取 sample_id 信息,后面画图要用到
cat ./9.pyclone/${case}_pyclone_analysis/tables/loci.tsv | cut -f 2 | sed '1d' | head -4 >./9.pyclone/${case}_pyclone_analysis/sample_id
done

之后得到的文件结构就是:

$ ls ./9.pyclone/*
./9.pyclone/case1_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

./9.pyclone/case2_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

./9.pyclone/case3_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

./9.pyclone/case4_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

./9.pyclone/case5_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

./9.pyclone/case6_pyclone_analysis:
cluster.txt  config.yaml  freq.txt  plots  sample_id  tables  trace  yaml

每个样本的 3 个文件是(以 case1 为例):

$ cat ./9.pyclone/case1_pyclone_analysis/freq.txt 
0.33636363636363636    0.41732283464566927 0.35135135135135137 0.12162162162162163
0.2191780821917808    0.28440366972477066 0.2125  0.11904761904761904
0.31137724550898205    0.2519480519480519  0.2246153846153846  0.13624678663239073
......

$ cat ./9.pyclone/case1_pyclone_analysis/cluster.txt 
1
0
0
0
......

$ cat ./9.pyclone/case1_pyclone_analysis/sample_id
case1_techrep_2
case1_biorep_A_techrep
case1_biorep_B
case1_biorep_C

然后就是运行,输入文件是 freq.txtcluster.txt,结果输出到 reasults.h5

for case in case{1..6}
do
run_citup_qip.py --submit local \
./9.pyclone/${case}_pyclone_analysis/freq.txt \
./9.pyclone/${case}_pyclone_analysis/cluster.txt \
./9.pyclone/${case}_pyclone_analysis/results.h5
done

接下来需要用 python 处理一下结果,因为得到的结果results.h5是一个hdf5格式,不是普通的文本文件,处理起来比较麻烦一些,用 R 语言或者 python 都可以处理,这里我尝试用 python 处理一下(代码有点丑):

然后就是从上述结果中取最佳拟合树、进化树节点,克隆组成等信息,我把处理的过程写成一个 python 脚本 citup.py

import sys
import h5py
hf=h5py.File(sys.argv[1],'r')
hf.keys()
opnum=hf["results/optimal/index"][0]
cellfreq=hf["trees/" + str(opnum) + "/clone_freq/block0_values"][:]
tree=hf["trees/" + str(opnum) + "/adjacency_list/block0_values"][:]
print cellfreq
print tree

然后就是批量处理:

for case in case{1..6}
do
python citup.py ./9.pyclone/${case}_pyclone_analysis/results.h5 | sed 's/^ \[//;s/\[//g;s/\]//g' | tr ' ' '\t'| grep '\.' >./9.pyclone/${case}_pyclone_analysis/cellfreq.txt
python citup.py ./9.pyclone/${case}_pyclone_analysis/results.h5 | sed 's/^ \[//;s/\[//g;s/\]//g' | tr ' ' '\t'| grep -v '\.' >./9.pyclone/${case}_pyclone_analysis/tree.txt
done

得到的 cellfreq.txttree.txt 这两个文件,前者行表示样本,列表示克隆 clusters,每个值代表克隆频率。后者表示的是进化树的克隆分支关系

$ cat ./9.pyclone/case1_pyclone_analysis/cellfreq.txt 
0.35833     0.076023    0.235088    0.330559        
0.348454    0.0699815   0.235574    0.34599         
0.401151    0.101147    0.217818    0.279883    
0.751012    0.00999059  0.120444    0.118553

$ cat ./9.pyclone/case1_pyclone_analysis/tree.txt 
0    1
1    2
0    3

Timescape 可视化

得到上面的结果,就可以进一步处理成 Timescape 的 input 了,处理的方法是:

The adjacency list can be written as a TSV with the column names source, target to be input into E-scape, and the clone frequencies should be reshaped such that each row represents a clonal frequency in a specific sample for a specific clone, with the columns representing the time or space ID, the clone ID, and the clonal prevalence.

也就是 tree.txt 不用做修改,cellfreq.txt 需要被处理为3列 time or space ID, the clone ID, and the clonal prevalence,R代码部分如下:

BiocManager::install("timescape")
library(timescape)
options(stringsAsFactors = F)
example("timescape")
browseVignettes("timescape") 
library(plotly)
library(htmlwidgets)
library(webshot)
for (i in 1:6) {
  # i = 1
  # tree 
  tree_edges = read.table(paste0("9.pyclone/case", i, "_pyclone_analysis/tree.txt"))
  colnames(tree_edges) = c("source","target")
  # clonal prevalences
  cellfreq = read.table(paste0("9.pyclone/case", i, "_pyclone_analysis/cellfreq.txt"))
  colnames(cellfreq) = 0:(length(cellfreq)-1)
  sample_id = read.table(paste0("9.pyclone/case", i, "_pyclone_analysis/sample_id"))
  cellfreq$timepoint = sample_id[ , 1]
  library(tidyr)
  clonal_prev = gather(cellfreq, key="clone_id", value = "clonal_prev", -timepoint)
  clonal_prev = clonal_prev[order(clonal_prev$timepoint),]
  # targeted mutations
  # mutations <- read.csv(system.file("extdata", "AML_mutations.csv", package = "timescape"))
  p = timescape(clonal_prev = clonal_prev, tree_edges = tree_edges,height=260)
  saveWidget(p, paste0("case", i,"_timescape", ".html"))
  }

因为 Timescape 的输出格式比较少见,上面我只能保存为网页文件,原文件是可交互的,找不到合适的转为图片的方法,所以下面我就直接放了截图。

可以看到,每个样本的克隆进化关系如图所示,当然结果并不是很理想,毕竟数据不是很完美,不过这并不影响我们对工具的学习:


生信技能树目前已经公开了三个生信知识库,记得来关注哦~

每周文献分享

https://www.yuque.com/biotrainee/weeklypaper

肿瘤外显子分析指南

https://www.yuque.com/biotrainee/wes

生物统计从理论到实践

https://www.yuque.com/biotrainee/biostat

本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:Nickier

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-04-27

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 每周文献分享第一季完结撒花暨第65期分享

    2019年3月11日是一个特别的日子,那一天我们发布了 生信菜鸟团每周文献分享第1期,也意味着生信技能树语雀知识库平台正式上线。

    生信菜鸟团
  • 一行代码的优雅| Python列表生成式

    列表是Python中非常常见的数据结构,在基础课中也占了不小的篇幅。今天的推送就列表相关的内容再整理。

    生信菜鸟团
  • ggplot2快速入门

    ggplot2是一个做科研都会用到的R包。其实它的使用并不难,这次推文我将会使用ggplot2自带的测试数据,和大家快速了解,入门ggplot2。

    生信菜鸟团
  • 字符串匹配 - KMP算法

    实现 strStr() 函数。给定一个 haystack 字符串和一个 needle 字符串,在 haystack 字符串中找出 needle 字符串出现的第一...

    憧憬博客
  • Spark HA集群搭建

    比如分别把这两个文件重命名为start-spark-all.sh和stop-spark-all.sh 原因: 如果集群中也配置HADOOP_HOME,那么在...

    CoderJed
  • gcRMA算法-聚类分析图/PCA图

    >source("https://bioconductor.org/biocLite.R")

    黑妹的小屋
  • 一起来学matlab-matlab学习笔记8 基本绘图命令_1 图形窗口简介

    本文为matlab自学笔记的一部分,之所以学习matlab是因为其真的是人工智能无论是神经网络还是智能计算中日常使用的,非常重要的软件。也许最近其带来的一...

    DrawSky
  • 设计师如何用原型中钢笔工具快速画图?

    其实只要学会使用摹客原型设计的钢笔工具,结合形状合成功能(布尔运算),就能自由绘制你想要的形状,让你的设计更加得心应手。

    奔跑的小鹿
  • docker下的spark集群,调整参数榨干硬件

    本文是《docker下,极速搭建spark集群(含hdfs集群)》的续篇,前文将spark集群搭建成功并进行了简单的验证,但是存在以下几个小问题:

    程序员欣宸
  • MySQL数据以全量和增量方式,向ES搜索引擎同步流程

    /usr/local/logstash/sync-config/cicadaes.conf

    知了一笑

扫码关注云+社区

领取腾讯云代金券