社区首页 >专栏 >细数绘制一张全景图所遇到的坑

细数绘制一张全景图所遇到的坑

生信技能树

发布于 2019-05-13 10:53:40

1.2K00

代码可运行

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

运行总次数：0

代码可运行

大家好，我是生信技能树学徒，前面我们带来了大量的表达数据挖掘实战演练，但是TCGA数据库之丰富程度，值得我们花费多年时间继续探索，现在带来的是突变全景图，如果你对之前的教程感兴趣，可以点击学习菜鸟团（周一数据挖掘专栏）成果展

就是上面这张全景，我重复出来的是下面这个样子。

数据准备

绘制全景图需要maf格式的突变信息文件以及临床信息文件。

还是从XENA上进行下载

需要注意，这里储存突变信息的文件需要是maf格式，和我们之前根据是否存在该基因的突变对样本进行分类的文件不同。

处理数据-R-maftools

1. 读取临床信息

tumor_type <- "LIHC"
Rdata_file <- paste('./data/', tumor_type, '.phenoData.Rdata', sep = '')
if (!file.exists( Rdata_file )) {
  phenoData <- read.table( destfile,
                           header = T,
                           sep = '\t',
                           quote = '' )
  rownames( phenoData ) <- phenoData[ , 1]
  colnames( phenoData )[1] <- "Tumor_Sample_Barcode"
  phenoData[1:5, 1:5]
  save( phenoData, file = Rdata_file )
}else{
  load( Rdata_file )
}

这里是遇到的第一个坑：我们看一下临床信息的“Tumor_Sample_Barcode”，是16位的短ID，但是后来在使用read.maf读取maf文件时，发现下载的maf文件的“Tumor_Sample_Barcode”是长ID，就存在了两个ID不匹配，从而导致临床信息被直接略过了。我去github上翻看了一下作者的代码，read.maf也可以接受数据框。所以就把maf文件先读取进来，处理一下ID。

2. 读取maf文件

maf <- data.table::as.data.table(read.csv(file = "./raw_data/TCGA.LIHC.mutect.DR-10.0.somatic.maf.gz", 
                                          header = TRUE, sep = '\t', 
                                          stringsAsFactors = FALSE, comment.char = "#"))
maf$Tumor_Sample_Barcode <- substr(maf$Tumor_Sample_Barcode, 1, 16)

require(maftools) 
## 作者用到了HBV和HCV的临床信息
phenoData$HBV <- ifelse(phenoData$hist_hepato_carc_fact == 'Hepatitis B', 'HBV', 'others')
phenoData$HCV <- ifelse(phenoData$hist_hepato_carc_fact == 'Hepatitis C', 'HCV', 'others')
phenoData[phenoData$neoplasm_histologic_grade == ""] <- 'no_reported'
## 这个函数不强求直接读取文本文件，也可以读取数据变量
laml <- read.maf(maf, clinicalData = phenoData)
laml 

laml@data <- laml@data[grepl('PASS', laml@data$FILTER), ]

接下来绘图遇到了第二个坑，关于factor的问题，以及颜色的对应关系的列表如何制作，绘图的函数怎么调用颜色信息。

3. 绘图

library(RColorBrewer)
png(paste0('oncoplot_top26_phone', tumor_type, '.png'), res = 150,
    width = 1500, height = 1080)

## 文章中这些driver gene是Mutsig挑选出来的，文章里面提供了，就直接使用了这个数据
genes = c("TP53", "CTNNB1", "ALB", "AXIN1", "BAP1", "KEAP1", "NFE2L2", "LZTR1", "RB1", "PIK3CA", "RPS6KA3", "AZIN1", "KRAS", "IL6ST", "RP1L1", "CDKN2A", "EEF1A1", "ARID2", "ARID1A", "GPATCH4", "ACVR2A", "APOB", "CREB3L3", "NRAS", "AHCTF1", "HIST1H1C")

## 为突变类型的分类数据设置颜色
variantClass <- names(table(laml@data$Variant_Classification))
col = c(RColorBrewer::brewer.pal(n = 4, name = 'Set1'),
        RColorBrewer::brewer.pal(n = 5, name = 'Set2'))
names(col) = variantClass
col

## 绘图的时候我们使用的数据是laml，临床信息在clinical.data里面
## 绘图函数要求这些设置颜色的数据是factor，所以我们要把加到图上的
## 临床信息转变为因子
laml@clinical.data$neoplasm_histologic_grade <- 
  as.factor(laml@clinical.data$neoplasm_histologic_grade)
gradecolors = RColorBrewer::brewer.pal(n = 4,name = 'Spectral')
names(gradecolors) = levels(laml@clinical.data$neoplasm_histologic_grade)

laml@clinical.data$race.demographic <- 
  as.factor(laml@clinical.data$race.demographic)
Racecolors = RColorBrewer::brewer.pal(n = 5,name = 'Spectral')
names(Racecolors) = levels(laml@clinical.data$race.demographic)

laml@clinical.data$gender.demographic <- 
  as.factor(laml@clinical.data$gender.demographic)
Gendercolors = c("#b3e2cd", "#fb9a99")
names(Gendercolors) = levels(laml@clinical.data$gender.demographic)

laml@clinical.data$HBV <- 
  as.factor(laml@clinical.data$HBV)
HBVcolors = c("#ffffb3", "#e31a1c")
names(HBVcolors) = levels(laml@clinical.data$HBV)

laml@clinical.data$HCV <- 
  as.factor(laml@clinical.data$HCV)
HCVcolors = c("#1b9e77", "#fc8d62")
names(HCVcolors) = levels(laml@clinical.data$HCV)

## 绘图函数需要一个list
phecolors = list(neoplasm_histologic_grade = gradecolors,
                 race.demographic = Racecolors,
                 gender.demographic = Gendercolors,
                 HBV = HBVcolors,
                 HCV = HCVcolors)
## clinicalFeatures是从laml@clinical.data里面挑取数据，所以
## 一定要是laml@clinical.data里面的列名
oncoplot(maf = laml, 
         colors = col, 
         bgCol = "#ebebeb", borderCol = "#ebebeb",
         genes = genes, GeneOrderSort = F, keepGeneOrder = T,
         fontSize = 7 , legendFontSize = 7,
         annotationFontSize = 7,
         annotationTitleFontSize = 7,
         sortByMutation = T,
         showTumorSampleBarcodes = F,
         annotationColor = phecolors,
         clinicalFeatures = c("neoplasm_histologic_grade",
                              "race.demographic",
                              "gender.demographic",
                              "HBV",
                              "HCV"))
dev.off()