泛癌分析，依然是比较火热的话题，尤其结合多组学的方式来进行研究，成为了热点中的热点，跟着时代的潮流，今天小云想为大家分享一下如何进行泛癌中肿瘤与正常组织甲基化分析，接下来我们开始今天的学习。

1. 为何做泛癌中肿瘤与正常组织甲基化分析

启动子区DNA甲基化标记在基因表达的转录调控中有着很重要的地位，参与了很多生物学过程，因此通过甲基化分析结合转录表达，可以更好的来解释生物学意义，了解相关的调控机制。

今天小云通过数据下载，提取感兴趣基因集在肿瘤组织和正常组织中甲基化数据，最后绘制好看的气泡图，展示出不同肿瘤组织中感兴趣基因的甲基化状态，一文拿捏，泛癌甲基化分析，下面就跟着小云的节奏开始今天的学习之旅吧！

2. 准备需要的R包

#安装需要的R包

BiocManager::install("ChAMP")

BiocManager::install("ChAMPdata")

BiocManager::install("ComplexHeatmap")

install.packages("ggpubr", "randomcoloR")

install.packages("ggplot2")

install.packages("data.table")

BiocManager::install("impute")

#加载需要的R包

library(ggplot2)

library(ChAMP)

library(ChAMPdata)

library(data.table)

library(randomcoloR)

library(ggpubr)

library(impute)

library(ComplexHeatmap)

3. 数据准备

在进行泛癌甲基化分析之前，小伙伴最关心的问题是我们需要准备哪些文件，以及这些文件该如何下载，不用担心，小云帮你解决数据下载问题，其实很简单，我们只需要泛癌样本信息文件，甲基化探针注释文件，泛癌甲基化文件，接下来小云给大家最一下展示奥！

#样本信息文件来自该网站https://gdc.cancer.gov/about-data/publications/pancanatlas，可以直接下载。

merged_sample_quality_annotations.tsv，只需要关注前三列信息就可以,其他列可以忽略。

#甲基化探针注释文件，来自ChAMP R包内置数据,记得需要加载一下就可以呀。

data("probe.features")

#泛癌甲基化文件，第一列为探针名，其他列为样本名，需要泛癌甲基化文件的小伙伴可以联系小云(ps:此处只列举了一种肿瘤组织的甲基化文件)。

4.泛癌甲基化分析

在准备好了数据后，小云带着大家进行今天的重头戏-泛癌甲基化分析，主要包括提取感兴趣基因集的启动子探针信息，通过提取的启动子探针信息，分别提取不同肿瘤组织中感兴趣基因集的甲基化数据等步骤，通过该分析小伙伴可以掌握很多数据格式处理技巧和方法，值的仔细阅读，小云强烈推荐哈。

# 获得同时有肿瘤和正常样本的肿瘤名

tumors <- c("BLCA","BRCA","CESC","CHOL","COAD",

            "ESCA","GBM","HNSC","KICH","KIRC",

            "KIRP","LIHC","LUAD","LUSC","PAAD",

            "PRAD","READ","STAD","THCA","UCEC")

# 获得感兴趣的基因集

frg<- c("CDKN1A","HSPA5","TTC35","SLC7A11","NFE2L2","MT1G","HSPB1","GPX4","FANCD2","CISD1","FDFT1","SLC1A5","SAT1",    "TFRC","RPL8","NCOA4","LPCAT3","GLS2","DPP4","CS","CARS","ATP5G3","ALOX15","ACSL4","EMC2")

# 读取样本信息文件

rawAnno <- read.delim("merged_sample_quality_annotations.tsv",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T) # 数据来自PanCanAtlas

#提取1-15个字符串

rawAnno$simple_barcode <- substr(rawAnno$aliquot_barcode,1,15)

#去除重复的样本名

samAnno <- rawAnno[!duplicated(rawAnno$simple_barcode),c("cancer type", "simple_barcode")]

#去除cancer type 为空的数据

samAnno <- samAnno[which(samAnno$`cancer type` != ""),]

write.table(samAnno,"simple_sample_annotation.txt",sep = "\t",row.names = F,col.names = T,quote = F)

# 提取匹配到FRG基因的启动子探针

promoter <- probe.features[which(probe.features$feature %in% c("TSS1500","TSS200")),] # 根据注释文件筛选启动子探针

# 取出感兴趣基因在感兴趣位点的探针

promoter <- promoter[which(promoter$gene %in% frg),]

write.table(promoter, "promoter_annotation_for_interested_genes.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

# 获取只有感兴趣探针/基因的甲基化子集

for (i in tumors) { # 比较慢请耐心

  message("--",i,"...")

  # 加载甲基化RData文件,保存该文件夹的路径

load(file.path("/media/desk16/wangd/pp85/tcga_methy450",paste0("TCGA-",i,"_methy450.Rdata")))

  #将数据类型转换为数据框

  methy450 <- as.data.frame(methy450)

  #将第一列设为行名

  rownames(methy450) <- methy450[,1]

  # 去掉第一列探针名

  methy450 <- methy450[,-1]

  # 获取行名和列名

  dimname <- dimnames(methy450)

  ## 将数据全部转为数值以防报错

  methy450 <- sapply(methy450, as.numeric)

  # 重新赋值行名和列名

  dimnames(methy450) <- dimname

  #将数据类型转换为数据框

  methy450 <- as.data.frame(methy450)

  # 获取甲基化数据和感兴趣探针的交集

  compb <- intersect(rownames(methy450),rownames(promoter))

  methy450 <- methy450[compb,] # 取出子集

   # 将探针名和基因名进行匹配

  methy450$gene <- promoter[compb,"gene"]

  # 如果基因匹配多个探针，则取中位数beta值

  methy450 <- apply(methy450[,setdiff(colnames(methy450), "gene")], 2, function(x) tapply(x, INDEX=factor(methy450$gene), FUN=median, na.rm=TRUE))

  methy450 <- as.data.frame(methy450)

  write.table(methy450, paste0("TCGA_",i,"_methy450_subset.txt"),sep = "\t",row.names = T,col.names = NA,quote = F)

}

deltaMeth <- NULL

for (i in tumors) {

  message("--",i,"...")

  # 获取只包含感兴趣探针/基因的甲基化数据

  meth_subset <- read.table(paste0("TCGA_",i,"_methy450_subset.txt"),sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)

  #提取样本名1-15个字符串

  colnames(meth_subset) <- substr(colnames(meth_subset),1,15)

  #根据TCGA命名规律提取肿瘤和正常组织样本

  tumsam <- colnames(meth_subset)[substr(colnames(meth_subset),14,14) == "0"] # 获取肿瘤样本

  norsam <- colnames(meth_subset)[substr(colnames(meth_subset),14,14) == "1"] # 获取正常样本

  outTab <- NULL

  for (j in rownames(meth_subset)) {

    if(i == "KICH") { # KICH没有正常的甲基化样本

      delta <- 0 # 如果在分析KICH则delta设置为0

      wt <- 1 # 如果在分析KICH则设置pvalue为1

      outTab <- rbind.data.frame(outTab,

                                 data.frame(gene = j,

                                            tumor = i,

                                            Delta = delta,

                                            Pvalue = wt,

                                            stringsAsFactors = F),

                                 stringsAsFactors = F)

    } else {

      tmp1 <- as.numeric(meth_subset[j,tumsam]) # 获取肿瘤样本的beta值

      tmp2 <- as.numeric(meth_subset[j,norsam]) # 获取正常样本的beta值

      # wilcox检验

      wt <- wilcox.test(tmp1,tmp2)$p.value

      avgt <- mean(tmp1) # 肿瘤样本的平均beta值

      avgn <- mean(tmp2) # 正常样本的平均beta值

      delta <- avgt - avgn # delta值由肿瘤样本减去正常样本计算得到

      #合并所有肿瘤组织的甲基化分析结果文件

      outTab <- rbind.data.frame(outTab,

                                 data.frame(gene = j,

                                            tumor = i,

                                            Delta = delta,

                                            Pvalue = wt,

                                            stringsAsFactors = F),

                                 stringsAsFactors = F)

    }

  }

  deltaMeth <- rbind.data.frame(deltaMeth,

                                outTab,

                                stringsAsFactors = F)

}

write.table(deltaMeth, "TCGA_pancan_delta_meth_subset.txt",sep = "\t",row.names = F,col.names = T,quote = F)

5.开始绘图

# 设置颜色

blue <- "#4577FF"

red <- "#C2151A"

orange <- "#E45737"

green <- "#6F8B35"

darkblue <- "#303B7F"

darkred <- "#D51113"

yellow <- "#EECA1F"

# 产生泡泡图

#确定基因的绘图顺序

deltaMeth$gene<-factor(deltaMeth$gene,                   levels=rev(c("CDKN1A","HSPA5","TTC35","SLC7A11","NFE2L2","MT1G","HSPB1","GPX4","FANCD2","CISD1",                                       "FDFT1","SLC1A5","SAT1","TFRC","RPL8","NCOA4","LPCAT3","GLS2","DPP4","CS","CARS","ATP5G3","ALOX15","ACSL4")))

#设置渐变色

my_palette <- colorRampPalette(c(blue,"white",orange), alpha=TRUE)(n=128)

#ggplot开始绘图

ggplot(deltaMeth, aes(x=tumor,y=gene)) +

  geom_point(aes(size=-log10(Pvalue),color=Delta)) +

  scale_color_gradientn('Delta(T-N)',

                        colors=my_palette) +

  #修改主题

  theme_bw() +

  theme(#panel.grid.minor = element_blank(),

        #panel.grid.major = element_blank(),

        axis.text.x = element_text(angle = 45, size = 12, hjust = 0.3, vjust = 0.5, color = "black"),

        axis.text.y = element_text(size = 12, color = rep(c(red,blue),c(14,10))),

        axis.title = element_blank(),

        panel.border = element_rect(size = 0.7, linetype = "solid", colour = "black"),

        legend.position = "right",

        plot.margin = unit(c(1,1,1,1), "lines"))

#保存图片

ggsave("tumor_and_normal_promoter_methylation_of_interested_genes_in_pancancer.pdf", width = 8,height = 6)

6.结果文件

(1）、 promoter_annotation_for_interested_genes.txt

该文件为提取感兴趣基因集的启动子探针信息文件(TSS1500,TS200)。

（2）、TCGA_BLCA_methy450_subset.txt

该文件为特定肿瘤组织中感兴趣基因集的甲基化数据文件。

（3）、.TCGA_pancan_delta_meth_subset.txt

该文件为计算的感兴趣基因集在泛癌中的甲基化delta值和Pvalue值文件

（4）、.simple_sample_annotation.txt

该文件为简化的泛癌样本信息文件。

（5）、.tumor_and_normal_promoter_methylation_of_interested_genes_in_pancancer.pdf

该图为感兴趣基因集在泛癌中的甲基化状态点图，横坐标为肿瘤组织名称，纵坐标为感兴趣基因名，点的颜色来表示甲基化Delta值(肿瘤组织-正常组织)，点的大小表示Pvalue值。

最终，小云一文完成了泛癌中肿瘤和正常组织中甲基化分析，泛癌分析的其他相关内容可以尝试使用本公司新开发的云平台生物信息分析小工具，零代码完成分析，云平台网址：http://www.biocloudservice.com/home.html，主要包括泛癌中特定通路基因的突变分析(http://www.biocloudservice.com/740/740.php)，泛癌中拷贝数变异分析(http://www.biocloudservice.com/700/700.php)，基因在肿瘤与正常组织中差异分析(http://www.biocloudservice.com/710/710.php)等相关泛癌分析小工具；该分析在服务器运行速度会很快(ps:甲基化文件比较大)，小云为大家带来一个福利奥，购买代码可以免费试用本公司服务器三天，服务器租赁用户所有付费代码可以折扣获取，有需要的可以联系小云哈！