Nothing
# plot the odds-ratio and the 95%-CI from Fishers excact test (two-sided)
# show A,B,C,D counts in pie-charts with size dependent on total nr. of genes
# take table with sums of scores per node from plot_anno_scores
# TODO: test plotted odds ratio with an example
plot_conti = function(aggrego){
# don't need root node for that
aggrego = aggrego[,1:5] # TODO: maybe avoid adding root-info in plot_anno_scores in the first place?
### perform fishers exact test for every GO
fish_odds = data.frame(aggrego, t(apply(aggrego, 1, fisher_conti)))
colnames(fish_odds)[c(2:9)] = c("A","B","C","D","odds_ratio","ci95_low","ci95_high","p")
out = fish_odds
# replace Inf or -Inf CI for plotting
inf_go = fish_odds[is.infinite(fish_odds$odds_ratio),"go_id"]
if(length(inf_go) > 0){
warning("The following GO-categories have an infinite odds ratio: ", paste(inf_go, collapse=", "))
}
zero_go = fish_odds[fish_odds$odds_ratio==0,"go_id"]
if(length(zero_go) > 0){
warning("The following GO-categories have an odds ratio of 0: ", paste(zero_go, collapse=", "))
}
############# TODO influenced by setting par(oma) before function call?
### plot the odds-ratios, with CI
op = par(no.readonly = TRUE)
# 3 panels
layout(matrix(c(1:2),ncol=1), widths=c(5),heights=c(3,3))
# odds-ratio
par(mar=c(0.5,4,3,2), bty="l") #, bty="n")
ymin = min(c(1,fish_odds$ci95_low[is.finite(fish_odds$ci95_low) & fish_odds$odds_ratio!=0]))
ymax = max(c(1,fish_odds$ci95_high[is.finite(fish_odds$ci95_high) & fish_odds$odds_ratio!=0]))
suppressWarnings(plot(fish_odds$odds_ratio, pch=19, ylab="", xaxt="n", xlab="", main="odds ratio (A/B) / (C/D)", xlim=c(0.5,nrow(fish_odds)+1), ylim=c(ymin,ymax), cex.axis=0.9, log="y", las=2,
panel.first={grid(0, NULL, lty=1, col=colors()[2])}))
# 95%-CI
suppressWarnings(arrows(c(1:nrow(fish_odds),1:nrow(fish_odds)),c(fish_odds$ci95_high,fish_odds$ci95_low), c(1:nrow(fish_odds),1:nrow(fish_odds)), c(fish_odds$odds_ratio, fish_odds$odds_ratio), angle=90, code=1, length=0.03))
# horizontal line at 1
abline(h=1, col="#F15A60")
axis(4,at=1,labels="1", col="#F15A60", col.axis="#F15A60", las=1, cex.axis=0.9)
# pie charts
pie_cols = colors()[c(124,132,59,137)]
# pie_cols = colors()[c(592,566,503,556)]
par(mar=c(5.5,4,3,2), bty="l")
plot(1,xlim=c(0.5,nrow(fish_odds)+1), ylim=c(0,2),type="n", main="annotated genes", xlab="", ylab="", xaxt="n", yaxt="n")
radi_units = 0.4/log(max(rowSums(fish_odds[,2:5]))+1)
for(i in seq_len(nrow(fish_odds))){
add.pie(z=as.numeric(fish_odds[i,2:3]), x=i, y=1.6, radius=log(sum(fish_odds[i,2:3])+1)*radi_units, labels="", col=pie_cols[1:2])
add.pie(z=as.numeric(fish_odds[i,4:5]), x=i, y=0.6, radius=log(sum(fish_odds[i,4:5])+1)*radi_units, labels="", col=pie_cols[3:4])
}
legend("right", fill=pie_cols, bty="n", legend=c("A","B","C","D"))
text(x=1:nrow(fish_odds), y=1.9, labels=paste0("n=",rowSums(fish_odds[,2:5])), xpd=TRUE, cex=0.8, pos=3, offset=0.7)
axis(1, at=1:nrow(fish_odds), labels=FALSE)
text(x=1:nrow(fish_odds), y=-0.35, labels=fish_odds$go_id, srt=45, adj=1, xpd=TRUE, cex=0.85)
par(op)
return(invisible(out))
}
fisher_conti = function(fish_line){
# contingency table
conti = matrix(as.numeric(fish_line[2:5]), ncol=2)
fish = fisher.test(conti)
out = c(fish$estimate, fish$conf.int, fish$p.value)
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.