dataesr/scientific-paragraphs-categorization huggingface dataset:
Error models and simulations We aligned RNA-Seq reads with permissive parameters (quality aware alignment, with no fixed mismatch cutoff ) using Bowtie  [19]  in order to estimate total mismatch profiles along the full length of the reads. As has been previously reported, we observed increased mismatch rates extending through the 3′ end of the read and a slight increase in mismatch rates in the first 5 bases  [4, 21, 22]  (Figure  2A ). We determined nucleotide mismatch frequencies by position in the read and by substitution type. These frequencies are used by the Spanki read simulator as weights in a weighted-random selection, to choose the total number of mismatches in a read, mismatch positions, and base substitution. We supplied these error models to the Spanki read simulator to generate a defined known input sample generated in silico from annotated transcript models, with error profiles matching these empirical models.
We then uniquely aligned the reads using TopHat  [23] , and compared alignment results to the known input to explore splice site detection parameters. This two-step process generated a simulated data set that mirrors the experimental dataset, except that the true input was known, providing us a platform for testing RNA-Seq junction alignment. Results that follow provide evaluation results for the TopHat aligner  [23] , although the same approach can be applied to any aligner output. We compared junction coverage with known input abundance for all junctions in the 10x transcript coverage pool (Figure  2B ). Since multiple transcripts at a locus may share a given junction, individual junction coverage was 1-400x (median 8x, 4. 2 million read pairs) reflecting both the random sampling of read positions and overlapping D.
The first step is junction quantification (C), where alignments are performed, junction alignments are curated, and junction coverages are calculated. Splicing event quantification (D), where a set of transcript models (from annotation or computed using a program such as Cufflinks  [11] ), are used to characterize pairwise splicing differences ("splicing events"). These events are merged with junction coverage data to quantify the mutually exclusive paths defined for each event. Splicing event comparison (E) uses these tabulated event-level quantifications to compare between replicates, and between pooled results for each sample, by Fisher's Exact Test on inclusion and exclusion junction counts. these simulated reads with a reference annotation (' Annotation guided') or without ('De novo'), and compared detection results with known input (Figure  2C ).
We did not observe overrepresentation of any splicing event type with increased depth. We generated normalized transcript abundance estimates in units of fragments per kilobase per million mapped reads (FPKM), and found that we obtained 6x coverage of 95% of the transcripts reliably detected at FPKM ≥ 1 in the full dataset (200 million mapped reads). We examined the false positive rate at multiple transcript coverage levels (Figure  2F ) and found that the rate increased with greater transcript coverage due to cumulative errors in alignment. These data indicate that greater read depth provides more opportunities to call false positives. However, the majority of the false positives can be filtered post hoc (Figure  2F ) as we explain later.
Studies in organisms with incompletely annotated genomes rely heavily on empirically detected junctions. Spanki's design allows the flexible application of these filters, which is critical to accommodate different sample types and alignment strategies. For example, aberrant splicing may be a feature of interest in mutant or cancer samples, rather than an artifact to filter out  [31]. Although we present results using the first generation TopHat aligner, other tools allow searching for fusion transcripts (TopHat2,  [28] ), or non-canonical splicing variants (MapSplice,  [32] ). In these cases, simulation allows for assessment of error rates and selective application of filters.
Differential splicing detection To identify cases of differential splicing between samples, it is essential to find where splicing patterns diverge to define a unit of comparison. Basic categories of alternative splicing incompletely describe complex splicing patterns, which can lead to under-reporting of differences. We applied standardized and exhaustive splicing event ontologies with AStalavista  [16]  to ensure that pairwise splicing differences are interrogated completely. Spanki parses AStalavista output to obtain sets of junctions that define mutually exclusive "paths" (Inclusion and Exclusion), to identify junctions that interrogate each path specifically. We use the detection of coverage over these junctions to calculate PSI. To analyze false positive calls of differential splicing, we generated a Spanki null model for splicing differences by simulating four read pools, each of which contained reads from all annotated transcripts in equal abundances of 300 Reads Per Kilobase (RPK).
We then analyzed this null dataset, applying junction filtering and Spanki junction-based event definitions, and also compared technical and biological replicates of real RNA-Seq data to each other. We performed a systematic categorization of all pairwise relationships using AStalavista  [16] , which constructs graphs from transcript models and outputs complete and non-redundant sets of splicing differences identified through graph alignment (see Methods). We used these pairwise definitions to make comparisons on the PSI metric calculated by Spanki. We found minor variation due to sampling alone and technical replication (Figure  4A, B ), but biological replication was a much greater source of variability (Figure  4C ). Low total abundance events showed the most disagreement between replicates (< 10 average coverage per site in either path).
Next we compared the number of differential splicing calls made in our simulated null dataset by Spanki and by several other methods. Spanki correctly called zero events differentially spliced in this dataset (Figure  4D ). We counted reads that map within exons using the script provided with DEXSeq  [12] , and performed an exon-level differential analysis. DEXSeq also called zero exons differentially expressed, however, with a reduced sensitivity to real alternative events in other data sets (not shown). Next we performed an isoform-centric analysis using MISO  [13] , which called differential splicing in transcripts of 222 genes. Analysis with Cuffdiff  [11] , with default parameters except for specifying upper quartile normalization, called 183 loci as differentially spliced, and 267 isoforms were called differentially expressed in our null input dataset with no true differences in splicing.
To analyze the false negative rate of differential splicing detection, we generated two simulated datasets with a known PSI splicing difference. To prevent cross-talk, we selected splicing events where both the inclusion and exclusion forms were composed of transcripts that were not part of any other splicing event (N = 1644 events), and simulated one pool so that PSI = 0. 25 for all events, and a second pool with 0. 75 PSI for all events, so that comparing the two would yield -0. 50ΔPSI. After processing these data with Spanki, the ΔPSI values clustered at -0. 50ΔPSI (median -0. 46) (Figure  4E ). Spanki failed to call significant differences in 315 events (19%) with FDR correction (Figure  4F ). Cufflinks with the Jensen-Shannon Divergence (JSD) metric had a 29% false negative rate, but performed better when comparing isoform abundances, failing to find differential isoform expression in only 5% of isoforms where an abundance difference was simulated.
At the exon level, DEXSeq had a false negative rate of 19%, but with false positive differences in 116 genes. MISO performed the best at this task, which did not detect a splicing difference in 115 genes (7% false negative rate). These results show that this simulation was a challenge for these tools and produced high false negative rates, but Spanki performed comparably to other tools at the same task. This junction-based calculation of splicing differences is a more balanced and less biased measure than exon counts. Although exon counts yield more data, ambiguity of assignment (Figure  1A ), coverage heterogeneity  [4] , and unprocessed transcripts  [33]  make these data unreliable. Exon counts are also a more imbalanced measurement of alternative forms. For example, in the case of skipped exons, the cassette inclusion form can be interrogated by reads within the cassette, but the exclusion form has zero exonic space that can be uniquely interrogated.
A B C Figure  5  Splicing event characterization in Drosophila heads. (A) Pairwise splicing events defined for all transcript models in Flybase 5. 39 annotation (black bars), and the subset of those detected as alternatively spliced in female and male Drosophila heads by Spanki (grey bars). Black bars indicate splicing event classes as a percent of all defined pairwise events. Grey bars indicate splicing events detected as a percent of all detected events. Cartoons of each event type are in the leftmost column, with the "inclusion" form indicated in green, and the "exclusion" form indicated in orange. A description of each type is adjacent to each cartoon. The "Unclassified" types includes diverse complex type with no concise verbal description. (B) Ratio vs average abundance scatterplots of gene expression in female and male heads.
We found that 7,894 splicing events in transcripts from 2,441 genes were alternatively spliced in head samples (5,450 internal events in 1,852 genes) (Figure  5B ). We found 182 events with significant differences between female and male heads (adjusted p-value < 0. 01, Benjamini and Hochberg, Additional file 3: Table  S3 ). To conservatively adjust the sex-biased expression calls for rare events and biological variability as outlined previously, we demanded >10 junction reads in each path, and that the unadjusted p-value for the between-sexes comparison was less than the unadjusted p-value between biological replicates. In light of the variance we observed in our null dataset (Figure  4A ) and between biological replicates (Figure  4C ), we also set a conservative threshold on the difference between samples at > 0.
6-99. 4% by Spanki (Figure  6A ). There are possible technical explanations for the non-sex-specific exon read results, including the sequencing of transcripts that are polyadenylated before splicing is complete (more prevalent on the 3-prime end in dsx). Inefficient poly-A selection due to oligo (dT) priming in unprocessed transcripts may also play a role  [33]. Our measurements using splice junctions better reflect processed transcriptional output. These results show that using junction coverage with Spanki results in more switch-like splicing difference calls. Similarly, we found highly sex-biased splicing of other members of the sex determination cascade, except for the tra locus where incomplete splicing to the female form is known to occur  [44]  (Figure  6B ). None of the other splice detection methods called sex-biased splicing for each of the sex determination gene transcripts (Figure  6C ), clearly highlighting the improvements made with filtered Spanki output.
Conclusions Junction-based splice calling is an important method for analysis of alternative splicing. Our results highlight many of the junction-read errors that can occur in these RNA-Seq datasets and outline simulation strategies for modeling these errors, either while developing an analysis pipeline, or tailored for each experiment. We have implemented tuneable filters in Spanki to remove false positives without sacrificing specificity, and clearly show that developing error models from real RNA-Seq data and applying these post hoc filters improves splicing detection in D. melanogaster heads. Failing to filter for sequencer and aligner error, to account for mapping ambiguity in transcriptome space, and to under-estimate the null distribution of splicing differences, results in inaccurate estimation of splicing differences in RNA-Seq studies.
PCR was performed in biological duplicates, with duplicate quantification of each biological duplicate. cDNA from 12. 5 ng of total RNA was amplified with Fast SYBR Green Master Mix (Applied Biosystems, Carlsbad, CA, USA) in a StepOne Real-Time PCR machine (Applied Biosystems, Carlsbad, CA, USA). Initial activation was performed at 95°C for 20 seconds followed by 40 cycles. Cycles were 95°C for 3 seconds followed by 60°C for 30 seconds. Then the melting curve was generated ranging from 60°C to 95°C with an increment of 0. 5°C each 5 seconds. Act5c (Actin 5C) was used as a control. Primers were designed with the web interface of the NCBI Primer-Blast software  [47]. All amplification products were analyzed by agarose gel electrophoresis and produced single fragments of predicted sizes. The relative transcript level was calculated using the cycle threshold value (Ct) by the method of 2 -ΔCt , where ΔCt = Ct transcript -Ct Act5c.
Read mapping We used reads that passed Chastity (score ≥ 0. 6) basecalling filtering (Illumina CASAVA pipeline 1. 6. 47. 1) and mapped using TopHat v1. 4. 1  [11] , with Bowtie v0. 12. 7  [19] , and samtools 0. 1. 12a [48] , and parameters "-g 1 -solexa1. 3-quals, -i 42. " We used D. melanogaster genome release 5  [49, 50] , as obtained from the UCSC genome browser (excluding "chrUextra")  [51] , for mapping. We also appended sequence for 96 exogenous controls to the genomic reference  [4]. A reference annotation (Ensembl release 67, corresponding to Flybase 5. 39) was also supplied in GTF format with the -G option. We made a minor modification in the annotation to remove the antisense transcripts of modifier of mdg4 (FBgn0002781), since these transcripts caused fatal errors in downstream analysis tools.
[11, 52]. We used Cufflinks  [11]  (v. 2. 0. 2) to generate abundance estimates of full-length isoforms, expressed in units of expected fragments per kilobase of transcript per million mapped reads (FPKM). We determined relative abundance differences using Cuffdiff v. 2. 0. 2  [11] , using upper quartile normalization, and setting "max-bundle-frags" very high (50E06), to ensure that very highly expressed features were not excluded. To provide alternative quantifications and comparisons, we used HTSeq  [52]  to generate counts of reads that fall within discrete features. The "htseq-count" program in HTseq v. 0. 5. 3, with the conservative "union" mode and default parameters, was used to generate counts. We used the R package DESeq (v. 1. 8. 3)  to test for differential expression  [52]. "Variance outliers" were identified as contrasts where the maximum residual variance is > 15.
Simulation and splicing analysis Simulations of junction reads, along with quantification of junctions and alternative splicing, was performed in the open-source python package: Splicing Analysis Kit (Spanki). A summary of Spanki's features is provided in Table  1. Spanki is available at http://www. cbcb. umd. edu/ software/spanki and http://github. com/dsturg/Spanki. Error models were built by performing quality-aware mapping on our reads with Bowtie v0. 12. 7  [19] , and supplying the map file to the program spankisim_models. These models are incorporated into the Spanki repository, so that they can be applied to simulations by choosing the "flyheads" error model. For generating simulated reads, we used the spanki-sim_transcripts command, supplying a reference annotation (Ensembl release 67, corresponding to Flybase 5.
39), using the default parameters for intron retention (-ir = 0. 20) and fragment size (-frag = 200 bp). The coverage to simulate was specified as either coverage values with the -cov parameter (for the simulations in Figure  2 ), or as reads per kilobase with the -rpk (for the simulations in Figure  4 ). We used the "flyheads" error model, which was built as described above, and is included in the program. The simulated reads in this study are available here: http://www. cbcb. umd. edu/software/spanki/ simulations. html. We also compared our detection performance to other tools, including Cuffdiff, MISO, and DEXSeq. We used the results reported in the "splicing. diff" file from Cuffdiff v. 2. 0. 2  [11]. We used MISO with event definitions for D. melanogaster from the MISO website, supplemented with additional custom definitions, and applied a Bayes factor cutoff > 40 for the false positive analysis, and > 10 for the false negative analysis  [13].
Spanki program design (Simulator) Spanki estimates error model parameters from a first pass alignment of real RNA-Seq reads using permissive quality aware mapping with Bowtie  [19]. The error modeling function within Spanki parses the alignments in Bowtie's map format, and produces probability weight matrices for mismatches by position in the read and by base substitution type, and for quality scores by position. The read simulator uses these models to introduce mismatches. Spanki's RNA-Seq simulator function generates simulated reads with errors incorporated. The user sets the transcripts to simulate (e. g. transcripts expressed in the biological sample under study), a depth of coverage (e. g. matching the experimental sequence depth), and mismatches are then introduced according to the specified model.
If a user does not want to build new error models, pre-built error models from D. melanogaster head RNA-Seq described in this work are included, along with models based on a sample from the modEN-CODE developmental timecourse  [29] , and a simple weighted-random model. The Spanki simulator takes transcript models in GTF format, and extracts transcript sequence from a genomic reference and chooses random positions in the designated transcript sequence to extract reads. For paired-end reads, fragment sizes are drawn from a normal distribution of mean 200 bp and standard deviation 20 bp. This is user tunable. For example, to simulate intron retention, Spanki generates a specified fraction of simulated reads from complete transcript sequence where introns are retained. Depth of coverage is specified in units of transcript coverage or reads-per-kilobase (RPK).
Reads can be generated for transcripts in fixed proportions, creating a null model for splicing differences between samples. Alternatively, Spanki accepts a text file where the user can list individual transcripts to simulate at different coverages, which allows simulating fixed quantitative splicing differences between alternative isoforms, or the user can specify a custom model built on the user's own data. Modeled error frequencies are applied as weights for mismatch number, position, and substitution. Weight matrices of quality scores are used to create a consensus quality value across all positions -one for matched positions, and one for mismatched positions, which are concatenated to create a quality string for the read. Spanki reports information that facilitates analysis of alignment and detection.
Coverage generated by the simulation for each splice junction is reported, along with read counts for each transcript. To enable the tracking of aligner errors, the genomic coordinates of origin for each read are incorporated into a unique read identifier. The true origin of simulated reads is also reported in a SAM file that represents a perfect alignment, which can be fed to an assembler such as Cufflinks  [11]  to allow the evaluation of error in transcript abundance estimates due to assembly separately from errors in alignment.
Spanki program design (Junction filtering quantification) For maximum flexibility, Spanki decouples the alignment and filtering steps, with a tool that applies posthoc analyses of alignment files. This allows alignments to be performed on multiple data sets, with consistent filtering applied later, and allows changing the filtering criteria without re-aligning. Spanki streams through a BAM file produced by any aligner, using the Pysam module (Andreas Hager, http://code. google. com/p/pysam/) and calculates junction coverage along with alignment diagnostic measurements. These measurements include the number of alignment offsets, alignment entropy  [29] , and Minimum Match on Either Side (MMES)  [53]. Qualitative diagnostic results are also reported, such as repetitiveness of exon anchor and intronic sequence.
Spanki program design (Splicing event definition and quantification) While there are only a few major splicing type events, the full diversity of possible splice forms is much more complex. Spanki provides utilities for parsing essentially any splicing event using definitions produced by AStalavista  [16]. The AStalavista algorithm begins by decomposing transcript models into "sites," which are exon boundaries. Graphs are built for each gene, where splice sites are nodes and intron or exon edges connect them. Splicing events are subgraphs with identical nodes on ends, but no common interior nodes. This process finds regions of the parent transcript where the donor/acceptor sites of two alternatives are present on a parent transcript, but utilized mutually exclusively in processed transcripts.
The Fisher's Exact Test (FET) is well suited to this task, since testing proportions accounts for differences in sample totals due to depth or transcription. Spanki constructs 2 × 2 contingency matrices from junction counts for each splicing event, to test the null hypothesis that the two samples have equal inclusion/exclusion proportions. The two cells of the first row of the matrix are the total read counts of the inclusion and exclusion junctions, respectively, for one sample. The second row contains the same data for the second sample. The test as constructed uses integer counts, as required for the Fisher's exact test. Each defined pairwise event is tested, and each gene may have multiple events. The test is performed using the fisher python package v. 0. 1. 4 (Brent Pederson, http://pypi.
python. org/pypi/fisher/). FDR correction is performed by the Benjamini-Hochberg method implemented in the StatsModels package (Skipper Seabold, Josef Perktold, http://statsmodels. sourceforge. net/). To help visualize splicing differences, Spanki includes R scripts to produce mosaic plots, where the relative size of each cell is proportional to real (non-normalized) cell counts. Code is also included to produce fourfold plots, which provide a visual test of the null hypothesis of the FET. This provides an effective simultaneous visualization of normalized proportions and significance. These plots are implemented in the "vcd" package for R  [56].
MISO results shown are based on isoform-centric analysis. (E) Volcano plot of ΔPSI calculated by Spanki for false negative analysis, comparing simulated datasets with PSI = 0. 25 and PSI = 0. 75 respectively, for 1644 events (expected ΔPSI = -0. 50, red line). (F) Counts of false negatives for differential splicing calls in simulated data with an input ΔPSI = -0. 50.
Acknowledgements We thank  Michael Sammeth  for sharing source code and  Trudy Mackay  for isogenic flies. Ryan Dale  provided technical advice. We also thank  Steve Mount , fellow members of our labs, and the modENCODE Drosophila transcription group ( Sue Celniker ,  Roger Hoskins ,  Steven Brenner ,  Brenton Graveley ,  Peter Cherbas , and  Thomas Gingeras ), for stimulating discussion, advice, feedback, and encouragement. This study utilized the highperformance computational capabilities of the Biowulf Linux cluster at the  National Institutes of Health, Bethesda, MD. (http://biowulf. nih. gov). This research was supported by the  Intramural Research Programs of the National Institutes of Health (NIH), National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK). DS was supported by a  pre-doctoral Intramural Research Training Award fellowship  through the  NIH Graduate Partnerships Program.