pages tagged r

Shrinking lists of gene names in R

Thu, 23 Mar 2017 13:45:44 -0700

I've been trying to finish a paper where I compare gene expression in 14 different placentas. One of the supplemental figures compares median expression in gene trees across all 14 species, but because tree ids like ENSGT00840000129673 aren't very expressive, and names like "COL11A2, COL5A3, COL4A1, COL1A1, COL2A1, COL1A2, COL4A6, COL4A5, COL7A1, COL27A1, COL11A1, COL4A4, COL4A3, COL3A1, COL4A2, COL5A2, COL5A1, COL24A1" take up too much space, I wanted a function which could collapse the gene names into something which uses bash glob syntax to more succinctly list the gene names, like: COL{11A{1,2},1A{1,2},24A1,27A1,2A1,3A1,4A{1,2,3,4,5,6},5A{1,2,3},7A1}.

Thus, a crazy function which uses lcprefix from Biostrings and some looping was born:

collapse.gene.names <- function(x,min.collapse=2) {
    ## longest common substring
    if (is.null(x) || length(x)==0) {
        return(as.character(NA))
    }
    x <- sort(unique(x))
    str_collapse <- function(y,len) {
        if (len == 1 || length(y) < 2) {
            return(y)
        }
        y.tree <-
            gsub(paste0("^(.{",len,"}).*$"),"\\1",y[1])
        y.rem <-
            gsub(paste0("^.{",len,"}"),"",y)
        y.rem.prefix <-
            sum(combn(y.rem,2,function(x){Biostrings::lcprefix(x[1],x[2])}) >= 2)
        if (length(y.rem) > 3 &&
            y.rem.prefix >= 2
            ) {
            y.rem <- 
                collapse.gene.names(y.rem,min.collapse=1)
        }
        paste0(y.tree,
               "{",paste(collapse=",",
                         y.rem),"}")
    }
    i <- 1
    ret <- NULL
    while (i <= length(x)) {
        col.pmin <-
            pmin(sapply(x,Biostrings::lcprefix,x[i]))
        collapseable <-
            which(col.pmin > min.collapse)
        if (length(collapseable) == 0) {
            ret <- c(ret,x[i])
            i <- i+1
        } else {
            ret <- c(ret,
                     str_collapse(x[collapseable],
                                  min(col.pmin[collapseable]))
                     )
            i <- max(collapseable)+1
        }
    }
    return(paste0(collapse=",",ret))
}

Adding a Table of Contents to PDFs from R

Tue, 13 Jan 2015 17:29:12 -0800

I routinely generate very large PDFs from R which have hundreds (or thousands) of pages, and navigating these pages can be very difficult. Unfortunately, neither R's pdf() nor its cairopdf() drivers support creating Table of Contents (or Index) while plots are being written out. In the case of cairo, the underlying library doesn't support it either, so this isn't something that can easily be added to R directly. I had been thinking about sitting down for months and writing the support into cairo and R's cairo package... but real life kept getting in the way.

Fast forward to a week ago, when I realized that pdftk does support dumping the table of contents and updating the table of contents using dump_data_utf8 and update_info_utf8! Armed with that knowledge, and a bit of hackery, we can save an index, and then update the pdf once it's been closed.

The R code then looks like the following:

 ..device.set.up <- FALSE
 ..current.page <<- 0

 save.bookmark <- function(text,bookmarks=list(),level=1,page=NULL) {
     if (!..device.set.up) {
         Cairo.onSave(device = dev.cur(),
                      onSave=function(device,page){
                          ..current.page <<- page
                      })
         ..device.set.up <<- TRUE
     }
     if (missing(page)|| is.null(page)) {
         page <- ..current.page+1
     }
     bookmarks[[length(bookmarks)+1]] <-
         list(text=text,
              level=level,
              page=page)
     return(bookmarks)
 }

 write.bookmarks <- function(pdf.file,bookmarks=list()) {
     pdf.bookmarks <- ""
     for (bookmark in 1:length(bookmarks)) {
         pdf.bookmarks <-
             paste0(pdf.bookmarks,
                    "BookmarkBegin\n",
                    "BookmarkTitle: ",bookmarks[[bookmark]]$text,"\n",
                    "BookmarkLevel: ",bookmarks[[bookmark]]$level,"\n",
                    "BookmarkPageNumber: ",bookmarks[[bookmark]]$page,"\n")
     }
     temp.pdf <- tempfile(pattern=basename(pdf.file))
     temp.pdf.info <- tempfile(pattern=paste0(basename(pdf.file),"info_utf8"))
     cat(file=temp.pdf.info,pdf.bookmarks)
     system2("pdftk",c(pdf.file,'update_info_utf8',temp.pdf.info,'output',temp.pdf))
     if (file.exists(temp.pdf)) {
         file.rename(temp.pdf,pdf.file)
     } else {
         warning("unable to properly create bookmarks")
     }
 }

and can be used like so:

 cairopdf(file="testing.pdf")
 bookmarks <- list()
 bookmarks <- save.bookmark("First plot",bookmarks)
 plot(1:5,6:10)
 bookmarks <- save.bookmark("Second plot",bookmarks)
 plot(6:10,1:5)
 dev.off()
 write.bookmarks("testing.pdf",bookmarks)

et voila. Bookmarks and a table of contents for PDFs.

This basic methodology can be extended to any language which writes PDFs and does not have a built-in method for generating a Table of Contents. Currently, the usage of Cairo.onSave is a horrible hack, and may conflict with anything else which uses the onSave hook, but hopefully R will report the current page number from Cairo in the future.

Plotting Embedded Bitmap in Vector Plot in R

Mon, 24 Feb 2014 21:02:00 -0800

Recently, one of my collaborators complained that one of my plots took forever to render on his machine. The particular plot in question had a few thousand points, many of which were overlapping. Ideally, R would be able to simplify the vector image which was drawn to avoid drawing points which were occluded by other points, but this is difficult to do properly, and R currently doesn't do it.

However, R is able to plot to a bitmap, and bitmap images have the nice property of automatically handling this for you. Furthermore, raster images have recently been made far less clunky in R, so it's pretty easy to shove an arbitrary bitmap image anywhere. With dev.capture in Cairo coupled with grid.raster in grid, we have everything we need to solve this problem:

require(grid)
require(Cairo)

start.rasterplot <- function(width=NULL,height=NULL) {
    x.y.ratio <- convertX(unit(1,"npc"),"mm",valueOnly=TRUE)/
        convertY(unit(1,"npc"),"mm",valueOnly=TRUE)
    width.points <- as.numeric(convertX(unit(1,"npc"),"points"))
    dpi <- as.numeric(convertX(unit(1,"inch"),"point"))
    if (is.null(width) && is.null(height)) {
        width <- 1024
    }
    if (is.null(width)) {
        width <- height*x.y.ratio
    }
    if (is.null(height)) {
        height <- width/x.y.ratio
    }
    Cairo(width=width,height=height,dpi=1024/width.points*dpi,file="/dev/null")
}

stop.rasterplot <- function(plot=TRUE) {
    raster.image <- dev.capture(native=TRUE)
    dev.off()
    if (plot) {
        grid.raster(raster.image,width=unit(1,"npc"),height=unit(1,"npc"))
        return(invisible())
    } else {
        return(raster.image)
    }
}

Now we can do the following:

pdf(file="raster.pdf")
start.rasterplot()
print(xyplot(y~x,
      data=data.frame(y=rnorm(1E8),x=rnorm(1E8))))
stop.rasterplot()
dev.off()

and our PDF will contain a raster image, and will load in seconds instead of taking forever to plot the file.

Bug Reporting Rate in Debian

Tue, 09 Oct 2012 18:04:53 -0700

Christian's most recent blog post got me wondering if the decline in the bug reporting rate in Debian was something new, or something which often happened during releases. So, lets try to figure that out. In the BTS, when a bug report is filed, the report is written to a file called bugnum.report, and then not touched from then on. Let's look at the modification date on that file to see when each bug was filed; and since we're going to plot this, lets only look at bugs ending in 00:

stat -c '%n %Y' /srv/bugs.debian.org/spool/{archive,db-h}/00/*.report > ~/reporting_rate.txt

Now, lets get the data into R and plot it. [For clarity, I'm not showing the R code, but it's available in the source code for this post.]

From the plot (Bugs reported per second over time with a red loess fit line), it looks like we do see a decline during certain periods. However, there's an even more alarming trend of a decrease in bug reporting in Debian which has been happening since 2006. (Note that I've truncated the y scale significantly; there are periods in Debian where the bug rate is astronomically high, usually corresponding to mass bug filings; I've also limited the plot to data from 2003 on, as I have to clean up that data significantly before I can plot it like this.)

Not sure exactly what that means, but it is troubling.

Plotting Ethnic Regions

Mon, 08 Oct 2012 11:41:19 -0700

I was asked over the weekend how to plot SNPs which are associated with specific ethnicities; the following code is a really quick stab at the problem using the grid graphics engine in R.

> require(grid)
> snp.position <- 1:10
> snp.integers <- sample.int(5,size=length(snp.position),replace=TRUE)
> ### the position is really the midpoint of the range
> snp.midpoint <- snp.position
> ### start of the SNP; we're assuming it should start at 0.
> snp.start <- c(0,snp.position[-1]-diff(snp.position)/2)
> ### stop of the snp; we're assuming it should stop at the last snp position
> snp.stop <- c(snp.position[-length(snp.position)]+diff(snp.position)/2,
>   snp.position[length(snp.position)])
> snp.width <- snp.stop - snp.start
> 
> ### these are the colors
> possible.colors <- c("red","blue","yellow","green","purple")
> snp.colors <- possible.colors[snp.integers]
> 
> ### this sets up the viewport that we'll plot into
> pushViewport(viewport(height=unit(1,"npc")-unit(7,"lines"),
> width=unit(1,"npc")-unit(7,"lines")
>   ))
> pushViewport(dataViewport(xscale=range(c(0,snp.stop)),yscale=c(0,1)))
> 
> ### this draws a rectangle around the graph
> grid.rect()
> ### this sets up the x axis
> grid.xaxis()
> ### this labels the X axis
> grid.text("Position on Chromosome",y=unit(-2.5,"lines"))
> 
> ### this draws all of the boxes corresponding to each SNP in the
> ### appropriate color
> grid.rect(x=unit(snp.start,"native"),
>         width=unit(snp.width,"native"),
>         y=unit(0.5,"native"),
>         height=unit(0.25,"native"),just=c("left","center"),
>         gp=gpar(col=snp.colors,fill=snp.colors))
> 
> ### this pops the viewport
> popViewport(2)

Unnecessary Use of data.frame is Slow

Mon, 16 Apr 2012 11:32:27 -0700

I've been working for a while on a reasonably large Genome-Wide Association Study dataset which has lead me through various interesting parts of handing large datasets in R. This dataset is approximately 320,000 rows by 5000 columns. After getting Rmpi working, and handling the dataset by row so I don't run out of memory, I've managed to get pretty decent performance. However, one small section of the code seemed to be taking forever to run.

It turns out that assigning data to a data.frame by row is incredibly slow in R. Thus, a section of my code which should have taken microseconds was taking tenths of seconds, and threatening to run all week. Using a matrix instead (which is basically what I want anyway) and converting to a data.frame at the very end makes the code multiple orders of magnitude faster.

Moral of the story? Don't use data.frame unnecessarily.

Introducing Sweavealike

Thu, 22 Mar 2012 14:06:37 -0700

I use R a lot. It's one of the primary tools I use in my day job as a scientist analyzing large datasets. If you use LaTeX with R (as I often do), you probably use Sweave to interleave R output and figures with your text describing those figures using the noweb method of literate programming.

Sweavealike is a plugin for IkiWiki that tries to do some of the useful things for IkiWiki that sweave does for R and LaTeX.

You use it like the following:

[[!sweavealike  echo=1 code="""
a <- 1
a <- a + 10
print(a)
"""]]

which produces this result when run:

> a <- 1
> a <- a + 10
> print(a)
[1] 11

You can also generate figures with it:

[[!sweavealike  fig=1 echo=1 results="hide" code="""
plot(1:10,(1:10)^2,xlab="x",ylab=expression(x^2),main="Example Figure")
"""]]

> plot(1:10,(1:10)^2,xlab="x",ylab=expression(x^2),main="Example Figure")

The plugin itself uses the neat Statistics::R perl module to handle all of the heavy lifting. I personally plan on using this plugin to help write some more entries in my learning R series of posts that I'm beginning to work on. Hopefully I'll find and fix most of the bugs as I embark on that process so anyone else who uses the plugin won't, but feel free to e-mail me if something isn't working as it should.

Finally, you shouldn't run this plugin on a publicly editable IkiWiki instance, because that would be a trivial local user exploit as R can run arbitrary code, read and write to arbitrary files, exhaust all memory, etc.

Debian R Archive

Fri, 02 Mar 2012 10:26:53 -0800

For those of you were were in the various Debian infrastructure channels, you might have noticed that I was playing around with wanna-build, dak, sbuild, and buildd and friends. [Thanks to everyone who answered questions, btw.] Over the past week, I've been building most of CRAN, Bioconductor, and omegahat for unstable, amd64. I plan to build the same set of packages for i386, and will start a build shortly for stable as well. This effort builds on top of Charles Blundell and Dirk Eddelbuettel's cran2deb, which does most of the heavy lifting.

If you're like me, and use lots of different R packages, or already use some of the R packages available on the previous build, you can simply point your sources.list to the [http://debian-r.debian.net] archive, load the appropriate GPG key, and away you go. I have a bit more information available here and I will try to keep that page updated as I build other architectures and build out for stable.