"The hunt for the genetic roots of common diseases has hit a blank wall."
...quoting the first sentence in Nicholas Wade's New York Times article reviewing this PLoS Biology research paper by David Goldstein and his colleagues at Duke University. Also be sure to see Richard Robinson's synopsis of this paper, both published this week in PLoS Bio.
The argument made by Goldstein and colleagues is that natural selection has done much better at eliminating disease-predisposing variants than we originally thought. Although we've found over 2000 disease-gene associations to common variants, it's likely that many of the common variants are tagging rare, ungenotyped variants. Yet most of the time we look at the gene nearest our top hits and try to spin a story about how that gene relates to our phenotype. As a proof of principle, the authors perform a simulation study using sickle-cell anemia as a model. Although the disease is caused by a variant in a single gene, the authors found statistically significant associations with 179 SNPs, spread across several megabases of DNA, containing many genes - and most of these SNPs were clearly pointing at the wrong thing.
New York Times: A New Way to Look for Diseases’ Genetic Roots
Synopsis (PLoS Biology): Common Disease, Multiple Rare (and Distant) Variants
Research article (PLoS Biology): Rare Variants Create Synthetic Genome-Wide Association
Abstract: Genome-wide association studies (GWAS) have now identified at least 2,000 common variants that appear associated with common diseases or related traits (http://www.genome.gov/gwastudies), hundreds of which have been convincingly replicated. It is generally thought that the associated markers reflect the effect of a nearby common (minor allele frequency >0.05) causal site, which is associated with the marker, leading to extensive resequencing efforts to find causal sites. We propose as an alternative explanation that variants much less common than the associated one may create “synthetic associations” by occurring, stochastically, more often in association with one of the alleles at the common site versus the other allele. Although synthetic associations are an obvious theoretical possibility, they have never been systematically explored as a possible explanation for GWAS findings. Here, we use simple computer simulations to show the conditions under which such synthetic associations will arise and how they may be recognized. We show that they are not only possible, but inevitable, and that under simple but reasonable genetic models, they are likely to account for or contribute to many of the recently identified signals reported in genome-wide association studies. We also illustrate the behavior of synthetic associations in real datasets by showing that rare causal mutations responsible for both hearing loss and sickle cell anemia create genome-wide significant synthetic associations, in the latter case extending over a 2.5-Mb interval encompassing scores of “blocks” of associated variants. In conclusion, uncommon or rare genetic variants can easily create synthetic associations that are credited to common variants, and this possibility requires careful consideration in the interpretation and follow up of GWAS signals.
Seminar: Translating Disease-Associate Genetic Variants into Pathways
Soumya Raychaudhuri is the lead author on the GRAIL paper in Nature Genetics. GRAIL (Gene Relationships Across Implicated Loci) is software available from the Broad Institute and looks like an interesting way to prioritize SNPs for followup. I'll cover it here in the future, but in the meantime, check out this seminar this Wednesday.
Title: Translating Disease-Associate Genetic Variants into Pathways
Location: 214 Light Hall
Date: January 27, 2010 12:00pm - 1:00pm
Speaker: Soumya Raychaudhuri, MD, PhD
Fellow, Division of Rheumatology, Allergy, and Immunology
Brigham and Women's Hospital
Boston, MA
Title: Translating Disease-Associate Genetic Variants into Pathways
Location: 214 Light Hall
Date: January 27, 2010 12:00pm - 1:00pm
Speaker: Soumya Raychaudhuri, MD, PhD
Fellow, Division of Rheumatology, Allergy, and Immunology
Brigham and Women's Hospital
Boston, MA
Fluctuation plot using ggplot2 in R
Found this nice way to visually summarize contingency tables using ggplot2 in R on Hadley Wickham's ggplot2 cheat sheet. Using the same data in my previous post on making scatterplots in small multiples, I'll demonstrate how to use ggfluctuation() to make a fluctuation plot. For a two-dimensional table, this simply puts different dimensions on the x and y axes, and the area of the rectangles is proportional to the density of observations in that cell of the table. Install ggplot2 as in the previous post, and run this code.
The table looks like this:
And the plot:
As a side note, I'm experimenting with a new way to post R code. If you hover your mouse over the code above you can click this button to copy all the text to your clipboard.
# Load ggplot2. Use install.packages("ggplot2") if you do not already have ggplot2 installed.
library(ggplot2)
# Load the diamonds dataset.
data(diamonds)
# Look at a few rows of the relevant columns
diamonds[1:20,2:3]
# Display the contingency table
table(diamonds$cut,diamonds$color)
# Make the fluctuation plot
ggfluctuation(table(diamonds$cut,diamonds$color))
The table looks like this:
D E F G H I J
Fair 163 224 312 314 303 175 119
Good 662 933 909 871 702 522 307
Very Good 1513 2400 2164 2299 1824 1204 678
Premium 1603 2337 2331 2924 2360 1428 808
Ideal 2834 3903 3826 4884 3115 2093 896
And the plot:
As a side note, I'm experimenting with a new way to post R code. If you hover your mouse over the code above you can click this button to copy all the text to your clipboard.
GWAS Manhattan plots and QQ plots using ggplot2 in R
*** Update April 25, 2011: This code has gone through a major revision. Please see the updated code and tutorial here. ***
Will posted earlier this week about how to produce manhattan plots of GWAS results using Stata, so I thought I'd share how I do this in R using ggplot2.
First, if you've never used ggplot2, you'll need to add it to your R installation by typing:
Once you've done that, copy and paste this command to download the functions I wrote necessary to produce these plots. If you'd like to see the source code yourself, copy the URL into your web browser.
Next, read in the PLINK results file to a data frame. Substitute plink.qassoc with your own results filename.
Finally, run this function which will produce and save in the current directory both a QQ-plot and a manhattan plot of your results:
QQ-Plot (these are simulated GWAS results):
Manhattan plot (click to enlarge):
A few notes: First, if you're doing this on a linux machine from your Windows computer, you'll need to be running the previously mentioned XMing server on your Windows computer for the plot to save correctly. Second, the qqman() function calls the manhattan() function, which is extremely slow and memory-intensive. It may take about 3 minutes to run for each dataset. The memory issue isn't a problem on 64-bit machines, but you may run out of memory on 32-bit machines if you're doing this with GWAS data. I'm going to try to improve this in the future. Finally, using that source command you also downloaded a function I wrote called qqmanall(), which does just what it sounds like - if you run it on a linux machine with no arguments it reads in ALL of the plink GWAS results stored in the current directory, and creates QQ and manhattan plots for all of them with a common upper limit for the y-axis corresponding to the most significant result. Enjoy.
...
Update Thursday, January 21, 2010: I neglected to mention yesterday the format of the plink.assoc or plink.qassoc files, in case you want to produce the same plots using results from another software other than plink. When you load your .assoc files in a data frame, the relevant columns are named "CHR", "BP", and "P". You can use this plotting function as long as you have these three columns in your data frame, regardless of whether you use PLINK or not.
The latest version of the code is reproduced below:
*** Update April 25, 2011: This code has gone through a major revision. Please see the updated code and tutorial here. ***
...
*** Update April 25, 2011: This code has gone through a major revision. Please see the updated code and tutorial here. ***
Will posted earlier this week about how to produce manhattan plots of GWAS results using Stata, so I thought I'd share how I do this in R using ggplot2.
First, if you've never used ggplot2, you'll need to add it to your R installation by typing:
install.packages("ggplot2")
Once you've done that, copy and paste this command to download the functions I wrote necessary to produce these plots. If you'd like to see the source code yourself, copy the URL into your web browser.
source("http://dl.dropbox.com/u/66281/0_Permanent/qqman.r")
Next, read in the PLINK results file to a data frame. Substitute plink.qassoc with your own results filename.
mydata=read.table("plink.qassoc", header=TRUE)
Finally, run this function which will produce and save in the current directory both a QQ-plot and a manhattan plot of your results:
qqman(mydata)
QQ-Plot (these are simulated GWAS results):
Manhattan plot (click to enlarge):
A few notes: First, if you're doing this on a linux machine from your Windows computer, you'll need to be running the previously mentioned XMing server on your Windows computer for the plot to save correctly. Second, the qqman() function calls the manhattan() function, which is extremely slow and memory-intensive. It may take about 3 minutes to run for each dataset. The memory issue isn't a problem on 64-bit machines, but you may run out of memory on 32-bit machines if you're doing this with GWAS data. I'm going to try to improve this in the future. Finally, using that source command you also downloaded a function I wrote called qqmanall(), which does just what it sounds like - if you run it on a linux machine with no arguments it reads in ALL of the plink GWAS results stored in the current directory, and creates QQ and manhattan plots for all of them with a common upper limit for the y-axis corresponding to the most significant result. Enjoy.
...
Update Thursday, January 21, 2010: I neglected to mention yesterday the format of the plink.assoc or plink.qassoc files, in case you want to produce the same plots using results from another software other than plink. When you load your .assoc files in a data frame, the relevant columns are named "CHR", "BP", and "P". You can use this plotting function as long as you have these three columns in your data frame, regardless of whether you use PLINK or not.
The latest version of the code is reproduced below:
*** Update April 25, 2011: This code has gone through a major revision. Please see the updated code and tutorial here. ***
...
*** Update April 25, 2011: This code has gone through a major revision. Please see the updated code and tutorial here. ***
Seminar: Cristen Willer - Mapping Genes for Complex Genetic Diseases: Lessons from Obesity and Lipid Levels
Christen J. Willer, (of the Willer et al Nature Genetics lipid genetics paper), postdoc at University of Michigan and departmental Faculty Candidate, will be giving a talk tomorrow (Wednesday Jan 20) at noon in 208 Light Hall. I'll be going to this one.
Topic: Mapping Genes for Complex Genetic Diseases: Lessons from Obesity and Lipid Levels
Speaker: Cristen J. Willer, PhD, Postdoctoral Research Fellow, Faculty Candidate
Institution: University of Michigan
Seminar Location: 208 Light Hall
Time: Noon-1pm
Sponsor: MPB Center for Human Genetic Research and Diabetes Research Training Center
Topic: Mapping Genes for Complex Genetic Diseases: Lessons from Obesity and Lipid Levels
Speaker: Cristen J. Willer, PhD, Postdoctoral Research Fellow, Faculty Candidate
Institution: University of Michigan
Seminar Location: 208 Light Hall
Time: Noon-1pm
Sponsor: MPB Center for Human Genetic Research and Diabetes Research Training Center
Genome-wide Manhattan Plots in STATA
I have used this chunk of code on numerous occasions to plot GWAS data, so I thought I'd share!
Enjoy!
twoway scatter logpval absposMB, msize(tiny) ysize(2) xsize(7) xaxis(1 2) yline(16.881438 2.9957323, lcolor(red))The variables needed are a log p-value (or some other statistic) and the absolute genomic position of each SNP (distance from the beginning of chromosome 1). If you need the offsets to compute this absolute position, they are listed in MB in the xline(---) portion of the plot. It makes something like this:
xline(245.522847 488.541076 688.046816 879.458034 1060.3159 1231.291599 1389.919738 1536.194564 1674.623832 1810.03746 1944.489844 2076.939655 2191.082635 2297.45122 2397.790135 2486.617389 2565.392131 2641.509284 2705.320935 2767.756899 2814.701222 2864.255932, lpattern(dash) lcolor(gs12))
xlabel(122 "1" 367 "2" 588 "3" 783 "4" 969 "5" 1145 "6" 1310 "7" 1463 "8" 1605 "9" 1742 "10" 1877 "11" 2010 "12" 2134
"13" 2244 "14" 2347 "15" 2442 "16" 2526 "17" 2603 "18" 2673 "19" 2736 "20" 2791 "21" 2839 "22" 2941 "X", axis(2))
xlabel( 0(50)3000, axis(1) alternate) subtitle(,bcolor(white) ring(30) orientation(vertical) position(9) nobexpand) ytitle("-Log p-value") xtitle("Physical Distance (Mb)", axis(1))
Enjoy!
Coming to R from SQL, Python, SAS, Matlab, or Lisp
Head over to Revolutions Blog for a list of PDF and powerpoint resources for making the transition to R from other programming or stats languages. All of these notes come from the New York R meetup. I enjoyed browsing the meetup's files - lots of powerpoints, PDFs, and example R data files for various topics, including several slideshows on ggplot2. Don't forget the ggplot2 tutorial I posted here earlier this week if you're completely new to ggplot2.
If you're coming from SPSS or SAS, I've read good things about Robert Muenchen's book R for SAS and SPSS Users, which can be found on Amazon (~$50), or downloaded for free from http://rforsasandspssusers.com/.
Revolutions: Coming to R from SQL, Python, SAS, Matlab, or Lisp
If you're coming from SPSS or SAS, I've read good things about Robert Muenchen's book R for SAS and SPSS Users, which can be found on Amazon (~$50), or downloaded for free from http://rforsasandspssusers.com/.
Revolutions: Coming to R from SQL, Python, SAS, Matlab, or Lisp
Subscribe to:
Posts (Atom)