Commit b72cdb5a authored by stevenn's avatar stevenn
Browse files

Diagnostic plots

parent 6ce64a03
GetDataFromBIOM <-function(dataBIOM)
{
GetDataFromBIOM <-function(dataBIOM)
{
counts = biom_data(dataBIOM)
taxo = observation_metadata(dataBIOM)
return(list(counts=counts,taxo=taxo))
}
counts = biom_data(dataBIOM)
taxo = observation_metadata(dataBIOM)
return(list(counts=counts,taxo=taxo))
}
GetDataFromCT <-function(dataC,dataT)
{
counts = dataC
taxo = dataT
return(list(counts=counts,taxo=taxo))
}
GetInteraction2 <- function(target)
{
res=c()
namesTarget = colnames(target)[2:ncol(target)]
k=1
for(i in 1:(length(namesTarget)-1))
GetDataFromCT <-function(dataC,dataT)
{
counts = dataC
taxo = dataT
return(list(counts=counts,taxo=taxo))
}
GetInteraction2 <- function(target)
{
for(j in (i+1):length(namesTarget))
res=c()
namesTarget = colnames(target)[2:ncol(target)]
k=1
for(i in 1:(length(namesTarget)-1))
{
res[k] = paste(namesTarget[i],":",namesTarget[j],sep="")
k = k+1
for(j in (i+1):length(namesTarget))
{
res[k] = paste(namesTarget[i],":",namesTarget[j],sep="")
k = k+1
}
}
return(res)
}
return(res)
}
PrintContrasts <- function (coefs, contrasts,contnames)
{
contrasts = as.matrix(contrasts)
out <-""
for (i in 1:ncol(contrasts))
{
contrast <- contrasts[,i]
contrast <- paste(ifelse(contrast > 0, "+ ", ""), contrast, sep = "")
contrast <- gsub("( 1)|(1)", "", contrast)
out = paste(out,paste("<b>",contnames[i], ":</b> <br/>", paste(contrast[contrast != 0], coefs[contrast != 0], collapse = " ", sep = " ")),"<br/>")
}
return(out)
}
## Get the counts for the selected taxonomy
GetCountsMerge <- function(dataInput,taxoSelect)
{
CT = dataInput$counts
taxo = dataInput$taxo
ordOTU = order(rownames(taxo))
indOTU_annot = which(rownames(CT)%in%rownames(taxo))
counts_annot = CT[indOTU_annot[ordOTU],]
taxoS = taxo[ordOTU,taxoSelect]
counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),sum)
rownames(counts)=counts[,1];counts=counts[,-1]
return(counts)
}
## Get the dds object of DESeq2
Get_dds_object <- function(input,counts,target,design)
{
dds <- DESeqDataSetFromMatrix(countData=counts, colData=target, design=design)
## Size factor estimation
dds <- estimateSizeFactors(dds,locfunc=eval(as.name(input$locfunc)))
dds <- estimateDispersions(dds, fitType=input$fitType)
dds <- nbinomWaldTest(dds)
return(list(dds = dds,counts=counts,target=target,design=design))
}
## Get the design according to the input
GetDesign <- function(input)
{
InterVar = input$InterestVar
Interaction = input$Interaction2
alltmp = c(InterVar,Interaction)
design = as.formula(paste("~", paste0(alltmp, collapse= "+")))
PrintContrasts <- function (coefs, contrasts)
{
contrast <- paste(ifelse(contrast > 0, "+ ", ""), contrast, sep = "")
contrast <- gsub("( 1)|(1)", "", contrast)
out <- paste(contrast[contrast != 0], coefs[contrast != 0], collapse = " ", sep = " ")
return(design)
}
return(out)
}
## Diagnostic Plots
Plot_diag <- function(input,resDiff)
{
colors = c("dodgerblue","firebrick1","MediumVioletRed","SpringGreen")
VarInt = input$VarInt
dds = resDiff$dds
counts = resDiff$counts
target = resDiff$target
group = as.data.frame(target[,VarInt])
if(input$DiagPlot=="barplotTot") barplotTot(counts,group = group, col=colors)
if(input$DiagPlot=="barplotNul") barPlotNul(counts, group = group, col=colors)
if(input$DiagPlot=="densityPlot") densityPlotTot(counts, group = group, col=colors)
if(input$DiagPlot=="MajTax") majTaxPlot(counts, group = group, col=colors)
if(input$DiagPlot=="SERE") SEREplot(counts)
if(input$DiagPlot=="Sfactors") diagSFactors(input,dds,frame=1)
if(input$DiagPlot=="SfactorsVStot") diagSFactors(input,dds,frame=2)
#if(input$DiagPlot=="pcaPlot") barplotTot(counts, group=target[,c(varInt1,varInt2)], col=colors)
}
## barplot total
barplotTot <- function(counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen"))
{
ncol1 <- ncol(group) == 1
barplot(colSums(counts), cex.names = cex.names, main = "Total read count per sample", ylab = "Total read count",
ylim = c(0, max(colSums(counts)) * 1.2), density = if (ncol1) {NULL}
else {15},
angle = if (ncol1) {NULL}
else {c(-45, 0, 45, 90)[as.integer(group[, 2])]}, col = col[as.integer(group[, 1])], las = 2)
legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
if (!ncol1) legend("topleft", levels(group[, 2]), density = 15,col = 1, angle = c(-45, 0, 45, 90)[1:nlevels(group[, 2])], bty = "n")
}
## barplot Nul
barPlotNul <-function (counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen"))
{
percentage <- apply(counts, 2, function(x) {sum(x == 0)}) * 100/nrow(counts)
percentage.allNull <- (nrow(counts) - nrow(removeNul(counts))) * 100/nrow(counts)
ncol1 <- ncol(group) == 1
barplot(percentage, las = 2, col = col[as.integer(group[,1])],
density = if (ncol1) {NULL}
else {15},
angle = if (ncol1) {NULL}
else {c(-45, 0, 45, 90)[as.integer(group[, 2])]},
cex.names = cex.names, ylab = "Proportion of null counts",
main = "Proportion of null counts per sample",
ylim = c(0, 1.2 * ifelse(max(percentage) == 0, 1, max(percentage))))
abline(h = percentage.allNull, lty = 2, lwd = 2)
legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
if (!ncol1) legend("topleft", levels(group[, 2]), density = 15, col = 1, angle = c(-45, 0, 45, 90)[1:nlevels(group[, 2])], bty = "n")
}
## Plot density
densityPlotTot <-function (counts, group, col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen"))
{
counts <- removeNulCounts(counts)
ncol1 <- ncol(group) == 1
plot(density(log2(counts[, 1] + 1)), las = 1, lwd = 2, main = "Density of counts distribution",
xlab = expression(log[2] ~ (raw ~ count + 1)),
ylim = c(0, max(apply(counts, 2, function(x) {max(density(log2(x + 1))$y)})) * 1.05),
lty = if (ncol1) {1}
else{c(1, 2, 3, 4)[as.integer(group[, 2])[1]]},
col = col[as.integer(group[, 1])[1]])
for (i in 2:ncol(counts))
{
lines(density(log2(counts[, i] + 1)), col = col[as.integer(group[,1])[i]], lwd = 2,
lty = if (ncol1) {1}
else {c(1, 2, 3, 4)[as.integer(group[, 2])[i]]})
}
legend("topright", levels(group[, 1]), lty = 1, col = col[1:nlevels(group[,1])], lwd = 2, bty = "n")
if (!ncol1) legend("topleft", levels(group[, 2]), lty = c(1, 2, 3, 4)[1:nlevels(group[, 2])], col = 1, lwd = 2, bty = "n")
}
## Table of maj. taxo
majTab <- function(counts,n)
{
seqnames <- apply(counts, 2, function(x) {
x <- sort(x, decreasing = TRUE)
names(x)[1:n]
})
seqnames <- unique(unlist(as.character(seqnames)))
sum <- apply(counts, 2, sum)
counts <- counts[seqnames, ]
sum <- matrix(sum, nrow(counts), ncol(counts), byrow = TRUE)
p <- round(100 * counts/sum, digits = 3)
return(p)
}
## Plot maj. taxo
majTaxPlot <-function (counts, n = 3, group, cex.names = 1, col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen"))
{
p = majTab(counts,n)
maj <- apply(p, 2, max)
seqname <- rownames(p)[apply(p, 2, which.max)]
ncol1 <- ncol(group) == 1
x <- barplot(maj, col = col[as.integer(group[, 1])], main = "Proportion of reads from\nmost expressed sequence",
ylim = c(0, max(maj) * 1.2), cex.main = 1,
cex.names = cex.names, las = 2, ylab = "Proportion of reads",
density = if (ncol1) {NULL}
else {15},
angle = if (ncol1) {NULL}
else {c(-45, 0, 45, 90)[as.integer(group[, 2])]})
legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
if (!ncol1) legend("topleft", levels(group[, 2]), density = 15, col = 1,
angle = c(-45, 0, 45, 90)[1:nlevels(group[, 2])], bty = "n")
for (i in 1:length(seqname)) text(x[i], maj[i]/2, seqname[i], cex = 0.8, srt = 90, adj = 0)
}
## plot SERE Coefs
SEREplot<-function (counts)
{
sere = SEREcoef(counts)
hc <- hclust(as.dist(sere), method = "ward.D")
plot(hc, las = 2, hang = -1, xlab = "SERE distance, Ward criterion",main = "Cluster dendrogram\non SERE values")
}
## Get the SERE COEF
SEREcoef<-function(counts)
{
sere <- matrix(NA, ncol = ncol(counts), nrow = ncol(counts))
for (i in 1:ncol(counts)) {
for (j in 1:ncol(counts)) {
sere[i, j] <- sigfun_Pearson_meta(counts[, c(i, j)])
}
}
colnames(sere) <- rownames(sere) <- colnames(counts)
sere <- round(sere, digits = 3)
return(sere)
}
## function for the SERE coef
sigfun_Pearson_meta <- function(observed) {
laneTotals <- colSums(observed)
total <- sum(laneTotals)
fullObserved <- observed[rowSums(observed) > 0, ]
fullLambda <- rowSums(fullObserved)/total
fullLhat <- fullLambda > 0
fullExpected <- outer(fullLambda, laneTotals)
fullKeep <- which(fullExpected > 0)
oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/fullExpected[fullKeep]
dfFull <- length(fullKeep) - sum(fullLhat != 0)
sqrt(sum(oeFull)/dfFull)
}
## Plots of size factors
diagSFactors<-function (input,dds,frame=1)
{
geomeans <- exp(rowMeans(log(counts(dds))))
samples <- colnames(counts(dds))
counts.trans <- log2(counts(dds)/geomeans)
xmin <- min(counts.trans[is.finite(counts.trans)], na.rm = TRUE)
xmax <- max(counts.trans[is.finite(counts.trans)], na.rm = TRUE)
if(!is.na(input$NbcolSfactors)) parCols = as.numeric(input$NbcolSfactors)
else parCols = ceiling(ncol(counts.trans)/3)
parRows = ceiling(ncol(counts.trans)/parCols)
if(frame==1)
{
par(mfrow=c(parRows,parCols))
for (j in 1:ncol(dds)) {
hist(log2(counts(dds)[, j]/geomeans), nclass = 100,
xlab = expression(log[2] ~ (counts/geometric ~ mean)), las = 1, xlim = c(xmin, xmax),
main = paste("Size factors diagnostic - Sample ",samples[j], sep = ""), col = "skyblue")
abline(v = log2(sizeFactors(dds)[j]), col = "red", lwd = 1.5)
}
}
if(frame==2)
{
plot(sizeFactors(dds), colSums(counts(dds)), pch = 19, las = 1,
ylab = "Total number of reads", xlab = "Size factors",
main = "Diagnostic: size factors vs total number of reads")
abline(lm(colSums(counts(dds)) ~ sizeFactors(dds) + 0), lty = 2, col = "grey")
}
}
############################################################
##
## CREATE THE CONTRAST DATABASE
##
############################################################
BaseContrast <- function(input,names,namesfile)
{
v_tmp = c()
filesize = file.info(namesfile)[,"size"]
for(i in 1:length(names))
{
Tinput = paste("input$",names[i],sep="")
expr=parse(text=Tinput)
val = eval(expr)
v_tmp[i] = as.numeric(val)
}
if(filesize!=0)
{
oldContrast = read.table(namesfile,header=TRUE)
colnamesTmp = c(colnames(oldContrast),input$ContrastName)
mat = cbind(oldContrast,v_tmp)
}
else{ colnamesTmp = input$ContrastName; mat = v_tmp}
write.table(mat,namesfile,row.names=FALSE,col.names = colnamesTmp)
}
## Remove nul counts
removeNulCounts <-function (counts)
{
return(counts[rowSums(counts) > 0, ])
}
......
......@@ -17,8 +17,11 @@ shinyServer(function(input, output,session) {
##
#####################################################
namesfile = "www/BaseContrast.txt"
#file.create(namesfile,showWarnings=FALSE)
## Create base for contrast
rand = floor(runif(1,0,1e9))
namesfile = paste("www/base/BaseContrast_",rand,".txt",sep="")
file.create(namesfile,showWarnings=FALSE)
## Counts file
dataInputCounts <-reactive({
......@@ -61,6 +64,14 @@ shinyServer(function(input, output,session) {
## Rownames
rownames(data)=data[,1];data=data[,-1]
## Add NA
data=as.matrix(data)
indNa = which(data=="")
data[indNa]=NA
return(as.data.frame(data))
})
......@@ -157,7 +168,7 @@ shinyServer(function(input, output,session) {
))
## Tab box for visualisation
## Tab box for data visualisation
output$TabBoxData <- renderUI({
data=dataInput()
......@@ -189,28 +200,11 @@ shinyServer(function(input, output,session) {
data = read.csv(inFile$datapath,sep="\t",header=TRUE)
return(as.data.frame(data))
rownames(target) <- as.character(target[, 1])
return((data))
})
# Infobox design
output$RowTarget <- renderInfoBox({
target = dataInputTarget()
InterVar = input$InterestVar
Interaction = input$Interaction2
alltmp = c(InterVar,Interaction)
if(!is.null(target))
{
#### Ajout fontion check target
infoBox(h6(strong("Target format")), subtitle = h6("Your target file is OK"), icon = icon("thumbs-o-up"),color = "green",width=NULL,fill=TRUE)
}
else infoBox(h6(strong("Warning")), subtitle = h6("Label of the target file must correspond to counts table column names") ,color = "orange",width=NULL,fill=TRUE, icon = icon("warning"))
})
## Interest Variables
output$SelectInterestVar <- renderUI({
......@@ -263,20 +257,58 @@ shinyServer(function(input, output,session) {
})
# Infobox design
output$RowTarget <- renderInfoBox({
target = dataInputTarget()
if(!is.null(target))
{
#### Ajout fontion check target
infoBox(h6(strong("Target file")), subtitle = h6("Your target file is OK"), icon = icon("thumbs-o-up"),color = "green",width=NULL,fill=TRUE)
}
else infoBox(h6(strong("Target file")), subtitle = h6("Label of the target file must correspond to counts table column names") ,color = "orange",width=NULL,fill=TRUE, icon = icon("warning"))
})
## taget table
output$DataTarget <- renderDataTable(
dataInputTarget(),
options = list(lengthMenu = list(c(10, 50, -1), c('10', '50', 'All')),
pageLength = 10,scrollX=TRUE
))
## Box for target visualisation
output$BoxTarget <- renderUI({
target = dataInputTarget()
if(!is.null(target) && nrow(target)>0)
{
box(title="Target file overview",width = NULL, status = "primary", solidHeader = TRUE,collapsible = TRUE,collapsed = TRUE,
dataTableOutput("DataTarget")
)
}
})
#####################################################
##
## DEFINE CONTRAST
##
#####################################################
output$contrastMat <- renderUI({
#dds=RunDESeq2()
#names = resultsNames(dds)
names=c('test1',"test2","test3")
resDiff = ResDiffAnal()
dds = resDiff$dds
names = resultsNames(dds)
Contrast=list()
for(i in 1:length(names)){Contrast[[i]] = textInput(names[i],names[i],0)}
......@@ -286,95 +318,248 @@ shinyServer(function(input, output,session) {
})
output$ContrastOverview <- renderUI({
output$ContrastOverview <- renderPrint({
resDiff = ResDiffAnal()
dds = resDiff$dds
names = resultsNames(dds)
#dds=RunDESeq2()
#names = resultsNames(dds)
names=c('test1',"test2","test3")
cont = input$ContrastList
ContrastBase = read.table(namesfile,header=TRUE)
filesize = file.info(namesfile)[,"size"]
res = PrintContrasts(names,ContrastBase[,cont])
return(res)
if(filesize!=0)
{
ContrastBase = read.table(namesfile,header=TRUE)
ind = which(colnames(ContrastBase)%in%cont)
div(HTML(PrintContrasts(names,sapply(ContrastBase[,ind],as.numeric),cont)))
}
})
## Add contrast function
AddCont <-eventReactive(input$AddContrast,{
resDiff = ResDiffAnal()
dds = resDiff$dds
names = resultsNames(dds)
BaseContrast(input,names,namesfile)
tmp = read.table(namesfile,header=TRUE)
Contrast = colnames(as.matrix(tmp))
updateSelectInput(session, "ContrastList","Contrasts",Contrast)
})
BaseContrast <- function(input,namesfile)
{
#dds=RunDESeq2()
#names = resultsNames(dds)
## Add contrast
observeEvent(input$AddContrast,{
oldContrast = read.table(namesfile,header=TRUE)
names=c('test1',"test2","test3")
v_tmp = c()
AddCont()
for(i in 1:length(names))
{
Tinput = paste("input$",names[i],sep="")
print(Tinput)
print(input$test1)
})
## Remove contrast function
RemoveCont <-eventReactive(input$RemoveContrast,{
## get the size of the contrast base file
filesize = file.info(namesfile)[,"size"]
if(filesize!=0)
{
tmp = read.table(namesfile,header=TRUE)
ind = which(colnames(tmp)%in%input$ContrastList)
matKept = as.matrix(tmp[,-ind])
ContrastKept = colnames(tmp)[-ind]
expr=parse(text=Tinput)
val = eval(expr)
print(val)
v_tmp[i] = as.numeric(val)
if(ncol(matKept)>0) write.table(matKept,namesfile,row.names=FALSE,col.names = ContrastKept)
else file.create(namesfile,showWarnings=FALSE)