Chapter 4 Differential gene analysis

Once you obtain gene counts from the preprocessing modules of RIMA, you can compare gene expression levels for different sample phenotypes. For example, you may wish to know which genes are differentially expressed between patients who respond and do not respond to a particular treatment. In addition, you may wish to know if particular gene sets are enriched in expression between different phenotypes. These gene sets may be associated with a known biological pathway, process or molecular function.

The differential gene analysis module of RIMA allows you to perform these comparisons.

4.1 Differential gene expression using DESeq2

DESeq2 is the most commonly used package that provides methods for comparison of gene expression levels. The package is described in:

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8

Details about how to run DeSeq2 are also included in the DESeq2 vignette.

4.1.1 Inputs of DESeq2

DESeq2 accepts four kinds of inputs: 1) transcript abundance files; 2) a read count matrix; 3) htseq-count files or 4) a SummarizedExperiment object.

Tximport is a tool that can import gene quantification data generated by Salmon (quant.sf) and convert the data into a transcript abundance format that DESeq2 can use. RIMA generated quant.sf files during the preprocessing module. vignette.

As a part of using tximport, Ensemble transcript IDs are converted into gene names. RIMA contains a tx2gene.csv file with transcript IDs and gene names as below:

ID,GENENAME
ENST00000000233.8,ARF5
ENST00000000412.6,M6PR
ENST00000000442.9,ESRRA
ENST00000001008.5,FKBP4
ENST00000001146.5,CYP26B1
ENST00000002125.7,NDUFAF7
ENST00000002165.9,FUCA2
ENST00000002501.9,DBNDD1
ENST00000002596.5,HS3ST1

Tximport uses the quant.sf files generated by Salmon and the tx2gene.csv file to generate matrices as described in DESeq2:

### filelist is a list of quant.sf files generated by Salmon gene quantification.
txi <- tximport(filelist, type=Salmon, tx2gene=tx2gene)

4.2 Different gene expression analysis

After creating the transcript abundance file, we are ready to find differentially expressed genes between your selected conditions. If you have identified a batch effect in your earlier analysis, DESeq2 can be run to account for this batch effect. We will first describe how to run DESeq2 without a batch effect. We will then describe how to run DESeq2 when a batch effect exists.

4.2.1 Running DESeq2 without batch effect

We will run DESeq2 to determine differential expression for a particular phenotype. In our case, we are interested in the genes that are differentially expressed between responders and non-responders. For this example, the Group column in the metasheet.csv has been set up with “R”, “NR” or “NA” as shown in Chapter 2.3 and the config.yaml file is set up with

  • the “batch” parameter set to false;
  • “design” set to Group
  • “Treatment” set to R
  • “Control” set to NR

as seen in Chapter 2.3.

This results in a DESeq call like the following:

ddsTxi <- DESeqDataSetFromTximport(txi,
                                       colData = samples,
                                       design = ~ Responder)

dds <- DESeq(ddsTxi)

4.2.2 Running DESeq2 with batch effect

DESeq2 can account for any batch effect by providing the batch variable as a covariate. If your batch effect analysis from the preprocessing module indicated that there is a batch effect in your samples, set the “batch” field in config.yaml to the appropriate column name in your metasheet.

#########Cohort level analysis parameters################
design: Group             # Condition on which to do comparsion (as set up in metasheet.csv)
Treatment: R              # Treatment use in DESeq2, corresponding to positive log fold change
Control: NR               # Control use in DESeq2, corresponding to negative log fold change
batch: syb_batch          # Options: 'false' or a column name from the metasheet.csv.

We use the information in “batch” for DESeq2 in the following manner, where “syn_batch” is the column designated in the config.yaml file:

ddsTxi <- DESeqDataSetFromTximport(txi,
                                       colData = samples,
                                       design = ~ syn_batch + Responder)
dds <- DESeq(ddsTxi)

4.2.3 Output from DESeq2

DESeq2 output includes values for baseMean, log2FoldChange, lfcSE, stat, pvalue, padj. These fields have the following meanings:

  • baseMean = the average of the normalized counts for all samples
  • Log2FoldChange = log base 2 fold changes for the condition tested. For example if the file is Responder_responder_VS_nonresponder_DESqe2_raw.txt, Log2FoldChange will be the estimated logarithmic fold change log2(responder/nonresponder).
  • lfcSE = Log2FoldChange Standard Error
  • stat= Wald statistic
  • pvalue = Wald test p value for the Log2FoldChange estimate
  • padj. = FDR adjusted p value (Benjamini-Hochberg adjusted p value)

{comparison condition}_DEseq2_raw.txt lists Gene names directly from the DESeq2 output. It also contains a value for -log10(padj).

4.2.4 Volcano Plot(s)

Volcano plots are often used to visualize the results of differential expression analysis. The x axis represents the logFC, and the y axis represent the -log10(P.value). Up-regulated genes have positive logFCs and down-regulated genes have negative logFCs.
Below is an example of a volcano plot from the previous differential expression analysis. Genes that meet the following criteria are labeled:

  • FDR < 0.05
  • Log2FoldChange > 1
Volcano plot

Figure 4.1: Volcano plot

4.3 Gene set enrichment analysis (GSEA)

GSEA was developed by UC San Diego and the Broad Institute. It is described in

Subramanian, Tamayo, et al. (2005, PNAS 102, 15545-15550) and Mootha, Lindgren, et al. (2003, Nat Genet 34, 267-273).

GSEA tools and software can be found on the GSEA website.

Gene Set Enrichment Analysis can be used to detect patterns in differential expression that affect particular gene pathways, molecular functions, cellular components or processes. RIMA uses the GSEA software to look for enrichment in the following gene sets:

  • KEGG canonical pathways
  • Gene Ontology sets – including Biological Process (BP), Cellular Component (CC) and Molecular Function (MF).

RIMA uses a ranked list of DESeq2 outputs. Ranking is based on the following:

sign(data$log2FoldChange) * (-log10(data$pvalue))

GSEA output is reported in 4 text files:

  • {comparison condition_phenotype1_VS_phenotype2}_KEGG_terms.txt
  • {comparison condition_phenotype1_VS_phenotype2}_GO_CC_terms.txt
  • {comparison condition_phenotype1_VS_phenotype2}_GO_BP_terms.txt
  • {comparison condition_phenotype1_VS_phenotype2}_GO_MF_terms.txt

Each of these tab separated files contains the following information:

  • ID = ID of the gene set
  • Description = description of the gene set
  • setSize = number of genes in the gene set after filtering out genes not in the expression data.
  • enrichmentScore = the degree to which the gene set is enriched in the top or bottom of the ranked set of genes.
  • NES = normalized gene enrichment score.
  • pvalue = nominal p value
  • p.adjust = family-wise error adjusted p value
  • qvalues = FDR adjusted p value
  • rank = the position in the ranked list at which the maximal enrichment score was achieved.
  • leading_edge = 3 statistics (tags, list and signal) on the subset of genes that most influence the enrichment score
  • core_enrichment group = the subset of genes that contribute most to the enrichment score.
  • significance = a representation of the p.adjust value using stars. One star is set for p.adjust <0.1, two stars <0.5 and three stars < 0.01.

Please see the GSEA documentation for more information about these values.

4.4 Single sample gene set enrichment analysis (ssGSEA)

ssGSEA is a different method of trying to find enriched pathways in your dataset. GSEA uses the differential expression patterns of all samples to generate enrichment scores based on phenotypic differences. In contrast, ssGSEA generates enrichment scores for each sample. Phenotypic differences can then be used to compare enriched pathways rather than differential gene expression levels. Please see ssGSEA for more information.

The input to ssGSEA is the tpm_convertID.batch file generated by limma. The file contains log2 transformed tpm values.

From the ssGSEA output, RIMA generates a heatmap of the gene sets that are most differentially enriched between the groups defined in the config.yaml file.