肿瘤外显子实战之探索一个肝癌样品的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
总共是23个肿瘤外显子测序，根据我们在生信菜鸟团的肿瘤外显子系列教程可以拿到全部肿瘤相关突变，教程列表：
虽然上面的教程演示的并不是吴教授文章里面的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"，可能是不同工具方法带来的差异结果。
你可能会需要：广州专场（全年无休）GEO数据挖掘课，带你飞（1.11-1.12）和 生信入门课全国巡讲2019收官--长沙站