Meta analysis in R

or 

the beneficial effect of teaching on research

I have been fascinated by meta analysis for a long time. It is so obviously the right way to approach the true effect of an intervention. Recently, an old binder presented itself in a pile of shite bunch papers I meant to read but found myself throwing in the bin tidying away. It contained the draft of a database of physiological data from the first years of my PhD. The idea was to compare all the baseline data from our kidney research group in Uppsala to look at the effect of the models as such and the interventions that were used repeatedly. With the insight of the intervening years it seems a lot less interesting now, but I still have the feeling that some areas of experimental research could benefit from meta analysis.

Which brings us to the the story I am about to tell. Four or five years ago, when I moved back to Uppsala, I got offered a lecture on physiological changes in the elderly. It was to be part of a final year course for the master in biomedicine programme. Just a single hour to show how all the physiology from the rest of the programme changed with age. To compare and contrast ageing as such with the accrued ailments of living for a long time, and distinguish these from the chronic and age related disease. It was not a huge success, but given the title I was not too disillusioned. The second year I was given two hours. Still a bit on the short side, but one hundred percent better than one.

It is not the most popular lecture, but I have had it for five years now and one of the things I teach is that some parts of ageing is caused by metabolism itself. The burning of oxygen singes the organism and with time it will break much like the paneling in an old sauna. As proof of this I used the idea of caloric restriction, which can prolong life in many strains of yeast, mice, and rats. Then, in 2012 an article was published on the effect of caloric restriction in the Rhesus monkey, a primate, and reasonably the closest relative to humans in which an experiment could be expected to be finished any time soon. It showed no effect. I happily included this in my lecture as a counter-point. Until in 2014, when updating the lecture for a new semester, I found that another experiment with caloric restriction in Rhesus monkeys had published their data and found a clear difference.

This made it hard to continue the lecture as I had done, I could just show both studies and say that we don't know. But the total number of animals included was quite large, and the effect measure very straight-forward. Death. So, I performed a meta analysis of mortality of these two studies, and a third smaller study published in 2003. This is the story of that analysis.

Quickly I installed the R package rmeta by Thomas Lumley and set to work. It is quite easy really, we start with setting up a table of results from the included studies. The table should include the total number of subjects in each group, and the number of deaths per group.

Hultström, M. Acta Physiol (Oxf). 2015 Feb 14. doi: 10.1111/apha.12468.

The we push this trough the rmeta function meta.MH(). To get a forest plot, we just run the plot() command, which has a default for handling the result of meta.MH() in the form of a forest plot. If you have a larger meta analysis there is also the funnelplot() that can be used to assess publication bias. Anyway, the result is quick and easily understood, which is really one of the major strengths of the forest plot.

Hultström, M. Acta Physiol (Oxf). 2015 Feb 14. doi: 10.1111/apha.12468.

There was no significant effect of caloric restriction on all cause mortality in Rhesus monkeys. Or, rather there was a small, clearly non-significant, effect. One of the reviewers asked what would be needed to show if this effect was true. That is, could I please perform a power analysis. So, I installed the pwr package and ran a 2p2n.test() using the most generous effect estimate, i.e. a hypothetical study that ran to completion where the whole control population had died giving an effect of 0.08. This resulted in a required population of 2806 subjects to reach 85% power. This is the power-level which is normally used as the basis for power calculation in clinical studies. However, the age-related mortality was a different story that you can find in the actual article.

The next thing that surprised me was how difficult it was to get this simple little analysis published. It appears that experimental journals don't publish meta analyses, and clinical journals that publish meta analyses, don't publish experimental results. Finally, I found a benevolent editor at Acta Physiologica who permitted it to be published as an editorial. So that is where it resides today, and finally I can give a fairly clear answer in my lecture on the effect of reducing metabolism by caloric restriction on ageing and on mortality. Only problem is, I now have to explain meta analysis and forest plots before I can show the actual data.

And, no I am not going to starve myself so that I can avoid some diseases we can treat in favour for a frailty for which the only known treatment is eating more.

Speech impediment - Locale in R under OSX

My R-installation developed a language problem ("L" is a statistical ploglamming language). The startup message, all errors and system messages started appearing in German. I do speak German but it was still quite annoying. This ended with a work-around that I published previously in a post about R-startup scripts. I have finally figured out what causes the error. My Language & Text settings for language were as follows:

1. British English 
2. Swedish
3. Norwegian
4. Danish
5. German 
6. French 

The problem lies in using the British English locale. I use it because I hate when all programs always correct my excellent spelling from Proper English™ to American. R does not have a locale for British English, nor for Swedish, Norwegian or Danish. So, to R, the locale priority list actually looks like this:

1.  
2. 
   
3. 
   
4. 
   
5. German 
6. French

So, understandably it goes for what it knows. The solution is to extend the list and put a standard English locale high up:

1. British English
2. English 
3. Swedish
4. Norwegian
5. Danish
6. German 
7. French 

Physiological illiteracy impedes progress in medicine!

A recent paper in PNAS has caused a bit of a ruckus discussion. The authors (Tim W. Fawcett and Andrew D. Higginson from U-Bristol) have compared citation rates with the number of equations in the text of papers, and found that more equations per page is associated with a lower number of citations in non-theoretical papers, but that there was no effect in theoretical papers. This carries over to all citations as well since the number of non-theoretical papers is much larger than the number of theoretical papers. They conclude that a higher density of equations leads to lower numbers of citations, and that this indicates a worse dissemination of the results to the wider scientific community. They go as far as using the title: "Heavy use of equations impedes communication among biologists."

This engendered a barrage of more or less agitated rebuttals with titles such as: "Do not throw equations out with the theory bathwater," "Mathematical illiteracy impedes progress in biology," "A suggestion on improving mathematically heavy papers," and "No evidence that equations cause impeded communication among biologists." Anyway, an alternative view is that theoretical papers (with lots of equations) are mostly cited in other theoretical papers, although this sounds obvious and would never have been published in PNAS.

I would wager that the same is true in all fields: Physiologists certainly cite other physiologists more often than they cite clinical scientists (and almost never cite mathematical biologists). It is further easy to extrapolate to the conclusion that detailed physiological data in a paper is associated with fewer citations in clinical medical journals. The obvious interpretation is that overzealous use of physiology in articles impairs dissemination of medical science, and that physiologists should use more common words and talk more about patient-survival and quality-of-life in the same way that mathematical biologists should relegate their equations to the appendix. The idea that they might be comparing apples to oranges (clearly a biological concept) is not even discussed in the paper, which makes the discussion a bit lopsided. More as if they were writing an opinion-piece than an article on empirical science.

Analysis of genome-scale data

I thought I should do a post on strategies for handling genome-scale data from large experiments. Like many happy novices we rushed in and performed a couple of huge microarray experiments. Now we are trying to compare these against eachother and against strain-differences in in the form of SNPs. To say that this is not trivial is an understatement. In the end we have come up with a couple of strategies that might be useful to others as well. They might seem simple, but they will help you wade through your data, without getting stuck there for too long.

Don't look at untested data for anything but quality control. In most analysis programs it is very easy to just plot a gene by its intensity, either over all samples or by crude difference between group-averages. This is fun, and can give you some idea of what you are looking for. However, there are anything from 20k-30k genes in your dataset and you will not have time to look at all of them and make a reasoned decision. That means that if you pick genes this way you will indeed pick some hits, but not in a systematic way, and that leaves you open to errors. Either of picking a gene where the variation is too large, or that may only be the n+100th gene in your set and not the top one. The point is: You just don't know. The suggestion is to test your genes first, and then use the list of significantly differentially expressed genes to pick your top hits.

Consider doing a single comparison across several groups. If you have run a large experiment, with several groups (where several is anything from four to whatever). It is a chore to keep track of which genes were significant in which comparison, and even more to then pick a small number of genes to validate that are somehow representative across groups. So, try to pool your data by phenotype. If you have one group where something happens and lots of different controls. Just pool the controls and do a single comparison. That will give you just one list and make your selection process much easier.

When comparing experiments, use rankings. So, now you have a list, or several if pooling didn't make sense. Now you have to find your top candidates for validation. It is then helpful to use rankings. The easiest (and often quite accurate) way is to rank by average difference between groups (that would be your fold-change or log-ratio). If you have just one list, then you are done. If you have two, or several lists, then you can order your lists by rank-sum or rank-product. The rank-sum is a kind method that gives you some hits even if a gene is poorly-ranked (high number) in one of your groups. The rank product will pronounce those that are highly-ranked in all lists. Then you have the very particular rank-difference and rank-ratio. These may be very interesting if you have two experiments with wildly different phenotypes, but that you can not compare directly. If a gene is ranked in the top 10 in one experiment and 7815th in the other, it may be a very interesting gene and will then be ranked highly by the rank-difference or ratio.

R-startup scripts

If there are some things you use all the time in R, you might not want to have to load them specifically every time. That is what ~/.Rprofile is for. As a start, have a look at R-bloggers: Customizing your .Rprofile and Customizing R: startup script. Then the question is what to put in your startup-script.

cat("Good morning Dr. Chandra.\n\n", sep = "") is an excellent choice.

As an alternative:
library(lubridate)
if (hour(Sys.time()) < 12)
  cat("\nGood morning Dr. Chandra!\n\n", sep = "")
if ((hour(Sys.time()) >= 12) & (hour(Sys.time()) < 18))
  cat("\nGood afternoon Dr. Chandra!\n\n", sep = "")
if (hour(Sys.time()) >= 18)
  cat("\nGood evening Dr. Chandra!\n\n", sep = "")
is more probable to be correct in a linguistic sense.

Then there are other important things. Like setting the language to English. For some reason R thinks I want it to speak German, even though nothing else on my system is in German. I have given up trying to find out how to fix it properly, so I use: Sys.setenv(LANG = "EN")

As they point out elsewhere though, you shouldn't put too many things in that you actually use for data analysis. That would make your code difficult to follow. It might be anyway, as mine is, but at least all libraries and such are loaded in the R-directory of the specific project. However, I do have a loading option for a .Rbasicfunctions.R:
if(file.exists("~/.Rbasicfunctions.R")){
    source("~/.RbasicFunctions.R")
    cat(".RbasicFunctions.R was loaded, providing the following functions:\n\n",sep="")
    print.functions()}

The actual file is empty though.

How many hits do you want?

Data mining is a huge business today and also an area of intense investigation. After an intense week of trying to compare gene-expressionmicroarray data from different experiments during my recent spout in Glasgow, I thought I would share some thoughts on analysis methods. Basically the choice comes down to the number of hits you might want, which in-turn is decided by how you intend to analyse and validate your results.

With modern high-throughput methods one of the main problems is choosing candidates for validation and further study. As long as you are just doing PCRs for validation, you can run quite a number, but proper validation will often require knock-out, and/or knock-in experiments, possibly both in-vivo and in-vitro. Picking the wrong candidate can cost you years of work. So, you had better be correct when you make your decision. This is a reason that it is common to filter high-throughput data by the size of change. Whether that is advisable depends on what you want to do with the data. Modern high-level data analysis often works much better if you have more to work with. Not to mention that many of the new methods of analysis are designed to look for systems of changes, each of which may be negligible, that together may cause significant physiological effects.

The calculation of a p-value comes in two steps. First you calculate a value and then you correct that value for multiple comparisons.

• Anova: Mathematically elegant but computationally heavy and quite strict, producing not so many significant hits. It also depends on there being similar variations and similar number of samples per group.
• Rank product: Basically ranks each gene in each sample multiply these across samples and compares the results to similar products from randomly generated datasets. Of all the heavy methods this is one of the heavier. i.e. long computations times.
• Empirical Bayesian methods: I won't say I understand what exactly this does, but it is described as "reducing the standard error toward a common mean." This means that you are in the end comparing the estimates of the means of your samples, i.e. a variant of Student's T-test, but one that somehow takes the rest of your expression data into account. This is a quick method that is also well regarded.
• T-test: Your favourite statistical test, seldom used by itself in microarrays, se SAM for more details.
• SAM: Significance Analysis of Microarrays, using T-tests to test the contrast of interest and compare it to other permutations of the samples. If the contrast of interest shows a better significance (by a set level) than the random configurations then it is deemed significant. Depending on the number of permutations you use the analysis is more or less computationally heavy, but it does by necessity take quite a bit of processing.
• Signal to noise ratio: You basically divide the the mean (difference in mean) by the standard deviation. This in itself does not produce a P-value, but it does provide a way to rank genes that can be used to pick the most probably changed genes. If you compare the list you get to lists produced by random permutations of your data, then you can produce proper p-values, but that makes it more onerous.

Once you have your p-values, they have to be corrected for the number of comparisons. The thing is that if you use a probability of chance difference cut-off of 5%, as is common, fully 1/20th of your significant genes are probably false-positives.

When correcting your p-values the first important distinction (that many are not aware of) is post-hoc versus à-priori. Mostly these relate to how to test samples after an ANOVA. Post-hoc testing means that you do not know which comparisons you want to make, so you just make them all. 

For example: You want to assess the effect of a treatment on two populations, say men and women. 

This gives rise to four groups: untreated-men, treated-men, untreated-women and treated-women. Or, you may even have before and after measurements in each of these, giving you eight groups. Between four groups there are six (x*(x-1)/2) possible comparisons, and between eight groups there are 28. If you correct your p-value for all these you may introduce a lot of false-negative results, i.e. expression changes that really are there, but that does not make your corrected cut-off. However, many of these comparisons may be nonsensical.

The difference between untreated-men-before-treatment to treated-women-after-treatment, could be one such example. It is difficult to interpret without knowing what the baseline for the women was. What you can do in this kind of situation is decide à-priori (eng. beforehand) to only run some of the comparisons. 

For example: You compare men to women before treatment, before and after treatment only in the same group, the change between before and after only within each sex, and the difference between the changes in the treated and untreated men to the same difference in women*.

This is only eight comparisons out of the possible 28, and some of these would not be included in the 28 at all, i.e. the differences. By cutting the number of comparisons you also cut the correction you have to make to your p-values. In addition, you can argue that the absolute differences are independent from the change caused by treatment and correct these comparisons separately, giving you only a four-fold correction to your p-value instead of a 28-fold. 

Anyway, that was only corrections between groups for one parameter. In a typical microarray experiment you will have 20000 parameters across your groups, so that the total number of comparisons would be 28*20000, or a lot.

Just multiplying your p-values by the number of comparisons was the original and most conservative correction. It is called Bonferroni, because that's the guy who described it. This is generally too conservative and will produce false-negatives even in experiments with a small number of comparisons. However, it is useful for narrowing the number of hits if you think you have too many, or if you need to be really certain of the hits you pick.

The family-wise error rate correction is what ANOVA uses across all samples and it is what post-hoc tests like Tukey's honest significant difference uses for post-hoc testing between individual samples. It is quite strict in large, relatively low signal-to-noise data, like microarrays, but quite relaxed in smaller data-sets. This is difficult to use in à-priori-type designs since these make the estimation of the family-wise variation difficult.

False discovery rate, or FDR, is an iterative correction across a ranked list. Basically, the first hit is uncorrected, the second somewhat corrected and so on. What it is actually saying is that; out of all the genes, from the top of the list to the point you are looking at, a given percentage will be false discoveries. So if you get 3000 with an FDR <0.05 you have 5% or 150 false positives in that sample. This is one of the least strict corrections. It is often used in microarray analysis for that very reason.

The final question is if you correct your FDR values by the number of group-wise comparisons you are doing. Or, maybe you are comparing similar, but independent, experiments, which you could argue strengthens your p-values by as much as p1*p2.

As they say: If it was easy, anyone could do it.

* in this kind of design you should really use a multi-way ANOVA instead of a plain one, and that choice often makes it obvious that you should use à-priori contrasts instead of a post-hoc.

Using the OS X clipboard in R

One annoyance with using a command line based statistical system, i.e. R, is that if you use a spreadsheet program to handle your data (as most of you do) you have to incessantly export to csv and import to R. For big data sheets and complex calculations this generally a Good Idea(tm). It improves reproducibility and readability when going back to the project at some later time. However, the complexity of that is enough that you often, at least I, end up doing statistics in Excel. This is generally a Bad Idea, and often entails much more work in the long run. In order to make importing and exporting data easier R can handle the clipboard directly, which is what this post is about.

In OS X the clipboard is easily accessed through the command line programs pbcopy and pbpaste (pb means "pasteboard" but everyone knows it as the clipboard so that's what I'm going to call it). These are easily scriptable programs that can be used with any unix command in the shell or scripts (enter "man pbpaste" in the terminal to learn more). "pbpaste" pastes the current clipboard to std-out (that is where you call it from), while "pbcopy" copies std-in (what you put there) to the clipboard. An important quirk is that they only handle plain text, RTF and EPS, that is formats that are transferred as plain text. This is because unix is a text based system. What it means is that you can't copy stuff from formatted text, e.g. word, but with spreadsheets you are OK.

As command line programs they can be handled through the R command: pipe(). Pipe does not produce the current clipboard contents, just a description of the pipe as such. You have to use a read command to get the actual contents.

> readLines(pipe("pbpaste"))

which will produce one string per line (completely useless as data, but often useful to see how a given copy is formatted). To format it for data you use the same as you would any import, using one of the higher level read-commands: read.table(), read.csv(), etc.

> x <- read.table(pipe("pbpaste"))

This will put your copied excel table directly into the table "x" in R. Often you will be forced to tweak the function a little to get a usable table.

> x <- read.table(pipe("pbpaste"), header=TRUE, sep="\t")

That should be able to read your copied Excel data ("header" specifies that the first row is the header, and "sep" specifies the separator which defaults to tab or "\t" in excel).

To write to the clipboard you use "pbcopy". It is a bit more involved, but not much. First you connect a variable to the pipe from pbcopy.

> osxclipboard <- pipe("pbcopy", "w")

Now you can write to the variable "osxclipboard" and it will turn up on the clipboard so that you can paste it into another application. You have to remember that it is a pipe to an external file, which means that you have to use "write()" to put data in the clipboard. If you just assign some data to "osxclipboard" it will be changed from a pipe to that data. Thus:

> write(x, file="osxclipboard")

where "x" is your data, will to the trick. Then you just have to clean up using:

> close(osxclipboard)

In the same way as with copying from the clipboard, this will not format your data optimally. so you have to add some formatting command for Excel to read it properly.

> write(x, file="osxclipboard", sep = "\t")

There we are, but that is rather alot of code. Not really that much easier than a regular export from excel and then import to R. So here are two simple functions to do it all at once, with the most used arguments (I'm sorry for the code, I haven't figured out how to get indenting to work properly in blogger).

Simple readlines from the clipboard:

clipboard.readlines <- function (

clipboardname = "pbpaste",

n = -1L) {

readLines(con = pipe(clipboardname),

n = n)}

Read table from the clipboard. I have only included the most commonly used arguments to read.table(). You could easily add all the rest if you need them:

clipboard.readtable <- function (

clipboardname = "pbpaste",

header = FALSE,

sep = "",

quote = "\"'",

dec = ".") {

read.table(file = pipe(clipboardname),

header = header,

sep = sep,

quote = quote,

dec = dec)}

Write to the clipboard:

clipboard.write <- function (x,

clipboardname = "pbcopy",

sep = " ") {

osxclipboard <- pipe(clipboardname, "w")

write(x = x,

file = clipboard,

sep = sep)

Close(osxclipboard)}

That was it for now. I hope someone finds it helpful. I'm probably not putting any of this into a formal extension, so you will just have to run the code to use it. Consider the code GPL-licensed, although the text is under regular copyright.

Don't log your data twice

I am sitting here with some RNA expression microarray data from an Illumina Beadarray experiment. In general the data consist of a long [23k genes] list of spot intensities for each array in the experiment. If you want to compare any two of these lists in R using Bioconductor and the beadarray package you can run a small piece of code like this:

plotMAXY(exprs(BSData),
arrays = 1:2,
labels = targets$sampleID[1:2])

[plotMAXY is a plotting function; exprs extracts expressions from the BeadSummaryData (BSData); arrays tells the function which arrays to compare and labels tells it what to call them]

Then you get a graph like this:


Figure 1: On the lower left you see a XY-plot, where the intensities of all genes on array C1 are plotted against the intensities on array C2. The dashed lines show a two fold change in expression between arrays. To the upper right is the MA-plot (M is the ratio of intensities on the two arrays plotted on the y-axis, and A is the average intensity of a gene plotted on the x-axis).

What should strike you is the fact that one of the arrays appears to have a generally higher intensity than the other and that the difference changes with intensity, giving a banana-shaped curve. Some very smart people have shown that most of this is just an artifact of the method and that the average intensity of all genes in comparable samples should be the same and not vary with the expression. So what you do is normalise the data of all arrays to the same quantile distribution.

BSData.log2.quantile = normaliseIllumina(BSData,
method = "quantile",
transform = "log2")

[BSData.log2.quantile is the output variable; normaliseIllumina is the function that normalises data from Illumina arrays using the method given by method and performs the transformation given by transform; BSData is the BeadSummaryData]

The result should be a graph such as:


Figure 2: plotMAXY result after normalisation, when all is well.

You can see that the average intensities of the two arrays are now equal so that the cloud lies around the unity line in the XY-plot and the 0-line in the MA-plot. There is no longer any banana shape to the curve. But, there are a numer of spots outside the two-fold lines showing differentially expressed genes.

What I got when I plotted the normalised data this time was this plot:


Figure 3: Example of a plotMAXY result that shouldn't exist.

It's like, D..N. It shouldn't look like that. S..T! It's impossible that there should be such a curve. The XY-plot shows one sample to be impossibly much stronger than the other, but only at high intensities. What the F..K are the fold change lines doing at that kind of angle any way? The MA curve at least looks like it should, but show no differently expressed genes at all. B....Y V......E! Even technical replicates aren't that similar.

When I can hear myself above the beating of my heart and can see the screen again through the torrent of cold sweat pouring into my eyes I look more closely at the XY-plot and see that the scales on the x- and y-axis are different. Actually no intensities are higher than 4.

What has happened is that the plotMAXY function in log2-transforms the data by default before plotting (that's why the 2-fold line is called 1), and I already did the log2-transformation when I quantile-normalised the data. What I should have done is added the "log = FALSE" line to the plotMAXY, like this:

plotMAXY(exprs(BSData.log2.quantile),
arrays = 1:2,
labels = targets$sampleID[1:2],
log = FALSE)

[log = FALSE tells the function plotMAXY not to log2 transform the expression data from BSData.log2.quantile]

Now i get a graph more like figure 2 and can continue with my analysis without further risk of apoplexy. Or maybe not, data analysis is a forest fraught with formidable foes and fearsome formulas.

I leave you now to venture again into the jungle of statistical tricks.

Wish me luck,

Michael

The lady tasting tea

Statistics is one of the cornerstones of modern science, the other three being Blood, Sweat and Tears.

First, it is important to distinguish between the different things that are called statistics. There is mathematical statistics, or so I have heard, and there is applied statistics, which is what we are generally talking about. In applied statistics there is descriptive statistics, which we call plotting, and finally inferential statistics.

Inferential statistics is the area that most (junior) scientists mean when they say that they hate statistics. It is the area of hypothesis testing or finding of significant differences, and what people hate about it is that it is never clear which test you should use.

You will have learnt that you should predefine the tests you are looking for and design the experiment after that. The problem lies in having a too theoretical understanding of statistics and not having enough experience with what ever kind of data you possess. This means that you have not been able to predict what kind of analysis gives the most power with the data you end up with.

In the end, the first time you handle any type of dataset has to be considered an exploratory investigation. And, it is important to accept that and really get to know how the data responds to analysis and what kind of results you can get out of it. Then, next time, you can predetermine the correct type of analysis, and feel like you are on higher ground, theoretically.

So, how do you get to understand statistics, or at least start to understand it?

I have found no statistics text books that are worth while in the beginning. None. Many are useful once you have understood something.

The same goes for statistics courses. They are just not useful, especially not for understanding statistics. The most useful courses I have had in statistics were those that focused on a specific problem and a specific approach, i.e. "Using R and Bioconductor for the analysis of mRNA expression micro arrays". Even with those, the real help with understanding comes from working with your own data.

I suggest two things:

1. Read the wonderful book "The lady tasting tea" by David Salsburg describes the development of statistics as a science in a most delightful way. He manages to present many of the basic ideas of statistics without being technical. I can only recommend it, especially if you like the history of science.
2. Start using a statistics program that can help you understand what you do. I suggest R. It's free, it's what all the statisticians use (so you can ask them how to do something), it includes a manual, and it includes original references to many of the methods.

Yours,

Michael Hultström