如何通过计算相关系数来降低基因表达矩阵的维数？答案

【问题标题】：How to reduce dimension of gene expression matrix by calculating correlation coefficients?如何通过计算相关系数来降低基因表达矩阵的维数？
【发布时间】：2019-11-05 06:10:45
【问题描述】：

我有兴趣找到一系列基因之间的 Pearson 相关系数。基本上，我有 Affymetrix 基因水平表达矩阵（行中的基因和列中的样本 ID），并且我有微阵列实验观察的注释数据，其中行中的样本 ID 和列中的描述标识。

数据

> expr_mat[1:8, 1:3]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at                6.062215        6.125023        5.875502
10_at               3.796484        3.805305        3.450245
100_at              5.849338        6.191562        6.550525
1000_at             3.567779        3.452524        3.316134
10000_at            6.166815        5.678373        6.185059
100009613_at        4.443027        4.773199        4.393488
100009676_at        5.836522        6.143398        5.898364
10001_at            6.330018        5.601745        6.137984

> anodat[1:8, 1:3]
               V1   V2    V3
1        SampleID   GA Batch
2 Tarca_001_P1A01   11     1
3 Tarca_013_P1B01 15.3     1
4 Tarca_025_P1C01 21.7     1
5 Tarca_037_P1D01 26.7     1
6 Tarca_049_P1E01 31.3     1
7 Tarca_061_P1F01 32.1     1
8 Tarca_051_P1E03 19.7     1

目标：

我打算看看每个样本中的基因如何与注释数据中对应样本的GA值相关，然后生成与目标观察数据保持高相关基因的子表达矩阵anodat$GA。

我的尝试：

gene_corrs <- function(expr_mat, anno_mat){
    stopifnot(ncol(expr_mat)==nrow(anno_mat))
    res <- list()
    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]))
                res <- cbind(res, cor=cor_mat[out])
                res <- as.dara.frame(res)
            }
        })
    })
    return(res)
}

但是，我的上述实现并没有返回我的预期，我需要通过查找与anodat$GA 具有强相关性的基因来过滤掉这些基因。

另一次尝试：

我读过一些关于类似问题的帖子，有些人讨论过使用limma 包。这是我使用limma 的尝试。这里我使用anodat$GA 作为协变量来拟合limma 线性模型：

library(limma)
fit <- limma::lmFit(expr_mat, design = model.matrix( ~ 0 + anodat$GA)
fit <- eBayes(fit)
topTable(fit, coef=2)

然后我希望从上面的代码中得到一个相关矩阵，并希望执行以下操作以获得过滤后的子表达式矩阵：

idx <- which( (abs(cor) > 0.8) & (upper.tri(cor)), arr.ind=TRUE)
idx <- unique(c(idx[, 1],idx[, 2])
correlated.genes <- matrix[idx, ]

但我仍然没有得到正确的答案。我对使用limma 方法很有信心，但我无法再次弄清楚上面的代码出了什么问题。谁能指出我如何使这项工作？有没有什么有效的方法可以做到这一点？

【问题讨论】：

标签： r matrix bioinformatics correlation limma

【解决方案1】：

不要让你的数据很难仔细检查，但在摘要中我会尝试这个：

library(matrixTests)
cors <- row_cor_pearson(expr_mat, anodat$GA)

which(cors$cor > 0.9)  # to get the indeces of genes with correlation > 0.9

【讨论】：

谢谢，但是当我将阈值设置为0.9时，就没有基因（行）了，这是不正确的。你确定这可以拟合基因表达数据吗？
你必须选择你想要的阈值，这个0.9只是一个例子。同样在执行此操作之前 - 确保样本名称在 expr_mat 和 anodat$GA 之间匹配。根据您的示例 - 它们可能不匹配。