Question: Variance unstable after fitting to linear model with limma-voom pipeline
0
gravatar for kentfung
13 days ago by
kentfung0
kentfung0 wrote:

So I was doing a differential expression analysis with RNA-Seq data using limma-voom pipeline. However, the fit for mean-variance is not like a flat line at 1 but a curve (top: plot after voom; bottom: plot after fitting to linear model). Is it normal? Can I still use the results? voom<em>plotSA eBayes</em>plotSA

Here are my codes:

library(edgeR)
library(limma)

#Import data
txi <- tximport(files = read.hd5s,
                     type = "kallisto",
                     tx2gene = tx2gene,
                     countsFromAbundance = "lengthScaledTPM",
                     ignoreTxVersion = TRUE)

# create DGEList object
raw_data = DGEList(txi$counts)

#Filtering low count reads (TPM < 2)
library(genefilter)
abundance = txi$abundance
abs_crit = pOverA(0.2, 2)
abs_filter <- filterfun(abs_crit)
keep <- genefilter(abundance, abs_filter)
raw_data <- raw_data[keep,,keep.lib.sizes=FALSE]
nrow(raw_data)

# calculate normalization factors to scale the raw library sizes
raw_data = calcNormFactors(raw_data)

# experiment design (cluster from previous clustering analysis)
design <- model.matrix(~condition)

# Batch effect correction
dat0 = as.matrix(raw_data$counts)
mod1 = design
mod0 = cbind(mod1[,1])
set.seed(12345)
svseq = svaseq(dat0,mod1,mod0)

# using limma to find differentially expressed genes
#model fit with batch
new_design <- model.matrix(~0+condition)
new_design = cbind(new_design,svseq$sv)

#voom transformation
voom_diff = voom(raw_data, new_design, normalize="none", plot = TRUE)

#contrast matrix
contr.matrix <- makeContrasts(conditionvscontrol = conditionTRUE - conditionFALSE,
                              levels = colnames(new_design))

#model fitting
vfit <- lmFit(voom_diff_ph, new_design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit,trend=TRUE)
plotSA(efit, main="Final model: Mean-variance trend")

First two columns of design matrix are like this:

> head(new_design[,1:2])
  conditionFALSE conditionTRUE
1            1           0
2            1           0
3            1           0
4            1           0
5            1           0
6            1           0

There are some 1 in TRUE and some 0 in FALSE but it's too long so I didn't include them.

limma voom • 39 views
ADD COMMENTlink modified 13 days ago • written 13 days ago by kentfung0

Sorry I know what is happening. It's normal if I choose trend=TRUE when fitting.

ADD REPLYlink written 12 days ago by kentfung0
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 16.09
Traffic: 187 users visited in the last hour