Entering edit mode
Morgan Mouchka
▴
30
@morgan-mouchka-5848
Last seen 10.5 years ago
Hello all,
Im currently using EdgeR to analyze RNASeq data and have very much
appreciated the software and its incredibly helpful user manual.
More specifically, I have two questions. The first is whether I have
chosen the correct method for analyzing my multi-factorial experiment.
I have 2 factors, state (apo, sym) and treatment (control, treatment).
Im interested if there is a differential response due to the
treatment that varies by state. In other words, do apos respond to
the treatment in the same way that syms respond?
To determine which genes were differentially expressed between the
control and treatment for each data set (i.e. apo and sym), I used a
GLM approach with contrasts to do pairwise comparisons between groups.
I have pasted my session output and preliminary code below.
Specifically, I could use advice on whether 1) my design matrix and 2)
the way I have called my contrasts are appropriate for my question.
Second, for genes that are differentially expressed in both the apo
and sym data sets, is it possible to test if the genes are
significantly differentially expressed from each other? For instance,
if a gene is 9-fold up regulated in apos, but only 3-fold up regulated
in syms, are these significantly different such that I could say that
this gene is significantly more up regulated in the apo state? It
seems that this could be some sort of interaction term, but Im
uncertain as to the best way to set this up.
Im using the 3.0.8 version of EdgeR and the 2.15.2 version of R. I
have 4 biological reps per treatment/state with the nomenclature of AC
for apo control, AT for apo treatment, SC for sym control, and ST for
sym treatment.
Thanks so much,
Morgan
Morgan Mouchka
PhD Candidate
E231 Corson Hall
Dept. Ecology and Evolutionary Biology
Cornell University
mep74@cornell.edu
http://www.eeb.cornell.edu/mouchka
> x <- read.delim("AC_AT_SC_ST_counts_cnidaln.txt",
row.names="Contig")
> head (x)
AC3 AC4
AC5 AC7 AT1 AT2
Locus_34979_Transcript_1/1_Confidence_1.000_Length_792 72 39
55 31 36 55
Locus_125599_Transcript_2/4_Confidence_1.000_Length_409 3 2
3 4 0 2
Locus_100348_Transcript_1/1_Confidence_0.000_Length_1519 1260 576
957 529 663 961
Locus_129112_Transcript_1/1_Confidence_0.000_Length_436 17 7
6 6 7 8
Locus_5661_Transcript_2/8_Confidence_1.000_Length_1090 0 0
0 0 0 0
Locus_79922_Transcript_1/1_Confidence_0.000_Length_6499 7774 3857
5588 3835 3633 4470
AT4 AT5
SC1 SC2 SC3 SC5
Locus_34979_Transcript_1/1_Confidence_1.000_Length_792 22 63
20 55 32 49
Locus_125599_Transcript_2/4_Confidence_1.000_Length_409 0 5
4 4 1 2
Locus_100348_Transcript_1/1_Confidence_0.000_Length_1519 516 903
407 1072 345 788
Locus_129112_Transcript_1/1_Confidence_0.000_Length_436 4 5
3 0 2 5
Locus_5661_Transcript_2/8_Confidence_1.000_Length_1090 0 0
0 0 0 0
Locus_79922_Transcript_1/1_Confidence_0.000_Length_6499 3668 5489
1092 2288 1381 2194
ST1 ST2
ST3 ST5
Locus_34979_Transcript_1/1_Confidence_1.000_Length_792 39 65
40 82
Locus_125599_Transcript_2/4_Confidence_1.000_Length_409 7 0
1 3
Locus_100348_Transcript_1/1_Confidence_0.000_Length_1519 958 1103
900 1488
Locus_129112_Transcript_1/1_Confidence_0.000_Length_436 7 11
6 8
Locus_5661_Transcript_2/8_Confidence_1.000_Length_1090 0 0
0 0
Locus_79922_Transcript_1/1_Confidence_0.000_Length_6499 3629 4837
3383 4221
> #assign samples to groups, A=AC, B=AT, C=SC, D=ST
> group <- factor(c("A","A","A","A","B","B","B","B","C","C","C","C","D
","D","D","D"))
> # Define data object (list-based, includes info for counts and info
about samples) and will calculate library size
> y <- DGEList (counts = x, group = group)
Calculating library sizes from column totals.
> #filter out very lowly expressed genes, based on number of lowest
number of reps (4 in this case)
> keep <- rowSums (cpm(y)>1) >= 4
> y2 <- y[keep,]
> dim(y2)
[1] 17739 16
> y2$samples$lib.size <- colSums(y2$counts)
> # Normalize for sample-specific effects
> y3 <- calcNormFactors (y2)
> y3$samples
group lib.size norm.factors
AC3 A 11095093 1.0237872
AC4 A 5604841 0.9912734
AC5 A 7517691 0.9824405
AC7 A 5918039 0.9919439
AT1 B 6004830 0.9609672
AT2 B 7692419 0.9543765
AT4 B 5197815 0.9575609
AT5 B 7759281 0.9643070
SC1 C 3542875 1.0693587
SC2 C 6720709 1.0488297
SC3 C 3881774 1.0022528
SC5 C 5109584 1.0429925
ST1 D 10429666 0.9815467
ST2 D 9662842 1.0103866
ST3 D 8267838 1.0198095
ST5 D 10698815 1.0069064
> #plot biological coefficient of variation (BVC)between samples
> plotMDS(y3)
> #check levels
> levels (y3$samples$group)
[1] "A" "B" "C" "D"
> #design matrix to describe treatment conditions
> design <- model.matrix(~0+group, data=y3$samples)
> colnames(design) <- c("A", "B", "C", "D")
> design
A B C D
AC3 1 0 0 0
AC4 1 0 0 0
AC5 1 0 0 0
AC7 1 0 0 0
AT1 0 1 0 0
AT2 0 1 0 0
AT4 0 1 0 0
AT5 0 1 0 0
SC1 0 0 1 0
SC2 0 0 1 0
SC3 0 0 1 0
SC5 0 0 1 0
ST1 0 0 0 1
ST2 0 0 0 1
ST3 0 0 0 1
ST5 0 0 0 1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
> y4 <- estimateGLMCommonDisp(y3, design, verbose=TRUE)
Disp = 0.02061 , BCV = 0.1436
> #estimate gene-wise dispersion
> y5 <- estimateGLMTrendedDisp(y4, design)
> y6 <- estimateGLMTagwiseDisp(y5,design)
> plotBCV(y6)
> fit <- glmFit(y6, design)
> # make contrasts
> my.contrasts <- makeContrasts(BvsA=B-A, DvsC=D-C, levels=design)
> #pairwise comparison of B and A (apo treatmetn and apo control)
> lrt.BvsA <- glmLRT(fit, contrast=my.contrasts[,"BvsA"])
> # summary of model fitting funtions
> summary(de<-decideTestsDGE(lrt.BvsA))
[,1]
-1 562
0 16091
1 1086
> lrt.DvsC <- glmLRT(fit, contrast=my.contrasts[,"DvsC"])
> summary(de<-decideTestsDGE(lrt.DvsC))
[,1]
-1 496
0 16710
1 533
[[alternative HTML version deleted]]
