TCGA体细胞突变系列教程--胃癌

生信技能树

发布于 2018-12-07 14:37:08

9.6K0

发布于 2018-12-07 14:37:08

文章被收录于专栏：生信技能树

jimmy

生信技能树联盟创始人

有这个想法很久了，我教了很多人如何批量下载TCGA数据，以及分析各个癌症的somatic突变信息以及TMB，还有突变的特征频谱。

下载TCGA所有癌症的maf文件计算TMB

下载TCGA所有癌症的maf文件做signature分析

TCGA计划的4个找somatic mutation的软件使用体验

但是限于时间和知识背景，虽然代码方面问题不大，但是即使把所有的数据全部走一波我的流程我也可能会看不懂，所以希望某些背景相关朋友能根据我教的知识来分析一波。

肿瘤突变分析越来越火了，一起来学习一下吧。今天和大家一起探索TCGA数据中胃癌突变的情况。

今天的探索分为两个部分：

1.Mutation

1）数据下载

目前TCGA突变分析的数据vcf格式数据是受限的，所以我们这应用maf文件进行分析。具体下载方式。直接去TCGA官网下载数据也不难，都很容易，并且工具都在添加一些新的功能，比如最近添加的CNV的分析，举一反三的看ICGC的应用方式几乎和TCGA的应用方式是一样的。

2）工具

主要用的R包，maftools,deconstructSigs以及他们的依赖包，后面的代码会详细的解释。安装过程有时候会非常痛苦，但是一般都能解决，安装R包的终结方式应该是“耐心”。如果网络安装不太好的话建议大家先下载下来，应用本地安装，安装后加载包时可能出现缺少依赖包，这时再耐心的安装依赖包即可。应用这种方式，目前没有遇到怎么也安装不上的包（当然了版本本身不支持的情况除外）。

3）代码和结果

source("https://bioconductor.org/biocLite.R")
chooseBioCmirror() #选择一个就近的镜像
biocLite("maftools")#安装
(maftools：http://bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html)
library(maftools)
Mandatory fields(必须字段): Hugo_Symbol, Chromosome, Start_Position, End_Position, Reference_Allele, Tumor_Seq_Allele2, Variant_Classification, Variant_Type and Tumor_Sample_Barcode(样本名，此字段沟通样本的maf文件和临床信息的关键，前者).

laml = read.maf(maf="STAD.mutectAdjustBarcode.maf.txt",clinicalData="clinical.STAD.tsv")#read.maf()函数有两个最关键的参数maf,clinicalData,这个两个数据框只需要共同的Tumor_Sample_Barcode，这个一点用起来非常方便，很多帖子并没有提到此处，详细信息可参考maftools的官网。注意：TCGA直接下载的maf文件第16列即为样品名(例如：TCGA-FP-A4BE-11A-11D-A24F-08)，但是直接下载的临床数据的样本名(例：TCGA-FP-A4BE)是不同的，此处需要整理成一致后读入。

getFields(laml) #查看有哪些字段，TCGA下载的数据120个字段

getClinicalData(laml) #查看临床信息
getSampleSummary(laml) #查看每个样品发生突变的情况，此处就可以计算tumor mutation load,TML=Missense_Mutation/外显子数。

plotmafSummary(maf = laml, rmOutlier = TRUE, 
               addStat = 'median', dashboard = TRUE,
               titvRaw=FALSE)#绘制整体的突变情况
#waterfall plot
#We will draw oncoplots for top ten mutated genes.
oncoplot(maf = laml, top = 20, fontSize = 12)
#绘制前20个突变基因的瀑布图。oncoplot()中参数gene=c()可以指定基因名，绘制感兴趣的基因的瀑布图

此图绘制出了胃癌中突变前20的基因。其中可以看到TP53等著名基因的突变。渴望探索的小伙伴可以去TCGA的官网试试一样可以绘制出此图，点选即可。

maftools的最核心的功能可能就介绍完了，但是maftools的功能远不止这些，其实临床信息方面我们还没有进行探索，同时maftools还可以进行某个基因突变情况的生存分析，这个下次带大家进行。网上有很多好的帖子介绍，最好的帖子还是官方的介绍。

2. Mutation Signature

首先了解mutation signature的概念，现在很多生信文章中提到Signature 1 ，Signature 3等，看到这个有点困惑，通读了这篇文章Signatures of mutational processes in human cancer[https://www.nature.com/articles/nature12477]。

文章中提到Different mutational processes often generate different combinations of mutation types, termed ‘signatures’.

(不同的突变过程可以产生不同的突变类型的组合)，称为“特征”。

这个概念给我的感觉是在混乱中用数学统计分析查找规律，并应用规律。

那肿瘤的突变特征如何计算呢？

source('http://bioconductor.org/biocLite.R');
chooseBioCmirror()
install.packages('deconstructSigs')
# dependencies 'BSgenome', 'BSgenome.Hsapiens.UCSC.hg19' 
BiocInstaller::biocLite('BSgenome')
BiocInstaller::biocLite('BSgenome.Hsapiens.UCSC.hg19')
BiocInstaller::biocLite('BSgenome.Hsapiens.UCSC.hg38')
BiocInstaller::biocLite('BSgenome.Hsapiens.NCBI.GRCh38')
#先把这些包装上，后面几个包比较大，建议大家先下载然后本地安装。都可以安装成功的。（BSgenome.Hsapiens.UCSC.hg38 vs BSgenome.Hsapiens.NCBI.GRCh38  https://www.biostars.org/p/340852/）

####这段代码简单可以看出，NCBI的chromosome直接用的数字，而UCSC用的是chrX，这点应该看到。
genome.ncbi <- BSgenome.Hsapiens.NCBI.GRCh38
seqlengths(genome.ncbi)

genome.ucsc<-BSgenome.Hsapiens.UCSC.hg38
seqlengths(genome.ucsc)
####


suppressPackageStartupMessages(library("deconstructSigs"))
suppressPackageStartupMessages(library("BSgenome"))
library("BSgenome.Hsapiens.NCBI.GRCh38")#此处加载这个注释包的原因为TCGA的maf文件目前注释GRCh38
options(stringsAsFactors = F)
str(laml) #查看数据结构，从maf文件中拿到我们需要的数据

mut = laml@data
head(mut)
getField(laml@data)

a=mut[,c(16,5,6,11,13)]
colnames(a)=c( "Sample","chr", "pos","ref",  "alt")
head(a)

a$Sample=as.character(a$Sample)
a=a[nchar(a$ref)==1 & nchar(a$alt)==1,]
#NCBI.GRCh38染色体号无“chr”,所以要替换掉

a2=cbind(a[,-2],chr=gsub("chr","",a$chr))
dim(a2)
head(a2)
tail(a2)
sigs.input <- mut.to.sigs.input(mut.ref = a2, 
                                sample.id = "Sample", 
                                chr = "chr", 
                                pos = "pos", 
                                ref = "ref", 
                                alt = "alt",
                                bsg = BSgenome.Hsapiens.NCBI.GRCh38 #此处一定要注意到版本的问题否则结果可能是错的。
                                )
class(sigs.input)
head(t(sigs.input))

#下面的语句，批量计算每个样本的signature。核心的函数就是whichSignatures()
w=lapply(unique(a2$Sample), function(i){
  ## signatures.cosmic signatures.nature2013
  sample_1 = whichSignatures(tumor.ref = sigs.input, 
                             signatures.ref = signatures.cosmic, 
                             sample.id =  i, 
                             contexts.needed = TRUE,
                             tri.counts.method = 'default')
  return(sample_1$weights)
})


w=do.call(rbind,w)

plotSignatures(plot_example, sub = 'example')
library(pheatmap)
#绘制热图
pheatmap(t(w),cluster_rows = F,cluster_cols = T)
pheatmap(w,cluster_rows = T,cluster_cols = F)