關(guān)于TCGAbiolinks包的學(xué)習(xí)前面一共介紹了5篇推文躬贡。
今天繼續(xù)學(xué)習(xí)如何使用TCGAbiolinks下載和整理MAF格式的突變數(shù)據(jù)。
之前的TCGA的MAF文件是可以下載的露懒，每個癌癥包含4種軟件得到的突變文件：
后來就改版了，不讓你隨便下載了升筏。但其實還是可以下載的涯雅，只不過沒有那么多選擇了鲜结！
現(xiàn)在的情況是每個樣本都是一個單獨的maf文件，需要下載后自己整理活逆，就像整理表達矩陣那樣精刷。
MAF文件的下載
但是現(xiàn)在我們有TCGAbiolinks，根本不需要自己動手蔗候，直接三步走即可得到我們需要的MAF文件怒允。
library(TCGAbiolinks)
query <- GDCquery(
project = "TCGA-COAD",
data.category = "Simple Nucleotide Variation",
data.type = "Masked Somatic Mutation",
access = "open"
)
GDCdownload(query)
GDCprepare(query, save = T,save.filename = "TCGA-COAD_SNP.Rdata")
這樣得到的這個Rdata文件其實是一個數(shù)據(jù)框，不過由于內(nèi)容和之前的MAF文件一模一樣锈遥，所以也是可以直接用maftools讀取使用的纫事。
maftools是一個非常強大的突變數(shù)據(jù)可視化和分析的R包勘畔，這個包在bioconductor上，需要的自行安裝丽惶。
無縫對接maftools
由于我們在第一步已經(jīng)下載過了炫七，所以這里就不用下載了，直接加載保存好的數(shù)據(jù)蚊夫。
我們以TCGA-COAD的數(shù)據(jù)作為演示诉字。
library(maftools)
load(file = "./TCGA-SNP/TCGA-COAD_SNP.Rdata")
maf.coad <- data
簡單看一下這個數(shù)據(jù)：
class(maf.coad)
## [1] "data.frame"
dim(maf.coad)
## [1] 252664 141
maf.coad[1:10,1:10]
## X1 Hugo_Symbol Entrez_Gene_Id Center NCBI_Build Chromosome Start_Position
## 1 1 AGRN 375790 BCM GRCh38 chr1 1046481
## 2 1 ACAP3 116983 BCM GRCh38 chr1 1295539
## 3 1 CALML6 163688 BCM GRCh38 chr1 1916980
## 4 1 PRKCZ 5590 BCM GRCh38 chr1 2150972
## 5 1 WRAP73 49856 BCM GRCh38 chr1 3635995
## 6 1 CHD5 26038 BCM GRCh38 chr1 6142440
## 7 1 CAMTA1 23261 BCM GRCh38 chr1 7663513
## 8 1 ERRFI1 54206 BCM GRCh38 chr1 8014193
## 9 1 SLC2A7 155184 BCM GRCh38 chr1 9022922
## 10 1 PGD 5226 BCM GRCh38 chr1 10411462
## End_Position Strand Variant_Classification
## 1 1046481 + Frame_Shift_Del
## 2 1295539 + Missense_Mutation
## 3 1916980 + Silent
## 4 2150972 + Silent
## 5 3635995 + Silent
## 6 6142440 + Missense_Mutation
## 7 7663513 + Silent
## 8 8014193 + Missense_Mutation
## 9 9022922 + Missense_Mutation
## 10 10411462 + Silent
可以看到是一個data.frame類型的文件。
這個文件一共有252664行知纷，141列壤圃，包含了gene symbol，突變類型琅轧，突變位置伍绳，導(dǎo)致的氨基酸變化等信息。
下面就直接用read.maf()函數(shù)讀取即可乍桂，沒有任何花里胡哨的操作冲杀！
maf <- read.maf(maf.coad)
## -Validating
## -Silent variants: 63597
## -Summarizing
## --Mutiple centers found
## BCM;WUGSC;BCM;WUGSC;BCM;BI--Possible FLAGS among top ten genes:
## TTN
## SYNE1
## MUC16
## -Processing clinical data
## --Missing clinical data
## -Finished in 6.000s elapsed (5.690s cpu)
然后就是進行各種可視化操作，毫無難度：
plotmafSummary(maf = maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE)
是不是非常簡單睹酌，雖然沒有直接提供單個癌癥的MAF文件权谁，但是使用TCGAbiolinks后，會直接幫我們整理好憋沿，沒有任何難度旺芽。
如果你由于各種原因不能使用這個包下載數(shù)據(jù)，那你可以直接用網(wǎng)頁下載辐啄，然后按照我之前的推文進行整理：TCGA下載和表達矩陣整理：最適合初學(xué)者的教程 - 簡書 (jianshu.com)
但是這個方法用在表達譜數(shù)據(jù)是沒有問題的采章，理論上用在其他類型的數(shù)據(jù)都是可以的，但是我并沒有嘗試過壶辜，歡迎大家使用后留言悯舟。
如果你在網(wǎng)絡(luò)上看見一個叫xxx.pl的文件，并且需要付費獲取砸民，建議你不要花這個冤枉錢抵怎，不值那個價，希望大家多多擦亮眼睛岭参！
如果你非要用手撕代碼的方式自己整理便贵，也是非常簡單的，比整理轉(zhuǎn)錄組數(shù)據(jù)的表達矩陣簡單100倍冗荸。
自己整理成MAF格式
首先你要去GDC TCGA的官網(wǎng)下載某個癌癥的所有的maf文件，還是以TCGA-COAD為例利耍，下載好之后是這樣的：
每個樣本第一個文件夾蚌本，每個文件夾下面有一個.gz結(jié)尾的壓縮文件盔粹，這個文件解壓縮之后就是大家熟悉的.maf文件大，但是只是一個樣本的~
把這個.maf文件用VScode打開后是這樣的：
不妨多解壓幾個打開看一看程癌，都是一樣的結(jié)構(gòu)舷嗡，所以就很簡單了，把所有的文件讀取進來然后直接rbind()即可嵌莉。
但是在此之前我們可以先讀取一個試試看：
# 路徑必須正確
tmp <- read.table("G:/tcga/GDCdata/TCGA-COAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation/007c2ae4-bbd2-42c6-ab67-bf016fbddb51/982004b5-52e1-4a69-97d3-25bdcb77b026.wxs.aliquot_ensemble_masked.maf.gz",
skip = 7, # 前面7行都不要
sep = "\t", # 必須指定
header = T # 有行名
)
tmp[1:10,1:8]
## Hugo_Symbol Entrez_Gene_Id Center NCBI_Build Chromosome Start_Position
## 1 NOC2L 26155 BCM GRCh38 chr1 945088
## 2 SDF4 51150 BCM GRCh38 chr1 1217753
## 3 B3GALT6 126792 BCM GRCh38 chr1 1233752
## 4 MIB2 142678 BCM GRCh38 chr1 1629520
## 5 NADK 65220 BCM GRCh38 chr1 1765325
## 6 GNB1 2782 BCM GRCh38 chr1 1804562
## 7 PANK4 55229 BCM GRCh38 chr1 2510063
## 8 PRXL2B 127281 BCM GRCh38 chr1 2588386
## 9 PRDM16 63976 BCM GRCh38 chr1 3412316
## 10 WRAP73 49856 BCM GRCh38 chr1 3633442
## End_Position Strand
## 1 945088 +
## 2 1217753 +
## 3 1233752 +
## 4 1629520 +
## 5 1765325 +
## 6 1804562 +
## 7 2510063 +
## 8 2588386 +
## 9 3412316 +
## 10 3633442 +
非常順利进萄，和上面那個整理好的格式一模一樣，唯一不同就是這個只是一個樣本的锐峭。
下面我們就批量讀取并合并就好了中鼠！
# 確定文件路徑！
dir.path <- "G:/tcga/GDCdata/TCGA-COAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation"
# 獲取所有maf文件路徑
all.maf <- list.files(path = dir.path, pattern = ".gz",
full.names = T, recursive = T)
# 看看前3個
all.maf[1:3]
## [1] "G:/tcga/GDCdata/TCGA-COAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation/007c2ae4-bbd2-42c6-ab67-bf016fbddb51/982004b5-52e1-4a69-97d3-25bdcb77b026.wxs.aliquot_ensemble_masked.maf.gz"
## [2] "G:/tcga/GDCdata/TCGA-COAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation/010f9040-294d-4d14-a2b4-80d7a11625dd/5083b949-1bf3-4bc2-bf4f-f668f8a13792.wxs.aliquot_ensemble_masked.maf.gz"
## [3] "G:/tcga/GDCdata/TCGA-COAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation/0148fff1-b8af-4bf0-8bcd-de1ff9f750f3/2c16cfe2-bf6d-4a39-af3a-9dfd5ada3e17.wxs.aliquot_ensemble_masked.maf.gz"
然后直接讀取就行了沿癞，覺得慢可以用data.table::fread()加快速度援雇。
maf.list <- lapply(all.maf, read.table,
skip = 7,
sep = "\t",
header = T)
然后直接合并即可，如果不放心可以看看列數(shù)列名是不是一樣椎扬，100%一樣惫搏，我們就不看了。
# lapply(maf.list, dim)
maf.merge <- do.call(rbind,maf.list)
目前為止看似一切順利蚕涤，本以為即將結(jié)束筐赔，但是沒想到橫生枝節(jié)！
竟然讀取不了揖铜，而且我們得到的這個maf.merge竟然只有137665行茴丰！和252664行的差距實在是太大了！
# 讀取失斅弧较沪！
maf1 <- read.maf(maf.merge)
## -Validating
## --Removed 5 duplicated variants
## --Non MAF specific values in Variant_Classification column:
## 3UTR DEL T T - novel TCGA-A6-6781-01A-22D-1924-10 TCGA-A6-6781-10A-01D-1924-10
果然不檢查數(shù)據(jù)是不行的！
然后只能回過頭去看哪里出了問題失仁，通過仔細使用VScode直接打開maf文件和我們讀取的文件對比尸曼，發(fā)現(xiàn)了問題。
在Variant_Classification這一列中萄焦，有一些3'UTR / 5'UTR這樣的類型控轿，但是在使用read.table()讀取的時候竟然識別不出來！
小丑竟是我自己拂封！
所以直接導(dǎo)致遇到這個之后的所有行都是錯位的茬射，而且少了非常多行。
生氣懊扒在抛！
但是找到問題之后解決就非常簡單，換個函數(shù)就行了萧恕，我們直接用data.table::fread()讀雀账蟆肠阱！
maf.list <- lapply(all.maf, data.table::fread,
sep = "\t",
header = T,
skip = 7
)
maf.merge <- do.call(rbind,maf.list)
dim(maf.merge)
## [1] 252664 140
現(xiàn)在就和前面的數(shù)據(jù)一模一樣了，252664行朴读，140列（少了一列是表示來自于第幾個樣本屹徘，沒有用）。
# 讀取成功衅金！
maf1 <- read.maf(maf.merge)
## -Validating
## -Silent variants: 63597
## -Summarizing
## --Mutiple centers found
## BCM;WUGSC;BCM;BI;BCM;WUGSC--Possible FLAGS among top ten genes:
## TTN
## SYNE1
## MUC16
## -Processing clinical data
## --Missing clinical data
## -Finished in 3.970s elapsed (3.730s cpu)
plotmafSummary(maf = maf1, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE)
簡單噪伊！下次說說這個maftools的使用。
覺得有用請多多轉(zhuǎn)發(fā)~拒絕不必要的花錢氮唯！難道免費的不如付費的香鉴吹？？
