R/plotCommunity.r

In cmfield/phylloR: R package for phyllosphere analysis

#' Function to produce a bar or violin plot of an individual community.
#'
#' @param counts A count table in which the rows are strains or OTUs and the columns are samples. The table should include only a single sample group.
#' @param type Determines the plot style, either 'bar', 'violin', 'points', 'swarm', 'barswarm' or 'violinswarm'.
#' @param xlabels  An optional vector of strain names, the default is to use the row names of the count table.
#' @param xcols An optional vector of strain colours, the default is to use rainbow(). 
#' @param res The resolution of the histograms used to create the violins for that plot style.
#' @param cutoff A numeric integer, which if provided reduces the plot to only the given number of most abundant strains.
#' @details
#' For each strain, samples in which they were not observed are counted and separated from the main bar or violin.
#' As a result, the number of samples in each bar or violin varies, and the histogram for each violin is therefore normalised.
#' @keywords phylloR
#' @return A list of the statistics for each strain as generated by boxplot().
#' @export
#' @author Chris Field <fieldc@@ethz.ch>
#' @examples
#' None

plotCommunity <- function(counts,type="bar",xlabels=NULL,xcols=NULL,res=50,cutoff=NULL){
    if(!type%in%c("bar","violin","points","swarm","barswarm","violinswarm")){
        cat("Not a valid type of community plot\n")
        return()
    }
    if(is.null(xlabels)){
        xlabels <- rownames(counts)
    }
    if(is.null(xcols)){
        xcols <- rainbow(nrow(counts))
    }

    par(mar=0.1+c(12,4,1,1))

    counts <- t(counts)
    rowSums <- apply(counts,1,sum)
    ncts <- 100*counts/rowSums
    medians <- apply(ncts,2,function(x) median(x[x>0]))
    ncts <- ncts[,order(medians,decreasing=T)]
    xlabels <- xlabels[order(medians,decreasing=T)]
    xcols <- xcols[order(medians,decreasing=T)]

    if(!is.null(cutoff)){
        cutoff <- as.integer(cutoff)
        ncts <- ncts[,1:cutoff]
        xlabels <- xlabels[1:cutoff]
        xcols <- xcols[1:cutoff]
    }

    brks = 10^seq(-2.8,2,length=res)

    stats <- apply(ncts,2,function(x) boxplot(x[x>(10^-2.8)],plot=F))
    hists <- apply(ncts,2,function(x) hist(x[x>(10^-2.8)],breaks=brks,plot=F))
    zeros <- apply(ncts,2,function(x) sum(x==0))

    plot(1,type="n",xaxs="i",xlim=c(0.5,0.5+ncol(ncts)),ylim=c(1e-3,1e2),log="y",xlab="",ylab="Relative Abundance",xaxt="none",yaxt="none")
    abline(h=c(1e-2,1e-1,1e0,1e1,1e2),col="grey")
    axis(1,xaxs="i",at=1:ncol(ncts),labels=xlabels,las=2)
    axis(2,at=c(1e-3,1e-2,1e-1,1e0,1e1,1e2),labels=c("Undetected","0.01%","0.1%","1%","10%","100%"))

    for(i in 1:ncol(ncts)){
        s = stats[[i]]
        if(type%in%c("bar","barswarm")){
            if(type=="bar"){
                points(rep(i,length(s$out)),s$out,pch=20,col=xcols[i])
                rect(i-0.4,s$stats[2],i+0.4,s$stats[4],col=xcols[i])
            }else{
                points(rep(i,length(s$out)),s$out,pch=20,col=paste(xcols[i],"40",sep=""))
                rect(i-0.4,s$stats[2],i+0.4,s$stats[4],col=paste(xcols[i],"40",sep=""))
            }
            segments(rep(i,2),s$stats[c(1,4)],rep(i,2),s$stats[c(2,5)])
        }else if(type%in%c("violin","violinswarm")){
            h = hists[[i]]
            if(sum(h$counts)>0){
                ycoords = approx(h$mids,n=5*res)$y
                lo <- loess(h$counts/max(h$counts*2)~h$mids,span=0.25)
                pred <- predict(lo,ycoords)
                pred[pred<=0] <- NA
                predlist <- split(pred,cumsum(is.na(pred)))
                ycoordslist <- split(ycoords,cumsum(is.na(pred)))
                ycoordslist <- lapply(1:length(ycoordslist),function(x) c(ycoordslist[[x]][!is.na(predlist[[x]])],rev(ycoordslist[[x]][!is.na(predlist[[x]])])))
                predlist <- lapply(predlist,function(x) c(i-x[!is.na(x)],i+rev(x[!is.na(x)])))
                for(p in 1:length(predlist)){
                    pred = predlist[[p]]
                    ycoords = ycoordslist[[p]]
                    if(type=="violin"){
                        polygon(pred,ycoords,col=xcols[i])
                    }else{
                        polygon(pred,ycoords,col=paste(xcols[i],"40",sep=""))
                    }
                }
            }
        }else if(type=="points"){
            points(rep(i,nrow(ncts)),ncts[,i],col=xcols[i],pch=20)
        }else if(type=="swarm"){
            lines(c(i,i),c(2e-3,1e2),col=1,lwd=1)
        }
        points(i,1e-3,cex=4*zeros[i]/nrow(ncts),col=xcols[i],pch=20)
        text(i,1.3e-3,zeros[i])
    }
    if(type%in%c("swarm","barswarm","violinswarm")){
        beeswarm(x=as.data.frame(ncts),add=TRUE,col=xcols,pch=20,corral="wrap")
    }
    segments(1:ncol(ncts)-0.4,unlist(lapply(stats,function(x) x$stats[3])),1:ncol(ncts)+0.4,unlist(lapply(stats,function(x) x$stats[3])),lwd=2)
    return(list(propCounts=ncts,stats=stats))
}