肿瘤外显子实战之探索一个肝癌样品的23个部位的差异

刚才听了一个讲座，题目很吸引人，关于癌症研究和进化生物学的主题，主讲人是吳仲義教授，澳大濠江傑出訪問學者、芝加哥大學教授。2004年當選中央研究院院士及美國科學促進協會外籍會員。曾擔任芝加哥大學生態與進化系主任，在其三任主任期間帶領該系升至美國排名第一。2008年被任命為中科院北京基因組研究所（BIG）的所長，並領導BIG的重組工作直至2014年。吳教授是第一批千人計畫學者，領軍國家自然科學基金委的大型項目“微進化和多基因相互作用”的研究。

吴教授提到了他们的一个工作，让我很眼熟，才发现就是我讲解肿瘤外显子课程的时候，拿出去作为例子的那篇文章。

很久以前我就分享过这个文章，发在 Proc Natl Acad Sci U S A. 2015 Nov， 题目是：Extremely high genetic diversity in a single tumor points to prevalence of non-Darwinian cell evolution ， 研究者对一个肿瘤 hepatocellular carcinoma (

测序超过300个部位，其中23个进行WES测序，另外286个进行基因型鉴定，发现近 1 亿的编码区的突变（这个说的并不是somatic突变，注意区分哦），说明了肿瘤的突变并不符合达尔文进化模式。

肿瘤外显子分析流程

原始测序数据在， BioProject: PRJCA000091, 见教程：https://mp.weixin.qq.com/s/DqgQlBaSGt73e4SeEeTJoQ

虽然上面的教程演示的并不是吴教授文章里面的BioProject: PRJCA000091数据，但是整体的流程，大家是可以走一波的，拿到maf文件，记录23个进行WES测序样品的全部somatic突变后，根据文章感兴趣的基因绘制一个热图。

根据文献里面的基因提取突变频率表格

其中 2015-hcc.maf 文件需要自己走肿瘤外显子流程，至少耗时3~10天才能拿到，需要学习linux和WES课程：

然后focus_genes.txt的基因列表，看文章即可拿到。

library(reshape2)
library(pheatmap)

# read the maf and get the mutation freqence
mut <- read.delim('./2015-hcc.maf')
mut <- mut[c('Tumor_Sample_Barcode','Hugo_Symbol','t_depth','t_alt_count')]
mut$fre <- mut$t_alt_count/mut$t_depth

# dcast 整合成突变频率矩阵
mut <- dcast(mut, Tumor_Sample_Barcode~Hugo_Symbol,
             fun.aggregate = sum, fill = 0, drop = F,
             value.var = 'fre')
rownames(mut) <- mut$Tumor_Sample_Barcode;mut$Tumor_Sample_Barcode <- NULL

# subset 提取文章关注的 gene 突变频率
gene <- read.table('./focus_genes.txt',header = F,stringsAsFactors = F)
tmp.1 <- intersect(gene$V1,colnames(mut))
mydata <- subset(mut,select=tmp.1)

得到的mydata就可以绘制热图，是一个矩阵，所以很简单的绘制热图代码即可

# 热图
pheatmap(mydata,
         scale='none',cluster_rows = T,
         color=colorRampPalette(c('white','red','yellow'))(100),
         show_rownames = T,
         clustering_method = 'median',
         border = NA)

可以看到，跟文章的图表类似：

根据原文：

Only two clones, indicated by blue bars, are represented by more than one sample.

认为只有两个小克隆群，分别是 D16、 D29 以及 C31、C74、C3。

我们的结果里面，同样存在两个小克隆群，分别是 D16、 D29 以及 C31、C74。

当然 B4 和 B6 也有相近的聚类距离部分原文报道的突变，在我们这里没有表现出来 "LRP2" "ITGB2" "MAP4" "PRKAR1B" "A2M" "TMPRRSS13"，可能是不同工具方法带来的差异结果。