R/summary_sheets_pats.r

In hloefflerwirth/oposSOM: Comprehensive analysis of transcriptome data

pipeline.summarySheetsPATs <- function(env)
{
  util.info("Processing PAT-centered Analyses")
  
  dir.create("Summary Sheets - PATs", showWarnings=FALSE)
  
  pat.metadata <- do.call(cbind, by(t(env$metadata), env$pat.labels, colMeans))


  filename <- file.path("Summary Sheets - PATs","PAT report.pdf")

  util.info("Writing:", filename)
  pdf(filename, 29.7/2.54, 21/2.54, useDingbats=FALSE)

  
  ### PAT clustering

  if( ncol(pat.metadata) > 2 )
  {
    hc <- hclust(dist(t(pat.metadata)))
  
    par(mar=c(8,6,4,5))
    plot(hc, main="PAT clustering", xlab="", sub="", hang=0.05 )
    
    cex.sample.portraits <- 0.025
    
    par(new=TRUE,mar=c(1,par("mar")[-1]))
    frame()
    
    for (i in seq_along(hc$order))
    {
      pat <- colnames(pat.metadata)[ hc$order[i] ]
                      
      m <- matrix(pat.metadata[,pat],
                  env$preferences$dim.1stLvlSom, env$preferences$dim.1stLvlSom)
  
      if (max(m) - min(m) != 0)  m <- 1 + (m - min(m)) / (max(m) - min(m)) * 999
  
      m <- cbind(apply(m, 1, function(x){x}))[nrow(m):1,]
      pix <- pixmapIndexed(m , col = env$color.palette.portraits(1000), cellres=10)
  
      addlogo(pix, (i-1)/(ncol(pat.metadata)-1)+cex.sample.portraits*c(-1,1), 0.05+0.08*(i%%2)+cex.sample.portraits*c(-1.2,1.2))
    }
  }
  
  ### PAT portraits
  
  layout(matrix(1:16,ncol=4,byrow=TRUE))
  
  for( pat in names(sort(table(env$pat.labels),decreasing=TRUE)) )
  {
       
    par(mar=c(0,0,0,0))
    frame()
      text(0.24, 0.94, pat, cex=3, adj=0)
      
    par(new=TRUE,mar=c(2,6,3,2))
    image(matrix(pat.metadata[,pat], env$preferences$dim.1stLvlSom, env$preferences$dim.1stLvlSom),
        axes=FALSE, col=env$color.palette.portraits(1000) )
        box()
        
    par(mar=c(0.5,0,4.5,2))  
    pat.samples <- names(which(env$pat.labels==pat))
    pie.table <- table(env$group.labels[pat.samples])[unique(env$group.labels)]
    pie.table <- pie.table[which(!is.na(pie.table))]
    pie(pie.table,col=env$groupwise.group.colors[names(pie.table)])

    mtext(paste("# samples:",sum(env$pat.labels==pat)),3) 
  }

  
  if( length(unique(env$pat.labels)) > 1 && length(unique(env$group.labels)) > 1 )
  {  
    ### Group to PAT relation barplots
    ##
    
    pat.group.assoc <- tapply(env$group.labels, env$pat.labels, function(x)
    {
      ret <- table(x)[unique(env$group.labels)]
      names(ret) <- unique(env$group.labels)
      return(ret)
    })
    pat.group.assoc <- do.call(cbind,pat.group.assoc)
    pat.group.assoc[which(is.na(pat.group.assoc))] <- 0
    
    layout(matrix(c(1, 2), 1, 2), c(4, 1), 1)
    par(mar=c(5, 5, 4, 2))
    barplot( pat.group.assoc[,ncol(pat.group.assoc):1],
            main="Association of PATs and groups", cex.main=1.5, cex.axis=1.4, cex.names=1,
            col=env$groupwise.group.colors, horiz=TRUE, las=1)
    
    par(mar=c(5, 1, 4, 2))
    plot(0, main="", cex.main=1, type="n", axes=FALSE, xlab="",
         ylab="", xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i")
  
    legend("topleft",legend=unique(env$group.labels),
         col=env$groupwise.group.colors, pch=15, pt.cex=1.4, title="groups")
  
    ##
    
    pat.group.assoc <- tapply(env$pat.labels, env$group.labels, function(x)
    {
      ret <- table(x)[sort(unique(env$pat.labels))]
      names(ret) <- sort(unique(env$pat.labels))
      return(ret)
    })[unique(env$group.labels)] 
    pat.group.assoc <- do.call(cbind,pat.group.assoc)
    pat.group.assoc[which(is.na(pat.group.assoc))] <- 0
  
    par(mfrow=c(1,1), mar=c(5, 8, 4, 2)) 
    y.coords <- barplot( pat.group.assoc[,ncol(pat.group.assoc):1],
             main="Association of PATs and groups", cex.main=1.5, cex.axis=1.4, cex.names=1,
             horiz=TRUE, las=1) 
    for( i in seq_along(unique(env$group.labels)) )
    {
      print.labels <- which( pat.group.assoc[,i] > max(colSums(pat.group.assoc,na.rm=TRUE))*0.05 )
      for( ii in seq_along(print.labels) )
      {
        text( x=sum( pat.group.assoc[1:(print.labels[ii]),i] ) - pat.group.assoc[print.labels[ii],i]/2,
              y=rev(y.coords)[i],
              names(print.labels)[ii] )
      }
    }
    
  }
  dev.off()
}