F1000Research F1000Research 2046-1402 F1000 Research Limited London, UK 10.12688/f1000research.18952.1 Software Tool Article Articles Detection and mitigation of spurious antisense expression with RoSA

[version 1; peer review: 2 approved with reservations]

 Mourão Kira Conceptualization Data Curation Formal Analysis Investigation Methodology Software Validation Visualization Writing – Original Draft Preparation Writing – Review & Editing 1 Schurch Nicholas J. Conceptualization Investigation Methodology Project Administration Software Supervision Writing – Review & Editing 2 Lucoszek Radek Data Curation Investigation Writing – Review & Editing 3 Froussios Kimon Conceptualization Methodology Software Writing – Review & Editing 4 MacKinnon Katarzyna Data Curation Investigation Writing – Review & Editing 5 Duc Céline Data Curation Writing – Review & Editing 6 Simpson Gordon Funding Acquisition Writing – Review & Editing a 3 Barton Geoffrey J. Funding Acquisition Project Administration Resources Supervision Writing – Review & Editing https://orcid.org/0000-0002-9014-5355 b 7 Synpromics Ltd, Edinburgh, Midlothian, EH25 9RG, UK Biomathematics and Statistics Scotland, James Hutton Institute, Aberdeen, Scotland, AB15 8QH, UK Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, Dundee, Scotland, DD1 5EH, UK Research Institute of Molecular Pathology, Vienna, 1030, Austria Cell & Molecular Sciences, James Hutton Institute, Invergowie, Dundee, Scotland, DD2 5DA, UK Équipe Épigénétique, Unité Fonctionnalité et Ingénierie des Protéines (UFIP) Faculté des Sciences et Techniques, Université de Nantes, NANTES, 92208 F44322 CEDEX 3, France Computational Biology, School of Life Sciences, University of Dundee, Dundee, Scotland, DD1 5EH, UK G.G.Simpson@dundee.ac.uk g.j.barton@dundee.ac.uk

*To whom correspondence should be addressed

No competing interests were disclosed.

 7 6 2019 2019 8 819 30 5 2019 Copyright: © 2019 Mourão K et al. 2019 This is an open access article distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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. Strand-specific 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 spike-ins (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 RNA-Seq experiment.

Availability: RoSA is open source software available under the GPL licence via the Barton Group GitHub page https://github.com/bartongroup.

 RNA-seq antisense expression gene expression Arabidopsis thaliana ENCODE Biotechnology and Biological Sciences Research Council BB/M004155/1 BB/M010066/1 This work has been supported by the Biotechnology and Biological Sciences Research Council [BB/M004155/1, BB/M010066/1] to G.J.B. and G.G.S. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 1. Introduction

Antisense RNAs are transcribed from the strand opposite to that of the sense transcript of either protein-coding or non- protein-coding 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 non-authentic 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.

 2. 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.

2.1 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 pre-processing 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 gene-specific 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:

Full read counts by gene

Antisense counts by gene (via the make_annotation script)

At least one of:

Spike-in sense and antisense counts

Spliced sense and antisense counts (via RoSA's count_spliced script)

RoSA calculates and returns antisense:sense ratios for the spike-in 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.

 2.2 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.

An example of spurious antisense reads displaying the same splicing structure as the sense strand.

The reverse strand gene AT4G18970 is strongly expressed in all 3 replicates (bottom tracks). Spurious antisense can also be seen in all replicates (top tracks), with splice points in the antisense signal matching splice points in the sense signal. Furthermore, the level of spurious antisense varies noticeably between replicates. (Figure generated by IGB 32 ).

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 splice-matched 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:

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.

If the splicing of the antisense expression is the same as the sense strand, then RoSA will remove this completely.

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 anti-sense 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 .

 2.3 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.

 2.4 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.

2.4.1 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 featureCounts 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.

 2.5 <italic toggle="yes">Arabidopsis thaliana</italic> 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).

2.5.1 Arabidopsis sample preparation and sequencing. Briefly, the RNA-Seq data for Experiments 1 is wild-type (WT) Arabidopsis thaliana Colombia-0 (Col-0) biological replicates. WT A. thaliana Col-0 seeds were sown aseptically on MS10 plates. The seeds were stratified for 2 days at 4°C and then grown at a constant 21°C under a 16-h light/8-h dark cycle for a further 14 days, at the end of which the seedlings were harvested. Total RNA was isolated from the seedlings with the RNeasy Plant Mini Kit (Qiagen). In Experiment 1, DNAse digestion was performed on column, as a part of RNA isolation, and 8 μl of ERCC spike-ins (External RNA Controls Consortium 2005) at a 1:100 dilution was added to 4 μg of total RNA. Libraries were prepared according to the TruSeq ® Stranded mRNA Sample Preparation Guide Rev E. The libraries were sequenced on a HiSeq2000 at the Genomic Sequencing Unit of the University of Dundee. This preparation largely mirror the sample preparation of the datasets take from Froussios et al. (2017, 39 , Experiments 2–4) In Experiments 2–4, however, the sequencing libraries were prepared using the Illumina TruSeq ® Stranded Total RNA with Ribo-Zero Plant kit.

Experiment 1 has 3 replicates, processed as one batch, with a total of 4 × 10 8 150-bp paired-end reads. Experiments 2 and 3 have 7 biological WT replicates, while Experiment 4 has 3, for a total of 17 biological WT replicates and ~1.7 × 10 9 100-bp paired-end reads across the three experiments. The same lab sowed, grew and harvested the plants, and prepared the libraries. The sequencing was performed on the same machine by the same people at the same sequencing facility and all the samples include the ERCC spike-ins which can verify the WT samples are consistent and comparable across experiments.

2.5.2 Quality control, alignment and quantification. The quality of the data was quantified using FastQC v0.11.2 40 with all the replicates performing as expected for high quality RNA-Seq data with good median per-base quality (≥28) across >90% of the read length. The read data for all experiments were aligned to the TAIR10 41 Arabidopsis thaliana genome using the splice-aware aligner STAR v2.4.2a 42 for Experiment 1 and STAR v2.5.0 for Experiments 2–4. The index was built with --sjdbOverhang 149 (Experiment 1) or –-sjdbOverhang 99 (Experiments 2–4) and the alignment was run with parameters: --outSJfilterIntronMaxVsReadN 5000 10000 15000 --outSAMAttributes All --outFilterMultimapNmax 2 --outFilterMismatchNmax 5 --outFilterType BySJout.

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.

 2.6 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).

 3 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.

Antisense versus sense counts for the ERCC spike-ins for each replicate in Experiment 1.

Points represent antisense and sense read counts for individual spike-ins. Each line is the average antisense:sense ratio for one replicate. Here, antisense:sense ratios vary by an order of magnitude across the 3 replicates, with values of 0.004 (WT1), 0.0008 (WT2), and 0.011 (WT3).

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).

Comparison of antisense:sense ratios calculated from spliced reads or spike-ins, by replicate.

Ratios estimated from spike-ins show good agreement with ratios estimated from spliced reads. Outliers not shown.

 Normalised antisense versus sense counts by replicate.

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.

Two genes with differing antisense expression profiles and the read count corrections proposed by RoSA.

The reverse strand gene AT4G18970 (left) has antisense expression which clearly matches the splice sites of the sense strand. RoSA eliminates almost all of the antisense reads. The forward strand gene AT5G66570 (right) has both antisense expression matching the sense strand splice sites, and a peak at the 5’ end which is unlikely to have resulted from incorrect read assignment. RoSA only reduces the antisense counts by around 40%. (Figures generated by IGB 32 ).

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.

3.1 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.

Spurious antisense:sense ratios for spike-ins, by research group.

Data are from either from ENCODE (Gingeras, Graveley and Lecuyer) or our own group (Simpson). Each point represents the ratio for a replicate. Ratios range from 0.0111 to 0.00003.

 4 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.

 Data availability Underlying data

Arabidopsis col-0 WT strand-specific RNA-Seq data from poly-A pulldown, Accession number E-MTAB-7990: https://identifiers.org/arrayexpress/E-MTAB-7990

RNA-seq data of wild type Arabidopsis seedlings, Accession number E-MTAB-5446: https://identifiers.org/arrayexpress/E-MTAB-5446

Extended data

Zenodo: bartongroup/RoSA: Initial, http://doi.org/10.5281/zenodo.2661378 43 .

This project contains the following extended data:

Accession numbers for ENCODE data: https://github.com/bartongroup/RoSA/tree/master/extras/F1000_manuscript/RoSA_Extended_Data.docx

Accession details for ENCODE data: https://github.com/bartongroup/RoSA/blob/master/extras/F1000_manuscript/ENCODE_accessions.xlsx

Arabidopsis seedlings RNA-seq read count expression counts: https://github.com/bartongroup/RoSA/tree/master/extras/F1000_manuscript/expression_data.csv

License: GNU General Public License 3.0.

 Software availability

Source code available from: https://github.com/bartongroup/RoSA

Archived source code as at time of publication: https://doi.org/10.5281/zenodo.2661378 43 .

License: GNU General Public License 3.0.

 Acknowledgments

A previous version of this article is available on BioRxiv: https://doi.org/10.1101/425900.

 Pelechano V Steinmetz LM : Gene regulation by antisense transcription. Nat Rev Genet. 2013;14(12):880–893. 24217315 10.1038/nrg3594 Matsui A Iida K Tanaka M : Novel Stress-Inducible Antisense RNAs of Protein-Coding Loci Are Synthesized by RNA-Dependent RNA Polymerase. Plant Physiol. 2017;175(1):457–472. 28710133 10.1104/pp.17.00787 5580770 Lin S Zhang L Luo W : Characteristics of Antisense Transcript Promoters and the Regulation of Their Activity. Int J Mol Sci. 2015;17(1): pii: E9. 26703594 10.3390/ijms17010009 4730256 Chan WY Wu SM Ruszczyk L : The complexity of antisense transcription revealed by the study of developing male germ cells. Genomics. 2006;87(6):681–92. 16458478 10.1016/j.ygeno.2005.12.006 Swiezewski S Liu F Magusin A : Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target. Nature. 2009;462(7274):799–802. 20010688 10.1038/nature08618 Liu F Marquardt S Lister C : Targeted 3' processing of antisense transcripts triggers Arabidopsis FLC chromatin silencing. Science. 2010;327(5961):94–97. 19965720 10.1126/science.1180278 Heo JB Sung S : Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science. 2011;331(6013): 76–79. 21127216 10.1126/science.1197349 Ietswaart R Wu Z Dean C : Flowering time control: another window to the connection between antisense RNA and chromatin. Trends Genet. 2012;28(9):445–453. 22785023 10.1016/j.tig.2012.06.002 Hobson DJ Wei W Steinmetz LM : RNA polymerase II collision interrupts convergent transcription. Mol Cell. 2012;48(3):365–374. 23041286 10.1016/j.molcel.2012.08.027 3504299 Crampton N Bonass WA Kirkham J : Collision events between RNA polymerases in convergent transcription studied by atomic force microscopy. Nucleic Acids Res. 2006;34(19):5416–5425. 17012275 10.1093/nar/gkl668 1636470 Onodera CS Underwood JG Katzman S : Gene isoform specificity through enhancer-associated antisense transcription. PLoS One. 2012;7(8):e43511. 22937057 10.1371/journal.pone.0043511 3427357 Kawano M Aravind L Storz G : An antisense RNA controls synthesis of an SOS-induced toxin evolved from an antitoxin. Mol Microbiol. 2007;64(3):738–754. 17462020 10.1111/j.1365-2958.2007.05688.x 1891008 Tufarelli C Stanley JA Garrick D : Transcription of antisense RNA leading to gene silencing and methylation as a novel cause of human genetic disease. Nat Genet. 2003;34(2):157–165. 12730694 10.1038/ng1157 Sherstnev A Duc C Cole C : Direct sequencing of Arabidopsis thaliana RNA reveals patterns of cleavage and polyadenylation. Nat Struct Mol Biol. 2012;19(8):845–52. 22820990 10.1038/nsmb.2345 3533403 Beiter T Reich E Weigert C : Sense or antisense? False priming reverse transcription controls are required for determining sequence orientation by reverse transcription-PCR. Anal Biochem. 2007;369(2):258–261. 17698025 10.1016/j.ab.2007.06.044 Haddad F Qin AX Giger JM : Potential pitfalls in the accuracy of analysis of natural sense-antisense RNA pairs by reverse transcription-PCR. BMC Biotechnol. 2007;7:21. 17480233 10.1186/1472-6750-7-21 1876213 Perocchi F Xu Z Clauder-Münster S : Antisense artifacts in transcriptome microarray experiments are resolved by actinomycin D. Nucleic Acids Res. 2007;35(19):e128. 17897965 10.1093/nar/gkm683 2095812 Thomason MK Storz G : Bacterial antisense RNAs: how many are there, and what are they doing? Annu Rev Genet. 2010;44:167–88. 20707673 10.1146/annurev-genet-102209-163523 3030471 Tzadok S Caspin Y Hachmo Y : Directionality of noncoding human RNAs: how to avoid artifacts. Anal Biochem. 2013;439(1):23–29. 23583907 10.1016/j.ab.2013.03.031 Yu WH Høvik H Olsen I : Strand-specific transcriptome profiling with directly labeled RNA on genomic tiling microarrays. BMC Mol Biol. 2011;12(1):3. 21235785 10.1186/1471-2199-12-3 3031212 Houseley J Tollervey D : Apparent non-canonical trans-splicing is generated by reverse transcriptase in vitro. PLoS One. 2010;5(8):e12271. 20805885 10.1371/journal.pone.0012271 2923612 van Dijk E Jaszczyszyn Y Thermes C : Library preparation methods for next-generation sequencing: tone down the bias. Exp Cell Res. 2014;322(1):12–20. 24440557 10.1016/j.yexcr.2014.01.008 Levin JZ Yassour M Adiconis X : Comprehensive comparative analysis of strand-specific RNA sequencing methods. Nat Methods. 2010;7(9):709–15. 20711195 10.1038/nmeth.1491 3005310 Parkhomchuk D Borodina T Amstislavskiy V : Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Res. 2009;37(18):e123. 19620212 10.1093/nar/gkp596 2764448 Jiang L Schlesinger F Davis CA : Synthetic spike-in standards for RNA-seq experiments. Genome Res. 2011;21(9):1543–1551. 21816910 10.1101/gr.121095.111 3166838 Zeng W Mortazavi A : Technical considerations for functional sequencing assays. Nat Immunol. 2012;13(9):802–807. 22910383 10.1038/ni.2407 4138050 Garalde DR Snell EA Jachimowicz D : Highly parallel direct RNA sequencing on an array of nanopores. Nat Methods. 2018;15(3):201–206. 29334379 10.1038/nmeth.4577 ENCODE Project Consortium: An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57–74. 22955616 10.1038/nature11247 3439153 Tarasov A Vilella AJ Cuppen E : Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032–2034. 25697820 10.1093/bioinformatics/btv098 4765878 The SAM/BAM Format Specification Working Group: Sequence Alignment/Map Format Specification.2017. Reference Source Winters-Hilt S : RNA-Dependent RNA Polymerase encoding Artifacts in Eukaryotic Transcriptomes. Int J Mol Genet Gene Ther. 2017;2(1). 10.16966/2471-4968.108 Freese NH Norris DC Loraine AE : Integrated genome browser: visual analytics platform for genomics. Bioinformatics. 2016;32(14):2089–2095. 27153568 10.1093/bioinformatics/btw069 4937187 Iyer LM Koonin EV Aravind L : Evolutionary connection between the catalytic subunits of DNA-dependent RNA polymerases and eukaryotic RNA-dependent RNA polymerases and the origin of RNA polymerases. BMC Struct Biol. 2003;3(1):1. 12553882 10.1186/1472-6807-3-1 151600 Pinzón N Bertrand S Subirana L : Functional lability of RNA-dependent RNA polymerases in animals. bioRxiv. 2018. 10.1101/339820 Baker SC Bauer SR Beyer RP : The External RNA Controls Consortium: a progress report. Nat Methods. 2005;2(10):731–734. 16179916 10.1038/nmeth1005-731 ERCC: NIST standard reference material 2374.2017. Reference Source Liao Y Smyth GK Shi W : featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–930. 24227677 10.1093/bioinformatics/btt656 Anders S Pyl PT Huber W : HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–169. 25260700 10.1093/bioinformatics/btu638 4287950 Froussios K Schurch NJ Mackinnon K : How well do RNA-Seq differential gene expression tools perform in a eukaryote with a complex transcriptome? bioRxiv. 2017. 10.1101/090753 Andrews S : FastQC: A quality control tool for high throughput sequence data.2010. Reference Source Arabidopsis Genome Initiative: Analysis of the genome sequence of the flowering plant Arabidopsis thaliana. Nature. 2000;408(6814):796–815. 11130711 10.1038/35048692 Dobin A Davis CA Schlesinger F : STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. 23104886 10.1093/bioinformatics/bts635 3530905 Schurch N : bartongroup/RoSA: Initial (Version v1.0). Zenodo. 2019. http://www.doi.org/10.5281/zenodo.2661378 10.5256/f1000research.20776.r52108 Reviewer response for version 1 Werner Andreas 2 Referee https://orcid.org/0000-0002-2274-0662 Casement John 1 Co-referee Bioinformatics Support Unit, Faculty of Medical Sciences, Newcastle University, Newcastle, UK RNA Interest Group, Institute for Cell and Molecular Biosciences, Newcastle University, Newcastle upon Tyne, UK

Competing interests: No competing interests were disclosed.

 28 8 2019 Copyright: © 2019 Casement J and Werner A 2019 This is an open access peer review report distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. recommendation approve-with-reservations

The authors address an important question, i.e. detection of real and experimentally generated antisense transcripts in strand-specific RNAseq datasets. With their software tool RoSA, they present a valid strategy to solve this problem. However, there are a number of concerns regarding the manuscript and the software.

Manuscript:

Introduction: The authors discuss the mechanisms of how antisense transcripts can regulate/interfere with the expression of the sense gene. The sentence ‘Post-transcriptionally, antisense transcripts can compete with sense transcripts for binding sites’, may be applicable in bacteria, but is misleading when studying eukaryotic systems. Post-transcriptional mechanisms involve RNAi, RNA editing or RNA-protein complexes but hardly competition for binding sites with the sense transcript. Moreover, the example given in the text (research by Tufarelli et al. 1 ) is a prime example of transcriptional interference, not of post-transcriptional regulation.

Introduction: ‘viral-dependent‘ should read ‘virus-dependent’.

2.2 Using spliced reads: The fact that RoSA cannot distinguish between spurious antisense reads and RdRP generated transcripts are only relevant if datasets from plants and C.elegans are interrogated (that express significant levels of RdRP). The issue becomes critical if particular loci affected by RdRP amplification contribute disproportionately to a genome wide correction factor. In line with reviewer 1, it is unclear how an open reading frame affects RNA dependent RNA polymerization.

2.4 Mitigating spurious antisense: ‘Since the spurious antisense reads are miss-assigned 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:

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.

Installation of the Python package using the pip installer worked without difficulties, as did the R package installation.

Requesting help at the command-line (via the -h option) for the Python scripts 'make_annotation' and 'count_spliced' had the side-effect of creating an empty log file each time.

Running the 'count_spliced' script over a modestly sized alignment file (around 0.5 GB) resulted in over 2000 log files being produced, each detailing the locations at which spliced reads were being counted. This may be related to running the script without drmaa (in serial mode), as a message to that effect was repeatedly printed to the terminal during operation.

The Python scripts can be run as R functions, but the help documentation is not available within R.

A detailed R package vignette including sample data sets and examples of code would be very helpful, especially for the less experienced users.

General comment: Is there a reason why the authors used a storage costly strategy aligning reads followed by counting rather than using alignment free software such as Kallisto or Salmon?

Are the conclusions about the tool and its performance adequately supported by the findings presented in the article?

Partly

Is the rationale for developing the new software tool clearly explained?

Yes

Is the description of the software tool technically sound?

Yes

Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?

Partly

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool?

Partly

Reviewer Expertise:

Andreas Werner: Molecular biology, natural antisense transcripts. John Casement: Bioinformatics, RNAseq analysis.

We confirm that we have read this submission and believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however we have significant reservations, as outlined above.

References : Transcription of antisense RNA leading to gene silencing and methylation as a novel cause of human genetic disease. Nat Genet .2003;34(2) : 10.1038/ng1157 157-65 12730694 10.1038/ng1157 10.5256/f1000research.20776.r49670 Reviewer response for version 1 Sudbery Ian 1 Referee https://orcid.org/0000-0002-5038-0190 Sheffield Institute For Nucleic Acids (SInFoNiA), Department of Molecular Biology and Biotechnology, University of Sheffield, Sheffield, UK

Competing interests: No competing interests were disclosed.

 2 7 2019 Copyright: © 2019 Sudbery I 2019 This is an open access peer review report distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. recommendation approve-with-reservations

Mourao et al. describe a package for estimation of the rate at which stranding information in an RNA-seq sample might be wrong, thus leading to inaccurate quantification of anti-sense transcription, which they apply to their own data and the data of others. This is a fascinating and useful contribution to RNA-seq analysis. The software described should be useful to a range of scientists working the analysis of RNA biology from transcriptomics data. I will split my evaluation of this into two sections – the manuscript/science, and the software:

Manuscript:

The manuscript is generally well written and describes an interesting problem. I believe that the approach is generally sound. I have a few small issues that the authors might like to address:

It is assumed that anti-sense spliced reads are spurious. However, one possibility is that anti-sense spliced reads arise by the action of RNA-dependent RNA-polymerases. The authors acknowledge this point, but state “only 8-30% of eukaryotic gene regions have significant length ORFs on their opposite strands, providing an upper limit on the potential impact of this method of transcript on the RNA complement within a cell”.  I am unsure how the presence of an ORF on the anti-sense strand has an effect on the RNA complement of a cell. The authors should clarify this argument.

An underlying assumption of the approach is that the ratio of spurious anti-sense reads to sense reads is consistent across the length of the transcript, and thus modelled by the reverse junction reads. This assumption is probably necessary and reasonable, but perhaps the authors might make it explicit.

In regions where two annotated transcripts overlap, the authors truncate their antisense copies of these transcripts to only cover the region that is unique to the sense transcript in question. This is necessary to avoid double counting of reads. However, it does mean that the number of reads counted as antisense will be an underestimate. This becomes problematic because the full length transcript, not the truncated version is used when calculating how many spurious antisense reads are present (this is possible because this is calculated from spliced reads). Thus the correction factor will be overestimated compared to the counted number of antisense reads. This could be detected by finding genes where the correction factor is larger than the total number of antisense reads leading to a negative corrected value for the antisense count. Unfortunately the package truncates the corrected count at 0, so these counts are not evident in output from the package.

The authors show that there is a large variation in the proportion of anti-sense reads that are spurious in three replicates of the same experimental condition. It would be useful to know if this variation is bigger or smaller than that between conditions. Variation within replicates reduces power to detect differential expression between conditions, but systematic variation between conditions might lead to false positive calls of differential expression. The authors might consider looking at the within condition and between condition variance of samples in the Graveley data set they have analysed.

The authors introduce the concept of a “maximal transcript”. This is very similar to the concept of the superTranscript 1  the authors might like to acknowledge this and cite the relevant work.

The authors note that for some genes there is more anti-sense expression than sense expression and suggest this might be a case of incorrect annotation. Do these genes tend to be spliced or single exon. If spliced, does the annotation carry the canonical splice sequence, its reverse complement, or some other sequence?

Figure 5 is reproduced at a very low resolution, such that it is not possible to read all the text on the panel marked (Right).

 Software:

The software is provided as an R package with a number of accompanying python modules and is available on GitHub. An archived version of the source code used to produce the figures in the paper is available on Zenodo. Dependencies are intended to be installed using a conda environment. The authors are to be commended on their good practice in this respect.

However, I did have a number of issues in attempting to use the software:

I could not install the conda environment from the provided environment file. Running the suggested commands gave the following error:

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 spike-in 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.

Are the conclusions about the tool and its performance adequately supported by the findings presented in the article?

Yes

Is the rationale for developing the new software tool clearly explained?

Yes

Is the description of the software tool technically sound?

Yes

Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?

Partly

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool?

Partly

Reviewer Expertise:

I am a lecturer in bioinformatics. My research focuses on using the tools for computational biology and functional genomics to study RNA biology and regulation of gene expression, particularly at the post transcriptional level.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

References : SuperTranscripts: a data driven reference for analysis and visualisation of transcriptomes. Genome Biol .2017;18(1) : 10.1186/s13059-017-1284-1 148 28778180 10.1186/s13059-017-1284-1