Blog for Sean Davis

flexible bioinformatics pipelines with snakemake

2013-05-10T04:00:00Z

Building bioinformatics pipelines can be a somewhat tedious process that seemingly never ends and often results in significant refactoring when assumptions or workflows change. The snakemake project is a small framework that facilitates building flexible, reuseable pipelines.

I am not going to go through the installation process in this entry, but since snakemake requires python3, suffice it to say that virtualenvs almost certainly will be involved since most systems default to python2.X.

I had a small revelation last night about how snakemake can be used to run the GATK pipeline. Like many folks, I use filename “tags” to track individual transformations and ordering of the tags in the filename to capture the ordering of the transformations.

Tool	Added filename extension
mark duplicates	md
quality recalibration	recal
indel realignment	realigned

So, a file with the name “stuff.realigned.recal.md.bam” has been realigned, then recalibrated, then had duplicates marked. Different orderings of tags like “md”, “realigned”, and “recal” represent different processing orders. A simple snakefile with rules for these transformations might look like this:

rule qualrecal:
    input: bam="{base}.bam"
    output: bam="{base}.recal.bam"
    shell: "touch {output}"

rule indelrealigner:
    input: bam="{base}.bam",intervals="{base}.intervals"
    output: bam="{base}.realigned.bam"
    shell: "touch {output}"

rule realignertargetcreator:
    input: "{base}.bam"
    output: "{base}.intervals"
    shell: 'touch {output}'

rule indexbam:
    input: '{base}.bam'
    output: '{base}.bam.bai'
    shell: 'touch {output}'

rule markdups:
    input: '{base}.bam'
    output: bam='{base}.md.bam',metrics='{base}.dupmetrics'
    shell: "touch {output.bam}; touch {output.metrics}"

I have not put any shell commands in this file just to keep things simple. Obviously, this Snakefile will not produce any meaningful output, but it is enough to be instructive.

If I add a rule, “final” that contains the output filename(s) of interest and then run snakemake -n, we can see how the finename(s) dictate the run order for the rules. This allows one to change the pipeline by simply adjusting output filenames. Quite cool!

rule final:
    input: "stuff.realigned.md.bam", "stuff.realigned.md.bam.bai"

Running snakemake -n with this rule included as the first rule in the snakefile yields:

rule realignertargetcreator:
        input: stuff.bam
        output: stuff.intervals
rule indelrealigner:
        input: stuff.bam, stuff.intervals
        output: stuff.realigned.bam
rule markdups:
        input: stuff.realigned.bam
        output: stuff.realigned.md.bam, stuff.realigned.dupmetrics
rule indexbam:
        input: stuff.realigned.md.bam
        output: stuff.realigned.md.bam.bai
rule final:
        input: stuff.realigned.md.bam, stuff.realigned.md.bam.bai
Job counts:
        count	jobs
        1	final
        1	indelrealigner
        1	realignertargetcreator
        1	indexbam
        1	markdups
        5

Changing the rule “final” to this:

rule final:
    input: "stuff.md.realigned.bam", "stuff.md.realigned.bam.bai"

and running snakemake -n again gives a different ordering of steps.

rule markdups:
        input: stuff.bam
        output: stuff.md.bam, stuff.dupmetrics
rule realignertargetcreator:
        input: stuff.md.bam
        output: stuff.md.intervals
rule indelrealigner:
        input: stuff.md.bam, stuff.md.intervals
        output: stuff.md.realigned.bam
rule indexbam:
        input: stuff.md.realigned.bam
        output: stuff.md.realigned.bam.bai
rule final:
        input: stuff.md.realigned.bam, stuff.md.realigned.bam.bai
Job counts:
        count	jobs
        1	final
        1	realignertargetcreator
        1	markdups
        1	indelrealigner
        1	indexbam
        5

Adding recalibration to the process is as simple as including “.recal” in the filename in the correct order. This works with one bam file or 1000 bam files. Nice….

autoreload in ipython to ease development

2013-05-09T04:00:00Z

New versions of ipython (after 0.12) ship with a convenient autoreload functionality. This is just a note-to-self about usage.

From the ipython documentation:

In [1]: %load_ext autoreload

In [2]: %autoreload 2

In [3]: from foo import some_function

In [4]: some_function()
Out[4]: 42

In [5]: # open foo.py in an editor and change some_function to return 43

In [6]: some_function()
Out[6]: 43

And here are the usage notes:

%autoreload: Reload all modules (except those excluded by %aimport) automatically now.
%autoreload 0: Disable automatic reloading.
%autoreload 1: Reload all modules imported with %aimport every time before executing the Python code typed.
%autoreload 2: Reload all modules (except those excluded by %aimport) every time before executing the Python code typed.
%aimport: List modules which are to be automatically imported or not to be imported.
%aimport foo: Import module ‘foo’ and mark it to be autoreloaded for %autoreload 1
%aimport -foo: Mark module ‘foo’ to not be autoreloaded.

Testing for Non-Differentially-Expressed Genes

2013-03-22T04:00:00Z

There has been an ongoing discussion on the Bioconductor mailing list dealing with a commonly-asked question (at least by folks I work with) but uncommonly answered. "How can I test for genes that are statistically not changing?"

Albyn Jones mentions that a common approach for tests of equivalence is the two-one-sided-tests approach that was apparently described by DJ Shuirmann in A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability.

Simon Anders answered the question quite thoroughly in text like so. And I quote:

One of the new features of DESeq2 is that it provides a empirical-Bayes style shrinkage estimates of coefficients (i.e., log fold changes) and also estimates standard errors for these coefficients (taken from the the reciprocal curvature of the posterior). Your task is one of the application we had in mind for this. You want to know whether the fold change is zero or nearly zero, and you also want to be ascertained that this estimate of a nearly-zero fold change is a precise one. So, to find genes whose absolute log fold change is _reliably_ smaller than some threshold, take all those genes for which the estimated abs log fold change is below the threshold and the standard error is well below the threshold, too. The reason, why I write "below a threshold" rather than "equal to zero" is that it's biologically unreasonable to assume that a treatment has exactly zero effect on a gene's expression. Gene regulation is such an interconnected network that, at least in my view, every gene will react ever so slightly to every perturbance, but often, this reaction is too small to be measurable. This is why you should define a threshold for "biologically unlikely to be relevant". Let's say, we don't believe that an expression change of less than 15% (0.2 on a log2 scale) is worth bothering. So, if you might take all genes with abs log2 fold change below 0.2, and furthermore require that the standard error of this log2 fold change estimate is below, say, 0.05, than you will get gene with true values well below at least 0.3 or so. In the end, the thresholds won't matter too much but if you want real p values, this is possible, too: Simply divide the log2 fold change estimates by their standard errors to get z scores (this works because the sampling distribution of our fold change estimates is reasonably close to normal), and do two one-sided tests (TOST).

Simon goes on to give code to do this analysis using the DESeq2 bioconductor package. He introduces the code with this comment:

To recap, you want to test the null hypothesis that your log fold change beta is larger than some threshold theta, i.e., you want to find genes for which you can reject this null hypothesis and hence say with some certainty that |beta| < theta. The TOST scheme suggests, as its name implies to perform two pone sided tests, one with the null hypothesis H0_A: beta > theta, the other with H0_B: beta < -theta, and then use the larger of the two p values.

The code is so clear and well-commented that it speaks for itself.

# Let's use the example data from the pasilla package
library( DESeq2 )
library( pasilla )
data( "pasillaGenes" )

# Create a DESeq2 data object from the pasilla data
dse <- DESeqSummarizedExperimentFromMatrix(
   counts(pasillaGenes), pData(pasillaGenes)[,c("condition","type")],
   ~ type + condition )

# Perform a standard DESeq2 analysis
dse <- DESeq(dse)

# The log2 fold changes are found here
beta <- results(dse)$log2FoldChange

# Just to make the following clearer, I should point out that the
# "results" accessor is just a short-cut for this access to the rowData:
all( beta == mcols(rowData(dse))$conditionuntreated, na.rm=TRUE )
# (returns TRUE)

# The log fold change estimates all come with standard error information
# which we find in the rowData (maybe we should copy this to the
# 'results', too)
betaSE <- mcols(rowData(dse))$SE_conditionuntreated

# Internally, the Wald test is implemented as a simple two-sided
# z test of beta/betaSE. Two demonstrate this, we to the test
# manually and compare
pvalDE <- 2 * pnorm( abs( beta ), sd = betaSE, lower.tail=FALSE )
all( abs( pvalDE - results(dse)$pvalue ) < 1e-15, na.rm=TRUE )
# (returns TRUE)

# This was the test for DE, of course, i.e., small pvalDE means that
# the gene's expression change (the true value of beta) is not zero

# What we want is the opposite, namely find gene, for which abs(beta)
# is smaller than some threshold, theta
theta <- .3

# So, we do our two one-sided tests. For a one-sided z test, we
# simply use tail probabilities from the normal distribution.

# First, the test of H0_A: true_beta > thr
pA <- pnorm( beta, thr, betaSE, lower.tail=TRUE )

# Next, the test of H0_B: true_beta < -thr
pB <- pnorm( beta, -thr, betaSE, lower.tail=FALSE )

# The overall p value is the maximum, because we want to reject H0_A
# and H0_B simultaneously
pvalTOST <- pmax( pA, pB )


# Let's adjust our two p values with BH:
sigDE <- p.adjust( pvalDE, "BH" ) < .1
sigSmall <- p.adjust( pvalTOST, "BH" ) < .1

# And make an MA plot, with sigDE in red and sigSmall in green
plot(
   rowMeans( counts(dse,normalized=TRUE) ), beta,
   log="x", pch=20, cex=.2,
   col = 1 + sigDE + 2*sigSmall )
# Plot is attached.

The plot that is generated by the last few lines shows in green the genes that are, indeed, statistically likely to be unchanged between conditions.

sqlalchemy with python 3 and mysql

2013-03-21T04:00:00Z

I have a few sqlalchemy-based projects that I have wanted to move over to python 3. Interestingly, the usual python-mysql is not supported. There are a number of other options, but I settled on Oracle’s mysql-connector-python. The usual pip incantation gets you the library (pure python).

pip install mysql-connector-python

Using the library is a matter of using the right connection string:

"mysql+mysqlconnector://..."

Latex conversion notes

2012-02-28T05:00:00Z

I have been avoiding latex lately because of a desire to publish documents in HTML as well as PDF and my sense was that converters from latex to other formats were lacking. After watching several beautiful presentations given based on the latex beamer class and being immersed in a Bioconductor course where vignettes are the coin of the realm, I decided to look a bit more fully. Since I work on a Mac much of the time, I started poking around in /usr/texbin and found a command, htlatex that looked promising.

To convert a .tex file to html is as simple as:

htlatex filename.tex

This, by default, produces HTML4 output as file "filename.html". To request xhtml output is as simple as supplying an argument:

htlatex filename.tex "xhtml"

Noting that OpenOffice is XML-based, someone has written a configuration file for converting latex to a file that OpenOffice (or, my current flavor, LibreOffice) can read and display.

htlatex filename.tex "xhtml,ooffice" "ooffice/! -cmozhtf" "-coo -cvalidate"

The little incarnation above creates an .odt file that opens directly in OpenOffice. Of course, OpenOffice can then convert to many other formats. Conversion to MS Word is also possible.

htlatex filename "html,word" "symbol/!" "-cvalidate"