Signifigance and effect size calculations for Maize Flooding dataset

Load ALDEx2 and the count data

library(ALDEx2)
library(BiocParallel)
load("/Users/rivers/Documents/chip_quantseq/resultscount-matrix.txt.Rdata")

Order libraries 1 to 19 (10 is not in the dataset)

rownames(df)<-df[,1]
df <-df[,c(2,12,14:20,3:11,13)]

Herbivory

Run DE analysis on no herbivory vs herbivory

conds <- c(rep("NoHerb", 9), rep("Herb", 10))
x.all <- aldex(df, conds, mc.samples=128, test="t", effect=TRUE,
     include.sample.summary=FALSE, denom="all", verbose=FALSE)

aldex.clr: generating Monte-Carlo instances and clr values
operating in serial mode
computing center with all features
aldex.ttest: doing t-test
aldex.effect: calculating effect sizes

View(x.all)
View(y.all)
View(y.all)
View(y.all)

par(mfrow=c(1,2))
aldex.plot(x.all, type="MA", test="welch", xlab="Log-ratio abundance",
    ylab="Difference")
aldex.plot(x.all, type="MW", test="welch", xlab="Dispersion",
    ylab="Difference")

Filter the data

# Get the Genes with a Benjamini-Hochberg corrected p-value less than 0.05, 
# AND an effect size greater than 1

herb.down <- rownames(x.all)[which(x.all$we.eBH < 0.05 & x.all$effect < -1)]
herb.up <- rownames(x.all)[which(x.all$we.eBH < 0.05 & x.all$effect > 1)]

Flooding

Run DE analysis on flooding vs No flooding

conds <- c(rep("NoFlood", 5), rep("Flood", 4),rep("NoFlood", 5), rep("Flood", 5))
x.all <- aldex(df, conds, mc.samples=128, test="t", effect=TRUE,
     include.sample.summary=FALSE, denom="all", verbose=FALSE, useMC=TRUE)

aldex.clr: generating Monte-Carlo instances and clr values
operating in serial mode
computing center with all features
aldex.ttest: doing t-test
aldex.effect: calculating effect sizes

par(mfrow=c(1,2))
aldex.plot(x.all, type="MA", test="welch", xlab="Log-ratio abundance",
    ylab="Difference")
aldex.plot(x.all, type="MW", test="welch", xlab="Dispersion",
    ylab="Difference")

Filter the data

# Get the Genes with a Benjamini-Hochberg corrected p-value less than 0.05, 
# AND an effect size greater than 1

flood.down <- rownames(x.all)[which(x.all$we.eBH < 0.05 & x.all$effect < -1)]
flood.up <- rownames(x.all)[which(x.all$we.eBH < 0.05 & x.all$effect > 1)]

GLM Modeling of flooding and herbivory interactions

covariates <- data.frame("Flooding" =c(rep(0, 9),rep(1, 10)), "Herbivory" = c(rep(0, 5), rep(1, 4), rep(0, 5), rep(1, 5) ))
mm <- model.matrix(~ Flooding * Herbivory, covariates)
z <- aldex.clr(df, mm, mc.samples=128, denom="all", useMC=TRUE)

multicore environment is is OK -- using the BiocParallel package
computing center with all features
'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading'package:stats' may not be available when loading

## operating in serial mode
## Warning in aldex.clr.function(reads, conds, mc.samples, denom, verbose, :
## values are unreliable when estimated with so few MC smps
## computing center with all features
glm.test <- aldex.glm(z, mm)

the condition has length > 1 and only the first element will be usedrunning tests for each MC instance:
the condition has length > 1 and only the first element will be used

|--

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-(25%)

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

-

# write data out

write.table(x.all, "/Users/rivers/Documents/chip_quantseq/herbivory.txt", sep="\t")

write.table(y.all, "/Users/rivers/Documents/chip_quantseq/flooding.txt", sep="\t")

write.table(glm.test, "/Users/rivers/Documents/chip_quantseq/glm.txt", sep="\t")