变异数据处理的总体思路如下：
1.数据下载
TCGA的突变数据有4个软件得到的不同版本：
这个可以在gdc的官网上找到，case选择KIRC，文件类型选择maf即可获得。
选择mutect，就一个文件，直接点进去，download就行，下载下来只有一个tar.gz文件，解压放在工作目录下。
tar -xzvf file.tar.gz 解压tar.gz，即可得到一个maf.gz文件。
同样的筛选条件，参考https://www.jianshu.com/p/559d9604fcdf下载临床信息数据并整理。
2.数据读取
使用R包maftools读取。
rm(list=ls())
options(stringsAsFactors = F)
require(maftools)
require(dplyr)
project='TCGA_KIRC'
laml = read.maf(maf = 'TCGA.KIRC.mutect.2a8f2c83-8b5e-4987-8dbf-01f7ee24dc26.DR-10.0.somatic.maf.gz')
laml
maf_df = laml@data
save(laml,maf_df,file = "maf.Rdata")
length(unique(maf_df$Tumor_Sample_Barcode))
length(unique(maf_df$Hugo_Symbol))
因此，有336个病人,9444个突变基因信息。了解maf还可以用下面的几个函数：
getSampleSummary(laml)
getGeneSummary(laml)
getFields(laml)
3.突变数据的可视化
3.1 plotmafSummary
maftools 自带可视化函数plotmafSummary，可以比较直观的统计maf文件的数据。
#if (as.numeric(dev.cur()) != 1) graphics.off()
plotmafSummary(maf = laml, rmOutlier = TRUE,
#showBarcodes = FALSE,
addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
就是将maf_df 数据框做了统计，用barplot和boxplot做了可视化。
3.2 突变频谱图
代码其实就一句！
oncoplot(maf = laml, top = 30, fontSize = 1)
#上面这幅图是TOP30的基因，如果想要挑出几个特定的基因也是可以的
oncoplot(maf = laml, genes = unique(maf_df$Hugo_Symbol[50:55]), fontSize = 1)
我们可以看一下oncoplot这个绘图包的帮助文档里面的示例代码
image.png
可以看到最后加了临床信息，配色也有所改变
可以加上一些临床信息
# pd 是临床信息数据框，第一列是Tumor_Sample_Barcode，后面几列是各种分组，如gender，age等
load("clinical.Rdata")
pd = cl_df[,c(12,5,6)] #选列
#colnames(pd)[1] = "Tumor_Sample_Barcode"
#调整pd的病人id顺序和laml里的样本id顺序一致
sam_id = unique(laml@variants.per.sample$Tumor_Sample_Barcode)
library(stringr)
pd = pd[match(str_sub(sam_id,1,12),pd$bcr_patient_barcode),]
pd$Tumor_Sample_Barcode = sam_id
pd = pd[,c(4,2,3)]
pd = na.omit(pd)
laml = read.maf(maf = 'TCGA.KIRC.mutect.2a8f2c83-8b5e-4987-8dbf-01f7ee24dc26.DR-10.0.somatic.maf.gz',clinicalData = pd)
oncoplot(maf = laml,
sortByAnnotation = TRUE,
clinicalFeatures = c('vital_status',
'gender')
)
下面展开一下这个图的解读
主体热图
一行是一个基因，总共是9444个基因，从中截取了top30；一列是一个样本，总共是336个样本。不同颜色代表不同类型的突变。
右侧条形图
右侧的条形图是每个基因的突变样本数、突变类型和比例
验证一下突变样本数
count(maf_df,Hugo_Symbol,sort = T) #count出自dplyr包，表示统计数据框中某一列的重复值
结果显示VHL在169样本中突变，样本总数336，所以是49%，以此类推
条形图的颜色是突变类型，以VHL基因为例，他的突变类型分别是：
maf_df %>% filter(Hugo_Symbol=="VHL") %>%
count(Variant_Classification,sort = T)
顶部条形图
显示每个样本里突变的基因个数，可以看到最高的是那个一枝独秀的1600多。
laml@variants.per.sample %>% head()
#####################################分割线###############################
这个图怎么做出来？附上小洁老师的代码
0.R包
rm(list=ls())
options(stringsAsFactors = F)
project='TCGA_KIRC'
library(maftools)
library(dplyr)
library(pheatmap)
library(ComplexHeatmap)
library(stringr)
library(deconstructSigs)
library(BSgenome)
library(BSgenome.Hsapiens.UCSC.hg38)
1.读取突变数据maf
是和上一篇一样的数据，从gdc下载的maf文件和临床信息，见：
https://www.jianshu.com/p/2d6cb5bd8771
laml = read.maf(maf = "TCGA.KIRC.mutect.2a8f2c83-8b5e-4987-8dbf-01f7ee24dc26.DR-10.0.somatic.maf.gz")
laml
mut=laml@data
mut[1:4,1:2]
mut=mut[mut$Variant_Type=='SNP',] #把SNP的数据挑出来
#每行记录了一个突变，所以统计样本列的行数即为样本的突变数量
plot(table(mut[,16]),las = 2)
2.制作signatures的输入数据（96频谱）
关于什么是mutation
signature：http://www.bio-info-trainee.com/1619.html
关于96频谱：http://www.biotrainee.com/thread-832-1-1.html
a = mut[,c(16,5,6,12,13)]
colnames(a)=c( "Sample","chr", "pos","ref", "alt")
a$Sample=as.character(a$Sample)
sigs.input <- mut.to.sigs.input(mut.ref = a,
sample.id = "Sample",
chr = "chr",
pos = "pos",
ref = "ref",
alt = "alt",
bsg = BSgenome.Hsapiens.UCSC.hg38)
class(sigs.input)
#第一个样本的突变绘图
barplot(as.numeric(sigs.input[1,])~colnames(sigs.input))
sigs.input[1:4,1:4]
一顿操作猛如虎，生成signature与样本对应关系的矩阵
if(!file.exists(paste0(project,"w.Rdata"))){
w=lapply(unique(a$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')
print(c(i,which(unique(a$Sample)==i)))
return(sample_1$weights)
})
w=do.call(rbind,w)
save(w,file = paste0(project,"w.Rdata"))
}
load(paste0(project,"w.Rdata"))
mat = t(w)
3.可视化
矩阵的可视化--热图
3.1 简单版本的热图
pheatmap(mat,
cluster_rows = F,
cluster_cols = T,
show_colnames = F)
3.2 好看一点的热图
看到了这个图，氮素不知道用什么包画的，反正闲的没事干，不如自己模仿一波。如果你知道用的什么包，可以告诉我~感谢
热图上面放直方图，直方图内容是每个样本的突变数量。想到了ComplexHeatmap，可以将直方图作为注释添加上去。
难点是调整顺序，R语言神技能加持。
捋一捋：
直方图横纵坐标是样本和突变数量
热图输入数据mat横纵坐标是样本和signature
注释栏来源于clinical数据，与样本一一对应。
所以需要调整样本顺序，三者统一。
注释栏简化一下，用生死和性别；先排一下序，再按照排好的顺序给mat和直方图输入数据排序。
load("clinical.Rdata")
table( str_sub(colnames(mat),1,12)%in% cl_df$bcr_patient_barcode)
#一一对应，只要mat列名对应的临床数据。
cl_df = cl_df[cl_df$bcr_patient_barcode %in% str_sub(colnames(mat),1,12),]
#按照性别和生死排序
cl_df = arrange(cl_df,gender,vital_status)
#mat排序
x = match(cl_df$bcr_patient_barcode,str_sub(colnames(mat),1,12))
mat = mat[,x]
#直方图输入数据排序
for_bar = laml@variants.per.sample$Variants[match(cl_df$bcr_patient_barcode,str_sub(laml@variants.per.sample$Tumor_Sample_Barcode,1,12))]
#确认顺序是否正确
cl_df$bcr_patient_barcode[1] == str_sub(colnames(mat),1,12)[1]
cl_df$bcr_patient_barcode[1] == for_bar[1]
annotation = HeatmapAnnotation(mut = anno_barplot(for_bar),
df = data.frame(gender = cl_df$gender,
vital = cl_df$vital_status)
)
# top_annotation参数在顶部添加注释信息
ht_list = Heatmap(mat, name = "mat",
#top_annotation = column_ha,
#right_annotation = row_ha,
cluster_rows = F,
cluster_columns = F,
#show_row_names = F,
show_column_names = F,
top_annotation = annotation)
draw(ht_list, heatmap_legend_side = "left", annotation_legend_side = "bottom")
箱线图
0.输入数据
rm(list=ls())
load("for_boxplot.Rdata")
这里面的三个数据：
exp和meta是miRNA的表达矩阵和临床信息，由GDC下载整理得到。
mut是突变信息，读取maf得到的数据框。
1.比较任意miRNA在tumor和normal样本中的表达量
这个只需要表达矩阵，以hsa-mir-143为例画图，可替换为其他任意miRNA。
exp[1:4,1:4]
group_list=ifelse(as.numeric(substr(colnames(exp),14,15)) < 10,'tumor','normal')
table(group_list)
library(ggstatsplot)
dat = data.frame(gene = exp["hsa-mir-143",],
group = group_list)
ggbetweenstats(data = dat, x = group, y = gene,title = "hsa-mir-143")
这样就画出来了上面那个图
还有相关性的图：
呼。~~~
代码就不贴了，反正还没有搞懂。。。慢慢来吧！！！
文末要特别感谢数据挖掘班的小洁老师，以上所有代码都是来自小洁老师❤
