In avastermark19/doseR: doseR

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(mclust)
library(doseR)
library(edgeR)

xassetCountrySector <- read.table('~/doseR/extdata/format.out.Heliconius.LEG',  sep="\t", header=TRUE, as.is=TRUE)
GFF <- read.table('~/doseR/extdata/gff_parser.out.Heliconius.LEG',  sep="\t", header=TRUE, as.is=TRUE)
gene.chr.match <- GFF[match(xassetCountrySector$Gene, GFF$Gene), ]

counts.leg   <- as.matrix(round(xassetCountrySector[3:ncol(xassetCountrySector)] ))
Groupings <- rep("X", ncol(xassetCountrySector)-2 )
segLengths <- xassetCountrySector$Len
cd.leg <- new("countDat", data = counts.leg, replicates = factor(Groupings), annotation = gene.chr.match  )
cd.leg@rowObservables$seglens = segLengths
libsizes(cd.leg) <-   unname(getLibsizes2(cd.leg,   estimationType = "edgeR"))

cd.leg@replicates<- as.factor(c("F","M","M","M","F","F" ))

cd.leg@RPKM<- make_RPKM(cd.leg)

cd.leg@annotation$something <- (gene.chr.match$Chr == "Z")
cd.leg@annotation$something[cd.leg@annotation$something==TRUE] <- "Z"
cd.leg@annotation$something[cd.leg@annotation$something==FALSE | is.na(cd.leg@annotation$something)] <- "A"

## Factorized anntoation column input:
cd.leg@annotation$something <- factor(x = cd.leg@annotation$something, levels = c("A", "Z"))

plotExpr(cd.leg, col=c("black","red","black","red"), notch=T, outline=FALSE, cex.axis=0.8, mode_mean=FALSE, LOG2=TRUE, clusterby_grouping=FALSE, groupings="something")

# REMOVING 28%
f_cd.leg <- dafsFilter(cd.leg)
plotExpr(f_cd.leg, col=c("black","red","black","red"), notch=T, outline=FALSE, cex.axis=0.8, mode_mean=FALSE, LOG2=TRUE, clusterby_grouping=FALSE, groupings="something")

# REMOVING 56%
f_cd.leg <- quantFilter(cd.leg, lo.bound=0.2) 
plotExpr(f_cd.leg, col=c("black","red","black","red"), notch=T, outline=FALSE, cex.axis=0.8, mode_mean=FALSE, LOG2=TRUE, clusterby_grouping=FALSE, groupings="something")

# REMOVING 71%
f_cd.leg <- iqrxFilter(cd.leg, iqr_multi = 1.5)
plotExpr(f_cd.leg, col=c("black","red","black","red"), notch=T, outline=FALSE, cex.axis=0.8, mode_mean=FALSE, LOG2=TRUE, clusterby_grouping=FALSE, groupings="something")

outlist <- generateStats(f_cd.leg, groupings="something", mode_mean=TRUE)
outlist$kruskal
outlist$summary

outlist <- test_diffs(f_cd.leg, groupings="something", mode_mean=TRUE,  treatment1="M", treatment2="F")
outlist$kruskal
outlist$summary

plotRatioBoxes(f_cd.leg, treatment1="M", treatment2="F", groupings="something", cex.axis=0.8, outline=FALSE)

plotRatioDensity(f_cd.leg, treatment1="M", treatment2="F", mode_mean=TRUE, LOG2=TRUE, col =c("black","red"), lty = 1, type = "l", groupings="something")

# LIN MODS
dm<- cD.DM(cd.leg)

#gl.MF  <-  glSeq(dm, "-1 + replicate")              # TOO SLOW TO RUN EVERY TIME..
#glZ.MF <-  glSeq(dm, "-1 + something + replicate")  # TOO SLOW TO RUN EVERY TIME..
#gl.zD  <-  glSeq(dm, "-1 + replicate*something")    # TOO SLOW TO RUN EVERY TIME..