I've analyzed some data from GEO (GSE52202) using RNA-seq to study gene expression in motor neurons differentiated from induced pluripotent stem cells (iPSCs) derived from ALS patients carrying the C9ORF72 repeat expansion. I aligned the data, counted with featureCounts, and analyzed with DESeq2. I uploaded the results to this GitHub Gist.
Here's how you can use R to create a simple volcano plot. First, download the results file here and save it as a text file called results.txt.
After reading in the data from GitHub the next section creates a basic volcano plot. A few more lines color the points based on their fold change and statistical significance. Finally, if you have the calibrate package installed, the last line labels a few genes of interest.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Download the data from github (click the "raw" button, save as a text file called "results.txt"). | |
# https://gist.github.com/stephenturner/806e31fce55a8b7175af | |
res <- read.table("results.txt", header=TRUE) | |
head(res) | |
# Make a basic volcano plot | |
with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2))) | |
# Add colored points: red if padj<0.05, orange of log2FC>1, green if both) | |
with(subset(res, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red")) | |
with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="orange")) | |
with(subset(res, padj<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="green")) | |
# Label points with the textxy function from the calibrate plot | |
library(calibrate) | |
with(subset(res, padj<.05 & abs(log2FoldChange)>1), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=.8)) |
Hi,
ReplyDeletethanks for sharing this ! It is very informative.
I am nonetheless bothered by the overlapping gene labels. Do you know how to solve this issue ?
thanks for your help !
Julien
You might take a look at Rcharts http://rcharts.io/
ReplyDeleteWhen I ran this command in R:
ReplyDeleteres <- read.gist("806e31fce55a8b7175af", header=TRUE)
it says: no "read.gist" function. What should I do? Thanks!
Ah. I removed the read.gist function from the package for some compatibility issues. I updated the text in the blog post.
DeleteThanks. It works now.
DeleteCould you please help me? I would like to add to the graph cut-off lines.
ReplyDeleteThank you,en-joy
See ?abline
ReplyDeletehttps://stat.ethz.ch/R-manual/R-devel/library/graphics/html/abline.html
Hello,
ReplyDeleteThis is very informative and useful. However, get this error whenever I run,
>bwadeseqres <- read.table("examplefile.txt", header=TRUE)
Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, :
line 102 did not have 5 elements
> head(bwadeseqres)
with(bwadeseqres, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2)))
Error in plot(log2FoldChange, -log10(pvalue), pch = 20, main = "Volcano plot", :
object 'log2FoldChange' not found
Kindly help!
Hi, thanks for the post. It was very useful and works perfect. Do you know a way to increase the size of the dots?
ReplyDeleteThanks