R/plot.clustal.R

   1 #' plot.clustal
   2 #'
   3 library(grid)
   4
   5 plot.clustal <- function(clustal,tree=NULL,refseq=1,refseq.start=1) {
   6     string.width <-
   7         function(x) {
   8             sapply(x,function(y){convertX(unit(1,"strwidth",y),"mm",
   9                                           valueOnly=TRUE)})
  10         }
  11
  12     ## http://www.jalview.org/help/html/colourSchemes/clustal.html
  13     ## https://www.cgl.ucsf.edu/chimera/docs/ContributedSoftware/multalignviewer/colprot.par
  14     calculate.color.column <-
  15         function(x) {
  16             ratio.true <- function(x){sum(x)/length(x)}
  17             colors <- rep(NA,length.out=length(x))
  18             clustal.colors <-
  19                 list("red"    = rgb(0.9,  0.2,  0.1, 1),
  20                      "blue"   = rgb(0.1,  0.5,  0.9, 1),
  21                      "green"  = rgb(0.1,  0.8,  0.1, 1),
  22                      "cyan"   = rgb(0.1,  0.7,  0.7, 1),
  23                      "pink"   = rgb(0.9,  0.5,  0.5, 1),
  24                      "magenta"= rgb(0.8,  0.3,  0.8, 1),
  25                      "yellow" = rgb(0.69, 0.69, 0.0, 1),
  26                      "orange" = rgb(0.9,  0.6,  0.3, 1))
  27             colors[x=="G"] <- "orange"
  28             colors[x=="P"] <- "yellow"
  29             ## identify consensus
  30             consensus <- " "
  31             tests <-
  32                 list("%"=list(aa=c("W","L","V","I","M","A","F","C","Y","H","P"),
  33                          min=0.6),
  34                      "#"=list(aa=c("W","L","V","I","M","A","F","C","Y","H","P"),
  35                          min=0.8),
  36                      "-"=list(aa=c("E","D"),
  37                          min=0.5),
  38                      "+"=list(aa=c("K","R"),
  39                          min=0.6),
  40                      "g"=list(aa="G",
  41                          min=0.5),
  42                      "n"=list(aa="N",
  43                          min=0.5),
  44                      "q"=list(aa=c("Q","E"),
  45                          min=0.5),
  46                      "p"=list(aa="P",
  47                          min=0.5),
  48                      "t"=list(aa=c("T","S"),
  49                          min=0.5))
  50             for (symbol in c("A","C","D","E","F","G","H","I","K","L","M","N","P",
  51                          "Q","R","S","T","V","W","Y")) {
  52                 tests[[symbol]] <- list(aa=symbol,min=0.85)
  53             }
  54             for (test in names(tests)) {
  55                 if (ratio.true(x %in% tests[[test]][["aa"]]) >=
  56                     tests[[test]][["min"]]) {
  57                     consensus <- test
  58                 }
  59             }
  60             ## handle colors
  61             color.tests <-
  62                 list("G"=list(color="orange",
  63                          consensus=NULL),
  64                      "P"=list(color="yellow",
  65                          consensus=NULL),
  66                      "T"=list(color="green",
  67                          consensus=c("t","S","T","%","#")),
  68                      "S"=list(color="green",
  69                          consensus=c("t","S","T","#")),
  70                      "N"=list(color="green",
  71                          consensus=c("n","N","D")),
  72                      "Q"=list(color="green",
  73                          consensus=c("q","Q","E","+","K","R")),
  74                      "W"=list(color="blue",
  75                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  76                      "L"=list(color="blue",
  77                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  78                      "V"=list(color="blue",
  79                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  80                      "I"=list(color="blue",
  81                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  82                      "M"=list(color="blue",
  83                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  84                      "A"=list(color="blue",
  85                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p",
  86                              "T","S","s","G")),
  87                      "F"=list(color="blue",
  88                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  89                      "C"=list(color="blue",
  90                          consensus=c("%","#","A","F","H","I","L","M","V","W","Y","S","P","p")),
  91                      "C1"=list(color="pink",
  92                          consensus="C",
  93                          aa="C"),
  94                      "H"=list(color="cyan",
  95                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  96                      "Y"=list(color="cyan",
  97                          consensus=c("%","#","A","C","F","H","I","L","M","V","W","Y","P","p")),
  98                      "E"=list(color="magenta",
  99                          consensus=c("-","D","E","q","Q")),
 100                      "D"=list(color="magenta",
 101                          consensus=c("-","D","E","q","Q")),
 102                      "K"=list(color="magenta",
 103                          consensus=c("+","K","R","Q")),
 104                      "R"=list(color="magenta",
 105                          consensus=c("+","K","R","Q"))
 106                      )
 107             for (test in names(color.tests)) {
 108                 aa <- test
 109                 if (!is.null(color.tests[[test]]$aa))
 110                     aa <- color.tests[[test]]$aa
 111                 colors[x==aa & all(consensus %in% color.tests[[test]][["consensus"]])] <-
 112                     clustal.colors[[color.tests[[test]][["color"]]]]
 113             }
 114             return(colors)
 115         }
 116
 117     name.spacing <- string.width("  ")
 118     symbol.spacing <- string.width(" ")/4
 119     clustal.names <- rownames(clustal)
 120     clustal.names.width <-
 121         string.width(clustal.names)
 122     ## find max character width in each string column
 123     symbol.widths <- rep(max(apply(clustal,2,function(x){max(string.width(x))})),
 124                          ncol(clustal))
 125     full.width <- convertX(unit(sum(symbol.widths)+
 126                                     name.spacing+max(clustal.names.width),"mm"),
 127                            "npc",valueOnly=TRUE)
 128     symbol.pos.x <- cumsum(symbol.spacing+symbol.widths)
 129     symbol.pos.x <- c(0,symbol.pos.x[-length(symbol.pos.x)])+symbol.widths/2+symbol.spacing
 130     ## going to be wider than the screen
 131     if (full.width > 1) {
 132         clustal.names.x <-
 133             unit(0,"npc")
 134         symbol.offset.x <-
 135             unit(max(c(convertX(unit(1,"npc"),"mm",valueOnly=TRUE)-
 136                            sum(symbol.widths)+name.spacing,0)),"mm")
 137     } else {
 138         clustal.names.x <-
 139             unit((1-full.width)/2,"npc")+unit(max(clustal.names.width)-clustal.names.width,"mm")
 140         symbol.offset.x <-
 141             unit((1-full.width)/2,"npc")+unit(max(clustal.names.width)+name.spacing,"mm")
 142     }
 143
 144     ## plot the symbol names
 145     line.pos.y <-
 146         unit(seq((length(clustal.names)-1)/2,
 147                  -(length(clustal.names)-1)/2,
 148                  length.out=length(clustal.names)),"lines")+
 149                      unit(0.5,"npc")
 150
 151
 152     ## figure out X width; if it's bigger than the plotting area, jam
 153     ## the names in anyway
 154     grid.text(clustal.names,
 155               y=line.pos.y,
 156               x=clustal.names.x,hjust=0)
 157     grid.rect(x=unit(rep(symbol.pos.x,each=nrow(clustal)),
 158                   "mm")+symbol.offset.x,
 159               y=rep(line.pos.y,times=ncol(clustal)),
 160               width=unit(symbol.widths+symbol.spacing,"mm"),
 161               height=unit(1,"line"),
 162               gp=gpar(col=NA,
 163                       fill=as.vector(apply(clustal,2,calculate.color.column))))
 164     ## plot the protein characters
 165     grid.text(as.vector(clustal),
 166               y=rep(line.pos.y,times=ncol(clustal)),
 167               x=unit(rep(symbol.pos.x,each=nrow(clustal)),
 168                   "mm")+symbol.offset.x,hjust=0.5)
 169     ## plot the position strings
 170     position.strings <- seq(refseq.start,by=1,length.out=ncol(clustal))
 171     if(length(position.strings) < 10) {
 172         wanted.position.strings <- position.strings %% 2 == 0
 173     } else {
 174         wanted.position.strings <- position.strings %% 5 == 0
 175     }
 176     position.strings <- as.character(position.strings)
 177     position.strings[!wanted.position.strings] <- ""
 178     grid.text(position.strings,
 179               y=unit((length(clustal.names)-1)/2+1,"lines")+unit(0.5,"npc"),
 180               x=unit(symbol.pos.x,"mm")+symbol.offset.x,
 181               rot=90
 182               )
 183 }