TEffectR: an R package for studying the potential effects of transposable elements on gene expression with linear regression model

Modeling gene expression with linear regression model

RNA-seq method yields count-type data rather than continuous measures of gene expression. Hence, generalized linear models (GLM) are used for modeling and statistical analysis of RNA-seq data sets, which are assumed to follow Poisson distribution or negative binomial distribution. In order to test differential gene expression, a number of analytical methods, including edgeR (Robinson, McCarthy & Smyth, 2010), and DESeq2 (Oshlack, Robinson & Young, 2010) use GLM where expression level of each gene is modeled as response variable while biological conditions (e.g., control vs experimental groups) are considered as explanatory variables or predictors. However, after the transformation of RNA-seq count data to log2-counts per million (logCPM) with Limma’s voom (Law et al., 2014) function, gene expression profiles can be ready for linear modelling.

LM has been used so far for modeling the regulatory effects of genetic and epigenetic factors on gene expression (Gerstung et al., 2015; Li, Liang & Zhang, 2014). For example, Li et al. developed RACER (Li, Liang & Zhang, 2014), using regression analysis approach, for exploring potential links between gene expression levels and a number of predictors, including DNA methylation level, copy number variation, transcription factor occupancy and microRNA expression level. Using a similar approach, TEffectR considers given biological conditions or covariates and TE expression levels as predictors to explain significant differences in gene abundances on a gene-by-gene basis. TEffectR assumes that each TE has a potential to influence the expression of a proximal gene. Accordingly, the expression level of a given gene can be modelled as follows: ${\displaystyle Geneexpressio{n}_i={\beta }_0+{\beta }_{1}T{E}_{1i}+\cdots +{\beta }_{n}T{E}_{ni}+{\beta }_{m}Covarite{s}_{mi}+{\epsilon }_i,{\epsilon }_i\sim N\left(0,{\sigma }^2\right)}$

• where i denotes the genes, and Geneexpression_i represents the normalized log2(CPM) value of the i^th gene.

• TE_ni stands for the normalized log2(CPM) expression value of n ^th TE which is located within the upstream of the i^th gene.

• Covarites_mi indicates covariate effects in the model, such as tissue type, age, gender, etc.

Implementation of the TEffectR package

The TEffectR package was written in R language (v.3.5.3) and it uses the functions of diverse computational tools. To extract gene annotation data from Ensembl database and manipulation of RepeatMasker annotation files, TEffectR respectively utilizes biomaRt (Durinck et al., 2009) and biomartr (Drost & Paszkowski, 2017) packages. The GenomicRanges (Lawrence et al., 2013) tool is used to identify TE sequences that are located in the neighborhood of the gene list provided by users. For data and string manipulation steps of the TEffectR workflow, we made use of dplyr (https://dplyr.tidyverse.org/), rlist (https://renkun-ken.github.io/rlist/) and stringr (https://stringr.tidyverse.org/) packages. BEDtools (Quinlan, 2014) and Rsamtools (https://bioconductor.org/packages/release/bioc/html/Rsamtools.html) were employed for the quantification of TE-derived sequencing reads in a given list of indexed and genome aligned BAM files. Two popular differential gene expression analysis packages for RNA-seq data sets, edgeR and limma, were used for filtering, normalization and transformation of expression values of both genes and TEs. Statistical significance of each LM and covariate effect in the corresponding regression model were calculated with lm(), which is a built-in function of R.

Data collection and processing for case study

In order to demonstrate the usage of TEffectR package, we made use of a publicly available whole transcriptome sequencing dataset including normal and tumour tissue specimens obtained from 22 ER+/HER2- breast cancer patients (GEO Accession ID: GSE103001) (Wenric et al., 2017). These transcriptome libraries were particularly included as they were prepared without poly(A) selection method and thereby allowing the measurement of TE expression uniformly (Solovyov et al., 2018). We downloaded sequencing reads in FASTQ file format from Sequence Read Archive (Leinonen et al., 2011) (SRA Accession ID: SRP116023) using SRA Tool Kit v.2.9.0 with “fastq-dump –gzip –skip-technical –readids –dumpbase –clip –split-3” command. Next, sequencing reads were aligned to the human reference genome GRCh38 (Ensembl version 78) using the splice-aware aligner HISAT2 v2.1.0 (Kim, Langmead & Salzberg, 2015) with “hisat2 -p -dta -x {input.index} -1 {input_1.fq} -2 {input_2.fq} -S {out.sam}” parameters. Stringtie v1.3.5 (Pertea et al., 2015) were used with “stringtie -e -B -p -G {input.gtf} -A {output.tab} -o {output.gtf} -l {input.label}{input.bam}” parameters for expression quantification at the gene level. We considered only the uniquely mapped reads overlapping TE regions for the expression quantification of TEs. Multi-mapped reads could cause ambiguity when analysing the local effects of TEs on proximal genes as the repeats have many copies on the genome (Goerner-Potvin & Bourque, 2018; Treangen & Salzberg, 2011). To remove multi- and unmapped reads from BAM files, “samtools view -bq 60 -o {out.bam}{input.bam}” command was used. If the user would like to include the multi-mapped reads, they can skip this last command.

Overview of the TEffectR package pipeline

The TEffectR package includes a set of functions (Fig. 1) that allows the identification of significant associations between TEs and nearby genes for any species whose repeat annotation is publicly available at the RepeatMasker website (http://www.repeatmasker.org/genomicDatasets/RMGenomicDatasets.html). Currently, the complete annotation for over 60 species (from primates to nematodes) can be downloaded from this main repository and can be analyzed with our R package. The TEffectR package works initially by manipulating the repeat annotation file and make it ready for downstream analysis. In the following step, our tool takes a raw count matrix of RNA-seq dataset from the user where the first column includes the gene symbols or Ensembl IDs and the other columns contain count values of genes across samples. Then, TEffectR retrieves genomic position of each gene in the respective genome using the given count matrix. Afterwards, based on the user-defined parameters, TEffectR determines all TE species that are located within the upstream regions of each gene individually.

The TEffectR package contains a handy function for obtaining sequencing read counts, which are aligned to each TE region from a given list of sorted and indexed BAM files. Additionally, it can calculate the total read counts of each TE associated with a certain gene. In the following step, TEffectR merges all read count values of both genes and TEs into a single count matrix. This count matrix is then filtered, normalized with Trimmed Mean of M-values (TMM) method and transformed for linear modeling using voom() function of the limma package. In the final step, TEffectR fits a linear regression model with customized design matrix for each gene, and it calculates adjusted R-square values and significance of the model, and estimates the model parameters. The users can output the results of all calculations in tab separated values (tsv) file format to assess the contribution of each covariate (e.g., individual repeats) to the model.

Descriptions of the functions in the TEffectR package

The TEffectR package provides six unique functions for predicting the potential influence of TEs on the transcriptional activity of proximal genes in the respective genome:

TEffectR::rm_format: This function takes RepeatMasker annotation file as input and extracts the genomic location of each TE along with repeat class and family information. The output of rm_format() function is used while searching TEs that are located in the upstream region of the genes of interest.

TEffectR::get_intervals: This function is used to retrieve the genomic locations of all genes in a given read count matrix by the user. Row names of the expression matrix must be one of the following: (i) official gene symbol, (ii) Ensembl gene or (iii) transcript ID. The output of this function is utilized while determining distance between genes and TEs.

TEffectR::get_overlaps: Takes the genomic intervals of genes and TEs as input. Besides, the user also requires to provide three additional parameters: (i) the maximum distance allowed between the start sites of genes and TEs, (ii) whether genes and TEs must be located in same strand and (iii) TE family or subfamily name (e.g., SINE, LINE). This function helps to detect TEs that are localized upstream of the genes of interest. The “distance” parameter of this function could be determined by the user based on the interest of the TE localization. The user can either give a positive value, which would take the TEs localized in the upstream region of the gene; or a negative value that would correspond to the downstream of the gene. Moreover, the absolute distance is not limited; however, we used the value of “5000” to test our R package with the TEs located within the 5kb upstream regions of the genes, as previous studies confirmed that TEs located within 5kb upstream of genes provide binding sites for transcription factors and provide a possible link to nearby gene expression (Bourque et al., 2008; Nikitin et al., 2018). However, if the user aims to study long-range effects, then the distance could be set to a higher value.

TEffectR::count_repeats: This function returns a raw count matrix of the total number of reads originated from TE sequences. Only the reads exhibiting 100% overlap with given TE regions are considered and the user needs to specify individual path of each BAM file as input.

TEffectR:: summarize_repeat_counts: Takes the output of count_repeats() function as input. It is used to calculate the total number of sequencing reads derived from each TE that is located upstream of a certain gene.

TEffectR::apply_lm: This core function applies filtering (≥10 reads), TMM normalization, voom transformation and LM to the given raw count expression values, respectively. It takes four arguments: (i) raw gene counts, (ii) raw TE counts, (iii) a data frame containing user-defined covariates (e.g., tissue type, disease state), and (iv) the output of get_overlaps() function. When covariates are determined, one may include all the biological factors to see if they could explain the expression of the gene in conjunction with TE expression. However, one may as well only use the TE expression as the single predictor without the inclusion of further covariates.

The apply_lm() function returns three outputs: (i) a tsv file containing the p-value of each model, significance level of covariates and associated adjusted R squared values. The generated tab delimited file contains the list for the LM results of all genes that have at least one TE within the region of interest as given by this function. (ii) another tsv file containing log2(CPM) values of genes and TEs included in LM, and (iii) a group of diagnostic plots for each significant model (p < 0.05).

A case example of TEffectR analysis using the RNA-seq data obtained from healthy and tumor tissues of ER+/HER2- breast cancer patients

Breast cancer pathogenesis was associated with genomic instability (Kwei et al., 2010), which often presents itself with the aberrant expression of TEs (Aguilera & Garcia-Muse, 2013; Burns, 2017). TEs including LINE, SINE and LTR elements were already shown to be dysregulated in this disease (Yandım & Karakülah, 2019a; Bakshi et al., 2016; Bratthauer, Cardiff & Fanning, 1994; Johanning et al., 2017); with little or no information on the subtypes of these repeats. Also, there is a paucity of information on the impact of such dysregulatory events on the expressions of genes. To demonstrate the usage of the TEffectR package on a real case example, we ran the TEffectR package on the transcriptome RNA-seq data set of ER+/HER2- breast cancer patients. Table 1 summarizes a group of significant genes involved in breast cancer pathogenesis, diagnosis and prognosis (Abdel-Fatah et al., 2014; Chung et al., 2015; Dunning et al., 2016; Forero et al., 2016; Hammerich-Hille et al., 2010; Heinrich et al., 2010; Heo et al., 2013; Kasper et al., 2005; Kwok et al., 2015; Li et al., 2014; Storm et al., 1995; Tang et al., 2018; Tishchenko et al., 2016; Wei et al., 2011), where TEffectR presented a linear regression model that shows the associations between the expressions of these genes and that of the uniquely mapped TE sequences located within their upstream 5 kb flanking regions (Fig. 2). These genes were only given to present a contextual example and were selected based on breast cancer literature from the tab-delimited file that contains all genes with at least one TE in their upstream regions. The LM could explain the effect of these TEs on the variation in gene expression along with other covariates such as the type of tissue (i.e., healthy or tumor) or the individual patients. For example, two of the dependent variables; the LINE element “L2b” and “Tissue type” (p = 0.0013) could statistically significantly predict 70.30% of the expression of the “SLC39A6 (LIV-1)” gene whereas the covariate “patient” (p = 0.1024) could not explain the variation in this particular gene’s expression in this statistically significant model (p < 0.001). On the other hand, even though the expression of LTR5_Hs could predict 38.92% of “RMND1” expression statistically significantly (p < 0.001), neither the type of tissue (healthy or tumor; p = 0.4280) nor the individual patient (p = 0.1193) could explain the expression of this gene. From the perspective of a molecular biologist, these results may imply that SLC39A6 gene could potentially be involved in the tumorigenesis of the breast whereas this was not the case for RMND1, and the relevant repeat motif upstream of both genes could be suitable for further experimental investigation in terms of its potential to influence the expression of the proximal gene. These results may have implications on the biological roles of TEs (e.g., L2b) on breast cancer-related gene expressions (e.g., SLC39A6) and could indicate their potential roles in the carcinogenesis of the breast where the tissue type (healthy vs. tumor) p-value of the LM result is less than a significant threshold (i.e., p < 0.05).