Skip to content

GitLab

Unverified Commit 49ea5e4e authored by Amine Ghozlane's avatar Amine Ghozlane Committed by GitHub
Browse files

Merge pull request #2 from pworinger/master 

Add plots to analyse the differential abundances and new visualisations
parents 62b4a71a cdddde0a
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2732 additions and 1307 deletions
LoadPackages.R
View file @ 49ea5e4e
options(download.file.method = 'wget', bitmapType='cairo')
#options(download.file.method = 'wget', bitmapType='cairo')
options(bitmapType='cairo')
if (!require("Rcpp")){
install.packages("Rcpp", repos="http://cran.univ-paris1.fr/")
......@@ -65,11 +66,11 @@ if(!require(treeWeightD3)){
devtools::install_github('pierreLec/treeWeightD3')
library(treeWeightD3)
}
if (!require(BiocInstaller)){
source("https://bioconductor.org/biocLite.R")
biocLite("BiocInstaller")
library(BiocInstaller)
}
# if (!require(BiocInstaller)){
# source("https://bioconductor.org/biocLite.R")
# biocLite("BiocInstaller")
# library(BiocInstaller)
# }
if (!require(d3heatmap)) {
devtools::install_github('aghozlane/d3heatmap')
......@@ -83,7 +84,8 @@ library(biomformat)
}
if (!require(scatterD3)) {
devtools::install_github('aghozlane/scatterD3')
#devtools::install_github('aghozlane/scatterD3')
install.packages("scatterD3")
library(scatterD3)
}
......@@ -109,7 +111,7 @@ if (!require(shinyjs)) {
}
if(!require(d3vennR)){
install_github("timelyportfolio/d3vennR")
devtools::install_github("timelyportfolio/d3vennR")
library(d3vennR)
}
......@@ -129,8 +131,7 @@ if (!require(gplots)) {
}
if (!require(DESeq2)) {
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
devtools::install_github("aghozlane/DESeq2shaman")
library(DESeq2)
}
......@@ -223,4 +224,47 @@ if (!require("htmltools")){
# library(plotly)
# }
if (!require("rAmCharts")){
install.packages("rAmCharts")
library(rAmCharts)
}
if(!require("colourpicker")){
install.packages("colourpicker")
library(colourpicker)
}
if(!require("data.table")){
install.packages("data.table")
library(data.table)
}
if(!require("UpSetR")){
install.packages("UpSetR")
library(UpSetR)
}
if(!require("ggrepel")){
install.packages("ggrepel")
library(ggrepel)
}
# if(!require("networkD3")){
# install.packages("networkD3")
# library(networkD3)
# }
if(!require("igraph")){
install.packages("igraph")
library(igraph)
}
if(!require("visNetwork")){
install.packages("visNetwork")
library(visNetwork)
}
if(!require("SummarizedExperiment")){
BiocManager::install("SummarizedExperiment")
library(SummarizedExperiment)
}
\ No newline at end of file
Rfunctions/CompPlot.R
View file @ 49ea5e4e
......@@ -5,20 +5,22 @@
##############################
## HEATMAP
##############################
Plot_Visu_Heatmap_FC <- function(input,BaseContrast,resDiff,export=FALSE){
Plot_Visu_Heatmap_FC <- function(input, BaseContrast, resDiff, ContrastListDebounce, SelectTaxoPlotCompDebounce, export=FALSE){
res = NULL
SelContrast = input$ContrastList_table_FC
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
colnames(selcontrast_matrix) = SelContrast
log2FC = Get_log2FC_padj(input,selcontrast_matrix,resDiff, info = NULL)$log2FC
if(!is.null(log2FC) && length(SelContrast)>=2)
{
cont = which(colnames(log2FC)%in%SelContrast)
log2FC = log2FC[,SelContrast]
ind_taxo = input$selectTaxoPlotComp
log2FC = log2FC[,SelContrast]
###
ind_taxo = SelectTaxoPlotCompDebounce()
#ind_taxo = input$selectTaxoPlotComp
ind = rownames(log2FC)%in%ind_taxo
log2FC = as.matrix(log2FC[ind,])
......@@ -41,21 +43,260 @@ Plot_Visu_Heatmap_FC <- function(input,BaseContrast,resDiff,export=FALSE){
return(res)
}
##############################
## P VALUE DENSITY PLOT
##############################
Plot_pValue_Density <- function(input, BaseContrast, resDiff, ContrastListDebounce, alphaVal){
res = NULL
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
if(length(SelContrast)>=1){
selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
colnames(selcontrast_matrix) = SelContrast
padj = data.frame(na.omit(Get_log2FC_padj(input,selcontrast_matrix,resDiff)$padj))
data <- data.frame()
for (cont in SelContrast){
data_cont <- data.frame(padj[,cont])
data_cont$contrast <- cont
colnames(data_cont) <- c("padj","contrast")
data <- rbind(data, data_cont)
}
data$contrast <- factor(data$contrast, levels = SelContrast)
# uncomment to zoom (X axis between 0 and input$alphaVal) (use "+ xlim(0,as.numeric(alphaVal))" instead, to redraw geom_density with only data under threshold)
p <- ggplot(data, aes(x = padj, color = contrast, fill = contrast)) + geom_density(alpha = input$fillOpacity, size = input$lineWidth) + theme_minimal() #+ coord_cartesian(xlim = c(0,as.numeric(alphaVal)))
p <- p + theme(axis.title = element_text(size = input$FontSizepValueDensity),
axis.text = element_text(size = input$FontSizepValueDensity),
legend.title = element_text(size = input$FontSizepValueDensity),
legend.text = element_text(size = input$FontSizepValueDensity))
res = p
}
return(res)
}
##############################
## LOGIT PLOT
##############################
Plot_Comp_Logit <- function(input, BaseContrast, resDiff, SelectTaxoPlotCompDebounce, export = FALSE)
{
plot = NULL
if(input$Contrast1 != input$Contrast2 && input$Contrast1 != "..." && input$Contrast2 != "..."){
SelContrast = c(input$Contrast1, input$Contrast2)
selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
colnames(selcontrast_matrix) = SelContrast
padj = na.omit(Get_log2FC_padj(input,selcontrast_matrix,resDiff)$padj)
#data to plot
x <- padj[,input$Contrast1]
logit_x <- log(x/(1-x))
y <- padj[,input$Contrast2]
logit_y <- log(y/(1-y))
#classification for colours
alpha_logit <- log(as.numeric(input$AlphaVal)/(1-as.numeric(input$AlphaVal)))
x_signif <- logit_x < alpha_logit
y_signif <- logit_y < alpha_logit
matrix <- cbind(x_signif, y_signif)
color_var <-
apply(matrix, 1, function(point)
if (point[[1]]) {
if (point[[2]]) {"Significant for both contrasts"}
else {"Significant for contrast 1"}
}
else{
if (point[[2]]) {"Significant for contrast 2"}
else{"Not significant"}
}
)
#labels
names <- row.names(padj)
###
points_to_label <- SelectTaxoPlotCompDebounce()
#points_to_label <- input$selectTaxoPlotComp
labels <-
sapply(names, function(name)
if (is.element(name, points_to_label)) {
points_to_label[match(name, points_to_label)]
} else{
""
})
if(!export){
plot <- scatterD3(x = logit_x,
y = logit_y,
xlab = paste("logit p value", input$Contrast1),
ylab = paste("logit p value", input$Contrast2),
hover_opacity = 1,
tooltip_text = row.names(padj),
# xlim = xlimits,
# ylim = ylimits,
lines = if (input$showSignifThresholdsLogitPlot) {
if(input$showDiagonal){
data.frame(
slope = c(Inf, 0, 1),
intercept = c(alpha_logit,alpha_logit,0),
stroke = "black",
stroke_width = input$linesWidthLogitPlot,
stroke_dasharray = 5
)
}
else{
data.frame(
slope = c(Inf, 0),
intercept = c(alpha_logit,alpha_logit),
stroke = "black",
stroke_width = input$linesWidthLogitPlot,
stroke_dasharray = 5
)}
}
else{if (input$showDiagonal){
data.frame(
slope = c(1),
intercept = c(0),
stroke = "black",
stroke_width = input$linesWidthLogitPlot,
stroke_dasharray = 5
)}
else{NULL}},
fixed = input$fixed11,
col_var = color_var,
col_lab = if(input$showLegendLogit){"Legend"}else{NA}
,
colors = c(
"Not significant" = input$colour01,
"Significant for contrast 1" = input$colour02,
"Significant for contrast 2" = input$colour03,
"Significant for both contrasts" = input$colour04
),
lab = labels,
point_opacity = input$pointOpacityLogit,
point_size = input$pointSizeLogit,
axes_font_size = paste(input$axisFontSizeLogit, "0%", sep = ""),
labels_size = input$labelsSizeLogit,
legend_font_size = paste(input$legendFontSizeLogit, "0%", sep = ""),
legend_width = input$legendWidthLogit,
# dom_id_reset_zoom = "scatterD3-reset-zoomLogit",
# dom_id_svg_export = "scatterD3-svg-exportLogit",
# menu = FALSE,
# height = input$heightVolcanoPlot,
# width=if(input$modifwidthVolcano){input$widthVolcanoPlot},
# disable_wheel = TRUE
)
}
if(export){
data <- data.frame(logit_x, logit_y, color_var)
plot <- ggplot(data, aes(x = logit_x, y = logit_y)) + theme_minimal()
plot <- plot + geom_point(aes(color = color_var), size = (input$pointSizeLogit %/% 30), alpha = input$pointOpacityLogit)
plot <- plot + scale_x_continuous(paste("logit p value", input$Contrast1))
plot <- plot + scale_y_continuous(paste("logit p value", input$Contrast2))
plot <- plot + theme(axis.title = element_text(size = input$axisFontSizeLogit),
axis.text = element_text(size = input$axisFontSizeLogit),
legend.title = element_text(size = input$legendFontSizeLogit),
legend.text = element_text(size = input$legendFontSizeLogit))
plot <- plot + geom_text(aes(label = labels, color = color_var), size = input$labelsSizeLogit %/% 3, hjust = 1, vjust = 1, show.legend = FALSE)
plot <- plot + scale_color_manual(name = "Legend title", values = c("Not significant" = input$colour01, "Significant for contrast 1" = input$colour02,
"Significant for contrast 2" = input$colour03, "Significant for both contrasts" = input$colour04))
if (input$fixed11) {
plot <- plot + coord_fixed()
}
if (input$showSignifThresholdsLogitPlot) {
plot <- plot + geom_hline(yintercept = alpha_logit, linetype = "dashed", size = input$linesWidthLogitPlot/2)
plot <- plot + geom_vline(xintercept = alpha_logit, linetype = "dashed", size = input$linesWidthLogitPlot/2)
}
if(input$showDiagonal){
plot <- plot + geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = input$linesWidthLogitPlot/2)
}
}
return(plot)
}}
##############################
## VENN DIAGRAM
##############################
Plot_Visu_Venn <- function(input,BaseContrast,resDiff,export=FALSE){
Plot_Visu_Venn <- function(input,BaseContrast,resDiff, ContrastListVennDebounce, export=FALSE){
res = NULL
SelContrast = input$ContrastList_table_FC
#SelContrast = input$ContrastList_table_FCVenn
SelContrast = ContrastListVennDebounce()
if(length(SelContrast)>=2 & length(SelContrast)<=4){
data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res
res = venn_tooltip(d3vennR(data=data))
}
return(res)
}
##############################
## UPSET
##############################
Plot_UpSet <- function(input,BaseContrast, resDiff, ContrastListDebounce, export=FALSE){
plot = NULL
df = NULL
#SelContrast = input$ContrastList_table_FC
SelContrast = ContrastListDebounce()
if(length(SelContrast)>=2){
data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$df.tot
listInput <- list()
n <- ncol(data)
for(i in 1:n){
new <- list(na.omit(data[,i]))
listInput <- c(listInput, new)
}
names(listInput) <- colnames(data)
# for plot
plot <- upset(fromList(listInput), nsets = n, nintersects = NA, set_size.show = input$showNumbers,
show.numbers = if(input$showNumbers){"yes"}else{"no"}, sets = SelContrast, keep.order = !(input$orderBySize),
point.size = input$pointSizeUpSet, line.size = 1.5, order.by = input$orderByUpset,
mainbar.y.label = "Number of differential features in common", sets.x.label = "Number of differential features",
text.scale = c(2, 1.6, 1.6, 1.6, 2, 1.5), decreasing = input$decreasingUpset)
# for table
combos <- Reduce(c,lapply(rev(1:length(listInput)), function(x) combn(rev(colnames(data)),x,simplify=FALSE) ))
listSignif <- Reduce(function(a,b) {unique(c(as.character(a),as.character(b)))}, listInput)
n_col = length(combos)
n_row = length(listSignif)
df <- matrix(rep(listSignif, times = n_col), ncol = n_col)
lst <- list()
for(j in 1:n_col){for(i in 1:n_row){group_contrasts <- unlist(combos[j])
if (is.element(df[i,j], lst)){df[i,j] <- NA}
else{if (!is.element(df[i,j], Reduce(intersect, listInput[group_contrasts]))){df[i,j] <- NA}
else{lst <- c(lst,df[i,j])}
}
}
}
colnames(df) <- lapply(combos, function(lstCont) paste(lstCont, collapse = " & "))
df = as.data.frame(apply(df,2,Go_data_top))
maxRow = max(apply(df,2,FUN=function(x) length(which(!is.na(x)))))
df = df[1:max(maxRow,1),]
df = df[,which(apply(!is.na(df),2,any))]
}
return(list(plot=plot,table=df))
}
##############################
## MULTIPLE VENN COMPAIR
##############################
Plot_MultipleVenn <- function(input,BaseContrast, resDiff, ContrastListDebounce){
plot = NULL
SelContrast = ContrastListDebounce()
if(length(SelContrast)>=2){
data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res_multiple_venn
plot <- ggplot(data,aes(x=x,y=y)) + geom_point() + theme_bw() +
geom_label_repel(aes(label=name), size = input$labelSizemultipleVenn) + xlim(c(0,1)) + ylim(c(0,1)) +
xlab(bquote(Contrast1 *intersect(Contrast2) ~ "/ Contrast1")) +
ylab(bquote(Contrast1 *intersect(Contrast2) ~ "/ Contrast2"))
plot <- plot + theme(axis.title = element_text(size = input$FontSizeMultipleVenn),
axis.text = element_text(size = input$FontSizeMultipleVenn - 2))
}
return(plot)
}
##############################################################
......@@ -77,6 +318,7 @@ Go_data_top <- function(vect)
## Get the data to plot the heatmap with the FC
## and for the logit plot
Get_log2FC_padj <-function(input,BaseContrast,resDiff, info = NULL)
{
log2FC = NULL
......@@ -100,17 +342,20 @@ Get_log2FC_padj <-function(input,BaseContrast,resDiff, info = NULL)
independentFiltering=input$IndFiltering,alpha=alpha)
}
log2FC = as.matrix(round(result[[SelContrast[1]]][, "log2FoldChange"], 3))
padj = as.matrix(round(result[[SelContrast[1]]][, "padj"], 3))
#padj = as.matrix(round(result[[SelContrast[1]]][, "padj"], 3))
padj = as.matrix((result[[SelContrast[1]]][, "padj"]))
if(nbCont>1)
{
for(i in 2:nbCont)
{
log2FC = cbind(log2FC,round(result[[SelContrast[i]]][, "log2FoldChange"], 3))
padj = cbind(padj,round(result[[SelContrast[i]]][, "padj"], 7))
#padj = cbind(padj,round(result[[SelContrast[i]]][, "padj"], 7))
padj = cbind(padj,(result[[SelContrast[i]]][, "padj"]))
}
}#
colnames(log2FC) = names(result)
colnames(padj) = names(result)
}
#}
rownames(log2FC) = rownames(result[[SelContrast[1]]])
rownames(padj) = rownames(result[[SelContrast[1]]])
......@@ -188,6 +433,7 @@ venn_tooltip <- function(venn){
GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
{
res = list()
res_multiple_venn = NULL
df.tot = NULL
VarInt = input$VarInt
dds = resDiff$dds
......@@ -208,11 +454,12 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
independentFiltering=input$IndFiltering,alpha=alpha)
}
padj = round(result[[SelContrast[1]]][, "padj"], 3)
# save(result,padj,SelContrast,file = "test1.RData")
df = as.matrix(rownames(result[[1]]))
if(length(which(padj>alpha))>0) df[which(padj>alpha),]=NA
if(any(is.na(padj))) df[which(is.na(padj)),]=NA
if(length(which(padj>alpha))>0) df[which(padj>alpha),]=NA
if(any(is.na(padj))) df[which(is.na(padj)),]=NA
if(nbCont>1)
{
for(i in 2:nbCont)
......@@ -226,32 +473,55 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
colnames(df) = SelContrast
df = as.data.frame(df)
}
## Keep the entire dataframe
df.tot = as.data.frame(apply(df,2,Go_data_top))
maxRow = max(apply(as.data.frame(apply(df,2,Go_data_top)),2,FUN=function(x) length(which(!is.na(x)))))
df.tot = df.tot[1:max(maxRow,1),]
## Remove col with only NA
df = df[,which(apply(!is.na(df),2,any))]
ncont = ncol(as.data.frame(df))
names.df = names(df)
rownames_multiple_venn = c()
cmp=1
if(ncont>1 && !is.null(ncont))
{
for(i in 1:(ncont))
{
for(j in i:ncont)
{
if(i!=j) res[[cmp]] = list(sets=list(names.df[i],names.df[j]),size= length(which(!is.na(intersect(df[,i],df[,j])))))
if(i==j) res[[cmp]] = list(sets=list(names.df[i]),size= length(which(!is.na(intersect(df[,i],df[,i])))))
cmp=cmp+1
# For Venn diagram (until intersection of 4 contrasts)
{ if(i==j) {res[[cmp]] = list(sets=list(names.df[i]),size= length(which(!is.na(intersect(df[,i],df[,i])))))
cmp=cmp+1}
if(i!=j) {res[[cmp]] = list(sets=list(names.df[i],names.df[j]),size= length(which(!is.na(intersect(df[,i],df[,j])))))
cmp=cmp+1
for(k in j:ncont){
if(k!=j){
res[[cmp]] = list(sets=list(names.df[i],names.df[j],names.df[k]), size = length(which(!is.na(intersect(intersect(df[,i],df[,j]),df[,k])))))
cmp=cmp+1
for(l in k:ncont){
if(l!=k){
res[[cmp]] = list(sets=list(names.df[i],names.df[j],names.df[k],names.df[l]), size = length(which(!is.na(intersect(intersect(intersect(df[,i],df[,j]),df[,k]),df[,l])))))
cmp=cmp+1}
}
}
}
}
# For 'Multiple Venn Compair'
if(i!=j) {
if (is.null(res_multiple_venn))
{res_multiple_venn <- data.frame(name = c(rownames_multiple_venn, paste(names.df[i],names.df[j], sep = " vs ")),
x = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,i],df[,i]))))),
y = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,j],df[,j]))))),
stringsAsFactors=FALSE)}
else{new_row <- data.frame(name = c(rownames_multiple_venn, paste(names.df[i],names.df[j], sep = " vs ")),
x = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,i],df[,i]))))),
y = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,j],df[,j]))))))
res_multiple_venn <- rbind(res_multiple_venn, new_row)}
}
}
}
}
}
return(list(res=res,df.tot=df.tot))
return(list(res=res,df.tot=df.tot,res_multiple_venn=res_multiple_venn))
}
Rfunctions/DiagPlot.R
View file @ 49ea5e4e
......@@ -138,7 +138,7 @@ HCPlot <- function (input,dds,group,type.trans=NULL,counts=NULL,CT,tree,col = c(
## Get the type of dendrogram
type <- input$typeHculst
dend <- set(dend, "labels_cex", input$cexLabelDiag)
#dend <- set(dend, "labels_cex", input$cexLabelDiag)
if(input$colorHC) labels_colors(dend)<-col[as.integer(as.factor(group))][order.dendrogram(dend)]
if(type=="hori")
{
......@@ -148,7 +148,8 @@ HCPlot <- function (input,dds,group,type.trans=NULL,counts=NULL,CT,tree,col = c(
else
{
par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
res = circlize_dendrogram(dend, labels_track_height = 0.2, dend_track_height = .3, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "))
res = circlize_dendrogram(dend, labels_track_height = 0.2, dend_track_height = .3, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "),
labels_cex = input$cexLabelDiag)
}
}
}
......@@ -554,8 +555,8 @@ PCAPlot_meta <-function(input,dds, group_init, n = min(500, nrow(counts(dds))),
{
type.trans <- type.trans[1]
if (type.trans == "VST") counts.trans <- assay(varianceStabilizingTransformation(dds))
else counts.trans <- assay(rlogTransformation(dds))
if (type.trans == "VST") counts.trans <- SummarizedExperiment::assay(varianceStabilizingTransformation(dds))
else counts.trans <- SummarizedExperiment::assay(rlogTransformation(dds))
counts.trans = counts.trans[,ind_kept]
rv = apply(counts.trans, 1, var, na.rm = TRUE)
......@@ -754,8 +755,8 @@ Get_pca_table <-function(input,dds, group_init, n = min(500, nrow(counts(dds))),
{
type.trans <- type.trans[1]
if (type.trans == "VST") counts.trans <- assay(varianceStabilizingTransformation(dds))
else counts.trans <- assay(rlogTransformation(dds))
if (type.trans == "VST") counts.trans <- SummarizedExperiment::assay(varianceStabilizingTransformation(dds))
else counts.trans <- SummarizedExperiment::assay(rlogTransformation(dds))
counts.trans = counts.trans[,ind_kept]
rv = apply(counts.trans, 1, var, na.rm = TRUE)
......
Rfunctions/Stat_Model.R
View file @ 49ea5e4e
......@@ -47,22 +47,22 @@ Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm,C
# }
# countsNorm = counts(dds, normalized = TRUE)
if(as.numeric(R.Version()$major)>=3 && as.numeric(R.Version()$minor) >=1.3){
if(nrow(counts)*nrow(target)>=50000)
{
try(DESeq(dds,fitType=input$fitType,parallel = TRUE,minReplicatesForReplace = Inf) ->ddstmp,silent = TRUE)
if(!is.null(ddstmp)) dds = ddstmp
}
if(nrow(counts)*nrow(target)<50000 || is.null(ddstmp))
{
dds <- estimateDispersions(dds, fitType=input$fitType)
dds <- nbinomWaldTest(dds)
}
}else{
# if(as.numeric(R.Version()$major)>=3 && as.numeric(R.Version()$minor) >=1.3){
# if(nrow(counts)*nrow(target)>=50000)
# {
# try(DESeq(dds,fitType=input$fitType,parallel = TRUE,minReplicatesForReplace = Inf) ->ddstmp,silent = TRUE)
# if(!is.null(ddstmp)) dds = ddstmp
# }
# if(nrow(counts)*nrow(target)<50000 || is.null(ddstmp))
# {
# dds <- estimateDispersions(dds, fitType=input$fitType)
# dds <- nbinomWaldTest(dds)
# }
# }else{
dds <- estimateDispersions(dds, fitType=input$fitType)
dds <- nbinomWaldTest(dds,modelMatrixType = "expanded")