新版TCGA的突變SNP數(shù)據(jù)下載和整理

關(guān)于TCGAbiolinks包的學(xué)習(xí)前面一共介紹了5篇推文。

今天繼續(xù)學(xué)習(xí)如何使用TCGAbiolinks下載和整理MAF格式的突變數(shù)據(jù)。

之前的TCGA的MAF文件是可以下載的，每個(gè)癌癥包含4種軟件得到的突變文件：

之前能下載，現(xiàn)在不能了

后來就改版了，不讓你隨便下載了。但其實(shí)還是可以下載的，只不過沒有那么多選擇了！

現(xiàn)在的情況是每個(gè)樣本都是一個(gè)單獨(dú)的maf文件，需要下載后自己整理，就像整理表達(dá)矩陣那樣。

MAF文件的下載

但是現(xiàn)在我們有TCGAbiolinks，根本不需要自己動(dòng)手，直接三步走即可得到我們需要的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")

這樣得到的這個(gè)Rdata文件其實(shí)是一個(gè)數(shù)據(jù)框，不過由于內(nèi)容和之前的MAF文件一模一樣，所以也是可以直接用maftools讀取使用的。

maftools是一個(gè)非常強(qiáng)大的突變數(shù)據(jù)可視化和分析的R包，這個(gè)包在bioconductor上，需要的自行安裝。

無縫對(duì)接maftools

由于我們?cè)诘谝徊揭呀?jīng)下載過了，所以這里就不用下載了，直接加載保存好的數(shù)據(jù)。

我們以TCGA-COAD的數(shù)據(jù)作為演示。

library(maftools)

load(file = "./TCGA-SNP/TCGA-COAD_SNP.Rdata")

maf.coad <- data

簡(jiǎn)單看一下這個(gè)數(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

可以看到是一個(gè)data.frame類型的文件。

這個(gè)文件一共有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)

然后就是進(jìn)行各種可視化操作，毫無難度：

plotmafSummary(maf = maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE)
image.png

是不是非常簡(jiǎn)單，雖然沒有直接提供單個(gè)癌癥的MAF文件，但是使用TCGAbiolinks后，會(huì)直接幫我們整理好，沒有任何難度。

如果你由于各種原因不能使用這個(gè)包下載數(shù)據(jù)，那你可以直接用網(wǎng)頁下載，然后按照我之前的推文進(jìn)行整理：TCGA下載和表達(dá)矩陣整理：最適合初學(xué)者的教程 - 簡(jiǎn)書 (jianshu.com)

但是這個(gè)方法用在表達(dá)譜數(shù)據(jù)是沒有問題的，理論上用在其他類型的數(shù)據(jù)都是可以的，但是我并沒有嘗試過，歡迎大家使用后留言。

如果你在網(wǎng)絡(luò)上看見一個(gè)叫xxx.pl的文件，并且需要付費(fèi)獲取，建議你不要花這個(gè)冤枉錢，不值那個(gè)價(jià)，希望大家多多擦亮眼睛！

如果你非要用手撕代碼的方式自己整理，也是非常簡(jiǎn)單的，比整理轉(zhuǎn)錄組數(shù)據(jù)的表達(dá)矩陣簡(jiǎn)單100倍。

自己整理成MAF格式

首先你要去GDC TCGA的官網(wǎng)下載某個(gè)癌癥的所有的maf文件，還是以TCGA-COAD為例，下載好之后是這樣的：

image.png

每個(gè)樣本第一個(gè)文件夾，每個(gè)文件夾下面有一個(gè).gz結(jié)尾的壓縮文件，這個(gè)文件解壓縮之后就是大家熟悉的.maf文件大，但是只是一個(gè)樣本的~

image.png

把這個(gè).maf文件用VScode打開后是這樣的：

image.png

不妨多解壓幾個(gè)打開看一看，都是一樣的結(jié)構(gòu)，所以就很簡(jiǎn)單了，把所有的文件讀取進(jìn)來然后直接rbind()即可。

但是在此之前我們可以先讀取一個(gè)試試看：

# 路徑必須正確
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      +

非常順利，和上面那個(gè)整理好的格式一模一樣，唯一不同就是這個(gè)只是一個(gè)樣本的。

下面我們就批量讀取并合并就好了！

# 確定文件路徑！
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個(gè)
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é)！

竟然讀取不了，而且我們得到的這個(gè)maf.merge竟然只有137665行！和252664行的差距實(shí)在是太大了！

# 讀取失??！
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ù)是不行的！

然后只能回過頭去看哪里出了問題，通過仔細(xì)使用VScode直接打開maf文件和我們讀取的文件對(duì)比，發(fā)現(xiàn)了問題。

在Variant_Classification這一列中，有一些3'UTR / 5'UTR這樣的類型，但是在使用read.table()讀取的時(shí)候竟然識(shí)別不出來！

小丑竟是我自己！

image.png

所以直接導(dǎo)致遇到這個(gè)之后的所有行都是錯(cuò)位的，而且少了非常多行。

生氣??！

但是找到問題之后解決就非常簡(jiǎn)單，換個(gè)函數(shù)就行了，我們直接用data.table::fread()讀?。?/p> 

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列（少了一列是表示來自于第幾個(gè)樣本，沒有用）。

# 讀取成功！
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)
image.png

簡(jiǎn)單！下次說說這個(gè)maftools的使用。

覺得有用請(qǐng)多多轉(zhuǎn)發(fā)~拒絕不必要的花錢！難道免費(fèi)的不如付費(fèi)的香？？

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容