Contents

1 Introduction

TEKRABber is used to estimate the correlations between genes and transposable elements (TEs) from RNA-seq data comparing between: (1) Two Species (2) Control vs. Experiment. In the following sections, we will use built-in data to demonstrate how to implement TEKRABber on you own analysis.

2 Installation

To use TEKRABber from your R environment, you need to install it using BiocManager:

install.packages("BiocManager")
BiocManager::install("TEKRABber")
library(TEKRABber)
library(SummarizedExperiment) # load it if you are running this tutorial

3 Examples

3.1 Comparing between two species, human and chimpanzee as an example

Gene and TE expression data are generated from randomly picked brain regions FASTQ files from 10 humans and 10 chimpanzees (Khrameeva E et al., Genome Research, 2020). The values for the first column of gene and TE count table must be Ensembl gene ID and TE name:

# load built-in data
data(speciesCounts)
hmGene <- speciesCounts$hmGene
hmTE <- speciesCounts$hmTE
chimpGene <- speciesCounts$chimpGene
chimpTE <- speciesCounts$chimpTE
# the first column must be Ensembl gene ID for gene, and TE name for TE
head(hmGene)
##            Geneid SRR8750453 SRR8750454 SRR8750455 SRR8750456 SRR8750457
## 1 ENSG00000000003        250        267        227        286        128
## 2 ENSG00000000005         13          2         15          9          5
## 3 ENSG00000000419        260        311        159        259        272
## 4 ENSG00000000457         86        131        100         94         80
## 5 ENSG00000000460         21         17         42         33         55
## 6 ENSG00000000938        162         75         95        252        195
##   SRR8750458 SRR8750459 SRR8750460 SRR8750461 SRR8750462
## 1        394        268        102        370        244
## 2          0          1          8          0          2
## 3        408        371        126        211        374
## 4        158        119         46         77         81
## 5         29         50         11         18         20
## 6        137         93        108        197         69

3.1.1 Query ortholog information and estimate scaling factor

In the first step, we use orthologScale() to get orthology information and calculate the scaling factor between two species for normalizing orthologous genes. The species name needs to be the abbreviation of scientific species name used in Ensembl. (Note: (1)This step queries information using biomaRt and it might need some time or try different mirrors due to the connections to Ensembl (2)It might take some time to calculate scaling factor based on your data size). For normalizing TEs, you need to provide a repeatmasker annotation table including four columns, (1) the name of TE (2) the class of TE (3) the average gene length of TE from your reference species (4) the average gene length from the species you want to compare. A way to download repeatmasker annotations is to query from UCSC Genome Table Browser and select the RepeatMasker track.

# You can use the code below to search for species name
ensembl <- biomaRt::useEnsembl(biomart = "genes")
biomaRt::listDatasets(ensembl)
# In order to save time, we have save the data for this tutorial.
data(fetchDataHmChimp)
fetchData <- fetchDataHmChimp

# Query the data and calculate scaling factor using orthologScale():
#' data(speciesCounts)
#' data(hg38_panTro6_rmsk)
#' hmGene <- speciesCounts$hmGene
#' chimpGene <- speciesCounts$chimpGene
#' hmTE <- speciesCounts$hmTE
#' chimpTE <- speciesCounts$chimpTE
#' 
#' ## For demonstration, here we only select 1000 rows to save time
#' set.seed(1234)
#' hmGeneSample <- hmGene[sample(nrow(hmGene), 1000), ]
#' chimpGeneSample <- chimpGene[sample(nrow(chimpGene), 1000), ]
#' 
#' fetchData <- orthologScale(
#'     speciesRef = "hsapiens",
#'     speciesCompare = "ptroglodytes",
#'     geneCountRef = hmGeneSample,
#'     geneCountCompare = chimpGeneSample,
#'     teCountRef = hmTE,
#'     teCountCompare = chimpTE,
#'     rmsk = hg38_panTro6_rmsk
#' )

3.1.2 Create inputs for differentially expressed analysis and correlation estimation

We use DECorrInputs() to return input files for downstream analysis.

inputBundle <- DECorrInputs(fetchData)

3.1.3 Differentially expressed analysis (DE analysis)

In this step, we need to generate a metadata contain species name (i.e., human and chimpanzee). The row names need to be same as the DE input table and the column name must be species (see the example below). Then we use DEgeneTE() to perform DE analysis. When you are comparing samples between two species, the parameter expDesign should be TRUE (as default).

meta <- data.frame(
    species = c(rep("human", ncol(hmGene) - 1), 
    rep("chimpanzee", ncol(chimpGene) - 1))
)

meta$species <- factor(meta$species, levels = c("human", "chimpanzee"))
rownames(meta) <- colnames(inputBundle$geneInputDESeq2)
hmchimpDE <- DEgeneTE(
    geneTable = inputBundle$geneInputDESeq2,
    teTable = inputBundle$teInputDESeq2,
    metadata = meta,
    expDesign = TRUE
)

3.1.4 Correlation analysis

Here we use corrOrthologTE() to perform correlation estimation comparing each ortholog and TE. This is the most time-consuming step if you have large data. For a quick demonstration, we use a relatively small data. You can specify the correlation method and adjusted p-value method. The default methods are Pearson’s correlation and FDR. Note: For more efficient and specific analysis, you can subset your data in this step to focus on only the orthologs and TEs that you are interested in.

# load built-in data
data(speciesCorr)
hmGeneCorrInput <- assay_tekcorrset(speciesCorr, "gene", "human")
hmTECorrInput <- assay_tekcorrset(speciesCorr, "te", "human")
chimpGeneCorrInput <- assay_tekcorrset(speciesCorr, "gene", "chimpanzee")
chimpTECorrInput <- assay_tekcorrset(speciesCorr, "te", "chimpanzee")

hmCorrResult <- corrOrthologTE(
    geneInput = hmGeneCorrInput,
    teInput = hmTECorrInput,
    corrMethod = "pearson",
    padjMethod = "fdr"
)

chimpCorrResult <- corrOrthologTE(
    geneInput = chimpGeneCorrInput,
    teInput = chimpTECorrInput,
    corrMethod = "pearson",
    padjMethod = "fdr"
)

head(hmCorrResult)
##          geneName      teName      pvalue        coef      padj
## 1 ENSG00000143125        L1MD 0.271964872  0.38497828 0.9990235
## 2 ENSG00000143125        MSTA 0.335873091  0.34036703 0.9990235
## 3 ENSG00000143125  MLT1N2-int 0.966658172  0.01524552 0.9990235
## 4 ENSG00000143125       LTR57 0.067870603  0.59794954 0.9990235
## 5 ENSG00000143125 HERVK11-int 0.001028058  0.87118988 0.3210294
## 6 ENSG00000143125        LTR5 0.855235258 -0.06647109 0.9990235

3.1.5 Explore your result using appTEKRABber():

TEKRABber provides an app function for you to quickly view your result. First, you will need to assign the differentially expressed orthologs/TEs results, correlation results and metadata as global variables: appDE, appRef, appCompare and appMeta. See the following example.

#create global variables for app-use
appDE <- hmchimpDE
appRef <- hmCorrResult
appCompare <- chimpCorrResult
appMeta <- meta # this is the same one in DE analysis

appTEKRABber()

In the Expression tab page, (1) you can specify your input gene and TE. The result will show in box plots with data points in normalized log2 expression level (2) DE analysis result will show in table including statistical information (3) Correlation result will indicate if these selected pairs are significantly correlated and the value of correlation coefficients. In the Correlation tab page (above figure), you can select your data in scatter plots in three ways. (1) Specify the data point and it will turn red in the distribution of your results (2) Show the distribution of all the data from your correlation results based on their correlation coefficients and adjusted p-value. The blue vertical dashed line indicates the boundary of adjusted p-value is 0.05, and the orange one is for adjusted p-value 0.01 (3) you can click the data points which you are interested in, and it will be listed in the table. You can also drag a certain area to show data points in it.

3.2 Comparing control and treatment samples within the same species

If you want to compare selected genes and TEs (1) from different tissue in same species or (2) control and drug treatment in same tissue in same species, please generate all the input files following the input format. Here we show an example data of prepared input files including expression counts from 10 control and 10 treatment samples. The format of input data: row names should be gene name or id, and column name is your sample id (please see details below).

# load built-in data
data(ctInputDE)
geneInputDE <- ctInputDE$gene
teInputDE <- ctInputDE$te

# you need to follow the input format as below
head(geneInputDE)
##                 control_1 control_2 control_3 control_4 control_5 treatment_6
## ENSG00000180263      1470      2072      1864      2238      2246        2599
## ENSG00000185985      1599      1045       946      1642      2199         665
## ENSG00000144355       517       380      1211       812        48         388
## ENSG00000234003         4         4        14        10         5           9
## ENSG00000257342         1         1         1         2         3           3
## ENSG00000223953       259       830       133       258       850         504
##                 treatment_7 treatment_8 treatment_9 treatment_10
## ENSG00000180263        2679        2562        2532         2682
## ENSG00000185985        1023        2477        1423         1731
## ENSG00000144355         275         633          59          248
## ENSG00000234003           4          18          13           22
## ENSG00000257342           0           6           1            5
## ENSG00000223953        1143        1500         498          864

3.2.1 DE analysis

For DE analysis in the same species, you also use DEgeneTE() function, however, you need to set the parameter expDesign to FALSE. You also need to provide a metadata which this time the column name must be experiment. See demonstration below:

metaExp <- data.frame(experiment = c(rep("control", 5), rep("treatment", 5)))
rownames(metaExp) <- colnames(geneInputDE)
metaExp$experiment <- factor(
    metaExp$experiment, 
    levels = c("control", "treatment")
)

resultDE <- DEgeneTE(
    geneTable = geneInputDE,
    teTable = teInputDE,
    metadata = metaExp,
    expDesign = FALSE
)

3.2.2 Correlation analysis

For a quick demonstration to perform correlation of genes and TEs in control and treatment sample, we use relatively small input tables which only include 10 genes and 10 TEs.

# load built-in data
data(ctCorr)
geneConCorrInput <- assay_tekcorrset(ctCorr, "gene", "control")
teConCorrInput <- assay_tekcorrset(ctCorr, "te", "control")
geneTreatCorrInput <- assay_tekcorrset(ctCorr, "gene", "treatment")
teTreatCorrInput <- assay_tekcorrset(ctCorr, "te", "treatment")

# you need to follow the input format as below
head(geneConCorrInput)
##                 control_1 control_2 control_3 control_4 control_5
## ENSG00000180263      1470      2072      1864      2238      2246
## ENSG00000185985      1599      1045       946      1642      2199
## ENSG00000144355       517       380      1211       812        48
## ENSG00000234003         4         4        14        10         5
## ENSG00000257342         1         1         1         2         3
## ENSG00000223953       259       830       133       258       850
controlCorr <- corrOrthologTE(
    geneInput = geneConCorrInput,
    teInput = teConCorrInput,
    corrMethod = "pearson",
    padjMethod = "fdr"
)

treatmentCorr <- corrOrthologTE(
    geneInput = geneTreatCorrInput,
    teInput = teTreatCorrInput,
    corrMethod = "pearson",
    padjMethod = "fdr"
)

head(treatmentCorr)
##          geneName teName    pvalue       coef      padj
## 1 ENSG00000180263    ALU 0.7350708 -0.2096204 0.8750842
## 2 ENSG00000180263  AluJb 0.6463988 -0.2814802 0.8729353
## 3 ENSG00000180263  AluJo 0.5879945 -0.3296607 0.8399922
## 4 ENSG00000180263  AluJr 0.6080013 -0.3130670 0.8563399
## 5 ENSG00000180263 AluJr4 0.7158070 -0.2251209 0.8729353
## 6 ENSG00000180263  AluSc 0.7130150 -0.2273721 0.8729353

3.2.3 Explore your result using appTEKRABber():

appDE <- resultDE
appRef <- controlCorr
appCompare <- treatmentCorr
appMeta <- metaExp

appTEKRABber()
sessionInfo()
## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## 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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] SummarizedExperiment_1.30.0 Biobase_2.60.0             
##  [3] GenomicRanges_1.52.0        GenomeInfoDb_1.36.0        
##  [5] IRanges_2.34.0              S4Vectors_0.38.0           
##  [7] BiocGenerics_0.46.0         MatrixGenerics_1.12.0      
##  [9] matrixStats_0.63.0          TEKRABber_1.4.0            
## [11] BiocStyle_2.28.0           
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.3            xfun_0.39               bslib_0.4.2            
##  [4] ggplot2_3.4.2           lattice_0.21-8          numDeriv_2016.8-1.1    
##  [7] vctrs_0.6.2             tools_4.3.0             bitops_1.0-7           
## [10] generics_0.1.3          parallel_4.3.0          tibble_3.2.1           
## [13] fansi_1.0.4             pkgconfig_2.0.3         Matrix_1.5-4           
## [16] lifecycle_1.0.3         GenomeInfoDbData_1.2.10 compiler_4.3.0         
## [19] munsell_0.5.0           DESeq2_1.40.0           codetools_0.2-19       
## [22] htmltools_0.5.5         sass_0.4.5              RCurl_1.98-1.12        
## [25] yaml_2.3.7              pillar_1.9.0            jquerylib_0.1.4        
## [28] MASS_7.3-59             BiocParallel_1.34.0     DelayedArray_0.26.0    
## [31] cachem_1.0.7            emdbook_1.3.12          tidyselect_1.2.0       
## [34] locfit_1.5-9.7          bdsmatrix_1.3-6         digest_0.6.31          
## [37] mvtnorm_1.1-3           dplyr_1.1.2             bookdown_0.33          
## [40] fastmap_1.1.1           grid_4.3.0              colorspace_2.1-0       
## [43] cli_3.6.1               magrittr_2.0.3          utf8_1.2.3             
## [46] scales_1.2.1            rmarkdown_2.21          XVector_0.40.0         
## [49] coda_0.19-4             evaluate_0.20           knitr_1.42             
## [52] bbmle_1.0.25            rlang_1.1.0             Rcpp_1.0.10            
## [55] glue_1.6.2              BiocManager_1.30.20     apeglm_1.22.0          
## [58] jsonlite_1.8.4          R6_2.5.1                plyr_1.8.8             
## [61] zlibbioc_1.46.0