I needed to install poretools, but the setup instructions assume you have sudo priviliges. After combing through 879868070879 posts on stack overflow, the R-help mailing list, serverfault, etc., I finally arrived at solution using the nuclear option. There may be a better way, but this worked for me.
First, nuke and rebuild your own R from source. First, remove your existing installation (as usual, rm -rf with caution...)
rm -rf ~/R
| ## RNA-seq analysis with DESeq2 | |
| ## Stephen Turner, @genetics_blog | |
| # RNA-seq data from GSE52202 | |
| # http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse52202. All patients with | |
| # ALS, 4 with C9 expansion ("exp"), 4 controls without expansion ("ctl") | |
| # Import & pre-process ---------------------------------------------------- | |
| # Import data from featureCounts |
| # Basic work flows: from SYMBOL to gene model and sequence | |
| # Orig from http://master.bioconductor.org/help/course-materials/2014/ISMB2014/brca1.R | |
| library(org.Hs.eg.db) | |
| eid <- select(org.Hs.eg.db, "BRCA1", "ENTREZID", "SYMBOL")[["ENTREZID"]] | |
| library(TxDb.Hsapiens.UCSC.hg19.knownGene) | |
| txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene | |
| txid <- select(txdb, eid, "TXNAME", "GENEID")[["TXNAME"]] | |
| cds <- cdsBy(txdb, by="tx", use.names=TRUE) |
| // ==UserScript== | |
| // @name Shorten Amazon Links (with affiliate ID) for Google Chrome | |
| // @namespace http://userscripts.org/users/427315 | |
| // @description Shorten all links to Amazon products for easy emailing. (Strips Amazon associate referral tags and replaces with gettgenedone-20.) | |
| // This is just a modification of Callum Locke's script at http://userscripts.org/scripts/show/55185 | |
| // ...but adjusted so it will work for google chrome | |
| // On Chrome, manually set the associateID below, as greasemonkey (GM_) functions won't work. | |
| // On chrome, must start with --enable-easy-off-store-extension-install | |
| // Then drag/drop to chrome://extensions |
June 19, 2014
| # 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")) |
| Gene log2FoldChange pvalue padj | |
| DOK6 0.51 1.861e-08 0.0003053 | |
| TBX5 -2.129 5.655e-08 0.0004191 | |
| SLC32A1 0.9003 7.664e-08 0.0004191 | |
| IFITM1 -1.687 3.735e-06 0.006809 | |
| NUP93 0.3659 3.373e-06 0.006809 | |
| EMILIN2 1.534 2.976e-06 0.006809 |
| Gene log2FoldChange pvalue padj | |
| DOK6 0.51 1.861e-08 0.0003053 | |
| TBX5 -2.129 5.655e-08 0.0004191 | |
| SLC32A1 0.9003 7.664e-08 0.0004191 | |
| IFITM1 -1.687 3.735e-06 0.006809 | |
| NUP93 0.3659 3.373e-06 0.006809 | |
| EMILIN2 1.534 2.976e-06 0.006809 | |
| TPX2 -0.9974 2.097e-06 0.006809 | |
| LAMA2 -1.425 2.39e-06 0.006809 | |
| CAV2 -1.052 3.213e-06 0.006809 |
| # Set up to run in GNU parallel. | |
| find *bam | parallel 'bedtools coverage -hist -abam {} -b target_regions.bed | grep ^all > {}.hist.all.txt' | |
| # What this actually runs: | |
| bedtools coverage -hist -abam samp.01.bam -b target_regions.bed | grep ^all > samp.01.bam.hist.all.txt | |
| bedtools coverage -hist -abam samp.02.bam -b target_regions.bed | grep ^all > samp.02.bam.hist.all.txt | |
| bedtools coverage -hist -abam samp.03.bam -b target_regions.bed | grep ^all > samp.03.bam.hist.all.txt | |
| bedtools coverage -hist -abam samp.05.bam -b target_regions.bed | grep ^all > samp.05.bam.hist.all.txt | |
| bedtools coverage -hist -abam samp.06.bam -b target_regions.bed | grep ^all > samp.06.bam.hist.all.txt | |
| bedtools coverage -hist -abam samp.10.bam -b target_regions.bed | grep ^all > samp.10.bam.hist.all.txt |
| # Assumes you've already run coverageBed -hist, and grep'd '^all'. E.g. something like: | |
| # find *.bam | parallel 'bedtools -abam {} -b capture.bed -hist | grep ^all > {}.all.txt' | |
| # Get a list of the bedtools output files you'd like to read in | |
| print(files <- list.files(pattern="all.txt$")) | |
| # Optional, create short sample names from the filenames. | |
| # For example, in this experiment, my sample filenames might look like this: | |
| # prefixToTrash-01.pe.on.pos.dedup.realigned.recalibrated.bam | |
| # prefixToTrash-02.pe.on.pos.dedup.realigned.recalibrated.bam |