ToPASeq 1.14.1
Note: the ToPASeq package currently undergoes a major rework due to the change of the package maintainer. It is recommended to use the topology-based methods implemented in the EnrichmentBrowser or the graphite package instead.
We start by loading the package.
library(ToPASeq)For RNA-seq data, we consider transcriptome profiles of four primary human airway smooth muscle cell lines in two conditions: control and treatment with dexamethasone Himes et al., 2014.
We load the airway dataset
library(airway)
data(airway)For further analysis, we only keep genes that are annotated to an ENSEMBL gene ID.
airSE <- airway[grep("^ENSG", rownames(airway)),]
dim(airSE)## [1] 63677 8assay(airSE)[1:4,1:4]## SRR1039508 SRR1039509 SRR1039512 SRR1039513
## ENSG00000000003 679 448 873 408
## ENSG00000000005 0 0 0 0
## ENSG00000000419 467 515 621 365
## ENSG00000000457 260 211 263 164The EnrichmentBrowser package incorporates established
functionality from the limma package for differential expression
analysis.
This involves the voom transformation when applied to RNA-seq data.
Alternatively, differential expression analysis for RNA-seq data can also be
carried out based on the negative binomial distribution with edgeR
and DESeq2.
This can be performed using the function EnrichmentBrowser::deAna
and assumes some standardized variable names:
For more information on experimental design, see the limma user’s guide, chapter 9.
For the airway dataset, the GROUP variable indicates whether the cell lines have been treated with dexamethasone (1) or not (0).
airSE$GROUP <- ifelse(airway$dex == "trt", 1, 0)
table(airSE$GROUP)##
## 0 1
## 4 4Paired samples, or in general sample batches/blocks, can be defined via a
BLOCK column in the colData slot.
For the airway dataset, the sample blocks correspond to the four different cell
lines.
airSE$BLOCK <- airway$cell
table(airSE$BLOCK)##
## N052611 N061011 N080611 N61311
## 2 2 2 2For RNA-seq data, the deAna function can be used to carry out differential
expression analysis between the two groups either based on functionality from
limma (that includes the voom transformation), or
alternatively, the frequently used edgeR or DESeq2
package. Here, we use the analysis based on edgeR.
library(EnrichmentBrowser)
airSE <- deAna(airSE, de.method="edgeR")## Excluding 50740 genes not satisfying min.cpm thresholdrowData(airSE, use.names=TRUE)## DataFrame with 12937 rows and 3 columns
## FC ADJ.PVAL edgeR.STAT
## <numeric> <numeric> <numeric>
## ENSG00000000003 -0.404945626610932 0.00213458295385677 35.8743710016452
## ENSG00000000419 0.182985434777531 0.0915691945173217 5.90960619951562
## ENSG00000000457 0.0143477674070905 0.922279475398735 0.0233923316993606
## ENSG00000000460 -0.141173372957313 0.619013213521584 0.492929955080683
## ENSG00000000971 0.402240426474171 0.00403820532305421 27.8509962017613
## ... ... ... ...
## ENSG00000273270 -0.12979385333726 0.495892935815196 0.901598359265221
## ENSG00000273290 0.505580471641003 0.00639218387702814 23.0905678847793
## ENSG00000273311 0.00161557580855132 0.996356136959404 8.04821151395742e-05
## ENSG00000273329 -0.222817127090519 0.388294594068834 1.42723325850574
## ENSG00000273344 0.0151704005097405 0.962777106053456 0.005435032737617Pathways are typically represented as graphs, where the nodes are genes and edges between the nodes represent interaction between genes.
The graphite package provides pathway collections from major pathway databases such as KEGG, Biocarta, Reactome, and NCI.
Here, we retrieve human KEGG pathways.
library(graphite)
pwys <- pathways(species="hsapiens", database="kegg")
pwys## KEGG pathways for hsapiens
## 301 entries, retrieved on 27-04-2018As the airway dataset uses ENSEMBL gene IDs, but the nodes of the pathways are based on NCBI Entrez Gene IDs,
nodes(pwys[[1]])## [1] "ENTREZID:10327" "ENTREZID:124" "ENTREZID:125" "ENTREZID:126"
## [5] "ENTREZID:127" "ENTREZID:128" "ENTREZID:130" "ENTREZID:130589"
## [9] "ENTREZID:131" "ENTREZID:160287" "ENTREZID:1737" "ENTREZID:1738"
## [13] "ENTREZID:2023" "ENTREZID:2026" "ENTREZID:2027" "ENTREZID:217"
## [17] "ENTREZID:218" "ENTREZID:219" "ENTREZID:220" "ENTREZID:2203"
## [21] "ENTREZID:221" "ENTREZID:222" "ENTREZID:223" "ENTREZID:224"
## [25] "ENTREZID:226" "ENTREZID:229" "ENTREZID:230" "ENTREZID:2538"
## [29] "ENTREZID:2597" "ENTREZID:26330" "ENTREZID:2645" "ENTREZID:2821"
## [33] "ENTREZID:3098" "ENTREZID:3099" "ENTREZID:3101" "ENTREZID:387712"
## [37] "ENTREZID:3939" "ENTREZID:3945" "ENTREZID:3948" "ENTREZID:441531"
## [41] "ENTREZID:501" "ENTREZID:5105" "ENTREZID:5106" "ENTREZID:5160"
## [45] "ENTREZID:5161" "ENTREZID:5162" "ENTREZID:5211" "ENTREZID:5213"
## [49] "ENTREZID:5214" "ENTREZID:5223" "ENTREZID:5224" "ENTREZID:5230"
## [53] "ENTREZID:5232" "ENTREZID:5236" "ENTREZID:5313" "ENTREZID:5315"
## [57] "ENTREZID:55276" "ENTREZID:55902" "ENTREZID:57818" "ENTREZID:669"
## [61] "ENTREZID:7167" "ENTREZID:80201" "ENTREZID:83440" "ENTREZID:84532"
## [65] "ENTREZID:8789" "ENTREZID:92483" "ENTREZID:92579" "ENTREZID:9562"we first map the gene IDs in the airway dataset from ENSEMBL to ENTREZ IDs.
airSE <- idMap(airSE, org="hsa", from="ENSEMBL", to="ENTREZID")## Loading required package: org.Hs.eg.db## Loading required package: AnnotationDbi## ## 'select()' returned 1:many mapping between keys and columns## Excluded 1133 genes without a corresponding to.ID## Encountered 8 from.IDs with >1 corresponding to.ID (a single to.ID was chosen for each of them)Next, we define all genes with adjusted p-value below 0.01 as differentially expressed, and collect their log2 fold change for further analysis.
all <- names(airSE)
de.ind <- rowData(airSE)$ADJ.PVAL < 0.01
de <- rowData(airSE)$FC[de.ind]
names(de) <- all[de.ind]This results in 2,426 DE genes - out of 11,780 genes in total.
length(all)## [1] 11795length(de)## [1] 2426The Pathway Regulation Score (PRS) incorporates the pathway topology by weighting the indiviudal gene-level log2 fold changes by the number of downstream DE genes. The weighted absolute fold changes are summed across the pathway and statistical significance is assessed by permutation of genes. Ibrahim et al., 2012
res <- prs(de, all, pwys[1:100], nperm=100)## 9 pwys were filtered outhead(res)## nPRS p.value
## Arachidonic acid metabolism 12.553137 0.00990099
## cGMP-PKG signaling pathway 10.836744 0.00990099
## Histidine metabolism 10.452106 0.00990099
## Linoleic acid metabolism 9.905204 0.00990099
## cAMP signaling pathway 8.612164 0.00990099
## Glycine, serine and threonine metabolism 7.888651 0.00990099Corresponding gene weights (number of downstream DE genes) can be obtained for a pathway of choice via
ind <- grep("Ras signaling pathway", names(pwys))
weights <- prsWeights(pwys[[ind]], de, all)
weights## 572 5898 3479 10928 3082 9846 4254
## 0 6 0 4 5 0 0
## 5601 5291 387 11186 10235 10681 2255
## 0 0 0 1 21 0 5
## 998 5566 9462 5337 5063 2260 2782
## 0 0 22 0 1 0 0
## 8315 25 5594 6655 6789 2254 5595
## 0 0 0 0 0 0 0
## 3551 208 8605 10156 4233 5599 8036
## 0 0 1 0 0 0 0
## 5878 4790 5602 2549 5869 2784 55770
## 0 0 1 1 0 0 0
## 7422 5924 2246 5159 59345 6654 5906
## 5 22 0 0 23 0 0
## 5321 10000 8503 5228 7010 5290 81579
## 1 0 4 0 0 4 0
## 5335 6237 2002 5605 5908 2791 91860
## 0 21 1 0 0 0 0
## 5338 25759 4301 10298 5894 3845 22800
## 0 0 0 0 2 0 21
## 5156 2113 5879 51196 2250 2247 2252
## 2 0 2 21 0 0 0
## 3480 8844 207 1969 5567 27 23179
## 2 0 0 0 23 1 7
## 805 5899 5868 56034 5295 5921 1956
## 0 0 0 5 4 0 0
## 4915 53358 5058 7424 284 5578 5922
## 0 1 0 5 5 0 0
## 8817 3815 2114 22808 808 5900 6464
## 0 0 0 0 0 7 0
## 2264 5966 2277 382 3481 5604 5881
## 2 0 0 0 0 1 2
## 100271927 80310 3643 598 5293 2783 55970
## 0 0 0 0 4 23 0
## 5970 7423 2787 3265 9610 11145 2788
## 0 0 23 0 2 0 0
## 627 2885 5781 5062 29110 1946 1435
## 0 0 1 0 0 0 0
## 8398 4303 4908 22821 54331 5320 4763
## 0 0 0 0 23 1 0
## 8831 5154 801 4893 1147 5863 1945
## 0 0 23 0 0 0 5Inspecting the genes with maximum number of downstream DE genes
weights[weights == max(weights)]## 59345 5567 2783 2787 54331 801
## 23 23 23 23 23 23reveals important upstream regulators including several G protein subunits such as subunit beta 2 (Entrez Gene ID 2783).