Вычислите глобальную нормализацию медианы в микрочипе

#r #statistics #bioinformatics #median

Вопрос:

Я хочу извлечь передний план/фон Cy5 для 4 массивов и вычесть фон из значений переднего плана, а затем преобразовать эти значения в log2. Затем я хочу рассчитать глобальную нормализацию медианы для этих 4 массивов, используя вычитаемые из фона значения Cy5. Для масштабирования будет использоваться медиана каждого массива. После нормализации все массивы должны иметь медиану 1.

Однако мой код ниже оценивает медиану как 0 вместо 1. Почему? И что я должен изменить, чтобы получить медиану 1?

 library(limma)
library(marray)


for(i in 1:4){
  name <- paste("sample", i, sep = ".")
  bg <- maRb(dat[,i])
  fg <- maRf(dat[,i])
  diff <- fg - bg
  diff[diff < 0] <- NA
  assign(name, log2(diff))
} 

data.prenorm <- cbind(sample.1, sample.2, sample.3, sample.4)
data.median  <- apply(data.prenorm, 2, median, na.rm = T)
data.norm    <- sweep(data.prenorm, 2, data.median)

colnames(data.norm) <- c("Array 1", "Array 2", "Array 3", "Array 4")

median(data.norm[ , 1], na.rm = T) 
median(data.norm[ , 2], na.rm = T)
median(data.norm[ , 3], na.rm = T)
median(data.norm[ , 4], na.rm = T)
 

данные:

 > dput(data.norm[1:4,1:4])
structure(c(0.335603031784438, 0.192645077942395, 0.280107919192734, 
4.59067615191555, 0, 0, -0.362570079384708, 6.14068778021722, 
-0.192645077942395, -0.263034405833793, -0.192645077942395, 3.4262647547021, 
-0.231325546106455, 0, -0.754887502163468, 6.13689620105484), .Dim = c(4L, 
4L), .Dimnames = list(NULL, c("Array 1", "Array 2", "Array 3", 
"Array 4")))
 

Ответ №1:

Я не являюсь экспертом в обработке микрочипов, поэтому могу лишь высказать некоторые предположения о том, чего вы пытаетесь достичь.

Видите ли, медиана нормализованного по медиане массива равна 0. Чтобы увидеть это, представьте, что у вас есть отсортированный массив данных x и массив равной длины, заполненный медианой x, назовите его y. Медиана x находится в его центральном положении, и медиана (x-y) также будет находиться в центральном положении, поскольку вычитание константы для каждого элемента x не изменяет их относительный размер. Центральное положение (x-y) равно медиане(x)-медиане(x) по определению, и это равно 0.

Здесь также есть несколько хороших дискуссий по перекрестной проверке по этому вопросу.

А теперь возвращаемся к вашей точке зрения. Я подозреваю, что то, что вы хотели обработать, должно быть в пространстве журнала. Это даст вам медиану нормализованных данных 0. Затем после преобразования нормализованных данных обратно в исходное пространство медиана будет равна 1.

Вот быстрая проверка с использованием некоторых фиктивных данных:

 dat = matrix(1:16, nrow = 4)
ldat = log2(dat)
ldat_norm = sweep(ldat,2,apply(ldat,2,median))
dat_norm = 2^ldat

median = apply(ldat,2,median)

> median
[1] 1.020621 1.002972 1.001136 1.000595