否则今天的文章你可能会有一种感觉：
本文主要解决的问题
1、下载目标基因/探针甲基化数据-Xena 文中抛出的问题
2、在线TCGA基因表达与甲基化相关性分析-cBioPortal 文中抛出的问题
3、使用R语言完成基因表达与甲基化相关性分析
首先解决第一个问题，在Xena中实现基因表达和甲基化的相关性分析：
不清楚如何到这步的请看友情提示...显示如下：
从这个界面中提取几个重点：
1、因为Xena无法获得基因水平甲基化，故无法直接完成基因水平的甲基化与表达相关性分析；
2、故，只能是基因上某个探针（X轴）与目标基因表达（Y轴）相关性分析；
3、探针甲基化定量方式是beta值（TCGA level3），基因表达为log2(norm_count+1)，即log2(RSEM+1)
4、此处得到的相关系数为负值，意味着CFTR基因表达与cg17616554探针甲基化beta值为负相关，pearson相关系数为-0.1076....
不知道最近几期持续关注生信控发布的甲基化分析文章的小伙伴，有没有突然想到一个事情~
提示一下：cBioPortal、负相关性最强、探针
还没来感觉的
一定是调皮的没好好看友情提示~
现在我们来解决第二、三个问题。由上文可知，我们要用到来自Firehose的探针甲基化数据，文件比较大(6.9GB+)，故可以使用linux系统从中提取目标基因CFTR所有探针的甲基化beta值：
le LIHC.methylation__humanmethylation450__jhu_usc_edu__Level_3__within_bioassay_data_set_function__data.data.txt | grep 'Hybridization\|Composite\|CFTR' > firehose.LIHC.CFTR.beta.txt
根据提取出的beta值可知，cBioPortal最终挑选用于代表基因的探针是cg22533025，故可以直接使用该探针甲基化与CFTR基因表达值进行相关性分析：
p_load(stringr)
#### firehose
data_firehose = read.table('firehose.LIHC.CFTR.beta.txt', h = F, sep = '\t', row.names = 1)
data_firehose = data_firehose[,which(data_firehose[2,] == 'Beta_value')]
colnames(data_firehose) = str_match(as.character(as.matrix(data_firehose[1,])), '(TCGA-[^-]*-[^-]*-\\d*)')[,2]
data_firehose = data_firehose[3:nrow(data_firehose),]
data_firehose[1:6,1:6]
#### cbioportal
data_cbio = read.table('plot_cBio.txt', h = T, sep = '\t', fill = T)
rownames(data_cbio) = as.character(data_cbio[,1])
#### firehose的探针甲基化 + cbioportal的基因表达
common = intersect(colnames(data_firehose), rownames(data_cbio))
rna_cbio = data_cbio[common, 3]
methy_firehose = as.matrix(data_firehose[complete.cases(data_firehose), common])
res_pea = cor.test(as.numeric(methy_firehose['cg22533025',]), rna_cbio)
res_spe = cor.test(as.numeric(methy_firehose['cg22533025',]), rna_cbio, method = 'spearman')
res = c(res_pea$p.value, res_pea$estimate, res_spe$p.value, res_spe$estimate)
names(res) = c('pearson_p.value', 'pearson_cor', 'spearman_p.value', 'spearman_cor')
与cBioPortal在线分析的结果一致：
发现问题没有？
cBioPortal挑选的负相关性最强的探针cg22533025，与基因表达的pearson相关系数为-0.08，而前文我们用Xena测试的cg17616554探针pearson相关系数为-0.1076....
是不是意味着cBioPortal的结果...有问题？
CFTR所有探针分析一遍：
tmp = t(sapply(1:nrow(methy_firehose), function(index){
res_pe = cor.test(as.numeric(methy_firehose[index,]), rna_cbio)
res_sp = cor.test(as.numeric(methy_firehose[index,]), rna_cbio, method = 'spearman')
c(res_pe$p.value, res_pe$estimate, res_sp$p.value, res_sp$estimate)
}))
rownames(tmp) = rownames(methy_firehose)
colnames(tmp) = c('pearson_p.value', 'pearson_cor', 'spearman_p.value', 'spearman_cor')
可见，无论是基于Spearman或是Pearson相关系数，cg22533025在所有探针的分析结果中均不是cBioPortal所谓的与基因表达负相关性最强：
For genes with multiple probes (usually from the Infinium arrays), we only include methylation data from the probe with the strongest negative correlation between the methylation signal and the gene's expression in the study.
当然，你可能觉得我们常用log2(RSEM+1)转化的基因表达量，毕竟cBioPortal中也提供了这个选项，会不会是基于log2(RSEM+1)转化的表达量做的相关性分析呢？实际是log2(RSEM+0.01)，结果如下：
与cBioPortal结果同样一致，也同样并非负相关性最强....反倒奔着近似最强正相关去了~~
更有意思的是， 作为cBioPortal数据来源的GDAC Firehose也提供了一个文件，同样是说，选择了与mRNAseq表达相关性最小的探针代表基因：
picking the probe with the min correlation with mRNAseq data for each gene
根据目标基因对应的Beta值能知道，这里pick的却是另外一个探针cg17204129！
不管这些'大佬'出于何种考虑，做了哪些玄幻的处理，还是希望大家对各类在线分析工具保持一定的怀疑态度，不仅仅是因为无法参透他们的骚操作，而且本身也存在一定的缺陷，例如：
1、cBioPortal最终pick的探针，并不一定就是你前期通过 甲基化差异表达分析 筛选得到的探针，即前后不一致导致分析结果的偏倚
2、cBioPortal中不包含正常对照样本
做个总结
如果想分析目标基因表达与甲基化相关性：
1、在线工具有cBioPortal、Xena
2、如果在线工具的分析结果不理想：
A:已有目标探针，则下载该探针的甲基化值，并使用文中的R脚本进行分析；
B:无目标探针，则下载该基因所有探针的甲基化值，并使用文中的R脚本进行分析。
有意思的问题
既然Xena和cBioPortal/Firehose均能得到目标基因表达log2(RSEM+1)及探针甲基化beta值，那么基于此数据，两者相关性分析结果是否一致？
解答
按理说应该一致，但实际并不一致！
学完本文，足以自行做个测试，多点思考🤔🤔🤔会有更多有意思的发现~
联系客服
