Analysing molecular data

In this practical session we will go through basics of preprocessing and visualizing molecular data (gene expression data). The data we use is single-cell RNA-seq data from healthy human brain samples. link The same dataset is used in the deconvolution project work and because of that, the main point of this practical session is for you to get familiar with handling the data, applying minimal preprocessing, visualizing the data, and finding differentially expressed genes.

Once gene expression has been quantified it is summarized as an expression matrix where each row corresponds to a gene (or transcript) and each column corresponds to a single cell. This matrix should be examined to remove poor quality cells which were not detected in either read QC or mapping QC steps. Failure to remove low quality cells at this stage may add technical noise which has the potential to obscure the biological signals of interest in the downstream analysis. Since there is currently no standard method for performing scRNASeq the expected values for the various QC measures that will be presented here can vary substantially from experiment to experiment. Thus, to perform QC we will be looking for cells and genes which are outliers with respect to the rest of the dataset rather than comparing to independent quality standards. Consequently, care should be taken when comparing quality metrics across datasets collected using different protocols.

Load data

We load the downloaded data, which is separated to csv-files per cell. The individual files are combined to form a data matrix where rows correspond to genes and columns correspond to individual cells. For convinience we modify the column names to match how the sample names are written in the annotation file. Also at this point we remove last three rows which correspond to some basic stats from the alignment step (not actual genes).

Load annotations

Annotations (metadata) corresponding the said data is loaded and ordered to match the ordering of samples in the data matrix.

Filtering cells/genes

Total number of RNA molecules detected per sample (library size) is import to check and remove cells with low library size. Wells with few reads/molecules are likely to have been broken or failed to capture a cell, and should thus be removed. Also, cells expressing unusually high number of reads might consist of multiple cells (multiplets) and should be removed. Note: Be careful when excluding cells in this manner, you might remove good cells (rare subtypes?)

In addition to ensuring sufficient sequencing depth for each sample, we also want to make sure that the reads are distributed across the transcriptome. Thus, we count the total number of unique genes detected in each sample and remove genes given a threshold.

In some experiments, control RNAs, e.g. ERCC spike-in RNAs, may be used. Spike-in RNAs are good for (at least) two things:

1. Cell quality can be measured as the ratio between ERCC spike-in RNAs and endogenous RNAs.
   • This ratio can be used to estimate the total amount of RNA in the captured cells. Cells with a high level of spike-in RNAs had low starting amounts of RNA, likely due to the cell being dead or stressed which may result in the RNA being degraded.
2. Normalization of molecule counts between cells. Unfortunately, the dataset we are using does not have any spike-ins, so running this step is meaningless.

Another useful QC metric is the ratio of reads aligned onto mitochrondrial genes per cell. High percentages may indicate broken cells and thus should be removed. However, no mitochondrial genes are expressed are found in the data and this step is skipped.

In addition to removing cells with poor quality, it is usually a good idea to exclude genes where we suspect that technical artefacts may have skewed the results. Moreover, inspection of the gene expression profiles may provide insights about how the experimental procedures could be improved.

It is typically a good idea to remove genes whose expression level is considered “undetectable”. We define a gene as detectable if at least two cells contain more than 5 transcripts from the gene. However, the threshold strongly depends on the sequencing depth. It is important to keep in mind that genes must be filtered after cell filtering since some genes may only be detected in poor quality cells.

Save non-normalized data for later usage.

Normalization

Library sizes vary because scRNA-seq data is often sequenced on highly multiplexed platforms the total reads which are derived from each cell may differ substantially. Some quantification methods (eg. Cufflinks, RSEM) incorporated library size when determining gene expression estimates thus do not require this normalization. However, if another quantification method was used then library size must be corrected for by multiplying or dividing each column of the expression matrix by a "normalization factor" which is an estimate of the library size relative to the other cells. Many methods to correct for library size have been developped for bulk RNA-seq and can be equally applied to scRNA-seq (eg. UQ, SF, CPM, RPKM, FPKM, TPM).

The simplest way to normalize this data is to convert it to counts per million (CPM) by dividing each column by its total then multiplying by 1,000,000. Note that spike-ins, if used, should be excluded from the calculation of total expression in order to correct for total cell RNA content. Also note that we are adding 1 to the data before normalization because we will log-transform the data eventually!

One potential drawback of CPM is if your sample contains genes that are both very highly expressed and differentially expressed across the cells. In this case, the total molecules in the cell may depend of whether such genes are on/off in the cell and normalizing by total molecules may hide the differential expression of those genes and/or falsely create differential expression for the remaining genes.

Log-transformation is pretty much a standard procedure in molecular data. Molecule count distribution is often skewed because cells not expressing a gene have low read counts but cells expressing the same genes have exponentially higher read counts. E.g. take a look at the distribution of FBXW7 before and after log-transformation. Note that the bimodal distribution seen in the histogram suggests that there are subpopulations of cells expressing FBXW7.

Feature selection

After preprocessing we usually want to see an overview of the data by visualizing it. But before visualization we want to find differentially expressing genes that help us discriminate different cell types from the data. One simple way is to order genes by their dispersion and use n top genes. High dispersion indicates the gene is expressed in different magnitudes among the cell population and could indicate existence of different cell types. Downside using this approach is that generally highly expressed genes are favored over lowly expressed. Quick fix to this is divide the dispersion by mean expression, i.e. calculating coefficient of variation.

Principal Component Analysis (PCA)

The easiest way to overview the data is by transforming it using the principal component analysis and then visualize the first two principal components. Principal component analysis (PCA) is a statistical procedure that uses a transformation to convert a set of observations into a set of values of linearly uncorrelated variables called principal components (PCs). The number of principal components is less than or equal to the number of original variables.

t-distributed Stochastic Neighbor Embedding (t-SNE)

An alternative to PCA for visualizing scRNASeq data is a tSNE plot. tSNE (t-Distributed Stochastic Neighbor Embedding) combines dimensionality reduction (e.g. PCA) with random walks on the nearest-neighbour network to map high dimensional data to a 2- or 3-dimensional space while preserving local distances between cells. In contrast with PCA, tSNE is a stochastic algorithm which means running the method multiple times on the same dataset will result in different plots. Due to the non-linear and stochastic nature of the algorithm, it is more difficult to intuitively interpret a tSNE visualization To ensure reproducibility, we fix the “seed” of the random-number generator in the code below so that we always get the same plot.

Dimensionality reduction evaluation

While visual inspection is useful, evaluation using quality metrics is also important. Here we assess the quality of the embedding with trustworthiness and continuity (unsupervised metrics) and with nearest neighbor error (supervised metric).

Modelling dropouts and identifying differentially expressed (DE) genes

Single-cell RNA sequencing is able to quantify the whole transcriptome from the small amount of RNA present in individual cells. However, a consequence of reverse-transcribing and amplifying small quantities of RNA is a large number of dropouts, genes with zero expression in particular cells. The frequency of dropout events is strongly non-linearly related to the measured expression levels of the respective genes. M3Drop posits that these dropouts are due to failures of reverse transcription, a simple enzyme reaction, thus should be modelled using the Michaelis-Menten equation as follows: P_i = 1 - S_i / (S_i + K) where P_i is the proportion of cells where gene i dropouts out, S_i is the mean expression of gene i and K is the Michaelis constant. The first step is to clean the data by remove cells with too few detected genes, genes that with very low expression, and to normalize the data. This can be done using any method but M3Drop includes a simple function that will clean the expression matrix and convert raw counts to counts per million (CPM). If alternative normalization methods are used the input expression matrix must not be log-transformed, nor contain negative values. If normalization adjusts zero values then M3Drop will use the minimal expression value in the entire matrix as the value for dropouts.

Next, we can compare the fits of different possible functions relating the proportion of dropouts to the average expression for each gene.

Visual inspection of the resulting plot shows that the Michaelis-Menten equation is the best fit to the data. However, we can also examine some statistics to confirm this:

Here we actually see that statistics favor the logistic curve (smaller value better)

Since the Michaelis-Menten equation is concave, averaging across a mixed population forces differentially expressed genes to be shifted to the right of the Michaelis-Menten curve. DE genes are identified by comparing the local K calculated for a specific gene to the global K fitted across all genes using a Z-test followed by multiple-testing correction

Note that this function runs directly from the expression matrix, hence one could skip straight to identifying DE genes without comparing models and any external normalisation method can be applied to the raw counts prior to DE analysis.

To check that the identified genes are truly differentially expressed we can plot the normalised expression of the genes across cells.

To get certain clustering out from the dendrogram of the heatmap:

We can then test all genes as marker genes for our cell type labels (or the provided clusters).

Marker genes are ranked by the area-under the ROC curve (AUC) for predicting the population with the highest expression of the gene from the other groups. Significance is calculated using a Wilcox-rank-sum test. Now we can examine the marker genes of different cell types more closely.