{'fraction_non_alphanumeric': 0.03888888888888889, 'fraction_numerical': 0.0426179604261796, 'mean_word_length': 4.804328621908128, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SQuIRE: Software for Quantifying Interspersed Repeat Elements

Transposable elements are interspersed repeat sequences that make up much of the human genome. Conventional approaches to RNA-seq analysis often exclude these sequences, fail to optimally adjudicate read alignments, or align reads to interspersed repeat consensus sequences without considering these transcripts in their genomic contexts. As a result, repetitive sequence contributions to transcriptomes are not well understood. Here, we present Software for Quantifying Interspersed Repeat Expression (SQuIRE), an RNA-seq analysis pipeline that integrates repeat and genome annotation (RepeatMasker), read alignment (STAR), gene expression (StringTie) and differential expression (DESeq2). SQuIRE uniquely provides a locus-specific picture of interspersed repeat-encoded RNA expression. SQuIRE can be downloaded at (github.com/wyang17/SQuIRE).

3 Due to the repetitive nature of TEs, short-read RNA sequences that originate from one locus can 46 ambiguously align to multiple copies of the same subfamily dispersed throughout the genome. This problem is most significant for younger TEs; older elements have accumulated nucleotide substitutions 48 over millions of years that can differentiate them and give rise to uniquely aligning TE reads (Giordano et 49 al. 2007). Because of these barriers, conventional RNA-seq analyses of TEs have either discarded multi-50 mapping alignments (Chuong et al. 2013) or combined TE expression to the subfamily level (Criscione et

119
In order to accurately quantify TE RNA expression at locus resolution, Count builds on these 120 previous methods by leveraging unique read alignments to each TE to assign fractions of multi-mapping 121 reads (Fig. 2). First, Count identifies reads that map to TEs (by at least 50% of the read length) and labels 122 them as "unique reads" or "multi-mapped reads." Second, Count assigns fractions of a read to each TE as 123 a function of the probability that the TE gave rise to that read. Uniquely aligning reads are considered 124 certain (e.g., probability = 100%, count = 1). Count initially assigns fractions of multi-mapping reads to 125 TEs in proportion to their relative expression as indicated by unique read alignments. 152 After the total counts (unique and multi-mapped) of each TE are re-calculated, multi-mapped reads 153 can be re-assigned in subsequent iterations of expectation (assigning multi-mapped read fractions to TEs) 154 and maximization (summation of unique and multi-mapped fraction counts). These iterations can be 155 repeated until a given iteration number set by the user or until the TE counts converge ("auto", when all 156 of the TEs with ≥ 10 counts change by < 1%). An example of Count output is provided in Supplemental 9 Figure 2. Schematic representation of the SQuIRE Count algorithm. First, Count labels reads as unique (filled arrows) or multi-mapping (empty arrows). Second, Count assigns fractions of multi-mapping reads in proportion to the normalized unique read expression of each TE. The partially filled arrows reflect the proportion of the read assigned to the TE of the corresponding color. Then, Count runs an Expectation-Maximization loop that estimates transcript length and reassigns multi-mapping reads for each TE (E-step), then re-estimates total read counts (M-step) until convergence.

160
To test the performance of Count, we simulated RNA-seq data from 100,000 randomly selected TEs 161 from the human GRCh38/hg38 (hg38) RepeatMasker annotation (see Methods). TEs were simulated with 162 read coverages of ranging from 2-4000X and simulated counts ranging from 2-4588.We first evaluated 163 accuracy by how closely SQuIRE Count output corresponded to the simulated read counts (i.e., % 164 Observed/Expected). However, using this calculation is not meaningful for TEs with low simulated 165 counts: a TE with 0 counts gives an infinite value, and a reported count of 1 for a TE with 2 simulated 166 reads gives a low 50% Observed/Expected. Thus, we were primarily interested in 'expressed' simulated

170
To test how well SQuIRE performed leveraging only uniquely aligning read information, we first 171 evaluated the % Observed/Expected of TE counts with 0 E-M iterations. We found that SQuIRE 172 accurately assigned read counts to most TEs, with a mean % Observed/Expected of 98.79% 173 (Supplemental Fig. S1). We predicted that this accuracy would be lower for TEs with less uniquely 174 alignable sequence. Indeed, SQuIRE was less accurate for elements with less than 10% divergence (mean 175 of 77.35 % Observed/Expected). The most frequently retrotranspositionally active TEs (i.e., AluYa5,

176
AluYa8, AluYb8, AluYb9, and L1HS) had counts ranging from 48-70% Observed/Expected, with a range 177 of 79-92% Observed/Expected at the subfamily level (Supplemental Table S2). This illustrates that even 178 without the EM-algorithm, SQuIRE is sensitive for highly homologous subfamilies at the subfamily level.

179
Given the low recovery of simulated counts for younger elements when relying solely on uniquely 180 aligning reads, we next evaluated how much adding the EM-algorithm improved Count's performance.

181
We anticipated that the counts for most TEs would not change, but that younger elements with less 182 divergence would have improved recovery of simulated reads. Indeed, the overall % Observed/Expected 183 counts of TE loci increased only slightly by 0.14% to a total of 98.93%. However, the change in % 11 Observed/Expected of TEs was much greater for the most homologous active elements, improving by 185 20.47% for young Alu elements and by 21.1% for L1HS loci (Fig. 3). At the subfamily level, the % 186 Observed/Expected of active TEs was improved by 8.1% for young Alu elements and by 2.2% for L1HS 187 (Supplemental Table S2). Using updated transcript information in the EM-algorithm is thus particularly 188 useful for TE biologists interested in younger elements that have previously been problematic to quantify 189 by RNA-seq.

190
We also wanted to evaluate SQuIRE's ability to distinguish whether a TE is expressed or not. To  Table S3). However, using the EM algorithm improved the TPR for Alus (TPR=85.22%-

267
We found that SQuIRE (99.86% ±1.46 %), TETools (100.14 ± 2.21%), and TEtranscripts (95.89 ± 268 16.41%) had comparable % Observed/Expected rates (Supplemental Fig. S3). In contrast, RepEnrich 269 (108.77 ± 40.67%) was less accurate in terms of % Observed/Expected. This is likely attributable to 270 RepEnrich's recommended use of Bowtie 1, which discards discordant reads and limits the number of 271 attempts to align both paired-end mates to repetitive regions. To support this, we compared how often 272 each aligner mapped a uniquely aligning simulated read to the correct location. We indeed found that 273 Bowtie 1 failed to report unique reads more often in a paired-end library compared to single-end 274 (Supplemental Table S6).

275
To compare SQuIRE to other TE analysis tools with biological data, we ran each pipeline on 276 publically available adult C57Bl/6 mouse tissue RNA-seq data (Brawand et al. 2011) using 277 GRCm38/mm10 (mm10) TE annotation. We compared the expression of subfamilies in testis compared 278 to pooled data from brain, heart, kidney, and liver tissues. To independently evaluate the fold-changes of 279 TE RNA between testis and somatic tissues, we also used our previously published adult C57Bl/6 mouse 280 Nanostring results (Gnanakkan et al. 2013). Unlike RNA-seq analysis, which infers transcript levels by 281 counting reads, Nanostring uses uniquely mapping probes to capture and count RNA molecules. We 282 compared the Nanostring log 2 fold changes (log 2 FC) of TE subfamily expression in testis and pooled 283 somatic tissue to the log 2 FC values found by SQuIRE, RepEnrich, TEtranscripts, and TETools 284 (Supplemental Fig. S4). We first looked at how often the direction of fold change corresponded between 285 each tool and Nanostring. For the 16 subfamilies queried, SQuIRE and TETools shared the same direction 18 of fold change as Nanostring more often than the other tools (SQuIRE: 12, TETools: 12, TEtranscripts: 9,

290
With SQuIRE, we can closely examine the mouse RNA-seq data at the locus level. For the 16 291 subfamilies analyzed by Nanostring and the TE analysis tools, we found that the reported subfamily-level 292 expression could be attributed to fewer than 7% of each subfamily's loci (Supplemental Fig. S5). This 293 suggests that regulation of TE transcription is not necessarily shared across all TEs from the same 294 subfamily. On the other hand, whereas the other subfamilies studied by Nanostring have only 1-4 295 significantly differentially expressed loci (log2FC >1, padj < 0.05), the IAPLTR3 subfamily has 11 loci 296 that are all differentially expressed in testis compared to somatic tissues (Fig. 5A). To test whether this 297 was an enrichment relative to the representation of IAPLTR3 in the mouse genome, we performed a 298 Fisher's exact test and found that IAPLTR3 loci were 10-fold more likely than expected to be 299 differentially expressed in testis (OR: 10.56, 95%CI: 5.25-18.97, padj < 1.61 e-08). This suggests that a 300 subset of TE locus expression may still be impacted by subfamily-specific regulation.

301
To further investigate the interplay between genomic context and TE subfamily, we identified the 302 closest genes to each differentially expressed locus and clustered the loci by their expression levels, as 303 shown in Figure 5A. We found a cluster of 3 loci exhibiting broad expression across somatic tissues from 304 the IAP1, MERVL, and MURVY LTR retrotransposon subfamilies. When we examined the genomic 305 context of these 3 loci, we found that all were located within genes with known broad tissue expression can vary across loci from the same subfamily. We also examined ERVB4-1B and IAPLTR3, the two LTR 314 retrotransposon subfamilies that exhibited the highest fold change by Nanostring. These subfamilies were 315 represented in the high-expressing, tissue-restricted loci cluster (Fig. 5A). While the transcription of the 316 ERVB4-1B locus on chr13 did not extend beyond annotations for that subfamily (Fig. 5C), the IAPLTR3 317 loci on chr14 (Fig. 5D) and chr18 are part of longer transcripts that initiate outside of the annotated TE.  Differentially expressed TE loci belonging to subfamilies previously analyzed by Nanostring a. The X-axis represents replicates of somatic and testis tissue samples from adult C57Bl/6 mouse. The Y-axis represents differentially expressed TE loci. The heatmap colors represent the log2 of total read counts +1 for each TE locus. b-d. Examples of intergenic TE loci differentially expressed in testis compared to somatic tissues. Tracks from brain, heart, kidney and liver replicates were collapsed into a single track. The scales of count expression are shown in brackets. The RefSeq track represents annotated genes. The RepeatMasker track represents transposable elements annotated in the reference genome. Transposable elements colored in red belong to the subfamily indicated; dark red indicates that that RepeatMasker entry meets significant differential expression thresholds (log2FC > 2, padj < 0.05).

339
Our efforts at making SQuIRE easy to use has resulted in a simple installation process in which the 340 user can copy and paste lines of code to install all prerequisite software and set up SQuIRE (Table 1). In 341 addition, SQuIRE is the only program that downloads reference annotation for assembled genomes 342 available on UCSC, allowing it to be easily adaptable to a variety of species. For genomes from non-343 model organisms or organism strains with high divergence from the reference annotation, SQuIRE can 344 also use RepeatMasker software output for even wider compatibility. To ensure that the pipeline is 345 streamlined and that the outputs are reproducible, SQuIRE also implements alignment and differential 346 expression for the user. In making SQuIRE as user-friendly as possible, we intend to improve the 347 reproducibility of bioinformatics in the TE field.

436
We randomly selected 100,000 TEs from the hg38 Repeatmasker annotation downloaded by Fetch.

437
We limited our list of potential TEs to those included in TEtranscripts (Jin et al. 2015) and RepEnrich 438 (Criscione et al. 2014) to enable comparisons between these different programs.  website. We aligned the data to the genome with STAR using "--winAnchorMultimapNmax 100","--496 outFilterMultimapNmax 100" parameters for multi-mapping. We then ran TEtranscripts with the "--mode 497 multi" setting to utilize its expectation-maximization algorithm for assigning multi-reads for the resulting 498 SAM file. Since TEtranscripts analyzes TE and gene expression together, we used refGene annotation 499 obtained by SQuIRE Fetch for the required gtf file. We used the parameters "--format SAM", "--mode 500 multi", "--stranded yes" for simulated human data, and "--format SAM", "--mode multi", "--stranded no" 501 for mouse data.