#下面的方法不能安装
source("https://bioconductor.org/biocLite.R")
biocLite("cqn")

#Install MDSeq from local source with
#install.packages("MDSeq_1.0.5.tar.gz", repos=NULL, type="source")

#Install MDSeq from GitHub with

library(devtools)
install_github("zjdaye/MDSeq")

data(sampleData)
# raw count table 
dat <- sample.exprs 
dim(dat)
dat[1:6, 1:6] 
# covariates 
X <- sample.pheno[,c("X1","X2")] 
head(X)
# group information 
group <- sample.pheno$group 
#键入summary(lmfit)将显示分析结果的统计概要，
summary(factor(group)) 
# lowly expressed genes were filtered by mean cpm value across all samples 
dat.filtered <- filter.counts(dat, mean.cpm.cutoff = 0.1)
dim(dat.filtered)
# Using normalized counts 
dat.normalized <- normalize.counts(dat.filtered, group=group, method="TMM") 
# including normalization factor as an offset in the mean-dispersion GLMs 
cnf <- calcNormFactors(dat.filtered, method="TMM") 
libsize <- colSums(dat.filtered) #normalization factor 
rellibsize <- libsize/exp(mean(log(libsize))) #relative library size 
nf <- cnf * rellibsize #normalization factor including library size
get.model.matrix(factor(rep(1:3,each=4))) 
get.model.matrix(factor(rep(1:3,each=4)), contrast.type = "contr.treatment")
# treatment group assignment 
group <- factor(sample.pheno$group, labels = c("Control", "Case")) 
table(group) 
# make design matrix with proper contrast setting 
groups <- get.model.matrix(group) 
groups 

# Check outliers using parallel process with 4 threads 
# For the sake of simplicity, we demonstrate outlier detection 
# by using the first 1,000 genes.
dat.checked <- remove.outlier(dat.normalized[1:1000, ], X=X, U=X, contrast = groups, mc.cores = 4) 
head(dat.checked$outliers) 

# status of outlier checking 
table(dat.checked$outliers$status) 
# frequency distribtuion of outliers 
table(dat.checked$outliers$num.outliers) 

# remove genes with status flag other than 0 
counts <- dat.checked$count[dat.checked$outliers$status==0,] 
dim(counts) 

# using parallel process with 4 threads 
fit <- MDSeq(counts, X=X, U=X, contrast = groups, mc.cores = 4) 

# simultaneously test zero-inflation 
head(fit$Dat[c("s","ZI.pval", "ZI")], 20)

# given log2-fold-change threshold = 1 
result <- extract.ZIMD(fit, compare = list(A="Case", B="Control"), log2FC.threshold = 1) 
head(result) 

sessionInfo()