Introduction

In this vignette, we describe usage of a suite of tools, SEESAW, Statistical Estimation of allelic Expression using Salmon and Swish.

Running SEESAW involves generation of a diploid transcriptome (e.g. using g2gtools, construction of a diploid Salmon index (specifying --keepDuplicates), followed by Salmon quantification with a number of bootstrap inferential replicates (we recommend 30 bootstrap replicates). These three steps (diploid reference preparation, indexing, quantification with bootstraps) provide the input data for the following statistical analyses in R/Bioconductor. The steps shown in this vignette leverage Bioconductor infrastructure including SummarizedExperiment for storage of input data and results, tximport for data import, and GRanges and Gviz for plotting.

In short the SEESAW steps are as listed, and diagrammed below:

  1. g2gtools (diploid reference preparation)
  2. Salmon indexing with --keepDuplicates
  3. Salmon quantification with bootstraps
  4. makeTx2Tss() aggregates data to TSS-level (optional)
  5. importAllelicCounts() creates a SummarizedExperiment
  6. Swish analysis: labelKeep() and swish() (skip scaling)
  7. Plotting

SEESAW allows for testing global allelic imbalance across all samples (pairwise testing within each individual), as well as differential or dynamic allelic imbalance (pairwise allelic fold changes estimated within individual, followed by testing across or along an additional covariate). Each of these allelic imbalance (AI) analyses takes into account the potentially heterogeneous amount of inferential uncertainty per sample, per feature (transcript, transcript-group, or gene), and per allele.

Below we demonstrate an analysis where transcripts are grouped by their transcription start site (TSS), although gene-level or transcript-level analysis is also possible. New plotting functions added to fishpond facilitate visualization of allelic and isoform changes at different resolutions, alongside gene models. In the first example, we perform global AI testing, and in the second example we perform dynamic AI testing, in both cases on simulated data associated with human genes.

Linking transcripts to TSS

We begin assuming steps 1-3 have been completed. We can use the makeTx2Tss function to generate a GRanges object t2g that connects transcripts to transcript groups.

## DataFrame with 216741 rows and 2 columns
##                           tx_id               group_id
##                     <character>            <character>
## ENST00000456328 ENST00000456328  ENSG00000223972-11869
## ENST00000450305 ENST00000450305  ENSG00000223972-12010
## ENST00000488147 ENST00000488147  ENSG00000227232-29570
## ENST00000619216 ENST00000619216  ENSG00000278267-17436
## ENST00000473358 ENST00000473358  ENSG00000243485-29554
## ...                         ...                    ...
## ENST00000420810 ENST00000420810 ENSG00000224240-2654..
## ENST00000456738 ENST00000456738 ENSG00000227629-2659..
## ENST00000435945 ENST00000435945 ENSG00000237917-2663..
## ENST00000435741 ENST00000435741 ENSG00000231514-2662..
## ENST00000431853 ENST00000431853 ENSG00000235857-5685..

Alternatively for gene-level analysis, one could either prepare a t2g data.frame with tx_id and gene_id columns, or a t2g GRanges object with a column group_id that is equal to gene_id.

Importing allelic counts

Here we will use simulated data, but we can import allelic counts with the importAllelicCounts() function. It is best to read over the manual page for this function. For TSS-level analysis, the t2g GRanges generated above should be passed to the tx2gene argument. This will summarize transcript-level counts to the TSS level, and will attach rowRanges that provide the genomic locations of the grouped transcripts.

Filtering features that have no information

Because we use --keepDuplicates in the step when we build the Salmon index, there will be a number of features in which there is no information about the allelic expression in the reads. We can find these features in bootstrap data by examining when the inferential replicates are nearly identical for the two alleles, as this is how the EM will split the reads. Removing these features avoids downstream problems during differential testing. Code for this filtering follows:

Testing for allelic imbalance across samples

We begin by generating a simulated data object that resembles what one would obtain with importAllelicCounts(). The import function arranges the a2 (non-effect) allelic counts first, followed by the a1 (effect) allelic counts. Allelic ratios are calculated at a1/a2.

## DataFrame with 20 rows and 2 columns
##          allele   sample
##        <factor> <factor>
## s1-a2        a2  sample1
## s2-a2        a2  sample2
## s3-a2        a2  sample3
## s4-a2        a2  sample4
## s5-a2        a2  sample5
## ...         ...      ...
## s6-a1        a1 sample6 
## s7-a1        a1 sample7 
## s8-a1        a1 sample8 
## s9-a1        a1 sample9 
## s10-a1       a1 sample10
## [1] "a2" "a1"

A hidden code chunk is used to add ranges from the EnsDb to the simulated dataset. For a real dataset, the ranges would be added either by importAllelicCounts (if using tx2gene) or could be added manually for transcript- or gene-level analysis, using the rowRanges<- setter function. The ranges are only needed for the plotAllelicGene plotting function below.

<hidden code chunk>

We can already plot a heatmap of allelic ratios, before performing statistical testing. We can see in the first gene, ADSS, there appear to be two groups of transcripts with opposing allelic fold change. SEESAW makes use of pheatmap for plotting a heatmap of allelic ratios.

The following two functions perform a Swish analysis, comparing the allelic counts within sample, while accounting for uncertainty in the assignment of the reads. The underlying test statistic is a Wilcoxon signed-rank statistic.

Plotting results

We can return to the heatmap, and now add q-values, etc. For details on adding metadata to a pheatmap plot object, see ?pheatmap.

In order to visualize the inferential uncertainty, we can make use of plotInfReps():

SEESAW provides plotAllelicGene() in order to build visualization of Swish test statistics, allelic proportions, and isoform proportions, in a genomic context, making use of Gviz. The first three arguments are the SummarizedExperiment object, the name of a gene (should match gene_id column), and a TxDb or EnsDb to use for plotting the gene model at the top. The statistics and proportions are then plotted at the first position of the feature (start for + features and end for - features).

You can also specify the gene using symbol:

In the allelic proportion and isoform proportion tracks, a line is drawn through the mean proportion for a2 and a1 allele, and for the isoform proportion, across samples, at the start site for each transcript group. The line is meant only to help visualize the mean value as it may change across transcript groups, but the line has no meaning in the ranges in between features. That is, unlike continuous genomic features (methylation or accessibility), there is no meaning to the allelic proportion or isoform proportion outside of measured start sites of transcription.

We can further customize the plot, for example, changing the labels displayed on the gene model, and changing the labels for the alleles. An ideogram can be added with ideogram=TRUE, although this requires connecting to an external FTP site.

See importAllelicGene() manual page for more details.

We can also customize the display of the alleles in the plotInfReps() plots, by adding a new factor, while carefully noting the existing and new allele labels, to make sure the annotation is correct:

## [1] "a2" "a1"

Testing for dynamic allelic imbalance

Above, we tested for global AI, where the allelic fold change is consistent across all samples. We can also test for differential or dynamic AI, by adding specification of a cov (covariate) which can be either a two-group factor, or a continuous variable. For continuous variable, the user should specify a correlation test, either cor="pearson" or "spearman".

## DataFrame with 20 rows and 3 columns
##          allele   sample      time
##        <factor> <factor> <numeric>
## s1-a2        a2  sample1      0.00
## s2-a2        a2  sample2      0.11
## s3-a2        a2  sample3      0.22
## s4-a2        a2  sample4      0.33
## s5-a2        a2  sample5      0.44
## ...         ...      ...       ...
## s6-a1        a1 sample6       0.56
## s7-a1        a1 sample7       0.67
## s8-a1        a1 sample8       0.78
## s9-a1        a1 sample9       0.89
## s10-a1       a1 sample10      1.00

Again, a hidden code chunk adds ranges to our simulation data.

<hidden code chunk>

In the following, we test for changes in allelic imbalance within sample that correlate with a covariate time.

Note the first two features have small q-values and opposite test statistic; here the test statistic is the average Pearson correlation of the allelic log fold change with the time variable, averaging over bootstrap replicates.

## DataFrame with 2 rows and 2 columns
##                     stat    qvalue
##                <numeric> <numeric>
## ADSS-244452134  0.870969     0.005
## ADSS-244419273 -0.861573     0.005

For plotting inferential replicates over a continuous variable, we must first compute summary statistics of inferential mean and variance:

Now we can examine the allelic counts across the time variable:

With a little more code, we can add a lowess line for each series:

Visualizing the allelic proportion in a heatmap helps to see relationships with the time variable, while also showing data from multiple features at once:

Finally, by binning the time covariate into a few groups, we can again draw the allelic and isoform proportions in the genomic context, now facetting across time.

First we create the binned covariate using cut, and rename the labels for nicer labels in our plot:

## 
## time-1 time-2 time-3 
##      3      4      3

We can then make our facetted allelic proportion plot:

If we also want to visualize how isoform proportions may be changing, we can add covFacetIsoform=TRUE, which additionally facets the isoform proportion plot by the covariate:

For further questions about the SEESAW steps, please post to one of these locations:

Session info

## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] SummarizedExperiment_1.25.3 MatrixGenerics_1.7.0       
##  [3] matrixStats_0.62.0          fishpond_2.1.45            
##  [5] EnsDb.Hsapiens.v86_2.99.0   ensembldb_2.19.10          
##  [7] AnnotationFilter_1.19.0     GenomicFeatures_1.47.14    
##  [9] AnnotationDbi_1.57.1        Biobase_2.55.2             
## [11] GenomicRanges_1.47.6        GenomeInfoDb_1.31.10       
## [13] IRanges_2.29.1              S4Vectors_0.33.17          
## [15] BiocGenerics_0.41.2        
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3            rjson_0.2.21               
##   [3] ellipsis_0.3.2              htmlTable_2.4.0            
##   [5] biovizBase_1.43.1           qvalue_2.27.0              
##   [7] XVector_0.35.0              base64enc_0.1-3            
##   [9] dichromat_2.0-0             rstudioapi_0.13            
##  [11] farver_2.1.0                bit64_4.0.5                
##  [13] fansi_1.0.3                 xml2_1.3.3                 
##  [15] splines_4.2.0               cachem_1.0.6               
##  [17] knitr_1.38                  Formula_1.2-4              
##  [19] jsonlite_1.8.0              Rsamtools_2.11.0           
##  [21] cluster_2.1.3               dbplyr_2.1.1               
##  [23] png_0.1-7                   pheatmap_1.0.12            
##  [25] compiler_4.2.0              httr_1.4.2                 
##  [27] backports_1.4.1             assertthat_0.2.1           
##  [29] Matrix_1.4-1                fastmap_1.1.0              
##  [31] lazyeval_0.2.2              cli_3.2.0                  
##  [33] htmltools_0.5.2             prettyunits_1.1.1          
##  [35] tools_4.2.0                 gtable_0.3.0               
##  [37] glue_1.6.2                  GenomeInfoDbData_1.2.8     
##  [39] reshape2_1.4.4              dplyr_1.0.8                
##  [41] rappdirs_0.3.3              Rcpp_1.0.8.3               
##  [43] jquerylib_0.1.4             vctrs_0.4.1                
##  [45] Biostrings_2.63.3           rtracklayer_1.55.4         
##  [47] xfun_0.30                   stringr_1.4.0              
##  [49] lifecycle_1.0.1             restfulr_0.0.13            
##  [51] gtools_3.9.2                XML_3.99-0.9               
##  [53] zlibbioc_1.41.0             scales_1.2.0               
##  [55] BSgenome_1.63.7             VariantAnnotation_1.41.3   
##  [57] hms_1.1.1                   ProtGenerics_1.27.2        
##  [59] parallel_4.2.0              RColorBrewer_1.1-3         
##  [61] SingleCellExperiment_1.17.2 yaml_2.3.5                 
##  [63] curl_4.3.2                  gridExtra_2.3              
##  [65] memoise_2.0.1               ggplot2_3.3.5              
##  [67] sass_0.4.1                  rpart_4.1.16               
##  [69] biomaRt_2.51.4              latticeExtra_0.6-29        
##  [71] stringi_1.7.6               RSQLite_2.2.12             
##  [73] highr_0.9                   BiocIO_1.5.0               
##  [75] checkmate_2.1.0             filelock_1.0.2             
##  [77] BiocParallel_1.29.21        rlang_1.0.2                
##  [79] pkgconfig_2.0.3             bitops_1.0-7               
##  [81] evaluate_0.15               lattice_0.20-45            
##  [83] purrr_0.3.4                 htmlwidgets_1.5.4          
##  [85] GenomicAlignments_1.31.2    bit_4.0.4                  
##  [87] tidyselect_1.1.2            plyr_1.8.7                 
##  [89] magrittr_2.0.3              R6_2.5.1                   
##  [91] generics_0.1.2              Hmisc_4.7-0                
##  [93] DelayedArray_0.21.2         DBI_1.1.2                  
##  [95] pillar_1.7.0                foreign_0.8-82             
##  [97] svMisc_1.2.3                survival_3.3-1             
##  [99] KEGGREST_1.35.0             abind_1.4-5                
## [101] RCurl_1.98-1.6              nnet_7.3-17                
## [103] tibble_3.1.6                crayon_1.5.1               
## [105] utf8_1.2.2                  BiocFileCache_2.3.4        
## [107] rmarkdown_2.13              jpeg_0.1-9                 
## [109] progress_1.2.2              grid_4.2.0                 
## [111] data.table_1.14.2           blob_1.2.3                 
## [113] digest_0.6.29               munsell_0.5.0              
## [115] Gviz_1.39.6                 bslib_0.3.1