SQuIRE reveals locus-specific regulation of interspersed repeat expression

1Transposable elements (TEs) are self-propagating mobile genetic elements. Their insertions have resulted in a complex distribution of interspersed repeats comprising almost half of the human genome. However, most TEs have lost the capacity for generating new insertions over their evolutionary history and are now fixed in the human population. Nevertheless, even elements that have lost the potential to retrotranspose can still be transcribed from their locations in the genome. TEs are significant contributors of promoters and cis-regulatory elements. Transcription of TEs has been implicated in physiological processes in development and early embryonic pluripotency. Conversely, TE expression can also be subject to transcriptional silencing. Loss of these regulatory mechanisms resulting in aberrant TE expression has been associated with cancer, neurodegenerative diseases, and infertility. However, a deeper understanding of how TE transcription impacts these biological processes has been limited by difficulties analyzing TE transcription in RNA sequencing (RNA-seq) data.

2TEs propagate using either DNA (‘transposons’) or RNA intermediates (‘retrotransposons’). Retrotransposons are further classified into Orders, namely long terminal repeats (LTR), long interspersed elements (LINEs), and short interspersed elements (SINEs). Most elements in animal genomes have accumulated nucleotide substitutions over millions of years. However, a subset remain retrotranspositionally active and generate new polymorphic insertions. The lack of unique sequence, particularly in newer TE insertions, has presented a problem for short-read RNA sequencing. Due to the repetitive nature of TEs, RNA-seq reads that originate from one locus can ambiguously align to many TEs sharing similar sequence dispersed throughout the genome. Because of these barriers, conventional RNA-seq analyses of TEs have either discarded multi-mapping alignments or combined TE expression to the subfamily level. Other groups have studied active LINE-1s using tailored pipelines, leveraging internal sequence variation and 3′ transcription extensions into unique sequence. However, these targeted approaches do not provide a global picture of expression from all classes of TEs. TE studies done at the subfamily level are unable to distinguish TE-intrinsic expression via TE-derived regulatory sequences from TE-extrinsic expression due to TE inclusion in a longer transcript. Conversely, RNA-seq pipelines that rely solely on uniquely aligning reads can miss differential expression of highly repetitive TEs. There is thus a need for locus-specific analyses of TE expression to understand their regulation in normal and disease states.

3To analyze global TE expression in conventional RNA-seq experiments, we have developed the Software for Quantifying Interspersed Repeat Expression (SQuIRE). SQuIRE provides a suite of tools to ensure the pipeline is user-friendly, reproducible, and broadly applicable. Like previous TE expression software, SQuIRE quantifies expression at the subfamily level and performs differential expression analyses on TEs and genes. Unlike past approaches however, SQuIRE quantifies TE expression at the locus level. We benchmarked this pipeline using both simulated and experimental datasets and compared its performance against other software to quantify TE expression. We applied SQuIRE to mouse tissue RNA-seq data and identified examples of locus-specific differential expression in testis compared to somatic tissues. SQuIRE enabled us to assess the physical context of expressed TEs (i.e. the location of these TEs in the genome and within potentially larger RNA transcripts). This revealed many examples of extrinsic regulation of TE expression as part of long transcripts, which subfamily-level analyses would otherwise miss. We also identify specific differentially expressed endogenous retroviral Gypsy loci in a Drosophila model of amyotrophic lateral sclerosis (ALS). Our findings confirm that locus-specific analysis, attainable with SQuIRE, is essential to get a true picture of the TE transcriptome.

4SQuIRE was written in Python 2 and tested with the following specific versions of software: STAR 2.5.3a, BEDtools 2.25.0, SAMtools 1.8, StringTie 1.3.3b, DESeq2 1.16.1, R 3.4.1 and Python 2.7.9. SQuIRE was developed for UNIX environments. Briefly, the SQuIRE pipeline includes Fetch to obtain reference annotation files, Map to align RNA-seq data, Count to quantify gene and TE expression, and Call to perform differential analysis. The algorithm for quantifying TE expression is exclusive to SQuIRE and described below. Details of the software parameters implemented in the SQuIRE pipeline are described in Supplementary Methods. We provide step-by-step instructions on our README to use the package manager Conda (conda.io) to download the correct versions of prerequisite software for SQuIRE (e.g. Python, R, STAR, BEDTools, StringTie, SAMtools, DESeq2). The README also instructs users how to create a non-reference table with the exogenous or polymorphic TE sequences and coordinates that they would like to add to the reference genome. Bash scripts to run each tool in the SQuIRE pipeline are also available on the website. Users can fill in crucial experiment information (raw data, read length, paired, strandedness, genome build, sample name and experimental design) into the ‘arguments.sh’ file, which the other scripts reference to run each step with the correct parameters.

5To quantify TE expression, Count first identifies reads that map to TEs. If a TE-mapping read aligns to a single locus after a genome-wide scan, it is labeled as a ‘unique read’; if the read maps to multiple locations, it is labeled as a ‘multi-mapped read’. Count allows for 50% of the read to map to flanking sequence to increase the detection of uniquely aligning reads. For paired-end reads, each individual end is first assessed for unique alignment before identifying their mates. If one multi-mapping end is paired with a uniquely aligning mate, the pair is considered ‘unique’ and other alignments of the multi-mapping mate are discarded. If the RNA-seq data is stranded, the sense and anti-sense direction of a TE are treated as separate transcripts to which a read can align. Second, Count assigns fractions of a read to each TE as a function of the probability that the TE gave rise to that read. Uniquely aligning reads are considered certain (i.e. probability = 100%, count = 1). Count initially assigns fractions of multi-mapping reads to TEs in proportion to their relative expression as indicated by unique read alignments. In doing so, Count also considers that TEs have varying uniquely alignable sequence lengths. To mitigate bias against the n number of TEs without uniquely aligning reads, these TEs receive fractions inversely proportional to the number of loci (N) to which each read aligned. Then Count assigns the remainder to the TEs with unique reads. If both mates of a read pair are multi-mapping, but only map concordantly to a single TE location, the discordant alignments are discarded. The read pair contributes a full read count to the TE, but they are not considered ‘unique’ and their positions do not contribute to the TE’s uniquely alignable length. To account for TEs that have fewer unique counts due to having less unique sequence, Count normalizes each unique count to the number of individual unique read start positions, or each TE’s uniquely alignable length. Among all TEs to which a multi-mapping read aligned, the TEs with unique reads ( are compared with each other. A fraction of a read is assigned to each TE in proportion to the contribution of the normalized unique count ( to the combined normalized unique count of all of the TEs being compared ( (Equation 1). Thus, the sum of unique counts and multi-mapped read fractions for each TE provides an initial estimate of TE read abundance based on empirically obtained unique read counts and uniquely alignable sequence.

6At this point, multi-mapping reads are assigned to TEs with no unique reads based only on the numbers of valid alignments for each read. This can result in over- or under-estimations of TE expression. To combat this issue, Count next refines this initial assignment by redistributing multi-mapping read fractions in proportion to estimated TE expression with an expectation-maximization algorithm. To estimate expression, Count uses the a TE’s total read count ( = unique read counts + multi-mapped fractions from the previous step) normalized by the effective transcript length: The effective transcript length is calculated as the estimated transcript length subtracted by the average fragment length aligned to that TE + 1,, as described previously. All of the TEs to which a multi-mapping read aligned ( are compared with each other. A fraction of a read is assigned to each TE in proportion to the relative normalized total count compared to the combined normalized total count of all of the TEs being compared (, as shown in Equation (2). Count assumes this value is proportional to the probability that the TE gave rise to the multi-mapping read, and assigns that fraction of a read count to the TE. Because TEs with a count fraction of less than 1 have a low probability of giving rise to any read, those TEs are assigned a count fraction of 0.

7After the total counts (unique and multi-mapped) of each TE are re-calculated, multi-mapped reads can be re-assigned in subsequent iterations of expectation (assigning multi-mapped read fractions to TEs) and maximization (summation of unique and multi-mapped fraction counts). These iterations can be repeated until a given iteration number set by the user or until the TE counts converge (‘auto’, when all of the TEs with ≥10 counts change by <1%).

8TEs with few uniquely aligning reads may be prone to misrepresentation. Users who want to identify TEs that are more likely to be false positives can examine a ‘score’ value provided in the Count output. The score is defined as × 100, where represents the number of all reads aligned to the locus (unique and multi-mapping), and represents the final read count from SQuIRE. A low score indicates that relatively more reads assigned to the TE are potentially derived from other loci, while a high score conveys greater certainty in the read count.

9An example of Count output is provided in Supplementary Table S1. Further details of the Count algorithm are in Supplemental Methods.

10To evaluate SQuIRE with known TE expression levels, we tested SQuIRE with simulated RNA-seq data. We randomly selected 100 000 TEs from the GRCh38/hg38 (hg38) Repeatmasker annotation downloaded by Fetch. We limited our list of potential TEs to those included in TEtranscripts and RepEnrich to enable comparisons between these different programs. Using the selected TE coordinates we generated a BED file using Clean and obtained FASTA sequences using Seek. To mimic intrinsically regulated expression which can be more difficult to detect, we did not include flanking sequence in the TE coordinates for simulation. From these TE sequences, we used the Polyester package from Bioconductor (R version 3.4.1;) to simulate 100 bp, paired-end, stranded RNA-seq reads with normally distributed fragment lengths around a mean of 250 bp. We simulated a uniformly distributed sequencing error rate of 0.5%. TEs were simulated with a mean read coverage of 20×, with 250 TEs deviating from that mean between 2- and 100-fold.

11To evaluate SQuIRE with induced TE expression, we transfected a LINE-1 (L1, L1RP) expressing plasmid into a cell line and used SQuIRE to evaluate L1 expression. LINE expression constructs were cloned into the pCEP4 backbone (Thermo Fisher Scientific, Waltham, MA) modified to confer puromycin resistance. Plasmids encoded either L1RP (MT302) or had no insert. Tet-On HEK293TLD (293T) cells were grown at 37°C, 5% CO2 in DMEM with 10% Tet-Free FBS (Takara, Mountain View, CA) and passaged every 3–5 days as needed with regular tests for mycoplasma contamination. For transfection, 300 000 293T cells were plated in 2 ml volume. 24 h later, cells were transfected using a cocktail of 2 μg plasmid DNA and 6 μl Fugene HD (Promega), and puromycin was added 24 h later for a total of 3 days of selection. 500 000 cells were then plated in three wells each, and doxycycline was added 2 h later (final concentration of 1 ug/ml) to induce L1 expression. RNA was collected after 72 h of L1 expression using the Zymo Quick-RNA MiniPrep kit (Zymo Research, Tustin, CA, USA). The RNA libraries of transfected 293T cells were prepared using the Illumina TruSeq Stranded Total Library Prep Kit with Ribo-Zero Gold (San Diego, CA, USA) to provide stranded, ribosomal RNA depleted RNA. The libraries were sequenced on an Illumina HiSeq 2500, using six samples per lane across eight lanes with paired-end 100 bp reads. We generated a mean of 263 127 067 paired reads per sample. The raw sequencing data were deposited to the NCBI Genome Expression Omnibus (GEO) with accession number GSE113960.

12We ran SQuIRE on HEK293T cells transfected with an L1RP expression construct (DA1) and an empty vector (DA5). To incorporate the L1RP vector into the alignment, the vector sequence was added to a custom table and referenced using the ‘—extra’ option in Map. We used SAMtools to identify reads that align to the construct (Supplementary Table S2). To test the effect of ectopic L1RP expression on the false positive rate of endogenous L1 expression estimation, we took L1RP-aligning reads from sample DA1 for in silico ‘spike-in’ to sample DA5. To downsample these L1RP-aligning reads, we used the SAMtools ‘-s <INT.FRAC> ’ option. We used values of 0.01, 1.001 and 3.0004 as inputs, resulting in expression levels of 0.43, 0.78 and 7.90 fpkm. [The integer before the decimal indicates the seed value and the number after the decimal indicates the fraction of total alignments desired for subsampling.] We used the SAMtools ‘merge’ tool to combine these L1RP-aligning reads with one lane equivalent (29.8 million reads) from the empty vector (DA5) sample.

13Mouse RNA-seq data were obtained from GEO with accession number GSE30352. This study included biological replicates of brain, heart, kidney, liver and testis tissue from adult C57BL/6 mice. The RNA-seq data was paired-end, unstranded, with 76 bp length reads.

14Because SQuIRE is the first to quantify TE expression at the locus level, we restricted comparisons of SQuIRE’s performance with other TE analysis software to subfamily level analyses. All pipelines were run on a server with a maximum of 128 GB memory available and 8 threads (-p setting). For SQuIRE, we used the Fetch tool to obtain hg38 and GRCm38/mm10 (mm10, based on the C56BL/6 strain) genome FASTA sequences and RepeatMasker annotation from UCSC and generate a STAR alignment index. We ran Map, Count and Call with default settings, specifying the –build, –read_length and –strandedness parameters for the simulated human and mouse datasets. For RepEnrich, we obtained the hg38 annotation for RepeatMasker from the RepEnrich GitHub website and mm10 RepeatMasker annotation from the RepeatMasker website. We mapped the RNA-seq data using Bowtie 1 according to RepEnrich's instructions. The alignments were then used for the RepEnrich software with the ‘–pairedend TRUE’ parameter for simulated human data, and ‘–pairedend FALSE’ for mouse data. For TETools, we generated rosette files for hg38 and mm10 by taking the Repeatmasker annotation from Clean for the first column and the repeat taxonomy for the second column (subfamily:family:superfamily). We used the BED file from Clean with Seek to obtain TE FASTA sequences for generation of a pseudogenome for TETools. TETools was run with the ‘-bowtie2’, ‘–RNApair’ and ‘–insert 250’ parameters for simulated human data and ‘-bowtie2’, ‘-insert 76’ for mouse data. For TEtranscripts, we obtained hg38 and mm10 GTF annotation from the TEtranscripts website. We aligned the data to the genome with STAR using ‘–winAnchorMultimapNmax 100’, ‘–outFilterMultimapNmax 100’ parameters for multi-mapping. We then ran TEtranscripts with the ‘–mode multi’ setting to utilize its expectation-maximization algorithm for assigning multi-reads for the resulting SAM file. Since TEtranscripts analyzes TE and gene expression together, we used refGene annotation obtained by SQuIRE Fetch for the required GTF file. We used the parameters ‘–format SAM’, ‘–mode multi’, ‘–stranded yes’ for simulated human data, and ‘–format SAM’, ‘–mode multi’, ‘–stranded no’ for mouse data.

15To compare aligners used by TE analysis tools on their ability to correctly identify uniquely mapping reads, we ran the aligners Bowtie1, Bowtie2, and STAR on the simulated TE RNA-seq data described above. We set each aligner to output a maximum of two valid alignments to quickly identify uniquely aligning reads with the parameter ‘-m2’ for Bowtie 1, ‘-k2’ for Bowtie 2 and ‘–outSAMmultNmax 2’ for STAR. We also ran STAR with the parameters ‘–outFilterScoreMinOverLread 0.4 –outFilterMatchNminOverLread 0.4 –chimSegmentMin 100’ to allow for discordant alignments, which STAR excludes by default. Bowtie2 reports discordant alignments by default, while Bowtie 1 can only report paired alignments. We used BEDTools to intersect the BAM outputs to RepeatMasker annotation to identify the TEs to which the aligners mapped the reads. Reads that only appeared once were labeled as ‘uniquely aligning’. We assessed whether the mapped TE matched the templating TE for the simulated read to determine if the uniquely aligning reads mapped to the correct location.

16Drosophila RNA-seq data were obtained from GEO (accession number GSE85398). This study included biological replicates of transgenic Drosophila expressing human TDP-43 protein (hTDP43). Expression of hTDP43 was activated by a pan-neuronal enhancer (Elav) or a pan-glial enhancer (Repo). RNA-seq was performed on paired-end, unstranded, total RNA libraries with 101 bp long reads. We ran SQuIRE on the resultant data using the dm6 genome assembly, with default settings specifying the –build, –read_length and –strandedness. We performed differential expression using Call comparing hTDP43-expressing Drosophila lines (Elav/TDP43 and Repo/TDP43) with control (TDP43/2U).

17The study that had previously analyzed this dataset used TEtranscripts to analyze TE RNA expression and further aggregated their subfamily level findings to the family level. They confirmed findings at the subfamily level by qRT-PCR, protein immunolabeling, and RNAi silencing studies based on the Gypsy consensus sequence (GenBank: M12927.1). The RepeatMasker track includes >9500 loci and 45 different subfamilies belonging to the Gypsy family. To corroborate past findings of Gypsy expression in this dataset, we focused our analysis on the subset of Gypsy loci that correspond with the Gypsy consensus sequence used in previous studies. We used BLAT to determine all consensus sequence-aligning loci belonging to ‘Gypsy_I’ subfamily annotations in the RepeatMasker track. Because the subfamily-level analyses on hTDP43-mediated Gypsy upregulation focused on Gypsy coding sequences, we used SQuIRE’s output for sequences corresponding to Gypsy-I’s internal sequence (‘Gypsy_I-int’), and excluded long terminal repeat entries which do not contain open reading frames (‘Gypsy_I-LTR’).

18Differential expression analysis of gene and TE expression was performed using DESeq2 via the SQuIRE Call tool (see Supplemental Methods). P-values were adjusted for multiple-comparisons with an FDR cutoff of 0.1. To determine if loci belonging to a TE subfamily was more likely to be differentially expressed in testis compared to other TE subfamily loci, a Fisher's exact test was performed. The Fisher's exact test was chosen due to the small percentage of TE loci that are expressed.

19SQuIRE provides a suite of tools for analyzing transposable element (TE) expression in RNA-seq data (Figure 1). SQuIRE’s tools can be organized into four stages: (i) Preparation, (ii) Quantification, (iii) Analysis and (iv) Follow-up. In the Preparation stage, Fetch downloads requisite annotation files for any species with assembled genomes available on University of California Santa Cruz (UCSC) Genome Browser. These annotation files include RefSeq gene information in BED and GTF format, and RepeatMasker TE information in a custom format. Fetch also creates an index for the aligner STAR from chromosome FASTA files. Clean reformats TE annotation information from RepeatMasker into a BED file for downstream analyses. The tools in the Preparation stage only need to be run once per genome build. The Quantification stage includes the alignment step Map and RNA-seq quantification step Count. Map aligns RNA-seq data using the STAR aligner with parameters tailored to TEs that allow for multi-mapping reads and discordant alignments. It produces a BAM file. Count quantifies TE expression using a SQuIRE-specific algorithm that incorporates both unique and multi-mapping reads. It outputs read counts and fragments per kilobase transcript per million reads (fpkm) for each TE locus, and aggregates TE counts and fpkm for TE subfamilies into a separate file. Count also quantifies annotated RefSeq gene expression with the transcript assembler StringTie to output annotated gene expression as fpkm in a GTF file, and as counts in a count table file. In the Analysis stage, Call performs differential expression analysis for TEs and RefSeq genes with the Bioconductor package DESeq2. To allow users to visualize alignments to TEs of interest visualized by the Integrative Genomics Viewer (IGV) or UCSC Genome Browser, the Follow-up stage tool Draw creates bedgraphs for each sample. Seek retrieves sequences for genomic coordinates supplied by the user in FASTA format. We describe further details of the SQuIRE pipeline in Supplemental Methods.

20SQuIRE’s Count algorithm addresses a fundamental issue with quantifying reads mapping to TEs: shared sequence identity between TEs from the same subfamily and even superfamily. When a read fragment originating from these non-unique regions is aligned back to the genome, the read may ambiguously map to multiple loci (‘multi-mapped reads’). This is not a major problem for older elements that have acquired relatively many nucleotide substitutions, and thus give rise to primarily uniquely aligning reads (‘unique reads’). However, TEs from recent genomic insertions that have high sequence similarity to other loci may have few distinguishing nucleotides. Among elements of approximately the same age, relatively shorter TEs also have fewer sequences unique to a locus. Thus, discarding or misattributing multi-mapped reads can result in underestimation of TE expression.

21Previous TE RNA-seq analysis pipelines have been able to quantify TE expression at subfamily-level resolution. The software RepEnrich ‘rescued’ multi-mapping reads by re-aligning them to pseudogenome assemblies of TE loci and assigning a fraction of a read inversely proportional to the number of subfamilies to which each read aligned. These multi-mapped fractions were combined with counts of unique reads aligned to each subfamily. This approach was an advance in that it used information from multi-mapped reads. However, this method results in assigning fractions that are proportional to the number of subfamilies that share the multi-mapped read's sequence, rather than each subfamily's approximate expression level. TEtranscripts expanded on this rescue method by assigning an initial fractional value inversely proportional to the number of TE loci (not subfamilies) to which each read aligned. This initial fractional value was then used in an expectation-maximization (EM) algorithm, which iteratively re-distributes fractions of a multi-mapping read among loci (E-step) in proportion to their relative multi-mapped read abundance estimated from a previous step (M-step). The total of multi-mapped reads and unique reads for each loci are then summed by subfamily. However, in excluding unique reads from the EM algorithm, TEtranscripts does not incorporate empirical high-confidence data to infer TE expression levels from unique TE alignments. Furthermore, in calculating the relative expression level of multi-mapped reads, TEtranscripts normalizes read counts based on annotated coordinates from RepeatMasker. This underestimates TE expression levels for transcripts shorter than the annotated genomic length. TEtranscripts then sums the unique and multi-mapping counts for each subfamily.

22In order to accurately quantify TE RNA expression at locus resolution, Count builds on these previous methods by leveraging unique read alignments to each TE to assign fractions of multi-mapping reads (Figure 2) and then iteratively improving those assignments. First, Count distinguishes reads that uniquely map to particular TE loci (‘unique reads’) from reads that ambiguously map to multiple locations (‘multi-mapped reads’). Second, the count of a multi-mapped read is divided into fractions allocated to different TE loci in proportion to each TE’s unique read count, normalized to uniquely alignable length. Third, the unique reads and multi-mapped read fractions are summed and normalized to each TE’s transcribed length. Finally, the normalized total read counts are used in an expectation-maximization (EM) loop to reallocate multi-mapped read fractions. SQuIRE Count is the only TE RNA-seq analysis tool to output the length and strandedness of each TE transcript based on aligned read positions. The outcomes of using only unique reads, ignoring unique reads in multi-mapping read assignment, or relying on annotated TE lengths are illustrated in Supplementary Figure S1. In using the empirically derived uniquely alignable length and transcribed length to normalize TE counts, SQuIRE improves multi-mapping read assignment to allow TE RNA quantification at the locus level.

23To test the performance of SQuIRE Count, we simulated RNA-seq data from 100 000 randomly selected TEs from the human GRCh38/hg38 (hg38) RepeatMasker annotation. To mimic intrinsically regulated TE expression with a wide range of expression levels, TEs were simulated with read coverages ranging from 2× to 4000× and simulated counts ranging from 2 to 4588. We first evaluated accuracy by how closely SQuIRE Count output corresponded to the simulated read counts (i.e. % Observed/Expected). However, using this calculation is not meaningful for TEs with low simulated counts: a TE with 0 counts gives an infinite value, and a reported count of 1 for a TE with two simulated reads gives a low 50% Observed/Expected. Thus, we were primarily interested in ‘expressed’ simulated TEs, considering only the 99 567 TEs with at least 10 simulated reads. Second, we evaluated SQuIRE by how often it correctly detected simulated TE expression (i.e. true positives) or misreported unexpressed TEs (i.e. false positives).

24To test how well SQuIRE performed leveraging only uniquely aligning read information, we first evaluated the % Observed/Expected of TE counts with 0 EM iterations. We found that SQuIRE accurately assigned read counts to most TEs, with a mean % Observed/Expected of 98.79% (Supplementary Figure S2). A subset of evolutionarily young subfamilies from the LINE-1 superfamily (i.e. L1PA1 or L1HS), the SINE Alu superfamily (e.g. AluYa5, AluYa8, AluYb8, AluYb9), as well as composite SVA (SINE-variable number tandem repeat (VNTR)-Alu) elements remain retrotranspositionally active in the human genome and generate new polymorphic insertions. These new insertions share sequence from their parent copy. We expected that SQuIRE’s accuracy would be lower for younger TEs with less uniquely alignable sequence. Indeed, SQuIRE was less accurate for elements with less than 10% divergence (mean of 77.35% Observed/Expected). The most frequently retrotranspositionally active TEs (i.e. AluYa5, AluYa8, AluYb8, AluYb9 and L1HS) had counts ranging from 48% to 70% Observed/Expected, with a range of 79–92% Observed/Expected at the subfamily level (Supplementary Table S3). This illustrates that even without the EM-algorithm, SQuIRE can distinguish expression from highly homologous TEs at the subfamily level.

25Given the low recovery of simulated counts for younger elements when relying solely on uniquely aligning reads, we next evaluated how much adding the EM-algorithm improved Count's performance. We anticipated that the counts for most TEs would not change, but that younger elements with less divergence would have improved recovery of simulated reads. Indeed, the overall % Observed/Expected counts of TE loci increased only slightly by 0.14% to a total of 98.93%. However, the change in % Observed/Expected of TEs was much greater for the most homologous active elements, improving by 20.47% for young Alu elements and by 21.1% for L1HS loci (Figure 3). At the subfamily level, the % Observed/Expected of active TEs was improved by 8.1% for young Alu elements and by 2.2% for L1HS (Supplementary Table S3). Using updated transcript information in the EM-algorithm is thus particularly useful for TE biologists interested in younger elements that have previously been problematic to quantify by RNA-seq.

26We also wanted to evaluate SQuIRE’s ability to distinguish whether a TE is expressed or not expressed. To examine how well Count detected expressed TEs, we calculated the true positive rate (TPR) as the percentage of TEs with at least 10 simulated reads that SQuIRE also reported to have ≥10 counts. Conversely, we evaluated how often SQuIRE falsely reports TE expression by calculating the positive predictive value (PPV) as the percentage of TEs with ≥10 reported counts that were in fact simulated to have ≥10 reads. The true negative rate, or how often SQuIRE correctly reports that a TE is not expressed, is less informative for evaluating TE estimation accuracy because the number of TEs in the hg38 genome is so high (>4 million TEs) that the true negative value would outweigh the false positive value. Overall, SQuIRE had both a high TPR of 98.5% and high PPV of 99.4%. These values were lower for frequently retrotranspositionally active Alu elements (TPR = 68.75–83.33%, PPV = 64.29–100%) and L1HS elements (TPR = 100%, PPV = 62.86%) using only unique reads for TE expression estimation (Supplementary Table S4). However, using the EM algorithm improved the TPR for Alu loci (TPR = 85.22–100%) by reducing false negative reports, and improved the PPV for L1HS loci (PPV = 78.57%) by reducing false positives. The inclusion of false positives in analysis can be further reduced by imposing a score threshold. A low score indicates that multi-mapping reads contribute significantly to the read count. When we plotted the TPR and PPV using various score thresholds, we found that using a score threshold of at least 50% maximized the combination of TPR and PPV for TEs in the hg38 genome build (Supplementary Figure S3).

27To evaluate how Count handles expression from young TEs, we transfected HEK293T cells with a plasmid containing an L1HS known as L1RP. Like endogenous TEs, RNA expression from the L1RP plasmid includes unique 5′ and 3′ sequence flanking the L1HS sequence. SQuIRE readily detected the ectopic transcript, which was 686-fold more highly expressed than the highest L1HS locus in control cells (301.97 fpkm versus 0.44 fpkm).

28To evaluate whether Count would favor the L1RP over endogenous loci due to its unique flanking sequence, we ‘spiked in’ L1RP plasmid-aligning reads to the RNA-seq alignment files of HEK293T cells transfected with empty plasmid. We used randomly downsampled reads at levels approximating 1×, 2× and 20× the highest expressing endogenous locus level. We then looked at the read counts of endogenous L1HS loci before and after ‘spike-in’. Before ‘spike-in’, 22 L1HS loci were detected with >10 counts, 5 of which initiated transcription at the L1HS promoter (Supplementary Figure S4). At all simulated L1RP expression levels, there were no L1HS loci with decreased read counts after L1RP spike in. This suggests that Count appropriately normalizes for each transcript's uniquely alignable sequence.

29Conversely, because L1RP has 99.9% sequence identity to the consensus sequence of all L1HS copies, we wanted to assess if ‘spiking-in’ multi-mapping L1RP reads would result in misattributed reads to low-expressing L1 loci. Of the L1RP-aligning reads that were spiked in, only 46–50% contributed to the total read count of the L1RP locus. To assess if the remaining reads affect estimates of expression at other L1 loci, we calculated the number of false positive L1 loci that became ‘expressed’ with >10 counts after the in silico ‘spike-in’ and how this affected the PPV. We focused on the three youngest L1 subfamilies that share the greatest homology with the L1RP sequence (i.e. L1HS or L1PA1, L1PA2 and L1PA3) and compared their false positive rates to older L1 loci (Figure 4). We found that ‘spiking-in’ L1RP-derived reads only resulted in 1 false positive L1HS locus for a PPV of 95.6%, while the older subfamilies had PPVs of 99.6–100% (Figure 4). The PPVs did not change with increasing L1RP expression. Thus, there is negligible misattribution of reads.

30Because some L1HS loci remain retrotranspositionally active, they can generate insertions that are polymorphic or novel compared to the reference human RepeatMasker annotation. To assess how expression from a novel L1HS locus can impact the quantitation of reference L1 loci, we removed the L1RP annotation from the Map alignment and re-ran Count on the ‘spiked-in’ data. We found that Count was still able to detect 21 out of 22 L1HS loci (95.4%). There was an increase in false positively reported loci, decreasing the PPV from 95.4% to 84.0% with increasing L1RP expression from 1× to 20× the highest expressing endogenous L1HS locus (Figure 4). The PPV remained high for the older L1 subfamilies, ranging from 97.5% to 100%. When we set a score threshold of >50, the PPV of L1HS returned to 95.2% for all ‘spike-in’ levels, with only 1 false positive L1HS locus reported. Using a score of >50 did not drastically reduce the detection of expressed L1HS loci, with 20 of 22 (90.9%) still meeting the threshold. Thus, while estimated expression of young L1HS elements may be affected by transcription from polymorphic insertions, accuracy can be improved by adding TE annotation or using a score threshold.

31Currently published TE analysis software include RepEnrich, TEtranscripts and TETools. TEtranscripts has previously illustrated the improvements of using TE-targeted software for quantifying TE expression compared to conventional pipelines. Because none of these programs is capable of reporting TE locus expression, we performed comparisons with SQuIRE with aggregated subfamily estimates. We used the simulated hg38 TE data described above to compare the recovery of simulated reads to the correct subfamily among TE quantification software (i.e., % Observed/Expected). For mapping, we ran each software's recommended aligner: STAR (used by SQuIRE and TEtranscripts), Bowtie 2 (used by TETools), and Bowtie 1 (used by RepEnrich). We found that SQuIRE (99.86 ± 1.46%), TETools (100.14 ± 2.21%), and TEtranscripts (95.89 ± 16.41%) had comparable % Observed/Expected rates (Supplementary Figure S5). In contrast, RepEnrich (108.77 ± 40.67%) was less accurate in terms of % Observed/Expected. This is likely attributable to RepEnrich's recommended settings for Bowtie 1, which discards discordant reads and limits the number of attempts to align both paired-end mates to repetitive regions. To support this, we compared how often each aligner mapped a uniquely aligning simulated read to the correct location. We indeed found that Bowtie 1 failed to report unique reads more often in a paired-end library compared to single-end (Supplementary Table S5).

32To compare SQuIRE to other TE analysis tools with biological data, we ran each pipeline on publicly available adult C57Bl/6 mouse tissue RNA-seq data using GRCm38/mm10 (mm10) TE annotation. We compared the expression of subfamilies in testis compared to pooled data from brain, heart, kidney and liver tissues. To independently evaluate the fold-changes of TE RNA between testis and somatic tissues, we also used our previously published adult C57Bl/6 mouse Nanostring results. Unlike RNA-seq analysis, which infers transcript levels by counting reads, Nanostring uses uniquely mapping probes to capture and count RNA molecules. It thus provides an orthogonal, alignment-independent approach with which to compare TE RNA-software pipelines. We compared the Nanostring log2 fold changes (log2FC) of TE subfamily expression in testis and pooled somatic tissue to the log2FC values found by SQuIRE, RepEnrich, TEtranscripts, and TETools (Supplementary Figure S6). Because the Nanostring probes were designed against TE consensus sequences, we do not expect exact correspondence with the RNA-seq analysis tools. We observe elements for which all TE RNA-seq tools report results opposing the Nanostring result (MMVL30, IAPLTR1a_Mm, RLTR13A1). Thus in addition to comparing each pipeline with Nanostring, we also evaluated when a result deviated from the other TE RNA-seq analysis pipelines. RepEnrich failed to detect differential expression of the L1_mus_musculus subfamily (L1_Mm), and reported a direction of log2FC for the MMETn subfamily that opposed Nanostring results. TEtranscripts similarly failed to detect differential expression of MMERVK10D3 subfamily that Nanostring and the other pipelines reported, and reported different log2FC from Nanostring, SQuIRE, and TETools for L1Mm. TETools deviated from Nanostring and the other RNA-seq pipelines for the MERVL subfamily, reporting decreased expression in testis while the other methods reported upregulation. Thus, SQuIRE is the only RNA-seq pipeline producing results that corresponded with at least two other methods.

33With SQuIRE, we can closely examine the mouse RNA-seq data at the locus level. For the 16 subfamilies analyzed by Nanostring and the TE analysis tools, using SQuIRE we found that the reported subfamily-level expression was due to expression from fewer than 7% of each subfamily's loci (Supplementary Figure S7). While most subfamilies studied by Nanostring have only 1–4 significantly differentially expressed loci (log2FC > 1, padj < 0.05), the IAPLTR3 subfamily has 11 loci that are all differentially expressed in testis compared to somatic tissues (Figure 5A). To test whether this was an enrichment relative to the representation of IAPLTR3 in the mouse genome, we performed a Fisher's exact test and found that IAPLTR3 loci were 10-fold more likely than expected to be differentially expressed in testis (OR: 10.56, 95% CI: 5.25–18.97, P-value < 1.61e–08). ERVB4-1B, another LTR retrotransposon that exhibited high fold change by Nanostring, was not similarly enriched among differentially expressed TE loci. In addition to a more careful analysis of which loci are transcribed, SQuIRE enables a closer look at TE transcript structure. In examining the TE loci with the greatest differential expression in testis, we found that the transcription of the ERVB4-1B locus on chr13 did not extend beyond annotations for that element (Figure 5B), suggesting intrinsically regulated expression. On the other hand, the IAPLTR3 loci on chr14 (Figure 5C) and chr18 are part of longer transcripts that initiate outside of the annotated TE. Altogether, this suggests while a subset of TEs may be regulated by shared TE sequence, most differential expression of TEs is locus-specific with varying transcript structures, a finding that was not evident until analysis at the locus level using SQuIRE.

34To further investigate the interplay between genomic context and TE subfamily, we identified the closest genes to differentially expressed TE loci. We found a cluster of three loci exhibiting broad expression across somatic tissues from the IAP1, MERVL, and MURVY LTR retrotransposon subfamilies. When we examined the genomic context of these 3 loci, we found that all were located within genes with known broad tissue expression (Gpbp1, Csnk2a1, Kyat1, respectively), with examples shown in Supplementary Figure S8. Another locus from the MURVY subfamily is in a cluster of TEs exhibiting high testis-restricted expression. In examining the transcript overlapping the MURVY locus, we see that the transcript initiates outside of the locus and find that the transcript is an alternative splicing isoform with splice donors from the third and fourth exons of a gene ∼5 kb away (Figure 5D). The gene, Gm11981, is a long noncoding RNA (lncRNA) known to exhibit testis-restricted expression. The different MURVY-containing transcript types illustrate how TE transcription can vary across loci from the same subfamily. Altogether, these findings would be lost without the use of SQuIRE to analyze TE transcription at the locus level.

35To further illustrate the importance of locus-level analysis of TE expression with SQuIRE, we applied the SQuIRE pipeline to a Drosophila melanogaster model of amyotrophic lateral sclerosis (ALS). Almost all (97%) of ALS patients develop cytoplasmic inclusions of the TDP-43 protein. TDP-43 is a DNA- and RNA-binding protein that is involved in the regulation of RNA splicing, microRNA biogenesis, transcriptional repression, and cell stress responses. Toxicity from cellular aggregation of TDP-43 has been previously shown to impact both neurons and glial cells. This Drosophila model conditionally expresses human TDP43 protein (hTDP43) in a Gal4/UAS system activated by Gal4 drivers. Overexpression of hTDP43 in neurons (Elav/TDP43) or glia (Repo/TDP43) replicates clinical and pathological features of ALS. This dataset had been previously analyzed using TEtranscripts at the family level, which found upregulation of several TE families with hTDP43 expression as compared to control samples with no hTDP43 expression (TDP43/2U). Pan-glial expressing hTDP43 Drosophila brains particularly feature increased TE expression, with upregulation of 23 among the 29 differentially expressed families. Among the 23 TE families increased in Repo/TDP43 Drosophila, Gypsy was of special interest. Gypsy is a retrotranspositionally active endogenous retrovirus that has previously been reported to generate new insertions in aging Drosophila brain. Pan-glial, but not pan-neuronal hTDP43-mediated Gypsy expression and its associated toxicity was confirmed by subfamily-level qRT-PCR, protein immunolabeling and RNAi silencing. Their findings suggested that upregulation of Gypsy in glial cells contribute to the decreased lifespan of Repo/TDP43 Drosophila.

36We investigated how Gypsy upregulation by hTDP43 expression reflected differential expression of individual Gypsy copies at the locus level. To corroborate previous subfamily-level findings, we excluded 44 out of 45 Gypsy subfamilies that did not correspond to the sequence used for Gypsy qPCR, antibody, and RNAi design. Because we were focused on examining the previously reported pan-glial hTDP43 mediated Gypsy upregulation, we excluded 46 Gypsy loci (of the 102) that exhibited negative log2 fold changes in Repo/TDP43 samples compared to controls or had infinite fold changes due to low expression across samples. To gauge the percentage that each remaining locus contributed to Gypsy differential upregulation (% DU contribution), we scaled the log2FC values by the total read count across all samples for each locus and divided by the sum of all scaled log2FC values. We identified 3 loci with the greatest (>10%) contribution to Gypsy differential upregulation in Repo/TDP43 samples (Figure 6A). All three loci were significantly differentially expressed in a pairwise comparison between Repo/TDP43 and control samples (chr2R:log2FC 1.31 padj 6.23e–24, chr3R: log2FC 2.58, padj 1.027e–13, chrX: log2FC 3.76, padj 1.43e–37). In examining the normalized counts at these loci, we were surprised to find that two loci located on chromosomes 3R and X exhibited even greater differential upregulation with pan-neuronal hTDP43 expression than with pan-glial hTDP43 expression, a result that was obscured by previous subfamily-level observation (Figure 6B, left). Only the Gypsy locus on chromosome 2R exhibited greater glial upregulation in Repo/TDP43 samples as previously reported. However, because the 2R Gypsy locus was expressed at higher levels, it made a greater contribution (49.1%) to the total measured Gypsy upregulation observed by TEtranscripts. Although the Gypsy element is still retrotranspositionally active and generating new insertions in the Drosophila genome, the alignments to the 2R locus include unique reads that are unlikely to come from a non-reference Gypsy element (Figure 6B, right). Thus, the expression pattern previously observed at the subfamily level is largely explained by a single locus. Our results demonstrate that subfamily-level analyses miss the locus-specific nature of Gypsy upregulation by hTDP43.

37To benchmark SQuIRE’s memory usage and running time for RNA-seq data of different sequencing depths, we subset the high-depth (mean 263 million reads across eight lanes) HEK293T cell line RNA-seq data into 1, 2 and 3-lane libraries with a mean sequencing depth of 32, 65 and 98 million reads. We evaluated the speed and memory performance of each Quantification and Analysis stage tool for each sequencing depth (Figure 7) using eight parallel threads and 64 Gb of available memory. We found that sequencing depth had the greatest effect on Count, taking 8.6 h to complete the 3-lane library compared to 2.4 h for the one lane library. The other tools took much less time and were less affected by sequencing depth. Map took 1–2 h for the different libraries. Call running time was also independent of library size, but it was greater when including all TE counts (10 min) compared to subfamily counts (2 min). We found that the memory usage of each tool was largely independent of sequencing depth, taking between 39–40 Gb of Memory for Map, 30–32 Gb for Count, and 7–8 Gb for Call.

38Our efforts at making SQuIRE easy to use has resulted in multiple features in addition to its ability to provide locus-level TE quantification (Table 1). To set up SQuIRE involves a simple installation process in which the user can copy and paste lines of code, which includes instructions for setting up prerequisite software. In addition, SQuIRE is the only program that downloads reference annotation for assembled genomes available on UCSC, allowing it to be easily adaptable to a variety of species. For genomes from non-model organisms or organism strains with high divergence from the reference annotation, SQuIRE can also use RepeatMasker software output for even wider compatibility. To ensure that the pipeline is streamlined and that the outputs are reproducible, SQuIRE also implements alignment and differential expression for the user. In making SQuIRE as user-friendly as possible, we intend to improve reproducibility of bioinformatics analyses in the TE field.

39We have developed SQuIRE to characterize TE expression using RNA-seq data. TEs are highly repeated in the genome, which can pose challenges for mapping reads unambiguously to specific transcribed loci. SQuIRE is the first RNA-seq analysis software that provides locus-specific TE expression quantification while also outputting subfamily-level expression estimates (Table 1). Our approach incorporates unambiguously mapping reads as well as ambiguously mapping reads, optimally adjudicating alignments of the latter using an EM algorithm. SQuIRE additionally provides empiric information on the structure of each TE transcript rather than relying on TE annotations, recognizing that TE transcripts can be shorter or longer, and sense or antisense compared to the genomic TE. We have shown that SQuIRE correctly attributes a high percentage of reads originating from TEs using simulated data. For older, retrotranspositionally inactive genomic repeats, SQuIRE very accurately assesses expression. These older elements represent the vast majority of TE loci in the human genome (>96.7%).

40Although the detection of reads is lower for frequently retrotranspositionally active, less divergent TEs (e.g. AluYa5, AluYa8, AluYb8, AluYb9, L1HS), we found that implementation of the EM algorithm improves accuracy and lowers both false positive and false negative calls of whether a TE locus is expressed. This finding also holds in biological settings, where SQuIRE is able to correctly detect L1HS expression when we express an ectopic sequence. It maintains a low false positive rate of misattributing these reads to endogenous L1HS loci. The ongoing activity of TEs also results in a significant number of mobile element insertion variants (MEI). Numerous commonly occurring structural variants owed to retrotransposition are missing in reference genome assemblies. Although these variants are not included in the default SQuIRE pipeline, SQuIRE provides users with two options to query transcription of these repeats. First, SQuIRE can detect transcription of polymorphic elements at the subfamily level. Secondly, SQuIRE can directly use sequences of known, non-reference TE insertion polymorphisms to detect locus-specific expression when these are supplied as a supplement to the reference build. For example, in the human genome, L1HS element sites and sequences can be obtained by targeted TE insertion mapping or whole genome sequencing. Polymorphic TE insertions have been reported to databases such as euL1db, dbRIP and by large studies like the 1000 Genomes Project. Using SQuIRE to detect expression of user-provided, non-reference TE sequences at these loci may be a useful feature for understanding functional consequences of these insertion variants. This confirms that SQuIRE can detect the expression of TEs in the reference genome that have in the past been problematic for global TE RNA expression analysis. For all TEs, SQuIRE provides the convenience of differential TE expression analysis with both locus-specific and subfamily-aggregated outputs.

41The SQuIRE algorithm builds on strategies used by previous TE analysis software. SQuIRE rescues multi-mapping reads aligned to TEs, which improves upon pipelines that only utilize uniquely aligning reads. A similar rescue strategy had been previously applied to multi-mapping CAGE-seq tags. SQuIRE reduces bias against TEs without unique sequence by first normalizing to uniquely alignable length. Although expectation-maximization algorithms have been previously used in TEtranscripts and the RNA-seq quantification software RSEM, SQuIRE differs from these by normalizing to transcribed TE length rather than annotated length. Without transcribed length information, repeated iterations can perpetuate a ‘poor gets poorer’ cycle in which an underestimation of partially transcribed TE expression levels worsens with each iteration. Furthermore, because TEs can be transcribed beyond their annotation, TE-analysis strategies that align to the transcriptome instead of the genome miss potential unique read alignments to flanking sequence and fail to capture the genomic context of TE expression. Here, we show that these additional features in SQuIRE’s Count algorithm improve on the accuracy of TE quantification, as assessed using both simulated reads and orthogonal approaches to measure log2 fold changes in mouse tissue comparisons. Our findings suggest that important biologic insights can be gained by examining TE transcription at the locus level.

42To date, locus-specific studies of TE expression and activity have mostly focused on identifying transcriptionally and retrotranspositionally active L1s in the human genome. While these targeted methods can enrich for expressed TEs, they require tailored approaches for sequencing library preparation and do not yet have accompanying software. SQuIRE provides a complementary tool that supplies a software package applicable to a broad array of conventional RNA-seq datasets. These focused studies have shown that rare, individual loci, widely distributed in the genome generate RNA transcripts. In applying SQuIRE to study locus-specific TE expression genome-wide in mouse tissues and a Drosophila disease model, we can see that this paradigm is not unique to L1s or humans. It seems a limited subset of TE loci are transcribed with complex patterns of tissue-specific expression. Furthermore, we found that the tissue expression patterns of TE loci reflect a variety of transcriptome contexts: broadly expressed mRNA transcripts, tissue-specific lncRNAs, and intrinsically regulated TE transcripts. How these TEs may affect gene regulation or biological processes remain open questions. Genome-wide analyses of TEs have indicated roles for cis-acting elements on transcriptional regulation, transcript splicing, and RNA function. By providing locus-level TE transcript estimations, we expect SQuIRE will enable studies that dissect the impacts of TE expression.

43SQuIRE is freely available for download through https://github.com/wyang17/SQuIRE under the GPL-3 license. The raw sequencing data and SQuIRE Count output for HEK293T cell transfection were deposited to the NCBI Genome Expression Omnibus with accession number GSE113960.