Question: Defining the design for a complicated multi-level experiment in Limma
0
gravatar for Colari19
12 months ago by
Colari1950
Colari1950 wrote:

Hello,

I've carried out a number of differential expression analyses using the Limma-Voom pipeline for a rather complicated data set. I'd just like to clarify that my approach detailed below is appropriate.

The design of the experiment is very much like those described in sections 9.6 and 9.7 of the Limma manual (time-course experiments and multi-level experiments).

Our data set comprises of RNA-Seq data from patients undertaking treatment with one of two drugs (drug A and drug B). For each subject in our data set, samples were collected at three time points from two different tissue types on the body. Samples from diseased tissue and healthy tissue were collected at week 0 (before commencing treatment), week 2 (two weeks after commencing treatment), and week 10 (ten weeks after commencing treatment).

In addition, each patient was genotyped for an allele of interest (we'll call it Allele.X), and designated as Allele.X-positive (1 or 2 copies of allele) or Allele.X-negative (0 copies of allele).

Therefore, a snippet of the data set looks something like this:

Subject      Sample     Drug    Time    Tissue      Allele.X
Subject.01  Sample.01   Drug.A  Week.0  Diseased    Positive
Subject.01  Sample.02   Drug.A  Week.2  Diseased    Positive
Subject.01  Sample.03   Drug.A  Week.10 Diseased    Positive
Subject.01  Sample.04   Drug.A  Week.0  Healthy     Positive
Subject.01  Sample.05   Drug.A  Week.2  Healthy     Positive
Subject.01  Sample.06   Drug.A  Week.10 Healthy     Positive
Subject.02  Sample.07   Drug.B  Week.0  Diseased    Negative
Subject.02  Sample.08   Drug.B  Week.2  Diseased    Negative
Subject.02  Sample.09   Drug.B  Week.10 Diseased    Negative
Subject.02  Sample.10   Drug.B  Week.0  Healthy     Negative
Subject.02  Sample.11   Drug.B  Week.2  Healthy     Negative
Subject.02  Sample.12   Drug.B  Week.10 Healthy     Negative

In total, 35 patients were treated with drug A and 31 were treated with drug B.

In the drug A group, 20 were Allele.X-negative and 15 were were Allele.X-positive. In the drug B group, 18 were Allele.X-negative and 13 were Allele.X-positive.

I would like to identify genes that change in expression over the course of treatment in both tissue types in both drug groups, irrespective of Allele.X status, i.e: week 2 and week 10 compared to week 0 in ALL of the drug A group and ALL of the drug B group.

I would also like to do the above comparisons, but also taking into account Allele.X status.

My approach do this has been to analyse each drug group separately, so when looking at drug A for example, this is done by subsetting the clinical data and expression data...

drug = "Drug.A"
clinical.data <- clinical.data[clinical.data$Drug == "Drug.A",]
exp.data <- exp.data[, clinical.data$Sample]

...and creating an interaction variable that combines tissue, time, and Allele.X status...

clinical.data$Tissue.Time.Allele.X = paste(clinical.data$Tissue, clinical.data$Time, clinical.data$Allele.X, sep = ".")

...and an interaction variable that combines subject and tissue...

clinical.data$Subject.Tissue = paste(clinical.data$Subject, clinical.data$Tissue, sep = ".")

The Tissue.Time.Allele.X variable is used to define the design...

limma_design = model.matrix(~ 0 + Tissue.Time.Allele.X, data = keep_samp)
voomed_counts = voom(edata, limma_design)

...and the Subject.Tissue variable is used to for sample blocking...

# Block correlation
my_cor = duplicateCorrelation(voomed_counts,  limma_design,  block = clinical.data$Subject.Tissue)
fit = lmFit(voomed_counts, limma_design, correlation = my_cor$cor, block = clinical.data$Subject.Tissue)

The makeContrasts function is used to extract the relevant contrasts...

# Extract contrasts of interest
my_contrasts <- makeContrasts(

Diseased_week2_vs_week0_All =
  ((15 * Diseased.Week.2.Positive + 20 * Diseased.Week.2.Negative) / 35) -
  ((15 * Diseased.Week.0.Positive + 20 * Diseased.Week.0.Negative) / 35),

Diseased_week10_vs_week0_All =
  ((15 * Diseased.Week.10.Positive + 20 * Diseased.Week.10.Negative) / 35) -
  ((15 * Diseased.Week.0.Positive + 20 * Diseased.Week.0.Negative) / 35),

Healthy_week2_vs_week0_All =
  ((15 * Diseased.Week.2.Positive + 20 * Diseased.Week.2.Negative) / 35) -
  ((15 * Diseased.Week.0.Positive + 20 * Diseased.Week.0.Negative) / 35),

Healthy_week10_vs_week0_All =
  ((15 * Healthy.Week.10.Positive + 20 * Healthy.Week.10.Negative) / 35) -
  ((15 * Healthy.Week.0.Positive + 20 * Healthy.Week.0.Negative) / 35),

Diseased_week2_vs_week0_Pos = Diseased.Week.2.Positive - Diseased.Week.0.Positive,
Diseased_week2_vs_week0_Neg = Diseased.Week.2.Negative - Diseased.Week.0.Negative,
Diseased_week10_vs_week0_Pos = Diseased.Week.10.Positive - Diseased.Week.10.Positive,
Diseased_week10_vs_week0_Neg = Diseased.Week.10.Negative - Diseased.Week.10.Negative,

Healthy_week2_vs_week0_Pos = Healthy.Week.2.Positive - Healthy.Week.0.Positive,
Healthy_week2_vs_week0_Neg = Healthy.Week.2.Negative - Healthy.Week.0.Negative,
Healthy_week10_vs_week0_Pos = Healthy.Week.10.Positive - Healthy.Week.10.Positive,
Healthy_week10_vs_week0_Neg = Healthy.Week.10.Negative - Healthy.Week.10.Negative

levels = limma_design
)

# Compute contrasts 
fit <- contrasts.fit(fit, my_contrasts)
fit <- eBayes(fit)

Hopefully it is clear that for the contrasts in all of the patients within the drug group, the average expression (weighted by sample number) of the Allele.X-positive and negative patients is being contrasted.

So my main question - is the above approach valid?

Thank you! And sorry for long post, but hopefully it is clear what I'm getting at.

voom rna-seq limma R • 465 views
ADD COMMENTlink written 12 months ago by Colari1950

Hello Colari19!

It appears that your post has been cross-posted to another site: https://support.bioconductor.org/p/125846/

This is typically not recommended as it runs the risk of annoying people in both communities.

ADD REPLYlink written 12 months ago by ATpoint40k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1652 users visited in the last hour