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

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 )
}

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), ]

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()