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

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.

• Welcome to the group. Can you specify what format you data is in and what format you need it in for the downstream GSEA. What tools have you tried or looked into? – Bioathlete Oct 17 '18 at 19:14
• Hi Jay, its not clear from your question what the biological question you are trying to address with a GSEA analysis. Do you want do do enrichment of the genes that are expresesd on average in the cells of this experiement? If this is the case, how would this differ from an analysis of bulk sequencing? Do you wish to know the pathways that differ between different scRNAseq samples? Do you wish to know the pathways that differ between clusters of cells within one sample? The biological question at hand will strongly influence the downstream analysis approach. – Ian Sudbery Oct 18 '18 at 13:08
• In addition to Ian question, you mention goseq. Are you aware of the differences between pathways and gene ontologies terms? Do you have in mind any pathway database or you don't care (there are several with different focus). Are we talking about human pathways or other organisms? For other organisms there are fewer databases of pathways or less consensus. Keep updating the question. I think it can become a great question for future readers – llrs Oct 19 '18 at 12:20

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.

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]

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.

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.

• DAVID is widely regarded as outdated. The databases that it uses are not kept up to date. This is not recommended. – Tom Kelly Nov 23 '18 at 7:36
• @TomKelly That is absolutely right. That's why I said it was good when it was updated. Not recommended until it starts to update again. – BC Wang Nov 29 '18 at 4:34