Script to allow gene set enrichment analysis of 10x genomics data in R

Question

I have 10x single cell RNA seq data. Which R package is best suited for analysis of the 10x data matrix. What is the script to prepare the data for downstream GSEA analysis.

I have already processed samples for single cell data (for multiple other purposes). my intention is to focus on subpopulations i.e. cell type x at different time points i.e. time 0, 1, 2 etc and compare the enriched gene pathways in this one cell type at different time points. there are three files that are used (at least in the R Seurat package) with the following file extensions: .mtx, .tsv (barcodes.tsv, genes.tsv and matrix.mtx). the R Seurat package merges these files to create a "Seurat object" which helps mainly with identifying individual diffl gene expression in each barcode cell

I found a package for R: goseq that at least has been used for bulkRNA seq but I was hoping someone has used this package specifically for single cell RNA seq data, and what script has been used to merge the 3 files to create an object for analysis. pardon my terminology, but I am a novice to R and have experience primarily through package use.

Answer 1

Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.

Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable. The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.

Answer 2

There is no purpose-built R package to perform gene set enrichment analysis on single-cell data but there does not need to be. You should be able to tools developed for bulk-RNA-Seq or microarray data, although you may not get as significant results from a sparse scRNA-Seq matrix as single-cell technologies have poor sensitivity and miss genes. What you need for a gene set enrichment analysis is:

• a database of pathways (or gene sets) to test enrichment for. Packages such as GO.db or reactome.db from Bioconductor support these in R.

• a gene expression matrix which you can get from the output of cellranger (as matrix.mtx). This should be normalised and corrected for batch effects. The Seurat R package is recommended to do this. Calling object@data will return the normalised seurat data matrix to pass it to other tools.

• an annotation to identify the genes (depending on the reference data you used for the cellranger "transcriptome").

• to define groups for differential expression analysis. You can import the clusters generated by cellranger in outs/analysis or export clusters of your choosing from the Loupe Browser as CSV files that can be imported into R. These should match the barcodes.tsv of columns of your data matrix to define the cluster for each cell barcode.

• to run a differential expression analysis to compare clusters. You can use clusters to define cell groups and perform analysis with tools you would use for bulk-RNA-Seq such as limma, edgeR, or deseq. Don't forget to adjust for multiple tests. You also need to correct for batch effects between replicates and log-transform you data to ensure the results reflect biological differences. This will give you a list of differentially expressed genes.

• to calculate gene set enrichment of differentially expressed genes in a gene set list such as those supported by pathway databases. You can do this like any other list of significant genes with tools used for bulk-RNA-Seq data, either with dedicated webtools or R packages such as safe or GSEA. You can also run the hypergeometric test in R yourself.

The cellranger output gives you all the data you need to use standard tools for which there are already guides on how to do differential analysis. These analyses do not differ from bulk-RNA-Seq analysis once the data has been normalised and the clusters have been defined.

[edit:corrected typo in package name]

Answer 3

I realized that your question is specifically gene set enrichment analysis. However, I personally have never tried any R package doing this. Although topGO may seem to be the one that may work for you? What I have been doing is that for gene enrichment analysis using Gene Ontology, simply paste and run, or DAVID when it was still updating. Sorry, I am not an expert in incorporating biological significance study in the computational result.

It seems that you are also not quite sure what tools you should use for downstream analysis. I think the following publication may be helpful.

https://www.ncbi.nlm.nih.gov/pubmed/30228881

The same primary author also published a paper titled, Cluster Headache: Comparing Clustering Tools for 10X Single Cell Sequencing Data (For some reason I could not find the link to it, will edit the answer if I find it. )

Hope that helps.