R-SAP: a multi-threading computational pipeline for the characterization of high-throughput RNA-sequencing data

The rapid expansion in the quantity and quality of RNA-Seq data requires the development of sophisticated high-performance bioinformatics tools capable of rapidly transforming this data into meaningful information that is easily interpretable by biologists. Currently available analysis tools are often not easily installed by the general biologist and most of them lack inherent parallel processing capabilities widely recognized as an essential feature of next-generation bioinformatics tools. We present here a user-friendly and fully automated RNA-Seq analysis pipeline (R-SAP) with built-in multi-threading capability to analyze and quantitate high-throughput RNA-Seq datasets. R-SAP follows a hierarchical decision making procedure to accurately characterize various classes of transcripts and achieves a near linear decrease in data processing time as a result of increased multi-threading. In addition, RNA expression level estimates obtained using R-SAP display high concordance with levels measured by microarrays.


TaqMan qRT-PCR data analysis:
Four replicates of TaqMan qRT-PCR measurements for MAQC Human Reference sample were obtained from GEO. Original dataset consisted of normalized expression values and their presence/absence calls for 1044 probes. Expression value for each probe was taken as the mean expression values across the four replicates. A probe was considered expressed if it had at least 75% presence call (present call in at least three replicates). 973 expressed probes were retained after the filtering. Expressed probes were then assigned to RefSeq transcripts (hg18) using TaqMan probe annotations provided under GEO's platform record for TaqMan (platform: GPL4097). 962 probes were successfully assigned to RefSeq transcripts and 727 of them were also detected (R-SAP RPKM > 0) in MAQC Human Reference RNA-Seq data.

b. Abundance quantification mode:
In order to estimate the expression values of RefSeq (hg18) transcripts, we ran Cufflinks v1.1.0 in its quantification mode. Aligned reads from TopHat and RefSeq transcripts in GTF format were used as input for Cufflinks. Cufflinks was run as: cufflinks -G Hg18RefSeq.gtf -u -F 0.0 --overhang-tolerance 3 <tophat-alignments> Parameter "--overhang-tolerance 3" was used to match R-SAP's default cutoff for "exonextension" (3 bp). Cufflinks reports abundance estimates as FPKM (fragments per kilobase of transcript per million fragments mapped) that are comparable to RPKM (reads per kilobase of transcript per million reads mapped) for single end read data.

Estimating expression:
rsem-calculate-expression --out-bam --seed-length 28 <bowtie-path> -bowtie-n 3 --bowtie-e 200 --bowtie-m 1 --phred33-quals --fragment-length-mean 200 --fragment-length-sd 80 <input-fastq-file> <reference-transcript-sequence-index> Cufflinks default values for fragment length distribution mean and standard deviation were used for RSEM run. RSEM estimated RSEM reports two measures of abundance estimates: Expected read count from each transcript and estimated fraction of transcripts made up by a given isoform or gene (τ value). τ value is generally converted to TPM (transcripts per million) by multiplying by 10 6 in order to get the expression value of transcripts. TPM value is not directly comparable to RPKM or FPKM value. We converted TPM values to comparable RPKM values using the conversion formula provided in (7): Here i corresponds to the i th RefSeq transcript and j varies from 1 to total number of RefSeq transcripts. l j is the length (in bp) of j th transcript.
Supplementary References:    Table 3: Intron-retention events detected in MAQC Reference Human dataset using R-SAP from high-scoring reads that were also characterized as internal-exon-extension.
Supplementary Table 4: Distribution of "multiple-annotations" reads that were detected in MAQC Reference Human dataset using R-SAP. Since more than one type of novel transcriptional event was detected in each of the "multiple-annotations" reads, reads representing characterization sub-categories here may be overlapping. "Exon-skipping" and "intronretention" events are already included in "exon-deletion" and "internal-exon-extension" events.  Table 5: Distribution of Trans-ABySS characterized reads that were also classified as "high-scoring" by R-SAP previously using MAQC Human Reference RNA-Seq data. RefSeq transcripts (hg18) were used as annotated set of transcripts. "novel-transcript" are those that could not be mapped to any of the known RefSeq exon by Trans-ABySS. Novel transcriptional event sub-categories (AS3, AS, AS53, novel_exon, novel_intron, novel_utr, retained_intron and skipped_exon) have 121 overlapping reads between them. Filtered-out reads are those that were not reported in any category by Trans-ABySS.