花花写于2020-01-17 小年快乐呀大家！
今天上班最后一天，接下来开启长达半个月的年假（距离下一场讲课还有20天，期待）。我今年第一年在广东过年，不回家了，我父母会从山东老家过来和我一起过年，后天就到咯。。。你呢
思维导图走起啦1.数据下载
1.1 突变数据
TCGA的突变数据有4个软件得到的不同版本：
这个可以在gdc的官网上找到，case选择KIRC，文件类型选择maf即可获得。
选择mutect，就一个文件，直接点进去，download就行，下载下来只有一个tar.gz文件，解压放在工作目录下。
tar -xzvf file.tar.gz 解压，即可得到一个maf.gz文件。
同样的筛选条件，参考https://www.jianshu.com/p/559d9604fcdf下载临床信息数据并整理。
1mkdir kirc-clinical
2 ./gdc-client.exe download -m gdc_manifest.2020-01-17\ \(1\).txt -d kirc-clinical
2.数据读取
2.1 突变数据
使用maftools读取。
1rm(list=ls())
2options(stringsAsFactors = F)
3require(maftools)
4require(dplyr)
5project='TCGA_KIRC'
6laml = read.maf(maf = 'TCGA.KIRC.mutect.somatic.maf.gz')
7#> -Reading
8#> -Validating
9#> -Silent variants: 8383
10#> -Summarizing
11#> --Mutiple centers found
12#> BCM;BI--Possible FLAGS among top ten genes:
13#> TTN
14#> MUC16
15#> HMCN1
16#> -Processing clinical data
17#> --Missing clinical data
18#> -Finished in 3.750s elapsed (3.430s cpu)
19laml
20#> An object of class MAF
21#> ID summary Mean Median
22#> 1: NCBI_Build GRCh38 NA NA
23#> 2: Center BCM;BI NA NA
24#> 3: Samples 336 NA NA
25#> 4: nGenes 9444 NA NA
26#> 5: Frame_Shift_Del 1732 5.155 4
27#> 6: Frame_Shift_Ins 1201 3.574 1
28#> 7: In_Frame_Del 238 0.708 0
29#> 8: In_Frame_Ins 350 1.042 0
30#> 9: Missense_Mutation 12997 38.682 36
31#> 10: Nonsense_Mutation 1259 3.747 2
32#> 11: Nonstop_Mutation 18 0.054 0
33#> 12: Splice_Site 490 1.458 1
34#> 13: Translation_Start_Site 25 0.074 0
35#> 14: total 18310 54.494 47
36maf_df = laml@data
37save(laml,maf_df,file = 'maf.Rdata')
38length(unique(maf_df$Tumor_Sample_Barcode))
39#> [1] 336
40length(unique(maf_df$Hugo_Symbol))
41#> [1] 9444
因此，有336个病人,9444个突变基因信息。了解maf还可以用下面的几个函数：
1getSampleSummary(laml)
2getGeneSummary(laml)
3getFields(laml)
2.2.临床信息
将下载好的临床信息xml文件整理成一个数据框。
1xmls = dir('kirc-clinical/',pattern = '*.xml$',recursive = T)
2library(XML)
3td = function(x){
4 result <- xmlParse(file.path('kirc-clinical/',x))
5 rootnode <- xmlRoot(result)
6 xmldataframe <- xmlToDataFrame(rootnode[2])
7 return(t(xmldataframe))
8}
9
10cl = lapply(xmls,td)
11cl_df <- as.data.frame(t(do.call(cbind,cl)))
12cl_df[1:3,1:3]
13#> additional_studies tumor_tissue_site histological_type
14#> 1 Kidney Kidney Clear Cell Renal Carcinoma
15#> 2 Kidney Kidney Clear Cell Renal Carcinoma
16#> 3 Kidney Kidney Clear Cell Renal Carcinoma
17save(cl_df,file = 'clinical.Rdata')
3.突变数据的可视化
3.1 plotmafSummary
maftools 自带可视化函数plotmafSummary，可以比较直观的统计maf文件的数据。
1dev.off()
2#> null device
3#> 1
4plotmafSummary(maf = laml, rmOutlier = TRUE,showBarcodes = FALSE,
5 addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
就是将maf_df 数据框做了统计，用barplot和boxplot做了可视化。
3.2 突变频谱图
代码其实就一句！
1oncoplot(maf = laml, top = 30, fontSize = 1)
下面展开一下这个图的解读
主体热图
一行是一个基因，总共是9444个基因，从中截取了top30；一列是一个样本，总共是336个样本。不同颜色代表不同类型的突变。
右侧条形图
右侧的条形图是每个基因的突变样本数、突变类型和比例
验证一下突变样本数
1count(maf_df,Hugo_Symbol,sort = T)
2#> # A tibble: 9,444 x 2
3#> Hugo_Symbol n
4#> <chr> <int>
5#> 1 VHL 169
6#> 2 PBRM1 148
7#> 3 TTN 77
8#> 4 SETD2 46
9#> 5 BAP1 37
10#> 6 MUC16 28
11#> 7 MTOR 23
12#> 8 KDM5C 21
13#> 9 HMCN1 20
14#> 10 ATM 19
15#> # … with 9,434 more rows
结果显示VHL在169样本中突变，样本总数336，所以是49%，以此类推
条形图的颜色是突变类型，以VHL基因为例，他的突变类型分别是：
1maf_df %>% filter(Hugo_Symbol=='VHL') %>%
2 count(Variant_Classification,sort = T)
3#> # A tibble: 7 x 2
4#> Variant_Classification n
5#> <fct> <int>
6#> 1 Missense_Mutation 60
7#> 2 Frame_Shift_Del 41
8#> 3 Nonsense_Mutation 27
9#> 4 Frame_Shift_Ins 22
10#> 5 Splice_Site 16
11#> 6 In_Frame_Del 2
12#> 7 Nonstop_Mutation 1
顶部条形图
显示每个样本里突变的基因个数，可以看到最高的是那个一枝独秀的1600多。
1laml@variants.per.sample %>% head()
2#> Tumor_Sample_Barcode Variants
3#> 1: TCGA-B8-4143-01A-01D-1806-10 1611
4#> 2: TCGA-B0-5098-01A-01D-1421-08 550
5#> 3: TCGA-A3-A8OV-01A-11D-A36X-10 120
6#> 4: TCGA-CJ-4920-01A-01D-1429-08 117
7#> 5: TCGA-CZ-5468-01A-01D-1501-10 102
8#> 6: TCGA-A3-A6NN-01A-12D-A33K-10 97
