Limma для сравнения Bulk RNA Seq с использованием makeContrasts и eBayes - PullRequest
Новая
       76

Limma для сравнения Bulk RNA Seq с использованием makeContrasts и eBayes

0 голосов
/ 28 февраля 2020

После дня поисков в интернете я решил, что лучше задать вопрос здесь.

Итак, эксперимент состоит в том, что у меня есть данные по объему РНК-последовательности 3 пациентов: A, B, C. И их данные RNA seq получены для предварительной обработки, цикла обработки 1, цикла обработки 2, цикла обработки 3.

Таким образом, в общей сложности у меня есть 12 образцов объемной последовательности РНК:

  • A.PreTreat -> A.Cycle1 -> A.Cycle2 -> A.Cycle3

  • B.PreTreat -> B.Cycle1 -> B.Cycle2 -> B .Cycle3

  • C .PreTreat -> C .Cycle1 -> C .Cycle2 -> C .Cycle3

Я хочу получить дифференциальный список генов между разными циклами (т.е. цикл 3 до предварительной обработки, цикл 3 до цикла 2), используя model.matrix(), lmFit(), makeContrasts(), contrasts.fit(), eBayes(), все из которых находятся в пакете лиммы.

Вот мой минимальный рабочий пример. 

library(limma)
# Already normalized expression set: rows are genes, columns are the 12 samples
normalized_expression <- matrix(data=sample(1:100), nrow=10, ncol=12)
colnames(normalized_expression) <- c("A.PreTreat", "A.Cycle1", "A.Cycle2", "A.Cycle3", "B.PreTreat", "B.Cycle1", "B.Cycle2", "B.Cycle3", "C.PreTreat", "C.Cycle1", "C.Cycle2", "C.Cycle3")

patient_and_treatment <- factor(colnames(normalized_expression), levels = colnames(normalized_expression))

design.matrix <- model.matrix(~0 + patient_and_treatment)
colnames(design.matrix) <- patient_and_treatment
fit <- lmFit(normalized_expression, design.matrix)

# I want to get a contrast matrix to get differential genes between cycle 3 treatment and pre-treatment in all patients
contrast.matrix <- makeContrasts("A.Cycle3+B.Cycle3+C.Cycle3-A.PreTreat-B.PreTreat-C.PreTreat",
                                 levels = levels(patient_and_treatment))

# Outputs Error of no residual degree of freedom
fit2 <- eBayes( contrasts.fit( fit, contrast.matrix ) )

# Want to run but cannot
summary(decideTests(fit2))

До сих пор я не застрял на остаточной ошибке степени свободы.

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

Любая помощь с благодарностью.

Спасибо!

1 Ответ

1 голос
/ 28 февраля 2020

Вы не можете иметь 1 наблюдение на группу, это делает регрессию бессмысленной, поскольку вы подгоняете каждую точку данных к себе.

Вкратце, вы ищете общие эффекты, наблюдаемые у всех пациентов, например, цикл 3 по сравнению с PreTreat и т. д. настройте модель следующим образом: 

library(limma)

metadata = data.frame(
Patient=gsub("[.][^ ]*","",colnames(normalized_expression)),
Treatment=gsub("^[A-Z][.]*","",colnames(normalized_expression))
)

   Patient    Treatment
1        A PreTreat
2        A   Cycle1
3        A   Cycle2
4        A   Cycle3
5        B PreTreat
6        B   Cycle1
7        B   Cycle2
8        B   Cycle3
9        C PreTreat
10       C   Cycle1
11       C   Cycle2
12       C   Cycle3

Теперь укажите матрицу модели. Термин «Пациент» должен учитывать различия в начальных уровнях между пациентами: 

design.matrix <- model.matrix(~0 + Treatment+Patient,data=metadata)
fit <- lmFit(normalized_expression, design.matrix)
contrast.matrix <- makeContrasts(TreatmentCycle3-TreatmentPreTreat,
TreatmentCycle1-TreatmentPreTreat,levels=design.matrix)
fit2 = contrasts.fit(fit, contrast.matrix)
fit2 = eBayes(fit2)

Вы можете проверить, что коэффициенты дают вам то, что вы хотели: 

fit2$coefficients
       Contrasts
        TreatmentCycle3 - TreatmentPreTreat
   [1,]                           -3.666667
   [2,]                          -13.666667
   [3,]                            1.666667
   [4,]                          -40.666667
   [5,]                           12.000000
   [6,]                          -46.000000
   [7,]                          -32.000000
   [8,]                            4.666667
   [9,]                           11.333333
  [10,]                            5.666667
       Contrasts
        TreatmentCycle1 - TreatmentPreTreat
   [1,]                           -11.33333
   [2,]                           -19.33333
   [3,]                           -27.33333
   [4,]                           -42.33333
   [5,]                            27.33333
   [6,]                           -32.66667
   [7,]                           -33.00000
   [8,]                           -30.66667
   [9,]                            46.00000
  [10,]                            17.33333
...