Working with a sparse matrix in R
Sparse matrices are necessary for dealing with large single-cell RNA-seq
datasets. They require less memory than dense matrices, and they allow some
computations to be more efficient. In this note, we'll discuss the internals of
dgCMatrix class with examples.
Harmony in motion: visualize an iterative algorithm for aligning multiple datasets
Harmony is a an algorithm for aligning multiple high-dimensional datasets, described by Ilya Korsunsky et al. in this paper. When analyzing multiple single-cell RNA-seq datasets, we often encounter the problem that each dataset is separated from the others in low dimensional space – even when we know that all of the datasets have similar cell types. To address this problem, the Harmony algorithm iteratively clusters and adjusts high dimensional datasets to integrate them on a single manifold. In this note, we will create animated visualizations to see how the algorithm works and develop an intuitive understanding.
Make tidy variant tables with MyVariant.info and Tabulator
Make tidy gene tables with MyGene.info and Tabulator
Extract data from a PDF file with Tabula
Kirkham et al. 2006 is a prospective 2-year study of 60 patients with rheumatoid arthritis (RA). It shows that “synovial membrane cytokine mRNA expression is predictive of joint damage progression in RA”. The PDF includes a few tables with data on cytokine measurements and correlations with joint damage. Here, we'll use Tabula to extract data from tables in the PDF file. Then we'll make figures with R.
Make heatmaps in R with pheatmap
Here are a few tips for making heatmaps with the pheatmap R package by Raivo Kolde. We'll use quantile color breaks, so each color represents an equal proportion of the data. We'll also cluster the data with neatly sorted dendrograms, so it's easy to see which samples are closely or distantly related.
Color points by density with ggplot2
Here, we use the 2D kernel density estimation function from the MASS R package to to color points by density in a plot created with ggplot2. This helps us to see where most of the data points lie in a busy plot with many overplotted points.
Build bioinformatics pipelines with Snakemake
Snakemake is a Pythonic variant of GNU Make. Recently, I learned how to use it to build and launch bioinformatics pipelines on an LSF cluster. However, I had trouble understanding the documentation for Snakemake. I like to learn by trying simple examples, so this post will walk you through a very simple pipeline step by step. If you already know how to use Snakemake, then you might be interested to copy my Snakefiles for RNA-seq data analysis here.
Determine if a transcription factor is bound to a genomic site with CENTIPEDE
I wrote a practical tutorial for how to use CENTIPEDE to determine if a transcription factor is bound to a site in the genome. The tutorial explains how to prepare appropriate input data and how to run the analysis. Please get in touch if you have any comments or suggestions.
Print bigWig data for each region in a BED file
I wrote a Bash script to call bigWigToBedGraph for each region in a BED file. You can quickly take a subset of bigWig data for regions of interest. In my particular case, I needed to get phastCons conservation scores for putative transcription factor binding sites.
Get human transcription factor target genes
I made a data package with human transcription factor target genes for use in R. It is a collection of data from three sources: TRED, ITFP, and ENCODE. I use them to test if the targets of a transcription factor are differentially expressed in my data. Also, I can test if a set of transcription factor target genes is enriched for some gene set of interest.
Quickly aggregate your data in R with data.table
In genomics data, we often have multiple measurements for each gene. Sometimes we want to aggregate those measurements with the mean, median, or sum. The data.table R package can do this quickly with large datasets.
In this note, we compute the average of multiple measurements for each gene in a gene expression matrix.
Make ribosomal RNA intervals for Picard CollectRnaSeqMetrics
Before you can use the CollectRnaSeqMetrics Picard tool, you must create a table of genomic intervals with the coordinates of all ribosomal genes in the genome. I wrote a bash script to prepare this ribosomal interval file from Gencode gene annotations.
Autocomplete gene names with mygene.info and typeahead.js
We can use mygene.info with typeahead.js to autocomplete gene names and retrieve every annotation you can think of (GO, Kegg, Ensembl, position, homologs, etc.). Try typing your gene name in the interactive example in this post, and download my code.
Create a quantile-quantile plot with ggplot2
After performing many tests for statistical significance, the next step is to check if any results are more extreme than we would expect by random chance. One way to do this is by comparing the distribution of p-values from our tests to the uniform distribution with a quantile-quantile (QQ) plot. Here's a function to create such a plot with ggplot2.
How to ssh to a remote server without typing your password
Here are a few tips to use
ssh more effectively. Login to your server using
public key encryption instead of typing a password. Use the
file to create short and memorable aliases for your servers. Also, use aliases
to connect through a login server into a work server.
Count the number of coding base pairs in each Gencode gene
We can use Python to count the coding base pairs in each Gencode gene. Here, we report the base pair count by gene rather than by transcript. When we encounter different transcripts for the same gene with overlapping exons, we only count those base pairs once rather than multiple times.
0-based and 1-based genomic intervals, overlap, and distance
Here, I describe two kinds of genomic intervals and include source code for testing overlap and calculating distance between intervals.