Я заинтересован в поиске коэффициентов корреляции Пирсона между списком генов.По сути, у меня есть матрица выражений на уровне генов Affymetrix (гены в строках и идентификатор образца в столбцах), и у меня есть данные аннотации наблюдения за экспериментом с микрочипами, где идентификатор образца в строках и идентификация описания в столбцах.
данные
> 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.
Еще одна попытка :
Iпрочитайте несколько постов о подобной проблеме, и некоторые люди обсуждали использование пакета 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, но я не мог понять, что пошло не так в коде снова.Кто-нибудь может указать мне, как заставить это работать?Есть ли эффективный способ сделать это?