TCGAbiolinks可以進(jìn)行甲基化分析，但是功能不如ChAMP強(qiáng)大，甲基化分析還是首推ChAMP包。

不過為了了解TCGAbiolinks包，里面關(guān)于甲基化分析的部分還是要學(xué)習(xí)一下。

主要是甲基化差異分析，甲基化的一些可視化，甲基化和轉(zhuǎn)錄組數(shù)據(jù)的聯(lián)合作圖。

加載數(shù)據(jù)

我們還是使用之前下載好的TCGA-COAD的甲基化β值矩陣。

數(shù)據(jù)下載見這篇：使用TCGAbiolinks批量下載最新版TCGA數(shù)據(jù)庫的各種組學(xué)數(shù)據(jù)！

library(TCGAbiolinks)
suppressMessages(library(SummarizedExperiment))
load(file = "G:/tcga/TCGA-dnaMethy/COAD_METHY_beta.Rdata")

查看一下數(shù)據(jù)，現(xiàn)在這個(gè)β值矩陣還是一個(gè)SummarizedExperiment對(duì)象：

beta.m <- data
beta.m
## class: RangedSummarizedExperiment 
## dim: 485577 352 
## metadata(1): data_release
## assays(1): ''
## rownames(485577): cg13869341 cg14008030 ... cg11478607 cg08417382
## rowData names(52): address_A address_B ... MASK_extBase MASK_general
## colnames(352): TCGA-D5-6530-01A-11D-1721-05
##   TCGA-AA-3660-11A-01D-1721-05 ... TCGA-A6-5664-01A-21D-1837-05
##   TCGA-D5-6533-01A-11D-1721-05
## colData names(107): barcode patient ... paper_vascular_invasion_present
##   paper_vital_status

甲基化差異分析

β矩陣不能含有缺失值，所以先去除缺失值，使用缺失值插補(bǔ)方法進(jìn)行插補(bǔ)也是可以的：

beta.m <- subset(beta.m, subset = (rowSums(is.na(assay(beta.m))) == 0))

如果你的甲基化矩陣是直接使用TCGAbiolinks包整理好的SummarizedExperiment對(duì)象，那么這個(gè)甲基化差異分析就非常簡(jiǎn)單。

我們需要確定誰和誰進(jìn)行相比，也就是要?jiǎng)?chuàng)建一個(gè)含有分組信息的列。

所有樣本的信息都在SummarizedExperiment對(duì)象中的colData中，包括分組信息、生存信息、分期信息等，我們這里只要用到分組信息即可，所以先把colData簡(jiǎn)化一下。

如果你這里都看不懂，可以翻看之前的推文：

TCGA轉(zhuǎn)錄組數(shù)據(jù)的表達(dá)矩陣提取 
使用TCGAbiolinks包進(jìn)行差異分析

# 只要兩列信息，其實(shí)只要sample_type一列也行
sample_info <- as.data.frame(colData(beta.m))[,c("barcode","sample_type")]

# sample_type改成normal和tumor
sample_info[sample_info == "Solid Tissue Normal"] <- "normal"
sample_info[sample_info != "normal"] <- "tumor"

table(sample_info$sample_type)
## 
## normal  tumor 
##     38    314

# 重新變成coldata
colData(beta.m) <- DataFrame(sample_info)

有了SummarizedExperiment對(duì)象和相應(yīng)的分組信息就可以進(jìn)行甲基化差異分析了。

# 非常慢
res <- TCGAanalyze_DMC(data = beta.m,
                       groupCol = "sample_type", # colData中這一列含有分組信息
                       group1 = "normal", # normal組
                       group2 = "tumor" # tumor組
                       )

## Group1:normal
## Group2:tumor
## Calculating the p-values of each probe...
## |=================================================================|100%                       Completed after 20 m 
## Saving volcano plot as: methylation_volcano.pdf

查看結(jié)果:

head(res)
##            mean.normal mean.tumor mean.normal.minus.mean.tumor
## cg21870274   0.7868787  0.5479701                   0.23890856
## cg16619049   0.2781233  0.2591982                   0.01892506
## cg18147296   0.7219909  0.6847735                   0.03721734
## cg13938959   0.6958773  0.4552416                   0.24063575
## cg12445832   0.4654454  0.2479300                   0.21751535
## cg23999112   0.7941180  0.5310289                   0.26308907
##            p.value.normal.tumor p.value.adj.normal.tumor
## cg21870274         9.732669e-18             2.668610e-16
## cg16619049         2.314279e-02             4.306785e-02
## cg18147296         1.024600e-01             1.566558e-01
## cg13938959         4.389698e-18             1.321596e-16
## cg12445832         2.608577e-18             8.391252e-17
## cg23999112         4.649043e-13             5.043637e-12
##                               status
## cg21870274 Hypermethylated in normal
## cg16619049           Not Significant
## cg18147296           Not Significant
## cg13938959 Hypermethylated in normal
## cg12445832 Hypermethylated in normal
## cg23999112 Hypermethylated in normal

可以看到和轉(zhuǎn)錄組的差異分析結(jié)果差不多，但是都是探針信息，基因信息需要注釋才行。

同時(shí)也會(huì)在當(dāng)前目錄下生產(chǎn)一個(gè)PDF格式的火山圖：

保存結(jié)果：

save(res, file = "tcga-coad_de_methy.rdata")

甲基化可視化

使用箱線圖可視化不同組別之間的甲基化值。

mdm <- TCGAvisualize_meanMethylation(data = beta.m,
                                     groupCol = "sample_type",
                                     print.pvalue = T
                                     )

會(huì)在目錄下生成一個(gè)PDF文件：

甲基化旭日?qǐng)D

可以在一張圖中查看轉(zhuǎn)錄組數(shù)據(jù)和甲基化數(shù)據(jù)的情況。

需要準(zhǔn)備轉(zhuǎn)錄組差異分析結(jié)果和甲基化差異分析結(jié)果。

轉(zhuǎn)錄組差異分析結(jié)果使用之前推文中得到的數(shù)據(jù)：xxxxxxxxxx

# 這個(gè)數(shù)據(jù)只是部分差異基因，你可以用所有基因的結(jié)果
load(file = "coadDEGs.Rdata")

甲基化差異分析結(jié)果就用本次的結(jié)果。

有了兩個(gè)差異分析的結(jié)果，就可以畫旭日?qǐng)D了：

starburst <- TCGAvisualize_starburst(
    met = res, 
    exp = coadDEGs,
    group1 = "normal",
    group2 = "tumor",
    met.platform = "Illumina Human Methylation 450",
    genome = "hg38",
    met.p.cut = 10^-5, 
    exp.p.cut = 10^-5,
    names = TRUE
)

會(huì)在當(dāng)前目錄下生成一個(gè)png格式的旭日?qǐng)D：