有时候,我们经常需要探讨,两组学数据如 基因与代谢物(转录组与代谢组),代谢物与微生物(代谢组与宏基因组),基因与蛋白质(转录组与蛋白组)等的表达关联,就需要用到两个组学的(差异)表达矩阵进行相关性分析。希望得出其中一组差异表达的物质对另一组差异表达的物质的的影响(如,正相关/负相关)。
- 数据准备
首先,需要获取两个组学各自的表达矩阵,这里以转录组和代谢组为例。
假设转录组表达矩阵如下:(fpkm)
df1 <- read.table(Mat1_PATH, sep = "\t", header = T, row.names = 1)
行名为基因ID,列名为样本名
ID | A01 | A02 | A03 | B01 | B02 | B03 |
---|---|---|---|---|---|---|
id1 | 1.48211 | 40.491178 | 131.64376 | 0.31306 | 0.39617 | 0.61075 |
id2 | 2.71537 | 2.16734 | 1.56108 | 1.70055 | 5.64507 | 1.194945 |
id3 | 136.59175 | 71.3823002 | 131.643764 | 158.059778 | 153.335882 | 132.485619 |
id4 | 0.374847 | 0.610750 | 0.137045 | 0.543599 | 0.279865 | 0.91199 |
id5 | 8.8307147 | 10.9066534 | 8.3346599 | 9.3874084 | 7.6636695 | 8.940899018 |
id6 | 8.22970427 | 10.05391844 | 20.135060 | 9.145283 | 19.45078117 | 21.5074579 |
代谢组表达矩阵如下:
df2 <- read.table(Mat2_PATH, sep = "\t", header = T, row.names = 1)
行名为代谢物ID,列名为样本名
ID | A01 | A02 | A03 | B01 | B02 | B03 |
---|---|---|---|---|---|---|
mid1 | 10709.86303 | 8011.184598 | 11131.11111 | 10886.6822 | 21090.99672 | 9360.52 |
mid2 | 12682.32 | 13572.24 | 11378.47 | 15024.66 | 17371.2 | 11949.45 |
mid3 | 136591.75 | 713823.002 | 11014.19899 | 6951.44 | 6825.125249 | 10599.68722 |
mid4 | 21267.584 | 21839.02115 | 15284.701 | 48790.3521 | 22752.39878 | 18467.26978 |
mid5 | 22859.54 | 45521.4 | 83815.55 | 53705.52 | 5053.8767 | 8557.899 |
mid6 | 82297.0427 | 10053.914 | 20135.06 | 91452.83 | 19450.78 | 21507.4 |
两个矩阵的样本名需要一一对应哦
- 数据预处理
当表达矩阵过于庞大的时候, 相关性计算非常耗费计算资源,我们可以仅对我们感兴趣的基因/代谢物进行分析。例如,可以对感兴趣的基因(如 有KEGG注释的基因,在各个样本中均广泛表达的基因,具有已知名称的代谢物等)
# 转置矩阵,变为,行为样本名,列为基因ID
df1 <- t(df1)
# 保留至少在5个样本中表达的基因
df1 <- df1[, colSums(df1[, , drop = F] != 0) >= 5, drop = F]
# 转置矩阵,变为,行为样本名,列为代谢物ID
df2 <- t(df2)
# 保留已命名的代谢物
nameID <- df_metaAnno$ID # 假设存在一个注释文件,保存着df2中已命名的代谢物
df2 <- df2[, nameID ]
- 相关性分析
相关性分析
分析连续变量之间的线性相关程度的强弱
常见的就是 Person 和 sperman 相关性分析
Ref:
Pearson相关分析:https://zhuanlan.zhihu.com/p/347467102
Person 相关系数分析要求,数据没有明显的异常值, 双变量正态分布存在
Sperman 相关系数,不要求正态分布
相关系数的绝对值越大,相关性越强,相关系数越接近于1或-1,负数表示负相关,蒸熟表示正相关,相关度越强,相关系数越接近于0,相关度越弱。
通常情况下通过以下取值范围判断变量的相关强度:
相关系数 0.8-1.0 极强相关
0.6-0.8 强相关
0.4-0.6 中等程度相关
0.2-0.4 弱相关
0.0-0.2 极弱相关或无相关
R中计算线性相关性的R包: Hmisc包中rcorr()函数计算p值和r值
(ref: https://zhuanlan.zhihu.com/p/671659920)
corr1 <- rcorr(df1, df2, type = c("pearson"))
- 数据后处理
4.1 提取 两组学矩阵中基因与代谢物之间的相关性
corr2 <- list(
corMatrix = corr1$r[1:ncol(df1), (ncol(df1) + 1):ncol(corr1$r)],
pMatrix = corr1$P[1:ncol(df1), (ncol(df1) + 1):ncol(corr1$P)]
)
解释,
corr1$r
中存放着相关性系数r
corr1$P
中存放着相关性P值
由于我们只关心 基因与代谢物之间的相关性,而不关心基因-基因,代谢物-代谢物之间的相关性,所有需要从结果表中仅提取 基因-代谢物的相关性结果
例如 corr1$r
矩阵的格式为:
cor | id1 | id2 | id3 | mid1 | mid2 | mid3 |
---|---|---|---|---|---|---|
id1 | 1 | 0.5 | 0.8 | 0.4 | 0.5 | 0.1 |
id2 | 0.5 | 1 | 0.55 | 0.77 | 0.69 | -0.83 |
id3 | 0.8 | 0.55 | 1 | -0.38 | 0.64 | -0.66 |
mid1 | 0.4 | 0.77 | -0.38 | 1 | 0.06 | 0.5 |
mid2 | 0.5 | 0.69 | 0.64 | 0.06 | 1 | 0.07 |
mid3 | 0.1 | -0.83 | -0.66 | 0.5 | 0.07 | 1 |
我们只需要提取除右上角的结果即可即
cor | mid1 | mid2 | mid3 |
---|---|---|---|
id1 | 0.4 | 0.5 | 0.1 |
id2 | 0.77 | 0.69 | -0.83 |
id3 | -0.38 | 0.64 | -0.66 |
corr1$P
矩阵也类似
4.2 合并相关系数r与相关性P值 -- 拉平表格
# flatten matrix of correlation coefficient
rDt <- data.frame(ID = rownames(corr2$corMatrix), corr2$corMatrix) %>%
setDT(.) %>%
data.table::melt(., id.vars = c("ID"), variable.name = "Meta", value.name = "corr")
ID | Meta | corr |
---|---|---|
id1 | mid1 | 0.4 |
id2 | mid1 | 0.77 |
id3 | mid1 | -0.38 |
id1 | mid2 | 0.5 |
id2 | mid2 | 0.69 |
id3 | mid2 | 0.64 |
id1 | mid3 | 0.1 |
id2 | mid3 | -0.83 |
id3 | mid3 | -0.66 |
相关性P值的矩阵处理也是如此
# flatten matrix of P-value
pDt <- data.frame(ID = rownames(corr2$pMatrix), corr2$pMatrix) %>%
setDT(.) %>%
data.table::melt(., id.vars = c("ID"), variable.name = "Meta", value.name = "p_value")
最后合并相关系数r和相关性P值到一个表格
# join data
rpDt <- rDt[pDt, on = .(ID = ID, Meta = Meta)]
例如:
ID | Meta | corr | p_value |
---|---|---|---|
id1 | mid1 | 0.4 | 0.04 |
id2 | mid1 | 0.77 | 0.001 |
id3 | mid1 | -0.38 | 0.05 |
id1 | mid2 | 0.5 | 0.03 |
id2 | mid2 | 0.69 | 0.0000006 |
id3 | mid2 | 0.64 | 0.0007 |
id1 | mid3 | 0.1 | 0.3 |
id2 | mid3 | -0.83 | 0.006 |
id3 | mid3 | -0.66 | 0.535 |
我们可以将结果保存为RDS或者excel的格式,用于后续分析是输入数据
saveRDS(corr2, file = "corr_result.rds"))
或者
fwrite(rpDt,
file = "corr_result.tsv.gz", quote = FALSE, sep = "\t",
col.names = TRUE, row.names = FALSE, compress = "gzip"
)
至此,相关性分析的内容就结束啦,后续可以基于该结果进行绘图,常见的表现形式有,网络图,桑基图,和弦图,热图等。