atena package provides methods to quantify the expression of transposable elements within R and Bioconductor.
Transposable elements (TEs) are autonomous mobile genetic elements. They are DNA sequences that have, or once had, the ability to mobilize within the genome either directly or through an RNA intermediate (Payer and Burns 2019). TEs can be categorized into two classes based on the intermediate substrate propagating insertions (RNA or DNA). Class I TEs, also called retrotransposons, first transcribe an RNA copy that is then reverse transcribed to cDNA before inserting in the genome. In turn, these can be divided into long terminal repeat (LTR) retrotransposons, which refer to endogenous retroviruses (ERVs), and non-LTR retrotransposons, which include long interspersed element class 1 (LINE-1 or L1) and short interspersed elements (SINEs). Class II TEs, also known as DNA transposons, directly excise themselves from one location before reinsertion. TEs are further split into families and subfamilies depending on various structural features (Goerner-Potvin and Bourque 2018; Guffanti et al. 2018).
Most TEs have lost the capacity for generating new insertions over their evolutionary history and are now fixed in the human population. Their insertions have resulted in a complex distribution of interspersed repeats comprising almost half (50%) of the human genome (Payer and Burns 2019).
TE expression has been observed in association with physiological processes in a wide range of species, including humans where it has been described to be important in early embryonic pluripotency and development. Moreover, aberrant TE expression has been associated with diseases such as cancer, neurodegenerative disorders, and infertility (Payer and Burns 2019).
The study of TE expression faces one main challenge: given their repetitive nature, the majority of TE-derived reads map to multiple regions of the genome and these multi-mapping reads are consequently discarded in standard RNA-seq data processing pipelines. For this reason, specific software packages for the quantification of TE expression have been developed (Goerner-Potvin and Bourque 2018), such as TEtranscripts (Jin et al. 2015), ERVmap (Tokuyama et al. 2018) and Telescope (Bendall et al. 2019). The main differences between these three methods are the following:
ERVmap (Tokuyama et al. 2018) is based on selective filtering of multi-mapping reads. It applies filters that consist in discarding reads when the ratio of sum of hard and soft clipping to the length of the read (base pair) is greater than or equal to 0.02, the ratio of the edit distance to the sequence read length (base pair) is greater or equal to 0.02 and/or the difference between the alignment score from BWA (field AS) and the suboptimal alignment score from BWA (field XS) is less than 5.
Telescope (Bendall et al. 2019) reassigns multi-mapping reads to TEs using their relative abundance, which like in TEtranscripts, is also estimated using an EM algorithm. The main differences with respect to TEtranscripts are: (1) Telescope works with an additional parameter for each TE that estimates the proportion of multi-mapping reads that need to be reassigned to that TE; (2) that reassignment parameter is optimized during the EM algorithm jointly with the TE relative abundances, using a Bayesian maximum a posteriori (MAP) estimate that allows one to use prior values on these two parameters; and (3) using the final estimates on these two parameters, multi-mapping reads can be flexibly reassigned to TEs using different strategies, where the default one is to assign a multi-mapping read to the TE with largest estimated abundance and discard those multi-mapping reads with ties on those largest abundances.
Because these tools were only available outside R and Bioconductor, the
atena package provides a complete re-implementation in R of these three methods to facilitate the integration of TE expression quantification into Bioconductor workflows for the analysis of RNA-seq data.
Another challenge in TE expression quantification is the lack of complete TE
annotations due to the difficulty to correctly place TEs in genome assemblies (Goerner-Potvin and Bourque 2018). The gold standard for TE annotations are
RepeatMasker annotations, available through the RepeatMasker track in UCSC
atena can fetch RepeatMasker annotations with the function
annotaTEs(). Moreover, this function can parse TE annotations by applying
parsefun. Examples of
parsefun included in
rmskidentity(): returns RepeatMasker annotations without any modification.
rmskbasicparser(): filters out non-TE repeats and elements without strand information from RepeatMasker annotations. Then assigns a unique id to each elements based on their repeat name.
rmskatenaparser(): attempts to reconstruct fragmented TEs by assembling together fragments from the same TE that are close enough. For LTR class TEs, tries to reconstruct full-length and partial TEs following the LTR - internal region - LTR structure.
OneCodeToFindThemAll(): implementation of the
OneCodeToFindThemAll.pl(Bailly-Bechet, Haudry, and Lerat 2014) tool for parsing RepeatMasker output files.
OneCodeToFindThemAll() functions try to address
the annotation fragmentation present in the output file of the RepeatMasker
software (i.e. presence of multiple hits (homology-based matches) corresponding
to a unique copy of an element). This is highly frequent for LTR class TEs,
where the consensus sequences are split into LTR and internal regions
separately, causing RepeatMasker to also report these two regions of a TE as
separate elements. These two functions attempt to identify these and other
cases of fragmented annotations and assemble them together into single
elements. To do so, the assembled elements must satisfy certain criteria.
Differences in these criteria, as
well as different approaches for finding equivalences between LTR and internal
regions to reconstruct LTR retrotransposons, is what differences these two
As an example, let’s retrieve TE annotations for Drosophila melanogaster
dm6 genome version. By setting
RepeatMasker annotations are not modified:
library(atena) library(GenomicRanges) rmskid <- annotaTEs(genome = "dm6", parsefun = rmskidentity) rmskid GRanges object with 137555 ranges and 11 metadata columns: seqnames ranges strand | swScore milliDiv milliDel <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric>  chr2L 2-154 + | 778 167 7  chr2L 313-408 + | 296 174 207  chr2L 457-612 + | 787 170 7  chr2L 771-866 + | 296 174 207  chr2L 915-1070 + | 787 170 7 ... ... ... ... . ... ... ...  chrUn_DS486004v1 99-466 - | 3224 14 0  chrUn_DS486005v1 1-1001 + | 930 48 0  chrUn_DS486008v1 1-488 + | 4554 0 0  chrUn_DS486008v1 489-717 - | 2107 9 0  chrUn_DS486008v1 717-1001 - | 2651 3 0 milliIns genoLeft repName repClass repFamily <numeric> <integer> <character> <character> <character>  20 -23513558 HETRP_DM Satellite Satellite  42 -23513304 HETRP_DM Satellite Satellite  19 -23513100 HETRP_DM Satellite Satellite  42 -23512846 HETRP_DM Satellite Satellite  19 -23512642 HETRP_DM Satellite Satellite ... ... ... ... ... ...  3 -535 ROVER-LTR_DM LTR Gypsy  1 0 (TATACATA)n Simple_repeat Simple_repeat  0 -513 NOMAD_LTR LTR Gypsy  0 -284 ACCORD_LTR LTR Gypsy  0 0 DMRT1A LINE R1 repStart repEnd repLeft <integer> <integer> <integer>  1519 1669 -203  1519 1634 -238  1516 1669 -203  1519 1634 -238  1516 1669 -203 ... ... ... ...  0 367 1  1 1000 0  31 518 0  -123 435 207  0 5183 4899 ------- seqinfo: 1870 sequences (1 circular) from dm6 genome
We get annotations with 137555 elements.
Now, let’s obtain the same annotations but processes them using the
rmskatenaparser function. We set parameter
strict = FALSE to avoid
applying a filter of minimum 80% identity with the consensus sequence and
minimum 80 bp length. The
insert parameter is set to 500 meaning that two
elements with the same name are merged if they are closer than 500 bp.
rmskat <- annotaTEs(genome = "dm6", parsefun = rmskatenaparser, strict = FALSE, insert = 500) loading from cache rmskat GRangesList object of length 1: $IDEFIX_LTR.1 GRanges object with 1 range and 11 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric>  chr2L 9726-9859 + | 285 235 64 15 genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer>  -23503853 IDEFIX_LTR LTR Gypsy 425 565 repLeft <integer>  29 ------- seqinfo: 1870 sequences (1 circular) from dm6 genome
How many elements are present in the annotations?
length(rmskat)  22848
As expected, we get a lower number of elements in the annotations because repeats that are not TEs have been removed. Furthermore, some fragmented regions of TEs have been assembled together.
rmskat object is of class
GRangesList. Each element of the
list represents an assembled TE containing a
GRanges object of length 1
(when the TE was not assembled with another element) or length > 1 (when
two or more fragments were assembled together into a single TE).
We can get more information of the parsed annotations by accessing the
metadata columns with
mcols(rmskat) DataFrame with 22848 rows and 3 columns status Rel_length Class <character> <numeric> <character> IDEFIX_LTR.1 LTR 0.225589 LTR DNAREP1_DM.2 noLTR 0.176768 DNA DNAREP1_DM.3 noLTR 0.151515 DNA DNAREP1_DM.4 noLTR 0.419192 DNA DNAREP1_DM.5 noLTR 0.409091 DNA ... ... ... ... Gypsy12_I-int.22844 noLTR 0.18211323 LTR PROTOP.22845 noLTR 0.20848214 DNA TAHRE.22846 noLTR 0.06518207 LINE TART.22847 noLTR 0.00780548 LINE FW_DM.22848 noLTR 0.08823529 LINE
There is information about the reconstruction status of the TE (status column), the relative length of the reconstructed TE (Rel_length) and the repeat class of the TE (Class). The relative length is computed by adding the length (in base pairs) of all fragments from the same assembled TE and dividing the sum by the length (in base pairs) of the consensus sequence. For full-length and partially reconstructed LTR TEs, the consensus sequence length used is the one resulting from adding twice the consensus sequence length of the long terminal repeat (LTR) and the one from the corresponding internal region. For solo-LTRs, the consensus sequence length of the long terminal repeat is used.
We can get an insight into the composition of the assembled annotations using the information from the status column. Let’s look at the absolute frequencies of the status and class of TEs in the annotations.
Here, full-lengthLTR are reconstructed LTR retrotransposons following the LTR - internal region (int) - LTR structure. Partially reconstructed LTR TEs are partialLTR_down (internal region followed by a downstream LTR) and partialLTR_up (LTR upstream of an internal region). int and LTR correspond to internal and solo-LTR regions, respectively. Finally, the noLTR refers to TEs of other classes (not LTR), as well as TEs of class LTR which could not be identified as either internal or long terminal repeat regions based on their name.
atena package provides useful functions to retrieve
TEs of a specific class, using a specific relative length threshold.
Those TEs with higher relative lengths are more likely to have intact open
reading frames, making them more interesting for expression quantification
and functional analyses. For example, to get LINEs with a minimum of 0.9
relative length, we can use the
getLINEs() function. We use the previously
obtained object (
rmskat) with annotations parsed with the
function, and set the
relLength to 0.9.
rmskLINE <- getLINEs(rmskat, relLength = 0.9) length(rmskLINE)  355 rmskLINE GRangesList object of length 1: $LINEJ1_DM.6 GRanges object with 1 range and 11 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric>  chr2L 47514-52519 + | 43674 5 0 0 genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer>  -23461193 LINEJ1_DM LINE Jockey 2 5007 repLeft <integer>  13 ------- seqinfo: 1870 sequences (1 circular) from dm6 genome
To get LTR retrotransposons, we can use the function
function also allows to get one or more specific types of reconstructed TEs.
To get full-length, partial LTRs and other fragments that could not be
reconstructed, we can:
rmskLTR <- getLTRs(rmskat, relLength = 0.8, full_length = TRUE, partial = TRUE, otherLTR = TRUE) length(rmskLTR)  1711 rmskLTR GRangesList object of length 1: $`BLOOD_I-int.20` GRanges object with 5 ranges and 11 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric>  chr2L 347941-348153 - | 3447 3 43 0  chr2L 348188-348355 - | 3447 3 43 0  chr2L 348356-354968 - | 60509 0 0 0  chr2L 354969-355181 - | 3434 5 43 0  chr2L 355216-355383 - | 3434 5 43 0 genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer>  -23165559 BLOOD_LTR LTR Gypsy 0 399  -23165357 BLOOD_LTR LTR Gypsy 213 186  -23158744 BLOOD_I-int LTR Gypsy 0 6613  -23158531 BLOOD_LTR LTR Gypsy 0 399  -23158329 BLOOD_LTR LTR Gypsy 213 186 repLeft <integer>  187  2  1  187  2 ------- seqinfo: 1870 sequences (1 circular) from dm6 genome
To obtain DNA transposons and SINEs, one can use the
getSINEs() functions, respectively.
Quantification of TE expression with
atena consists in the following two
Building of a parameter object for one of the available quantification methods.
Calling the TE expression quantification method
qtex() using the previously
built parameter object.
The dataset that will be used to illustrate how to quantify TE expression with
atena is a published RNA-seq dataset of Drosophila melanogaster available
at the National Center for Biotechnology Information (NCBI) Gene Expression
Omnibus (accession no. GSE47006).
The two selected samples are: a piwi knockdown and a piwi control (GSM1142845
and GSM1142844). These files have been subsampled. The piwi-associated
silencing complex (piRISC) silences TEs in the Drosophila ovary, thus, the
knockdown of piwi causes the de-repression of TEs.
Here, the expression of full-length LTR retrotransposons present in
will be quantified.
To use the ERVmap method in
atena we should first build an object of the class
ERVmapParam using the function
singleEnd parameter is set to
TRUE since the example BAM files are single-end. The
ignoreStrand parameter works analogously to the same parameter in the function
summarizeOverlaps() from package GenomicAlignments and should be set to
TRUE whenever the RNA library preparation protocol was stranded.
One of the filters applied by the ERVmap method compares the alignment score of a given primary alignment, stored in the
AS tag of a SAM record, to the largest alignment score among every other secondary alignment, known as the suboptimal alignment score. The original ERVmap software assumes that input BAM files are generated using the Burrows-Wheeler Aligner (BWA) software (Li and Durbin 2009), which stores suboptimal alignment scores in the
XS tag. Although
AS is an optional tag, most short-read aligners provide this tag with alignment scores in BAM files. However, the suboptimal alignment score, stored in the
XS tag by BWA, is either stored in a different tag or not stored at all by other short-read aligner software, such as STAR (Dobin et al. 2013).
To enable using ERVmap on BAM files produced by short-read aligner software other than BWA,
atena allows the user to set the argument
suboptimalAlignmentTag to one of the following three possible values:
The name of a tag different to
XS that stores the suboptimal alignment score.
The value “none”, which will trigger the calculation of the suboptimal alignment score by searching for the largest value stored in the
AS tag among all available secondary alignments.
The value “auto” (default), by which
atena will first extract the name of the short-read aligner software from the BAM file and if that software is BWA, then suboptimal alignment scores will be obtained from the
XS tag. Otherwise, it will trigger the calculation previously explained for
Finally, this filter is applied by comparing the difference between alignment and suboptimal alignment scores to a cutoff value, which by default is 5 but can be modified using the parameter
suboptimalAlignmentCutoff. The default value 5 is the one employed in the original ERVmap software that assumes the BAM file was generated with BWA and for which lower values are interpreted as “equivalent to second best match has one or more mismatches than the best match” (Tokuyama et al. 2018, pg. 12571). From a different perspective, in BWA the mismatch penalty has a value of 4 and therefore, a
suboptimalAlignmentCutoff value of 5 only retains those reads where the suboptimal alignment has at least 1 mismatch more than the best match. Therefore, the
suboptimalAlignmentCutoff value is specific to the short-read mapper software and we recommend to set this value according to the mismatch penalty of that software. Another option is to set
suboptimalAlignmentCutoff=NA, which prevents the filtering of reads based on this criteria, as set in the following example.
bamfiles <- list.files(system.file("extdata", package="atena"), pattern="*.bam", full.names=TRUE) empar <- ERVmapParam(bamfiles, teFeatures = rmskLTR, singleEnd = TRUE, ignoreStrand = TRUE, suboptimalAlignmentCutoff=NA) empar ERVmapParam object # BAM files (2): control_KD.bam, piwi_KD.bam # features (1711): BLOOD_I-int.20, ..., NOMAD_LTR.22823 # single-end, unstranded
In the case of paired-end BAM files (
singleEnd=FALSE), two additional arguments can be specified,
strandMode defines the behavior of the strand getter when internally reading the BAM files with the
GAlignmentPairs() function. See the help page of
strandMode in the GenomicAlignments package for further details.
fragments controls how read filtering and counting criteria are applied to the read mates in a paired-end read. To use the original ERVmap algorithm (Tokuyama et al. 2018) one should set
fragments=TRUE (default when
singleEnd=FALSE), which filters and counts each mate of a paired-end read independently (i.e., two read mates overlapping the same feature count twice on that feature, treating paired-end reads as if they were single-end). On the other hand, when
fragments=FALSE, if the two read mates pass the filtering criteria and overlap the same feature, they count once on that feature. If either read mate fails to pass the filtering criteria, then both read mates are discarded.
An additional functionality with respect to the original ERVmap software is the integration of gene and TE expression quantification. The original ERVmap software doesn’t quantify TE and gene expression coordinately and this can potentially lead to counting twice reads that simultaneously overlap a gene and a TE. In
atena, gene expression is quantified based on the approach used in the TEtranscripts software (Jin et al. 2015): unique reads are preferably assigned to genes, whereas multi-mapping reads are preferably assigned to TEs.
In case that a unique read does not overlap a gene or a multi-mapping read does not overlap a TE,
atena searches for overlaps with TEs or genes, respectively. Given the different treatment of unique and multi-mapping reads,
atena requires the information regarding the unique or multi-mapping status of a read. This information is obtained from the presence of secondary alignments in the BAM file or, alternatively, from the
NH tag in the BAM file (number of reported alignments that contain the query in the current SAM record). Therefore, either secondary alignments or the
NH tag need to be present for gene expression quantification.
The original ERVmap approach does not discard any read overlapping gene annotations. However, this can be changed using the parameter
geneCountMode, which by default
geneCountMode="all" and follows the behavior in the original ERVmap method. On the contrary, by setting
atena also applies the filtering criteria employed to quantify TE expression to the reads overlapping gene annotations.
atena also allows one to aggregate TE expression quantifications. By default, the names of the input
GRangesList object given in the
teFeatures parameter are used to aggregate quantifications. However, the
aggregateby parameter can be used to specify other column names in the feature annotations to be used to aggregate TE counts, for example at the sub-family level.
To use the Telescope method for TE expression quantification, the
TelescopeParam() function is used to build a parameter object of the class
As in the case of
aggregateby argument, which should be a character vector of column names in the annotation, determines the columns to be used to aggregate TE expression quantifications. This way,
atena provides not only quantifications at the subfamily level, but also allows to quantify TEs at the desired level (family, class, etc.), including locus based quantifications. For such a use case, the object with the TE annotations should include a column with unique identifiers for each TE locus and the
aggregateby argument should specify the name of that column. When
aggregateby is not specified, the
names() of the object containing TE annotations are used to aggregate quantifications.
Here, TE quantifications will be aggregated according to the
names() of the
bamfiles <- list.files(system.file("extdata", package="atena"), pattern="*.bam", full.names=TRUE) tspar <- TelescopeParam(bfl=bamfiles, teFeatures=rmskLTR, singleEnd = TRUE, ignoreStrand=TRUE) tspar TelescopeParam object # BAM files (2): control_KD.bam, piwi_KD.bam # features (CompressedGRangesList length 1711): BLOOD_I-int.20, ..., NOMAD_LTR.22823 # aggregated by: CompressedGRangesList names # single-end; unstranded
In case of paired-end data (
singleEnd=FALSE), the argument usage is similar to that of
ERVmapParam(). In relation to the BAM file, Telescope follows the same approach as the ERVmap method: when
fragments=FALSE, only mated read pairs from opposite strands are considered, while when
fragments=TRUE, same-strand pairs, singletons, reads with unmapped pairs and other fragments are also considered by the algorithm. However, there is one important difference with respect to the counting approach followed by ERVmap: when
fragments=TRUE mated read pairs mapping to the same element are counted once, whereas in the ERVmap method they are counted twice.
As in the ERVmap method from
atena, the gene expression quantification method in Telescope is based on the approach of the TEtranscripts software (Jin et al. 2015). This way,
atena provides the possibility to integrate TE expression quantification by Telescope with gene expression quantification. As in the case of the ERVmap method from
atena, either secondary alignments or the
NH tag are required for gene expression quantification.
Finally, the third method available is TEtranscripts. First, the
TEtranscriptsParam() function is called to build a parameter object of the class
TEtranscriptsParam. The usage of the
aggregateby argument is the same as in
ERVmapParam(). Locus based quantifications in the TEtranscripts method from
atena is possible because the TEtranscripts algorithm actually computes TE quantifications at the locus level and then sums up all instances of each TE subfamily to provide expression at the subfamily level. By avoiding this last step,
atena can provide TE expression quantification at the locus level using the TEtranscripts method. For such a use case, the object with the TE annotations should include a column with unique identifiers for each TE and the
aggregateby argument should specify the name of that column.
In this example, the
aggregateby argument will be set to
aggregateby = "repName" in order to aggregate quantifications at the repeat
name level. Moreover, gene expression will also be quantified. To do so,
gene annotations are loaded from a TxDb object.
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene) txdb <- TxDb.Dmelanogaster.UCSC.dm6.ensGene txdb_genes <- exonsBy(txdb, by = "gene") length(txdb_genes)  17807
bamfiles <- list.files(system.file("extdata", package="atena"), pattern="*.bam", full.names=TRUE) ttpar <- TEtranscriptsParam(bamfiles, teFeatures = rmskLTR, geneFeatures = txdb_genes, singleEnd = TRUE, ignoreStrand=TRUE, aggregateby = c("repName")) ttpar TEtranscriptsParam object # BAM files (2): control_KD.bam, piwi_KD.bam # features (CompressedGRangesList length 86732): BLOOD_I-int.20, ..., FBgn0286941 # aggregated by: repName # single-end; unstranded
For paired-end data (
singleEnd=FALSE), the usage of the
fragments argument is the same as in
Regarding gene expression quantification,
atena has implemented the approach of the original TEtranscripts software (Jin et al. 2015). As in the case of the ERVmap and Telescope methods from
atena, either secondary alignments or the
NH tag are required.
Following the gene annotation processing present in the TEtranscripts algorithm, in case that
geneFeatures contains a metadata column named “type”, only the elements with “type” = “exon” are considered for the analysis. Then, exon counts are summarized to the gene level in a
GRangesList object. This also applies to the ERVmap and Telescope methods for
atena when gene features are present. Let’s see an example of this processing:
# Creating an example of gene annotations annot_gen <- GRanges(seqnames = rep("2L",8), ranges = IRanges(start = c(1,20,45,80,110,130,150,170), width = c(10,20,35,10,5,15,10,25)), strand = "*", type = rep("exon",8)) # Setting gene ids names(annot_gen) <- paste0("gene",c(rep(1,3),rep(2,4),rep(3,1))) annot_gen GRanges object with 8 ranges and 1 metadata column: seqnames ranges strand | type <Rle> <IRanges> <Rle> | <character> gene1 2L 1-10 * | exon gene1 2L 20-39 * | exon gene1 2L 45-79 * | exon gene2 2L 80-89 * | exon gene2 2L 110-114 * | exon gene2 2L 130-144 * | exon gene2 2L 150-159 * | exon gene3 2L 170-194 * | exon ------- seqinfo: 1 sequence from an unspecified genome; no seqlengths ttpar_gen <- TEtranscriptsParam(bamfiles, teFeatures = rmskLTR, geneFeatures = annot_gen, singleEnd = TRUE, ignoreStrand=TRUE) ttpar_gen TEtranscriptsParam object # BAM files (2): control_KD.bam, piwi_KD.bam # features (CompressedGRangesList length 1714): BLOOD_I-int.20, ..., gene3 # aggregated by: CompressedGRangesList names # single-end; unstranded
Let’s see the result of the gene annotation processing:
features_tt <- atena::features(ttpar_gen) features_tt[!attributes(features_tt)$isTE$isTE] GRangesList object of length 3: $gene1 GRanges object with 3 ranges and 13 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric> gene1 chr2L 1-10 * | <NA> NA NA NA gene1 chr2L 20-39 * | <NA> NA NA NA gene1 chr2L 45-79 * | <NA> NA NA NA genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer> gene1 <NA> <NA> <NA> <NA> <NA> <NA> gene1 <NA> <NA> <NA> <NA> <NA> <NA> gene1 <NA> <NA> <NA> <NA> <NA> <NA> repLeft isTE type <integer> <logical> <character> gene1 <NA> FALSE exon gene1 <NA> FALSE exon gene1 <NA> FALSE exon ------- seqinfo: 1870 sequences (1 circular) from dm6 genome $gene2 GRanges object with 4 ranges and 13 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric> gene2 chr2L 80-89 * | <NA> NA NA NA gene2 chr2L 110-114 * | <NA> NA NA NA gene2 chr2L 130-144 * | <NA> NA NA NA gene2 chr2L 150-159 * | <NA> NA NA NA genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer> gene2 <NA> <NA> <NA> <NA> <NA> <NA> gene2 <NA> <NA> <NA> <NA> <NA> <NA> gene2 <NA> <NA> <NA> <NA> <NA> <NA> gene2 <NA> <NA> <NA> <NA> <NA> <NA> repLeft isTE type <integer> <logical> <character> gene2 <NA> FALSE exon gene2 <NA> FALSE exon gene2 <NA> FALSE exon gene2 <NA> FALSE exon ------- seqinfo: 1870 sequences (1 circular) from dm6 genome $gene3 GRanges object with 1 range and 13 metadata columns: seqnames ranges strand | swScore milliDiv milliDel milliIns <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric> <numeric> gene3 chr2L 170-194 * | <NA> NA NA NA genoLeft repName repClass repFamily repStart repEnd <integer> <character> <character> <character> <integer> <integer> gene3 <NA> <NA> <NA> <NA> <NA> <NA> repLeft isTE type <integer> <logical> <character> gene3 <NA> FALSE exon ------- seqinfo: 1870 sequences (1 circular) from dm6 genome
Finally, to quantify TE expression we call the
qtex() method using one of the previously defined parameter objects (
TelescopeParam) according to the quantification method we want to use. The
qtex() method returns a
SummarizedExperiment object containing the resulting quantification of expression in an assay slot. Additionally, when a
DataFrame, object storing phenotypic data is passed to the
qtex() function through the
phenodata parameter, this will be included as column data in the resulting
SummarizedExperiment object and the row names of these phenotypic data will be set as column names in the output
In the current example, the call to quantify TE expression using the ERVmap method would be the following:
emq <- qtex(empar)
emq class: RangedSummarizedExperiment dim: 1711 2 metadata(0): assays(1): counts rownames(1711): BLOOD_I-int.20 ROO_LTR.23 ... MDG1_LTR.22807 NOMAD_LTR.22823 rowData names(4): status Rel_length Class isTE colnames(2): control_KD piwi_KD colData names(0): colSums(assay(emq)) control_KD piwi_KD 137 131
In the case of the Telescope method, the call would be as follows:
tsq <- qtex(tspar)
tsq class: RangedSummarizedExperiment dim: 1712 2 metadata(0): assays(1): counts rownames(1712): BLOOD_I-int.20 ROO_LTR.23 ... NOMAD_LTR.22823 no_feature rowData names(4): status Rel_length Class isTE colnames(2): control_KD piwi_KD colData names(0): colSums(assay(tsq)) control_KD piwi_KD 150 126
For the TEtranscripts method, TE expression is quantified by using the following call:
ttq <- qtex(ttpar)
ttq class: RangedSummarizedExperiment dim: 85131 2 metadata(0): assays(1): counts rownames(85131): ACCORD2_I-int ACCORD2_LTR ... FBgn0286941 FBgn0286941 rowData names(6): status Rel_length ... exon_name isTE colnames(2): control_KD piwi_KD colData names(0): colSums(assay(ttq)) control_KD piwi_KD 150 133
As mentioned, TE expression quantification is provided at the repeat name level.
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:  LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C  LC_TIME=en_GB LC_COLLATE=C  LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8  LC_PAPER=en_US.UTF-8 LC_NAME=C  LC_ADDRESS=C LC_TELEPHONE=C  LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C time zone: America/New_York tzcode source: system (glibc) attached base packages:  stats4 stats graphics grDevices utils datasets methods  base other attached packages:  TxDb.Dmelanogaster.UCSC.dm6.ensGene_3.12.0  GenomicFeatures_1.52.0  AnnotationDbi_1.62.0  RColorBrewer_1.1-3  atena_1.6.0  SummarizedExperiment_1.30.0  Biobase_2.60.0  GenomicRanges_1.52.0  GenomeInfoDb_1.36.0  IRanges_2.34.0  S4Vectors_0.38.0  BiocGenerics_0.46.0  MatrixGenerics_1.12.0  matrixStats_0.63.0  knitr_1.42  BiocStyle_2.28.0 loaded via a namespace (and not attached):  tidyselect_1.2.0 dplyr_1.1.2  blob_1.2.4 filelock_1.0.2  Biostrings_2.68.0 bitops_1.0-7  fastmap_1.1.1 RCurl_1.98-1.12  BiocFileCache_2.8.0 GenomicAlignments_1.36.0  promises_220.127.116.11 XML_3.99-0.14  digest_0.6.31 mime_0.12  lifecycle_1.0.3 ellipsis_0.3.2  KEGGREST_1.40.0 interactiveDisplayBase_1.38.0  RSQLite_2.3.1 magrittr_2.0.3  compiler_4.3.0 progress_1.2.2  rlang_1.1.0 sass_0.4.5  tools_4.3.0 utf8_1.2.3  yaml_2.3.7 rtracklayer_1.60.0  prettyunits_1.1.1 bit_4.0.5  curl_5.0.0 DelayedArray_0.26.0  xml2_1.3.3 BiocParallel_1.34.0  withr_2.5.0 purrr_1.0.1  grid_4.3.0 fansi_1.0.4  xtable_1.8-4 biomaRt_2.56.0  cli_3.6.1 rmarkdown_2.21  crayon_1.5.2 generics_0.1.3  rjson_0.2.21 httr_1.4.5  DBI_1.1.3 cachem_1.0.7  stringr_1.5.0 zlibbioc_1.46.0  parallel_4.3.0 restfulr_0.0.15  BiocManager_1.30.20 XVector_0.40.0  vctrs_0.6.2 Matrix_1.5-4  jsonlite_1.8.4 bookdown_0.33  hms_1.1.3 bit64_4.0.5  magick_2.7.4 jquerylib_0.1.4  glue_1.6.2 codetools_0.2-19  stringi_1.7.12 BiocVersion_3.17.1  later_1.3.0 BiocIO_1.10.0  tibble_3.2.1 pillar_1.9.0  rappdirs_0.3.3 htmltools_0.5.5  GenomeInfoDbData_1.2.10 R6_2.5.1  dbplyr_2.3.2 sparseMatrixStats_1.12.0  evaluate_0.20 shiny_1.7.4  lattice_0.21-8 highr_0.10  AnnotationHub_3.8.0 png_0.1-8  Rsamtools_2.16.0 memoise_2.0.1  SQUAREM_2021.1 httpuv_1.6.9  bslib_0.4.2 Rcpp_1.0.10  xfun_0.39 pkgconfig_2.0.3
Bailly-Bechet, Marc, Annabelle Haudry, and Emmanuelle Lerat. 2014. “‘One Code to Find Them All’: A Perl Tool to Conveniently Parse Repeatmasker Output Files.” Mobile DNA 5 (1): 1–15.
Bendall, Matthew L, Miguel De Mulder, Luis Pedro Iñiguez, Aarón Lecanda-Sánchez, Marcos Pérez-Losada, Mario A Ostrowski, R Brad Jones, et al. 2019. “Telescope: Characterization of the Retrotranscriptome by Accurate Estimation of Transposable Element Expression.” PLoS Computational Biology 15 (9): e1006453.
Dobin, Alexander, Carrie A Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R Gingeras. 2013. “STAR: Ultrafast Universal Rna-Seq Aligner.” Bioinformatics 29 (1): 15–21.
Goerner-Potvin, Patricia, and Guillaume Bourque. 2018. “Computational Tools to Unmask Transposable Elements.” Nature Reviews Genetics 19 (11): 688–704.
Guffanti, Guia, Andrew Bartlett, Torsten Klengel, Claudia Klengel, Richard Hunter, Gennadi Glinsky, and Fabio Macciardi. 2018. “Novel Bioinformatics Approach Identifies Transcriptional Profiles of Lineage-Specific Transposable Elements at Distinct Loci in the Human Dorsolateral Prefrontal Cortex.” Molecular Biology and Evolution 35 (10): 2435–53.
Jin, Ying, Oliver H Tam, Eric Paniagua, and Molly Hammell. 2015. “TEtranscripts: A Package for Including Transposable Elements in Differential Expression Analysis of Rna-Seq Datasets.” Bioinformatics 31 (22): 3593–9.
Li, Heng, and Richard Durbin. 2009. “Fast and Accurate Short Read Alignment with Burrows–Wheeler Transform.” Bioinformatics 25 (14): 1754–60.
Payer, Lindsay M, and Kathleen H Burns. 2019. “Transposable Elements in Human Genetic Disease.” Nature Reviews Genetics 20 (12): 760–72.
Tokuyama, Maria, Yong Kong, Eric Song, Teshika Jayewickreme, Insoo Kang, and Akiko Iwasaki. 2018. “ERVmap Analysis Reveals Genome-Wide Transcription of Human Endogenous Retroviruses.” Proceedings of the National Academy of Sciences 115 (50): 12565–72.