Detection and mitigation of spurious antisense expression with RoSA [version 1; peer review: 2 approved with reservations]

Antisense transcription is known to have a range of impacts on sense gene expression, including (but not limited to) impeding transcription initiation, disrupting post-transcriptional processes, and enhancing, slowing, or even preventing transcription of the sense gene. Strandspecific RNA-Seq protocols preserve the strand information of the original RNA in the data, and so can be used to identify where antisense transcription may be implicated in regulating gene expression. However, our analysis of 199 strand-specific RNA-Seq experiments reveals that spurious antisense reads are often present in these datasets at levels greater than 1% of sense gene expression levels. Furthermore, these levels can vary substantially even between replicates in the same experiment, potentially disrupting any downstream analysis, if the incorrectly assigned antisense counts dominate the set of genes with high antisense transcription levels. Currently, no tools exist to detect or correct for this spurious antisense signal. Our tool, RoSA (Removal of Spurious Antisense), detects the presence of high levels of spurious antisense read alignments in strand-specific RNA-Seq datasets. It uses incorrectly spliced reads on the antisense strand and/or ERCC spikeins (if present in the data) to calculate both global and gene-specific antisense correction factors. We demonstrate the utility of our tool to filter out spurious antisense transcript counts in an Arabidopsis thaliana RNASeq experiment. Availability: RoSA is open source software available under the GPL Open Peer Review


Introduction
Antisense RNAs are transcribed from the strand opposite to that of the sense transcript of either protein-coding or non-proteincoding genes. They appear to be widespread in all kingdoms of life and can play distinct roles in regulating gene expression or function. Typically, antisense RNAs are non-coding and expressed at lower levels than sense gene transcripts. However, they can exhibit a range of sizes, and may or may not have 5' cap or 3' poly(A) tails depending on whether they arise from either their own promoters, from divergent promoters, or from copying of sense transcripts by RNA-dependent RNA polymerases (see 1 and references therein, 2-4 ). In Arabidopsis thaliana, for example, the transcription of the Flowering Locus C (FLC) gene is known to be affected by transcription of antisense ncRNAs: COOLAIR 5,6 , a set of ncRNAs antisense to FLC, and COLDAIR 7 , antisense to COOLAIR. Both COLDAIR and COOLAIR are associated with different changes in sense strand gene expression at the FLC locus 8 . Antisense transcription is known to affect sense gene expression through multiple mechanisms 1 . During transcription, RNA polymerases may physically interfere with each other if both sense and antisense transcription take place simultaneously. Interference can prevent or slow down transcription (e.g. through RNA polymerase collisions 9,10 ) or force particular isoforms to be produced preferentially 11 . Post-transcriptionally, antisense transcripts can compete with sense transcripts for binding sites 12 . For example, the transcription of the human haemoglobin gene HBA1 is affected when the LUC7L gene on the opposite strand does not terminate, due to a deletion. It produces an antisense transcript that overlaps with HBA1, and which methylates the HBA1 promoter, repressing its expression 13 . In addition, since regions of protein coding genes on opposite DNA strands can overlap, their expression effectively generates transcripts that are, to varying extents, antisense to each other. Such overlapping gene pairs are a common feature of genome organization. We and others have shown that in some eukaryotic genomes tail-tail overlap enables the use of pre-mRNA 3' processing signals in different registers for genes coded on either strand 14 .
Incorporating antisense RNAs into genome annotation and properly quantifying their expression patterns is thus crucial, but remains challenging. Transcriptome-wide identification of RNAs is currently dominated by RNA-Seq. In this widespread experimental approach RNA is rarely sequenced directly, but instead is fragmented and first copied into cDNA and then copied again, so that libraries of DNA are sequenced. However, the copying of RNA by viral-derived reverse transcriptases is problematic. First, these polymerases exhibit DNA dependent polymerase activity, which can result in copies of the cDNA that can be incorrectly interpreted as antisense transcription. Second, just as reverse transcriptases switch template strand in viral biology, they can similarly switch templates in RNA-Seq library preparation, resulting again, in the interpretation of nonauthentic antisense RNAs 15-21 . Historically, in microarray and RT-PCR experiments, this step is known to assign some transcripts to the wrong strand, creating spurious antisense transcripts. Preparing samples with actinomycin D can help to reduce the number of spurious antisense transcripts 17 but can have unwanted side effects 20 . Alternative approaches to make strand-specific RNA-Seq libraries have been developed to mitigate artefacts arising from reverse transcription, however most of these also use reverse transcription 22 and so have similar problems with incorrect assignment. For example, the highly-rated 22,23 and widely used dUTP protocol for stranded RNA-Seq 24 is known to generate low levels of spurious antisense reads ranging from 0.6-3% of the sense signal 22,25,26 . Ultimately, the direct sequencing of full-length RNA molecules 27 will overcome many of the problems of distinguishing authentic antisense RNAs. However, currently, reverse-transcriptase based approaches dominate and the extent of spurious antisense RNAs identified in RNA-Seq datasets is rarely exposed.
In this paper, we analyse spurious antisense reads in 199 RNA-Seq experiments, across multiple organisms from both ENCODE 28 and our own work. Our results show that spurious antisense reads are often present in experiments, and can manifest at levels greater than 1% of sense transcript levels. Furthermore, the number of spurious antisense reads can vary substantially between replicates within the same experiment. In some cases, this variation may be sufficient to disrupt downstream analysis of antisense gene expression, by causing spurious antisense counts to dominate the set of genes with high antisense transcription levels.
To detect and correct for wrongly assigned reads we developed a tool, RoSA (Removal of Spurious Antisense), which calculates an antisense correction factor by identifying subsets of reads where all antisense reads are spurious. We evaluate the effect of using RoSA on Arabidopsis thaliana experimental data where varying levels of spurious antisense were present in different replicates. RoSA reduces the overall dependence of antisense counts on sense counts, a key indicator of the presence of spurious antisense. For individual genes with different real and spurious antisense characteristics, RoSA reduces spurious antisense counts while retaining the antisense signal.

Methods
As noted by Jiang et al. (2011,25 ), spurious antisense read counts can be estimated by analysing either ERCC spike-in data or counts of sense and antisense reads around splice sites. Each approach has different advantages: using spike-ins is simpler and faster, while using spliced reads allows a gene-by-gene estimate to be made. RoSA implements both approaches, in conjunction with pre-processing scripts to generate specialised read counts required by the tool. Once RoSA has an estimate of the levels of spurious antisense, it can adjust the raw antisense counts to account for the incorrectly stranded reads.

RoSA: Removal of Spurious Antisense
Our scripts and analysis code are bundled as a tool, RoSA (Removal of Spurious Antisense), available from the Barton Group's github pages at https://github.com/bartongroup/RoSA. RoSA is an R package supported by two python pre-processing scripts, callable from R.
For genes with spliced transcripts which are expressed in the data, RoSA uses the subset of reads from either strand that map across the splice junctions. The antisense reads in this subset are almost certainly spurious, and so RoSA can use the read counts to calculate a gene-specific antisense correction factor (Section 2.2). For genes without spliced transcripts, RoSA uses ERCC spike-in data, if present. Here any antisense read mappings are, by definition, spurious and the ratio of sense to antisense reads mapping to the spike-ins thus provides a global, rather than gene-specific, antisense correction factor (Section 2.3). If ERCC spike-in data is not available, RoSA instead calculates a global estimate of the spurious antisense fraction from the set of spliced reads. Counting all, or spliced-only, antisense reads is not directly supported by existing tools. RoSA's preprocessing scripts perform these functions. The make_annotation script creates an antisense annotation (as gtf) from a standard annotation (as gff or gtf), which can then be used to generate antisense read counts via a standard read counting tool (Section 2.4.1). RoSA doesn't specify how the sense and antisense gene expression is counted leaving users free to apply whichever tool they feel will best represent the gene expression in their experiment. However, the accuracy of the corrections calculated by RoSA will be affected by this choice in the same way as the calculation of differential gene expression. If counting methods are used that only consider regions within sense features that do not overlap any antisense feature, the genespecific corrections calculated by RoSA may be less accurate where the overlap is large and/or the sense or antisense expression is low.
RoSA then adjusts these raw counts to produce corrected antisense counts (Section 2.4). The count_spliced script generates sense and antisense counts of reads at splice junctions, used when estimating spurious antisense from spliced reads. The script takes a standard annotation (as gtf/gff) and corresponding alignment (as bam) and outputs counts of spliced sense and antisense reads to a designated output file.
RoSA takes several datasets containing different read counts as its input, for each replicate: RoSA calculates and returns antisense:sense ratios for the spikein data, or spliced read data, or both, and, for each gene and replicate, outputs new read counts values corrected for spurious antisense. RoSA also plots antisense versus sense counts of the original and corrected data, by replicate.

Using spliced reads
RoSA's main approach to estimating spurious antisense is to use spliced reads within the main dataset. Reads which map antisense to a multi-exon gene, and that also show the same splicing pattern as spliced sense-mapping reads are almost certainly spurious, as the splicing motif (canonically GU-AG) will be incorrect on the opposite strand ( Figure 1). An estimate of spurious antisense can be calculated by considering only spliced reads whose splices match annotated splice sites (splice-matched reads), and, as with the spike-ins, calculating the ratio of antisense to sense reads.
Splice-matched reads are identified by first filtering all spliced reads in the data. In a bam file of aligned reads, spliced reads have a CIGAR string containing 'N', indicating a skipped region. SAM processing tools such as sambamba 29 support filtering on the CIGAR string and can extract spliced reads rapidly. A second filtering step pulls out only those reads whose splice locations match at least one intron in the annotation, by processing each read in turn, identifying the spliced positions (based on the read location and the CIGAR string) and checking the annotation for a matching intron. Finally, the strand of each spliced read can be determined from its flag field value 30 , and compared to the strand of the matching intron(s). Reads on the same strand as the intron(s) are counted as sense reads, and remaining reads as antisense reads. Since spurious antisense reads are misallocated sense reads, the number of antisense splice-matched reads assigned to a gene is strongly positively correlated with the number of its sense splice-matched reads (see Section 3). The ratio of antisense:sense counts on the splicematched reads thus gives a simple global estimate of the level of spurious antisense across the whole dataset. Using spliced reads has the advantage that an antisense:sense ratio can be calculated on a gene-by-gene basis, for any spliced gene. Genes without any spliced reads fall back on the global estimate, calculated either from the spike-ins (see Section 2.3) or the spliced reads.
In the case of real, unannotated, antisense expression at a gene locus, the behaviour of RoSA falls into three categories: 1. If the splicing of the true antisense transcript differs from the sense transcript (including no splicing) then RoSAs gene specific correction will remove any spurious antisense expression (identified by antisense matches to the sense splicing) and leave the true antisense expression unchanged.
2. If the splicing of the antisense expression is the same as the sense strand, then RoSA will remove this completely.
3. If the true antisense splicing is the same as the sense strand in some parts of the transcript, but not across the entire transcript, then RoSA will remove a fraction of the true antisense expression depending on how similar the splicing patters are.
We anticipate that occurrences of 2 & 3 will be uncommon in RNA-seq datasets. Point 2 highlights a minor potential limitation of the gene-specific splicing-based corrections calculated by RoSA, namely that it cannot distinguish between spurious antisense signal and potential biological transcription from RNA-dependent RNA polymerases (RdRPs). Although RdRPs are widespread in eukaryotic genomes, only 8-30% of eukaryotic gene regions have significant length ORFs on their opposite strands 31 , providing an upper limit on the potential impact of this method of transcription on the RNA complement within a cell. Eukaryotic RdRPs evolved independently from their viral counterparts and, in plants, are involved in siRNA transcriptional silencing 33 . This is not the case in animals however (except in C. elegans) where their function remains elusive 34 .

Using ERCC spike-ins
An alternative approach to estimating spurious antisense is to use ERCC spike-in data. Developed by the External RNA Control Consortium (ERCC) 35 , the ERCC spike-in controls are synthetic RNA transcripts that are added to RNA-Seq experiments to act as controls 36 . The 92 spike-ins are designed to mimic a range of eukaryotic mRNA characteristics, varying in length, GC-content and concentration, with a 20bp poly-A tail. They have minimal sequence similarity with known eukaryotic transcripts. Since the spike-ins are synthetic, they are unidirectional, and so any reads assigned as antisense to a spike-in can be assumed to be spurious. As the spike-ins are present at a wide range of concentrations, they are detected with a wide range of read counts, permitting an estimate of the ratio of antisense to sense read counts on the spike-ins to be calculated, which can then be used to estimate the contribution of spurious antisense transcripts across the full dataset. Obtaining sense and antisense counts for the spike-ins is straightforward. First we align the reads to the spike-ins (using the spike-in annotation file ERCC92.gtf, available at https://www.thermofisher.com/order/catalog/product/4456739) and then count reads, using a strand-aware read counting tool such as featureCounts 37 , HT-SeqCount 38 , etc. Now averaging the spurious antisense:sense ratio across all of the spike-ins gives a global estimate of the spurious antisense, in just the same way as for the spliced reads.

Mitigating spurious antisense
Having identified high or differing levels of spurious antisense in an RNA-Seq experiment, we also want to correct for the incorrectly assigned reads so that true differential expression calling can be performed. The ratio of spurious antisense:sense read counts can be used as a simple correction factor. Defining r as the ratio of spurious antisense:sense and S and A respectively as the number of sense and antisense counts for a gene, the number of spurious antisense read counts A S is estimated for each gene as: A S = r.S .
Then the antisense count can be corrected to account for the spurious antisense by taking A -A S . This correction simply adjusts read counts for each gene, and does not identify specific reads as incorrectly assigned, so pile-ups cannot be adjusted.
Since the spurious antisense reads are mis-assigned sense reads, RoSA then adds the spurious antisense count for each gene to its sense count.

Counting antisense reads.
In order to apply the antisense correction factor, counts of antisense reads for each gene are required. Counting antisense reads is not directly supported by read counting tools. However, it can be performed with feature-Counts 37 by setting parameters to indicate that reads are stranded in the opposite direction to which they are. Unfortunately, if there are overlapping genes then reads in the overlaps will be counted twice using this tactic. As reads in regions of gene overlap are necessarily ambiguous, they cannot be considered to be antisense, spurious or otherwise. RoSA avoids this issue by building a custom antisense annotation based on the input sense annotation but excluding regions where genes on opposite strands overlap. Different gene transcripts are accounted for by merging all transcripts for a gene into a single maximal transcript. Whenever exons of different transcripts overlap in the annotation, the exon in the maximal transcript is the maximum extent of both exons. Given a maximal transcript, the script creates an antisense feature on the opposite strand which runs for the full extent of the maximal transcript. If the maximal transcript of another gene overlaps with the antisense feature, then the antisense feature is truncated to avoid overlapping. Once an antisense annotation is available, a read counting tool can be used to count antisense reads, by providing the antisense annotation instead of the standard annotation.

Arabidopsis thaliana datasets with spurious antisense
A procedure to experimentally generate RNA-Seq data with specific levels of spurious antisense is not known. Our main experimental data (Experiment 1) is therefore drawn from the study which originally motivated our investigation into incorrectly assigned antisense reads. In this study, spurious antisense occurred by chance at varying orders of magnitude across different replicates. Additionally, we perform a meta-analysis using three other Arabidopsis thaliana datasets (Experiments 2-4, 39 ) and data from ENCODE (see Underlying data for the full list of the ENA and ENCODE accessions).
The read data were also aligned to the ERCC spike-ins annotation, using the same parameters. Read counts per gene were then quantified from these alignments with featureCounts v1.5.0-p1 using the publicly available TAIR10 annotation with the parameters: -s 2 -p -t exon --largestOverlap. After running RoSA's make_annotation script to build an antisense annotation, antisense read counts per gene were quantified in the same way, with the parameters: -s 2 -p -t antisense --largestOverlap. Finally, spliced sense and antisense reads were counted using RoSA's count_spliced script with the TAIR10 annotation.

Operation
A full description of RoSAs environment, dependencies, installation and basic operation can be found on the RoSA GitHub repository. Briefly, RoSA is a combination of an R package and python scripts for data preprocessing. Minimal system requirements for the package are R v3.5+, python 2 v2.7+ the LSD R package and the python packages scipy (v0.16.1 -0.17.1), numpy, pandas (not v0.20.1), six and, optionally, drmaa for cluster integration. The python scripts to find and count spliced antisense and sense reads also depends on sambamba. To facilitate ease-of-use, a conda environment that captures all the relevant dependencies is included as part of the RoSA codebase. RoSA's python scripts are provided as a python package ad are installed via pip, while the R package can be installed directly from within R using the devtools package.
RoSA operates on the total and spliced read counts from sense and antisense bam format read alignments of stranded RNA-Seq datasets, either with or without ERCC spike-in standards.
To facilitate easy generation of this read count data, RoSA includes helper pre-processing scripts to generate the antisense counterpart of the provided gtf/gff format sense-strand genome annotations (make_annotation), and to generate spliced-read gene count data from the bam format read alignments using both the sense-and anti-sense annotations (count_spliced). Both of these helper scripts can be called directly within R as part of the RoSA R package, Detailed help for the R RoSA functionality can be accessed within R with the command, help(rosa).

Results
We used RoSA to analyze our data from Experiment 1 for spurious antisense, using both the spike-in and spliced reads counts. RoSA calculated antisense:sense ratios for the spike-ins (Figure 2) showing that the 3 replicates have antisense:sense ratios on the spike-ins of 0.0008, 0.004 and 0.011. Although these ratios are small, if the replicates were being compared for differential expression, the differences are potentially substantial for highly expressed genes, and could lead to differential antisense expression being called erroneously.
For each replicate we calculated the spurious antisense: sense ratios for the spliced reads with RoSA, and compared them to the spike-ins. An overview of the results for all three replicates shows that the spurious antisense levels calculated from the spike-ins are in good agreement with the levels calculated from the spliced reads ( Figure 3 and Figure 4, Row 1).   Each column presents data for one replicate. Row 1: Antisense:sense ratios calculated from spike-ins (Black points & fit line) and spliced reads for each gene (density heatmap). The antisense:sense ratios for both the spike-ins and spliced reads are in good agreement. The strong correlation between the sense and antisense spliced counts, and the constant antisense:sense ratio across all genes, indicates that the majority of the antisense expression in the data is not a sequence-, or gene-specific, phenomenon. Rather, this is what would be expected from a systematic process affecting a constant fraction of the sequenced reads. Row 2: Antisense:sense ratios calculated for the full gene counts (spliced & unspliced). The correlation between the sense and antisense expression persists, however it is weaker than the correlation using just the spliced and spike-in sense and antisense expression. This reflects the inclusion of true biological antisense expression, unspliced genes where a global correction is less accurate, and low expression genes where the splicing correction is not well measured. Row 3: Corrected antisense:sense ratios calculated for the full gene counts (spliced & unspliced). The corrected antisense counts show much weaker correlation with the corresponding gene counts reflecting the removal of the systematic spurious antisense count signal. On all plots the dashed line marks y=x; points above this line correspond to genes where the antisense:sense ratio is > 1.
Finally, RoSA calculated a spurious antisense correction across the whole of each replicate. Every spliced gene was corrected with the antisense:sense ratio specific to the gene, and unspliced genes were corrected using the mean ratio calculated from the spike-ins. (RoSA also allows the unspliced correction to be calculated from the mean spliced reads ratio, for datasets without ERCC spike-ins). Overall, RoSA reduces the correlation between antisense and sense counts in the data (Figure 4, Rows 2 & 3), as would be expected with a reduction in incorrectly assigned reads. Two examples of corrections made by RoSA are shown in Figure 5, where the antisense signal appears to be almost entirely spurious, RoSA's correction factor reduces the antisense counts substantially, but where there also appears to be some real antisense signal, RoSA's correction factor leaves a higher proportion of counts.
As well as identifying instances of antisense expression, looking at antisense counts in this way can also be useful in identifying misannotated genes. For example, in our data there are many genes where the antisense:sense ratio is more than 1 (e.g. see points lying above x=y in Figure 4, Row 2), which may indicate an incorrect strand assignment in the annotation.

Comparing antisense:sense ratios
Calculating antisense:sense ratios allows comparisons of spurious antisense to be made between replicates and between experimental condition, and can reveal whether there are systematic differences which might confound experimental comparisons. For example, Figure 1 presents results from an RNA-Seq experiment where spurious antisense levels differed by an order of magnitude between replicates. In this experiment, the WT replicates had spurious antisense:sense ratios of 0.0031 (SD 0.0116), 0.0009 (SD 0.0070) and 0.0111 (SD 0.031).
To determine the extent of this problem for RNA-Seq datasets in general, we investigated the spurious antisense levels across a range of experiments and research groups. We analysed antisense reads assigned to the spike-ins from three other experiments in our lab (Experiments 2-4), as well as 195 publicly available human datasets from the ENCODE project that included the ERCC spike-ins 28 (see Underlying data for details of the sense, antisense and RoSA-corrected antisense expression for all A. thaliana genes in the datasets from Experiments 1-4). A separate antisense:sense ratio was calculated for each replicate in each experiment (Figure 6), showing that spurious antisense reads are present at varying levels and can range across several orders of magnitude. This presents a serious quality control issue for anyone investigating differential antisense expression: a real difference in antisense expression could be completely masked by a difference in spurious antisense.

Conclusions
Spurious antisense is common in strand-specific RNA-Seq datasets and can occur at varying levels across replicates in the same experiment. Differing levels of such incorrectly assigned reads are enough to disrupt differential expression analyses of antisense gene expression.
We have developed a new tool, RoSA, which can detect, quantify and correct for spurious antisense. RoSA provides an important quality control step for researchers analyzing antisense expression in their data.  sense reads, RoSA then adds the spurious antisense count for each gene to its sense count'. Adding the spurious antisense counts to the sense counts assumes that reverse transcriptase errors happen at the cost of 'real' RNA -> cDNA processing. Is this correct? Software:

Data availability
The package as it stands is not user friendly and has several bugs that need to be resolved before the software can be used smoothly.
Solving environment: failed ResolvePackageNotFound: -r-base==3.5.1=h4fe35fd_0 -libssh2==1.8.0=h5b517e9_2 -matplotlib==2.2.3=py27h8e2386c_0 To try and fix this problems I removed the pinning on the exact build, leaving only the pinning on version number. This returned the following error: Solving environment: failed UnsatisfiableError: The following specifications were found to be in conflict: -matplotlib=2.2.3 -sip=4.18.1 -tk=8.6.8 I was able to install by manually creating an environment using the dependencies listed in the readme, although using the pip install for rosa downgraded the pandas install from 0.24 to 0.19.
The ability to install all dependencies as a conda environment is very useful to the average user, and the authors should consider either fixing the supplied env file, or (ideally) creating a proper rosa conda package, so that rosa and all its dependencies (both R and python) can be installed with conda install rosa.
Once installed I found the documentation to be a little on the sparse side, and it times confusing. For example, the help for the `rosa` function states that the parameter "data" should be a list of dataframes with a particular list of columns. But the help for the "group" parameter says that it is a character vector specifying the replicates for each column in data. Thus it is unclear if replicates should be provided as a member of a list of data.frames, or as separate columns in one data.frame.

○
There are required arguments for matrices of spike-in counts, as well as an argument to provide the positions of spike-ins in the "data" matrix. It is unclear whether spike-ins should be provided using only one or both of these mechanisms, or what to do if there is no spikein data.

○
The help states "This function was written with the intention of obtaining count and group parameters from an edgeR DGEList object d, via d$counts and d$samples$group." But d$counts would not return a data.frame of the format specified in the help for the "data" argument.

○
The manuscript doesn't really mention TPMs, but the help for rosa states that TPMs are required to run the function. I believe it would not be clear to a less experienced user how to calculate a TPM for the anti-sense counts generated on the custom anti-sense annotation.
○ Not help page is provided for either the `count_spliced` nor `make_annotation` functions. Nor is much in the way of documentation provided for the python scripts these functions ○ call.
The package has no associated vignette. I think it would be of great use to users, particularly less experienced ones, to provide a walk-through of the analysis of a dataset from start to finish, with a simplified example dataset.

○
Overall, I think that this work is both novel, interesting, useful and basically sound. Before I can recommend approval, I would need to see minimally that the issues identified with installation and documentation have been dealt with. I would also be interested in the authors replies to my comments on the manuscript itself.