RNA-Seq datasets

TissueEnrich defines tissue-specific genes using RNA-Seq data from the HPA, GTEx, and mouse ENCODE. In order to make the tissue-specific gene calculations more robust, we only used tissues that had ≥2 biological replicates. The datasets used in the tool are:

• HPA Dataset: RNA-Seq data across 35 human tissues (Uhlén et al. 2015).
• GTEx Dataset: RNA-Seq data across 29 human tissues (Ardlie et al. 2015).
• Mouse ENCODE Dataset: RNA-Seq data across 17 mouse tissues (Shen et al. 2012).

When using teEnrichment, the user can specify the RNA-Seq dataset (rnaSeqDataset) to be used for the tissue-specific gene enrichment analysis.

• 1 for “Human Protein Atlas” (default)
• 2 for “GTEx”
• 3 for “Mouse ENCODE”

Note: The tissues isolated from embryonic stages are prefixed with ‘E’ followed by the timepoint. For example in Mouse ENCODE data, the placenta tissue isolated from embryonic day 14.5 is named as E14.5-Placenta. All the other tissues are isolated from adults.

Defining Tissue-specific Genes

Tissue-specific genes are defined using the algorithm from the HPA (Uhlén et al. 2015), and can be grouped as follows:

• Tissue Enriched: Genes with an expression level greater than 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a particular tissue compared to all other tissues.
• Group Enriched: Genes with an expression level greater than 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a group of 2-7 tissues compared to all other tissues, and that are not considered Tissue Enriched.
• Tissue Enhanced: Genes with an expression level greater than 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a particular tissue compared to the average levels in all other tissues, and that are not considered Tissue Enriched or Group Enriched.

In teEnrichment, the user can specify the type of tissue-specific genes (tissueSpecificGeneType) to be used for the tissue-specific gene enrichment analysis.

• 1 for “All” (default)
• 2 for “Tissue-Enriched”
• 3 for “Tissue-Enhanced”
• 4 for “Group-Enriched”

Hypergeometric test

The hypergeometric test is used to calculate tissue-specific gene enrichment. The p-value is calculated as:

\[ P(X \gt k) = \sum\limits_{i=k+1}^n \frac{{{K}\choose{i}} {{N-K}\choose{n-i}}}{{{N}\choose{n}}} \] and the fold-change is calculated as: \[ Fold-change =\left( \frac{k}{n} \right)/\left( \frac{K}{N} \right) \]

Where, N is the total number of genes, K is the total number of tissue-specific genes for a tissue, n is the number of genes in the input gene set, k is the number of tissue-specific genes in the input gene set. The p-values can be corrected for multiple hypothesis testing using the Benjamini & Hochberg correction by setting multiHypoCorrection = TRUE (It is TRUE by default).

Background genes

TissueEnrich now enables users to provide the background genes for carrying out tissue-specific gene enrichment. In this case, instead of using all the genes in the dataset, a background gene set is being used to carry out the enrichment analysis. It should be noted that the background gene set must have all the genes of the input gene set. The p-value is calculated as: \[ P(X \gt k) = \sum\limits_{i=k+1}^n \frac{{{K_b}\choose{i}} {{N_b-K_b}\choose{n-i}}}{{{N_b}\choose{n}}} \] and the fold-change is calculated as: \[ Fold-change =\left( \frac{k}{n} \right)/\left( \frac{K_b}{N_b} \right) \] Where, \(N_b\) is the total number of background genes, \(K_b\) is the total number of tissue-specific genes for a tissue in background genes, n is the number of genes in the input gene set, k is the number of tissue-specific genes in the input gene set. The p-values can be corrected for multiple hypothesis testing using the Benjamini & Hochberg correction by setting multiHypoCorrection = TRUE (It is TRUE by default). If the background gene set is not provided all the genes will be used as background.

Example: Tissue-specific gene enrichment

This example uses trophectoderm (TE) specific genes identified from single cell RNA-Seq analyses, performed on human blastocysts on days 5, 6, and 7 of preimplantation development (Petropoulos et al. 2016). The single cells are assigned to either the inner cell mass (epiblast plus emerging extraembryonic endoderm) or the TE using PCA. After that, a list of 100 TE-specific genes was generated using differential gene expression analysis (Petropoulos et al. 2016). We used those 100 genes as the input gene set and carried out tissue-specific gene enrichment using the tissue-specific genes defined by the HPA dataset.

Note: The input gene set can either contain Ensembl Ids (ENSEMBLIdentifier()) or Gene Symbols (SymbolIdentifier()) (specify this using the geneIdType parameter in the input GeneSet object).

library(TissueEnrich)
genes<-system.file("extdata", "inputGenes.txt", package = "TissueEnrich")
inputGenes<-scan(genes,character())
gs<-GeneSet(geneIds=inputGenes,organism="Homo Sapiens",geneIdType=SymbolIdentifier())
output<-teEnrichment(inputGenes = gs)

The output is a list object containing the enrichment results. These results are explained in the next section.

seEnrichmentOutput<-output[[1]]
enrichmentOutput<-setNames(data.frame(assay(seEnrichmentOutput),row.names = rowData(seEnrichmentOutput)[,1]), colData(seEnrichmentOutput)[,1])
enrichmentOutput$Tissue<-row.names(enrichmentOutput)
head(enrichmentOutput)
#>                 Log10PValue Tissue.Specific.Genes fold.change
#> Adipose Tissue    0.0000000                     1   1.2394054
#> Adrenal Gland     0.0000000                     0   0.0000000
#> Appendix          1.0631035                     4   4.4829558
#> Bone Marrow       0.5552268                     4   2.7907142
#> Breast            0.0000000                     0   0.0000000
#> Cerebral Cortex   0.0000000                     3   0.4191623
#>                          Tissue
#> Adipose Tissue   Adipose Tissue
#> Adrenal Gland     Adrenal Gland
#> Appendix               Appendix
#> Bone Marrow         Bone Marrow
#> Breast                   Breast
#> Cerebral Cortex Cerebral Cortex
ggplot(enrichmentOutput,aes(x=reorder(Tissue,-Log10PValue),y=Log10PValue,label = Tissue.Specific.Genes,fill = Tissue))+
      geom_bar(stat = 'identity')+
      labs(x='', y = '-LOG10(P-Adjusted)')+
      theme_bw()+
      theme(legend.position="none")+
      theme(plot.title = element_text(hjust = 0.5,size = 20),axis.title = element_text(size=15))+
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),panel.grid.major= element_blank(),panel.grid.minor = element_blank())

ggplot(enrichmentOutput,aes(x=reorder(Tissue,-fold.change),y=fold.change,label = Tissue.Specific.Genes,fill = Tissue))+
      geom_bar(stat = 'identity')+
      labs(x='', y = 'Fold change')+
      theme_bw()+
      theme(legend.position="none")+
      theme(plot.title = element_text(hjust = 0.5,size = 20),axis.title = element_text(size=15))+
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),panel.grid.major= element_blank(),panel.grid.minor = element_blank())

library(tidyr)
seExp<-output[[2]][["Placenta"]]
exp<-setNames(data.frame(assay(seExp), row.names = rowData(seExp)[,1]), colData(seExp)[,1])
exp$Gene<-row.names(exp)
exp<-exp %>% gather(Tissue=1:(ncol(exp)-1))

ggplot(exp, aes(key, Gene)) + geom_tile(aes(fill = value),
     colour = "white") + scale_fill_gradient(low = "white",
     high = "steelblue")+
     labs(x='', y = '')+
      theme_bw()+
      guides(fill = guide_legend(title = "Log2(TPM)"))+
      #theme(legend.position="none")+
      theme(plot.title = element_text(hjust = 0.5,size = 20),axis.title = element_text(size=15))+
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),panel.grid.major= element_blank(),panel.grid.minor = element_blank())

seGroupInf<-output[[3]][["Placenta"]]
groupInf<-data.frame(assay(seGroupInf))
print(head(groupInf))
#>      Gene           Group
#> 1     CGA Tissue-Enriched
#> 2    GCM1 Tissue-Enriched
#> 3 CYP19A1 Tissue-Enriched
#> 4   GPR32 Tissue-Enriched
#> 5  CLEC1A Tissue-Enriched
#> 6 SLC13A4 Tissue-Enriched

print(geneIds(output[[4]]))
#> [1] "C10ORF54" "CGB"      "GRAMD3"   "PVRL4"

Orthologous gene enrichment

The teEnrichment function can calculate mouse tissue-specific gene enrichment from a human gene list or vice versa. The user simply specifies if the input data is from mouse or human, and selects the tissue-specific gene data of interest, whether it is from mouse or human. The function will automatically carry out orthologous tissue-specific gene enrichment using one-to-one protein coding orthologous genes between human and mouse, downloaded from Ensembl V91 database (Aken et al. 2016).

library(TissueEnrich)
library(ggplot2)
genes<-system.file("extdata", "inputGenes.txt", package = "TissueEnrich")
inputGenes<-scan(genes,character())
gs<-GeneSet(geneIds=inputGenes,organism="Homo Sapiens",geneIdType=SymbolIdentifier())
output<-teEnrichment(inputGenes = gs,rnaSeqDataset = 3)
seEnrichmentOutput<-output[[1]]
enrichmentOutput<-setNames(data.frame(assay(seEnrichmentOutput), row.names = rowData(seEnrichmentOutput)[,1]), colData(seEnrichmentOutput)[,1])
enrichmentOutput$Tissue<-row.names(enrichmentOutput)
ggplot(enrichmentOutput,aes(x=reorder(Tissue,-Log10PValue),y=Log10PValue,label = Tissue.Specific.Genes,fill = Tissue))+
      geom_bar(stat = 'identity')+
      labs(x='', y = '-LOG10(P-Adjusted)')+
      theme_bw()+
      theme(legend.position="none")+
      theme(plot.title = element_text(hjust = 0.5,size = 20),axis.title = element_text(size=15))+
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),panel.grid.major= element_blank(),panel.grid.minor = element_blank())

Gene groups

The genes are divided into six groups based on their gene expression across the tissues. These groups are:

• Not Expressed: Genes with an expression level less than 1 (TPM or FPKM) across all the tissues.
• Tissue Enriched: Genes with an expression level greater than or equal to 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a particular tissue compared to all other tissues.
• Group Enriched: Genes with an expression level greater than or equal to 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a group of 2-7 tissues compared to all other tissues, and that are not considered Tissue Enriched.
• Tissue Enhanced: Genes with an expression level greater than or equal to 1 (TPM or FPKM) that also have at least five-fold higher expression levels in a particular tissue compared to the average levels in all other tissues, and that are not considered Tissue Enriched or Group Enriched.
• Expressed in all: Genes with an expression level greater than or equal to 1 (TPM or FPKM) across all of the tissues that are not in any of the above 4 groups.
• Mixed: Genes that are not assigned to any of the above 5 groups.

Genes from the Tissue Enriched, Group Enriched, and Tissue Enhanced groups are classified as tissue-specific genes.

Example: Tissue-specific gene retrieval

In the example below, we supplied a subset of mouse ENCODE data, consisting of expression data of 36 genes across 17 tissues.

library(TissueEnrich)
library(SummarizedExperiment)
data<-system.file("extdata", "test.expressiondata.txt", package = "TissueEnrich")
expressionData<-read.table(data,header=TRUE,row.names=1,sep='\t')
se<-SummarizedExperiment(assays = SimpleList(as.matrix(expressionData)),rowData = row.names(expressionData),colData = colnames(expressionData))
output<-teGeneRetrieval(se)
head(assay(output))
#>      Gene                 Tissue           Group             
#> [1,] "ENSMUSG00000003200" "All"            "Expressed-In-All"
#> [2,] "ENSMUSG00000003206" "Bone.Marrow"    "Tissue-Enhanced" 
#> [3,] "ENSMUSG00000003208" "All"            "Mixed"           
#> [4,] "ENSMUSG00000004530" "All"            "Expressed-In-All"
#> [5,] "ENSMUSG00000004535" "All"            "Expressed-In-All"
#> [6,] "ENSMUSG00000004540" "E14.5.Placenta" "Tissue-Enriched"

As seen above, the output consists of the tissue-specific genes information in a SummarizedExperiment object with columns for Gene name, Tissue name, and Tissue-Specific group.

Example: Tissue-specific gene enrichment in custom dataset

The example uses 10 genes, randomly selected from the 36 genes used in the Tissue-specific gene retrieval example. The tissue-specific genes identified from the custom gene expression are used to calculate tissue-specific gene enrichment in the input gene set.

library(TissueEnrich)
library(ggplot2)
genes<-system.file("extdata", "inputGenesEnsembl.txt", package = "TissueEnrich")
inputGenes<-scan(genes,character())
gs<-GeneSet(geneIds=inputGenes)
output2<-teEnrichmentCustom(gs,output)
enrichmentOutput<-setNames(data.frame(assay(output2[[1]]), row.names = rowData(output2[[1]])[,1]), colData(output2[[1]])[,1])
ggplot(enrichmentOutput,aes(x=reorder(Tissue,-Log10PValue),y=Log10PValue,label = Tissue.Specific.Genes,fill = Tissue))+
      geom_bar(stat = 'identity')+
      labs(x='', y = '-LOG10(P-Value)')+
      theme_bw()+
      theme(legend.position="none")+
      theme(plot.title = element_text(hjust = 0.5,size = 20),axis.title = element_text(size=15))+
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),panel.grid.major= element_blank(),panel.grid.minor = element_blank())

As seen above, the metadata of the output consists of the tissue-specific gene enrichment information in a list object containing the enrichment results.