Pretty printing progress bars from .Call in R

I wanted to do something that takes a while in C (called from R) and wanted to print a progress bar.

In the R function, I did

pBar <- txtProgressBar( min = 0, max = n, style = 3 )
ret <- .Call("myfunction", arg1, arg2, as.integer(type), pBar, PACKAGE="annotSnpStats")

In the C function, I had

SEXP myfunction(SEXP arg1, SEXP arg2, SEXP Rtype, SEXP pBar) {

  int nprotect=0;
  SEXP utilsPackage, percentComplete;
  PROTECT(utilsPackage = eval(lang2(install("getNamespace"), ScalarString(mkChar("utils"))), R_GlobalEnv));
  PROTECT(percentComplete = allocVector(INTSXP, 1));
  nprotect+=2;
  int *rPercentComplete = INTEGER(percentComplete);
 
  ...

  for(i=0; i<nx; i++) { // index rows of x
  *rPercentComplete = i; //this value increments
  eval(lang4(install("setTxtProgressBar"), pBar, percentComplete, R_NilValue), utilsPackage);
  }

  ... 

  UNPROTECT(nprotect);
  return(myreturn);

}

And … it … just … worked! Happy day.

To see the full R function, go to dups in snp-match.R and the C is in countdiffs in comparisons.c.

Coloc v2.2

Coloc v2.2 is up on CRAN. The last two changes are shown below. The important thing is that the arguments to coloc.abf() have changed. Please do revise your code!

I do try and avoid completely changing arguments to released functions, but in this case, the function was introduced relatively recently, and the change makes sense because it allows us to analyse either datasets for which coefficients and standard errors are available, or for which only p values and minor allele frequencies are available, in a single function.

2013-19-06  Chris Wallace  <chris.wallace@cimr.cam.ac.uk>

        * v2.2
        Merged coloc.abf and coloc.abf.imputed(), so that datasets for
        wheich beta, var(beta) are available can be matched to datasets
        with only p values and maf.2

        This means the arguments to coloc.abf() have been changed!  Please
        check ?coloc.abf for the new function.
2013-03-06  Chris Wallace  <chris.wallace@cimr.cam.ac.uk>

        * v2.1
        Bug fix for coloc.abf() function, which used p12 instead of
        log(p12) to calculate L4.

        New function coloc.abf.imputed() to make better use of fuller
        information on imputed data.

Major update to coloc package

I am a co-author on another paper about colocalisation posted on arXiv. It’s a novel approach, using Bayesian inference based on Approximate Bayes Factors derived from p values, making colocalisation testing much more practical when data is not often as open access as claimed. My co-author, Vincent Plagnol, has written a nice post about it on Haldane’s Sieve.

The software to conduct these tests has been included in my coloc package, and I think such a change deserves a bump in version, so coloc v2.0 is now available on CRAN, together with a new vignette, explaining the different methods of analysis available, but see the papers 1 2 for the nitty gritty.

Writing R vignettes in emacs org mode using ox-ravel

I want to be able to write R vignettes in org-mode. Crazy, perhaps, but it has become my default way to write everything, notes on current projects, even beamer slides and papers. Now I need to write a vignette for the coloc R package, and switching to Sweave reminds me what a pain that is. So I had a play with knitr before, and liked it, but again it means learning a new markup syntax.

I have found ox-ravel, which means I can now write .org files, and export part or all to Sweave, knitr, brew, etc. I cannot stress just how fantastic this is, in allowing me to write everything in one syntax!

But setting it up required updating org-mode, because it requires the new ox export engine found in org-8. In the process, I didn’t want to break org2blog which I use to write these blog posts. Knowing this could get complicated, I documented each step here.

Continue reading

snpStatsWriter v1.5-1 and knitr vignettes

A new version of snpStatsWriter is up on CRAN. It includes a new function, write.snphap(), for writing snphap files, and a brief vignette has been added, my first vignette written in markdown with knitr.

I’m not absolutely sure I want to learn another markdown syntax, I spend so much time in Org mode that translating from one to another is a pain, but I wanted to give knitr a go and markdown is friendlier than Sweave-like syntax. Perhaps next time I will try ravel to combine org-mode and knitr.

Anyway, it worked well, and you can see the markdown vignette source and resulting html.

The joys of Rscript and parsing commandArgs

I’ve had a hacky, but workable, setup using a mixture of m4 and R CMD BATCH to run jobs on our local grid varying the parameters each time. Today I switched to Rscript, and it’s been a breeze.

At the top of my .R file I put

#!/usr/bin/Rscript
args <- getArgs()

and I put the getArgs() function below in my ~/.Rprofile.

Then I can run the script using

./myfile.R --args myarg=something

and in the .R file refer to args$myargs. –args is used to separate arguments you want to pass to R from those you want to be available in your script.

Continue reading

git cheatsheet

There are lots of good docs on git and github around. Perhaps, sometimes, too many! I am not attempting to add to them, just collating the specific sets of commands that I find useful in one place. I don’t promise they are useful for you, but the references to more complete docs may be.

Continue reading

Posted in git

ggplot2 cheatsheet

These are some of the most common things I want to do with ggplot2, but always forget the syntax (not helped when things are still evolving with each release!)

Titles, axis labels

  • main title
+ ggtitle("...")
  • set axis title (but see scale_x_continuous below)
+ xlab("...")
  • remove axis title, labels and tick marks (respectively)
+ theme(axis.title.x=element_blank(), 
        axis.text.x=element_blank(), 
        axis.ticks=element_blank())
  • make axis labels at a 45 degree angle
+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
  • customize axis breaks and labels. title is optional, labels are optional and default to as.character(breaks)
scale_x_continuous("x axis title (optional)",
                   breaks=c(2000,5000),
                   labels=c("foo","bar"))
  • stretch axis to full width/length of plot
scale_x_continuous(expand=c(0,0))

Legend

  • no legend
+ theme(legend.position="none")

ggheatmap, version 2

It looks like my ggplot2 heatmap function gets most traffic on this blog. That’s a bit unfortunate, because it’s the first function I wrote in earnest using ggplot2 and ggplot2 itself has undergone some updates since then, meaning my code is clunky, outdated and, er, broken.

So, with a bit more knowledge of ggplot2 and grid gained over the last few months, I have updated it today, and it is working. I hope it’s useful!

Get the code from github, and run it with

## simulate data
library(mvtnorm) 
sigma=matrix(0,10,10)
sigma[1:4,1:4] <- 0.6
sigma[6:10,6:10] <- 0.8
diag(sigma) <- 1
X <- rmvnorm(n=100,mean=rep(0,10),sigma=sigma)

## make plot
p <- ggheatmap(X)

## display plot
ggheatmap.show(p)

The result should look something like

https://cwcode.wordpress.com/wp-content/uploads/2013/01/wpid-ggheatmap.jpg