如何通过相关值正确过滤掉基因表达矩阵？答案

【问题标题】：How to correctly filter out gene expression matrix by correlation value?如何通过相关值正确过滤掉基因表达矩阵？
【发布时间】：2019-11-03 16:30:57
【问题描述】：

我已经预处理了 Affymetrix 微阵列基因表达数据（行中的 32830 个探针组，列中的 735 个 RNA 样本）。这是我的表达式矩阵的样子：

> exprs_mat[1:6, 1:4]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04 Tarca_005_P1A05
1_at                6.062215        6.125023        5.875502        6.126131
10_at               3.796484        3.805305        3.450245        3.628411
100_at              5.849338        6.191562        6.550525        6.421877
1000_at             3.567779        3.452524        3.316134        3.432451
10000_at            6.166815        5.678373        6.185059        5.633757
100009613_at        4.443027        4.773199        4.393488        4.623783

我还有这个 Affymetrix 表达式的表型数据（行中的 RNA 样本标识符，列中的样本描述）：

 > pheno[1:6, 1:4]
                       SampleID   GA Batch     Set
Tarca_001_P1A01 Tarca_001_P1A01 11.0     1 PRB_HTA
Tarca_013_P1B01 Tarca_013_P1B01 15.3     1 PRB_HTA
Tarca_025_P1C01 Tarca_025_P1C01 21.7     1 PRB_HTA
Tarca_037_P1D01 Tarca_037_P1D01 26.7     1 PRB_HTA
Tarca_049_P1E01 Tarca_049_P1E01 31.3     1 PRB_HTA
Tarca_061_P1F01 Tarca_061_P1F01 32.1     1 PRB_HTA

由于在 phenodata 中，行中的样本标识符，我需要找到将 phenodata 中的 sampleID 与表达式矩阵 exprs_mat 中的 sampleID 匹配的方法。

目标：

我想通过测量每个基因与phenodata 中的目标谱数据之间的相关性来过滤出表达矩阵中的基因。这是我最初的尝试，但不太确定准确性：

更新：我在 R 中的实现：

我打算看看每个样本中的基因如何与注释数据中相应样本的 GA 值相关。这是我在 R 中找到这种相关性的简单函数：

getPCC <- function(expr_mat, anno_mat, verbose=FALSE){
stopifnot(class(expr_mat)=="matrix")
stopifnot(class(anno_mat)=="matrix")
stopifnot(ncol(expr_mat)==nrow(anno_mat))
final_df <- as.data.frame()
lapply(colnames(expr_mat), function(x){
    lapply(x, rownames(y){
        if(colnames(x) %in% rownames(anno_mat)){
            cor_mat <- stats::cor(y, anno_mat$GA, method = "pearson")
            ncor <- ncol(cor_mat)
            cmatt <- col(cor_mat)
            ord <- order(-cmat, cor_mat, decreasing = TRUE)- (ncor*cmatt - ncor)
            colnames(ord) <- colnames(cor_mat)
            res <- cbind(ID=c(cold(ord), ID2=c(ord)))
            res <- as.data.frame(cbind(out, cor=cor_mat[res]))
            final_df <- cbind(res, cor=cor_mat[out])
        }
    })
})
return(final_df)

}

但上面的脚本没有返回我期望的正确输出。有什么想法可以正确实现吗？有什么想法吗？

【问题讨论】：

如何过滤？高相关，低相关？另外，请注意expr_mat 中的列名与pheno 中的顺序（Sample_ID）不匹配（您可能需要先匹配它们）。
@PoGibas 我刚刚更新了我的帖子，在找到表达矩阵中的基因与表型数据中的基因之间的相关性之前，需要找到expr_mat 中的sampleID 和phenodata 中的sampleID 的匹配。你有什么主意吗？如何纠正上述做法？谢谢
@PoGibas 我用我的代码更新了我的帖子。有什么想法吗？
@Jerry 你能显示示例输出吗
@cephalopod 样本输出应该与expr_mat 矩阵具有相同的格式，其中应包含具有高相关值的过滤基因列表。

标签： r bioinformatics

【解决方案1】：

这样做有帮助吗：

library(tidyverse)

x <- data.frame(stringsAsFactors=FALSE,
     Levels = c("1_at", "10_at", "100_at", "1000_at", "10000_at", "100009613_at"),
     Tarca_001_P1A01 = c(6.062215, 3.796484, 5.849338, 3.567779, 6.166815,
                           4.443027),
     Tarca_003_P1A03 = c(6.125023, 3.805305, 6.191562, 3.452524, 5.678373,
                           4.773199),
     Tarca_004_P1A04 = c(5.875502, 3.450245, 6.550525, 3.316134, 6.185059,
                           4.393488),
     Tarca_005_P1A05 = c(6.126131, 3.628411, 6.421877, 3.432451, 5.633757,
                           4.623783)
     )


y <- data.frame(stringsAsFactors=FALSE,
     gene = c("Tarca_001_P1A01", "Tarca_013_P1B01", "Tarca_025_P1C01",
              "Tarca_037_P1D01", "Tarca_049_P1E01", "Tarca_061_P1F01"),
     SampleID = c("Tarca_001_P1A01", "Tarca_013_P1B01", "Tarca_025_P1C01",
                    "Tarca_037_P1D01", "Tarca_049_P1E01", "Tarca_061_P1F01"),
     GA = c(11, 15.3, 21.7, 26.7, 31.3, 32.1),
     Batch = c(1, 1, 1, 1, 1, 1),
     Set = c("PRB_HTA", "PRB_HTA", "PRB_HTA", "PRB_HTA", "PRB_HTA", "PRB_HTA")
     )



x %>% gather(SampleID, value, -Levels) %>% 
  left_join(., y, by = "SampleID") %>% 
  group_by(SampleID) %>% 
  filter(value == max(value)) %>% 
  spread(SampleID, value)

【讨论】：