1 Example data

In this example, we use the ?airway data set.

library("airway")
data("airway")

2 Annotating data

This section demonstrates one of many possible workflows for adding annotations to the data set. Those annotations are meant to make it easier for users to identify genes of interest, e.g. by displaying both gene symbols and ENSEMBL gene identifiers as tooltips in the interactive browser.

First, we make a copy of the Ensembl identifiers – currently stored in the rownames() – to a column in the rowData() component.

rowData(airway)[["ENSEMBL"]] <- rownames(airway)

Then, we use the org.Hs.eg.db package to map the Ensembl identifiers to gene symbols. We store those gene symbols as an additional column of the rowData() component.

library("org.Hs.eg.db")
rowData(airway)[["SYMBOL"]] <- mapIds(
    org.Hs.eg.db, rownames(airway),
    "SYMBOL", "ENSEMBL"
)

Next, we use the uniquifyFeatureNames() function of the scuttle package to replace the rownames() by a unique identifier that is generated as follows:

  • The gene symbol if it is unique.
  • A concatenate of the gene symbol and Ensembl gene identifier if the gene symbol is not unique.
  • The Ensembl identifier if the gene symbol is not available.
library("scuttle")
rownames(airway) <- uniquifyFeatureNames(
    ID = rowData(airway)[["ENSEMBL"]],
    names = rowData(airway)[["SYMBOL"]]
)
airway
#> class: RangedSummarizedExperiment 
#> dim: 63677 8 
#> metadata(1): ''
#> assays(1): counts
#> rownames(63677): TSPAN6 TNMD ... ENSG00000273492 ENSG00000273493
#> rowData names(12): gene_id gene_name ... ENSEMBL SYMBOL
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSample

3 Differential expression

To generate some example results, we first use edgeR::filterByExpr() to remove genes whose counts are too low to support a rigorous differential expression analysis. Then we run a standard Limma-Voom analysis using edgeR::voomLmFit(), limma::makeContrasts(), and limma::eBayes(); alternatively, we could have used limma::treat() instead of limma::eBayes().

The linear model includes the dex and cell covariates, indicating the treatment conditions and cell types, respectively. Here, we are interested in differences between treatments, adjusted by cell type, and define this comparison as the dextrt - dexuntrt contrast.

The final differential expression results are fetched using limma::topTable().

library("edgeR")

design <- model.matrix(~ 0 + dex + cell, data = colData(airway))

keep <- filterByExpr(airway, design)
fit <- voomLmFit(airway[keep, ], design, plot = FALSE)
contr <- makeContrasts("dextrt - dexuntrt", levels = design)
fit <- contrasts.fit(fit, contr)
fit <- eBayes(fit)
res_limma <- topTable(fit, sort.by = "P", n = Inf)
head(res_limma)
#>                 gene_id gene_name entrezid   gene_biotype gene_seq_start gene_seq_end seq_name
#> CACNB2  ENSG00000165995    CACNB2       NA protein_coding       18429606     18830798       10
#> DUSP1   ENSG00000120129     DUSP1       NA protein_coding      172195093    172198198        5
#> MAOA    ENSG00000189221      MAOA       NA protein_coding       43515467     43606068        X
#> SPARCL1 ENSG00000152583   SPARCL1       NA protein_coding       88394487     88452213        4
#> PRSS35  ENSG00000146250    PRSS35       NA protein_coding       84222194     84235423        6
#> STEAP2  ENSG00000157214    STEAP2       NA protein_coding       89796904     89867451        7
#>         seq_strand seq_coord_system  symbol         ENSEMBL  SYMBOL     logFC  AveExpr         t
#> CACNB2           1               NA  CACNB2 ENSG00000165995  CACNB2  3.205606 3.682244  36.52011
#> DUSP1           -1               NA   DUSP1 ENSG00000120129   DUSP1  2.864778 6.644455  28.94270
#> MAOA             1               NA    MAOA ENSG00000189221    MAOA  3.256085 5.950559  28.27912
#> SPARCL1         -1               NA SPARCL1 ENSG00000152583 SPARCL1  4.489009 4.166904  27.98254
#> PRSS35           1               NA  PRSS35 ENSG00000146250  PRSS35 -2.828184 3.224885 -27.49302
#> STEAP2           1               NA  STEAP2 ENSG00000157214  STEAP2  1.894559 6.790009  26.92657
#>              P.Value    adj.P.Val        B
#> CACNB2  2.463606e-11 4.153640e-07 15.90582
#> DUSP1   2.090804e-10 9.944584e-07 14.60579
#> MAOA    2.586665e-10 9.944584e-07 14.38248
#> SPARCL1 2.849300e-10 9.944584e-07 13.91713
#> PRSS35  3.349813e-10 9.944584e-07 13.71315
#> STEAP2  4.054069e-10 9.944584e-07 13.99293

Then, we embed this set of differential expression results in the airway object using the embedContrastResults() method and we use the function contrastResults() to display the contrast results stored in the airway object.

library(iSEEde)
#> Loading required package: iSEE
airway <- embedContrastResults(res_limma, airway,
    name = "dextrt - dexuntrt",
    class = "limma"
)
contrastResults(airway)
#> DataFrame with 63677 rows and 1 column
#>                  dextrt - dexuntrt
#>                 <iSEELimmaResults>
#> TSPAN6          <iSEELimmaResults>
#> TNMD            <iSEELimmaResults>
#> DPM1            <iSEELimmaResults>
#> SCYL3           <iSEELimmaResults>
#> FIRRM           <iSEELimmaResults>
#> ...                            ...
#> ENSG00000273489 <iSEELimmaResults>
#> ENSG00000273490 <iSEELimmaResults>
#> ENSG00000273491 <iSEELimmaResults>
#> ENSG00000273492 <iSEELimmaResults>
#> ENSG00000273493 <iSEELimmaResults>
contrastResults(airway, "dextrt - dexuntrt")
#> iSEELimmaResults with 63677 rows and 18 columns
#>                         gene_id   gene_name  entrezid   gene_biotype gene_seq_start gene_seq_end
#>                     <character> <character> <integer>    <character>      <integer>    <integer>
#> TSPAN6          ENSG00000000003      TSPAN6        NA protein_coding       99883667     99894988
#> TNMD                         NA          NA        NA             NA             NA           NA
#> DPM1            ENSG00000000419        DPM1        NA protein_coding       49551404     49575092
#> SCYL3           ENSG00000000457       SCYL3        NA protein_coding      169818772    169863408
#> FIRRM           ENSG00000000460    C1orf112        NA protein_coding      169631245    169823221
#> ...                         ...         ...       ...            ...            ...          ...
#> ENSG00000273489              NA          NA        NA             NA             NA           NA
#> ENSG00000273490              NA          NA        NA             NA             NA           NA
#> ENSG00000273491              NA          NA        NA             NA             NA           NA
#> ENSG00000273492              NA          NA        NA             NA             NA           NA
#> ENSG00000273493              NA          NA        NA             NA             NA           NA
#>                    seq_name seq_strand seq_coord_system      symbol         ENSEMBL      SYMBOL
#>                 <character>  <integer>        <integer> <character>     <character> <character>
#> TSPAN6                    X         -1               NA      TSPAN6 ENSG00000000003      TSPAN6
#> TNMD                     NA         NA               NA          NA              NA          NA
#> DPM1                     20         -1               NA        DPM1 ENSG00000000419        DPM1
#> SCYL3                     1         -1               NA       SCYL3 ENSG00000000457       SCYL3
#> FIRRM                     1          1               NA    C1orf112 ENSG00000000460       FIRRM
#> ...                     ...        ...              ...         ...             ...         ...
#> ENSG00000273489          NA         NA               NA          NA              NA          NA
#> ENSG00000273490          NA         NA               NA          NA              NA          NA
#> ENSG00000273491          NA         NA               NA          NA              NA          NA
#> ENSG00000273492          NA         NA               NA          NA              NA          NA
#> ENSG00000273493          NA         NA               NA          NA              NA          NA
#>                     logFC   AveExpr         t     P.Value  adj.P.Val         B
#>                 <numeric> <numeric> <numeric>   <numeric>  <numeric> <numeric>
#> TSPAN6          -0.464216   5.02559 -6.412594 0.000108256 0.00129908   1.04829
#> TNMD                   NA        NA        NA          NA         NA        NA
#> DPM1             0.125077   4.60191  1.599859 0.143110360 0.25794747  -6.30770
#> SCYL3           -0.042107   3.47269 -0.429777 0.677175139 0.78307084  -7.22549
#> FIRRM           -0.228517   1.40857 -0.983508 0.350328761 0.49314902  -6.19572
#> ...                   ...       ...       ...         ...        ...       ...
#> ENSG00000273489        NA        NA        NA          NA         NA        NA
#> ENSG00000273490        NA        NA        NA          NA         NA        NA
#> ENSG00000273491        NA        NA        NA          NA         NA        NA
#> ENSG00000273492        NA        NA        NA          NA         NA        NA
#> ENSG00000273493        NA        NA        NA          NA         NA        NA

4 Live app

In this example, we use iSEE::panelDefaults() to specify rowData() fields to show in the tooltip that is displayed when hovering a data point.

The application is then configured to display the volcano plot and MA plot for the same contrast.

Finally, the configured app is launched.

library(iSEE)

panelDefaults(
    TooltipRowData = c("SYMBOL", "ENSEMBL")
)

app <- iSEE(airway, initial = list(
    VolcanoPlot(ContrastName = "dextrt - dexuntrt", PanelWidth = 6L),
    MAPlot(ContrastName = "dextrt - dexuntrt", PanelWidth = 6L)
))

if (interactive()) {
    shiny::runApp(app)
}