Identification and quantification of circular RNA from RNA-Seq data

Here we run an analysis pipeline in R for the identification of circRNA. The pipeline attempts to recreate parts of the KNIFE pipeline described in the paper by Szabo et al. (2015). The code provided here uses the systemPipeR package with functions designed to run the pipleline on a local desktop. Consult the systemPipeR documentation to modify the code so it can be run on a cluster.

Prepartory steps

To install and resolve bioinformatics software dependecies, you can use the Bioconda package manager. First download Miniconda, then create an environment, activate it, add bioconda channels, and install the required software. For example, here we create a local environment named rnaseq and install bowtie2 and trim-galore.

conda create -n rnaseq python=2 anaconda
source activate rnaseq
conda config --add channels conda-forge
conda config --add channels defaults
conda config --add channels r
conda config --add channels bioconda
conda install bowtie2
conda install trim-galore

1. Create project directories

First start a project and create subdirectories under the project’s working directory

folders <- c(
  "parameters", "src",
  "data/genome_data/hg19/fasta",
  "data/genome_data/hg19/bowtie2_index",
  "data/biosamples",
  "results/trimmed",
  "results/aligned/genome",
  "results/aligned/ribosomal",
  "results/aligned/junctions_reg",
  "results/aligned/junctions_scrambled",
  "/results/analysis/reports"
)

for (i in 1:length(folders))  {
  dir.create(paste(getwd(), folders[i], sep = "/"), recursive = TRUE)
}

2. Download genome data

To get started, we need to download the hg19_fastas.tar.gz archive that contains fasta files for the ribosome, the linear junctions, the scrambled junctions, and the reference genome. The compressed tar file is available at this link

3. Build index files

We first load the systemPipeR package and other required packages. We also create external files for the command-line software. These are the param files that specify the parameters and the target files that contain the relative path locations of any required data files.

#source("http://bioconductor.org/biocLite.R")
#biocLite("tgirke/systemPipeR", dependencies = TRUE)
library(systemPipeR)
library(dplyr)

The index files are large and might take a long time to download on a regular connection. It could be faster to just run the Bowtie index builder to create the index files. These two lines build the indices. The param and target files are available here.

args_bowtie2_index <- systemArgs(sysma = "./parameters/bowtie2_index.param",
                                 mytargets = "./parameters/hg19_fasta_filepaths.txt")
runCommandline(args = args_bowtie2_index)

4. Download sample data

We use the same HER2 Positive Breast Tumor dataset used in the KNIFE as sample data. The data is made up of paired-end reads from a total RNA rRNA-depleted RNA-seq library. More details about the experiment can be found here. To ensure reproducibility and data provenance, we download the SRA data files directly from the GEO database.

download.file("ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX/SRX374/SRX374866/SRR1027187/SRR1027187.sra", "./data/biosamples/SRR1027187.sra")

Next we convert the SRA file to fastq using fastq-dump. The options provided specify the following: split the paired reads into two files, dump the biological reads only, and compress the output using gzip. We also keep the defaults by not adding read ids and leaving the quality encoding as ASCII+33.

system("fastq-dump --split-files --skip-technical --gzip ./data/biosamples/SRR1027187.sra --outdir ./data/biosamples/")

5. Trim data

Reads in the raw data can contain low quality bases and partial adapter sequences that needs to be removed. Trimmed reads that are too short should also be filtered in order to obtain good alignment.

args_trim <- systemArgs(sysma = "./parameters/trim_galore.param",
                           mytargets = "./parameters/data_filepaths.txt")
trim_command <- paste(sysargs(args_trim)[1],
                      unlist(strsplit(sysargs(args_trim)[2], " "))[8],
                      sep = " ")
(trim_command)

[1] "trim_galore --length 50 --paired --output_dir /hd2/KNIFE_project_01/results/trimmed/  /hd2/KNIFE_project_01/data/biosamples/SRR1027187_1.fastq.gz  /hd2/KNIFE_project_01/data/biosamples/SRR1027187_2.fastq.gz"

Although we are treating each paired end read independently, we add the the –paired option argument in the trim_galore command to peform a validation test that ensures that both sequence pairs have a certain minimum length specified by the option --length. To run trim-galore we use the system function

system(trim_command)

6. Run data quality control checks

Lastly, we generate and plot FASTQ quality summary statistics. The fastqc manual provides a good explantion of the fastqc plots. First we define a helper function

runFastqPlots <- function(file_list, report_name) {
  fqlist <- bplapply(seq_along(file_list),
                     function(x)
                       seeFastq(
                         fastq = file_list[x],
                         batchsize = 100000,
                         klength = 8
                       ),
                     BPPARAM = MulticoreParam(workers = 4))
  png(
    sprintf(
      "./results/analysis/reports/fastq_report_%s.png",
      report_name
    ),
    height = 800,
    width = 200* length(fqlist),
  )
  seeFastqPlot(unlist(fqlist, recursive = FALSE))
  dev.off()
}

then we run it on the raw data

runFastqPlots(file_list = infile1(args_trim),
              report_name = "SRR1027187")

and on the trimmed data

data_trimmed_filepaths <- paste0(outfile1(args_trim),
                                 SampleName(args_trim),
                                 c("_val_1.fq.gz", "_val_2.fq.gz"))
names(data_trimmed_filepaths) <- paste0(SampleName(args_trim),"_trimmed")

runFastqPlots(file_list = data_trimmed_filepaths,
              report_name = "SRR1027187_trimmed")

Next we plot the reports side by side

The plots for the trimmed data look much better. Although the mean quality distribution for the reads is bimodal, there are no low quality peaks.

Alignment Steps

Now we independently align each paired-end read dataset to four separate Bowtie2 indices: for the genome, ribosomal RNA, linear exon–exon junctions, and scrambled exon–exon junctions. First we load Bowtie2 parameter files and filepaths

args_align_genome <- systemArgs(sysma = "./parameters/bowtie2_align_genome.param",
                                mytargets = "./parameters/data_trimmed_filepaths.txt")
args_align_ribo <- systemArgs(sysma = "./parameters/bowtie2_align_ribosomal.param",
                              mytargets = "./parameters/data_trimmed_filepaths.txt")
args_align_juncts_reg <- systemArgs(sysma = "./parameters/bowtie2_align_junctions_reg.param",
                                    mytargets = "./parameters/data_trimmed_filepaths.txt")
args_align_juncts_scram <- systemArgs(sysma = "./parameters/bowtie2_align_junctions_scrambled.param",
                                      mytargets = "./parameters/data_trimmed_filepaths.txt")

then we run the Bowtie2 alignment commands through systemPipeR

# align to genome
genome_bam_files <- runCommandline(args = args_align_genome)
# align to ribosome
ribo_bam_files <- runCommandline(args = args_align_ribo)
# align to linear junctions
juncts_reg_bam_files <- runCommandline(args = args_align_juncts_reg)
# align to srambled junctions
juncts_scram_files <- runCommandline(args = args_align_juncts_scram)

We can open “submitargs01_log” for each alignment to display the alignment summary and the system command used. For the genome alignment, we can open the file inside R like this

readLines(file.path("./results/aligned/genome","submitargs01_log"))

 [1] "bowtie2 -p 8 --no-unal --score-min L,0,-0.24 --rdg 50,50 --rfg 50,50 -x /hd2/KNIFE_project_01/data/genome_data/hg19/bowtie2_index/hg19_genome -U /hd2/KNIFE_project_01/results/trimmed/SRR1027187_1_val_1.fq.gz  -S /hd2/KNIFE_project_01/results/aligned/genome/SRR1027187_1_val_1.sam"
 [2] "19308315 reads; of these:"
 [3] "  19308315 (100.00%) were unpaired; of these:"
 [4] "    6692805 (34.66%) aligned 0 times"
 [5] "    8628824 (44.69%) aligned exactly 1 time"
 [6] "    3986686 (20.65%) aligned >1 times"
 [7] "65.34% overall alignment rate"
 [8] "bowtie2 -p 8 --no-unal --score-min L,0,-0.24 --rdg 50,50 --rfg 50,50 -x /hd2/KNIFE_project_01/data/genome_data/hg19/bowtie2_index/hg19_genome -U /hd2/KNIFE_project_01/results/trimmed/SRR1027187_2_val_2.fq.gz  -S /hd2/KNIFE_project_01/results/aligned/genome/SRR1027187_2_val_2.sam"
 [9] "19308315 reads; of these:"
[10] "  19308315 (100.00%) were unpaired; of these:"
[11] "    6753734 (34.98%) aligned 0 times"
[12] "    8563704 (44.35%) aligned exactly 1 time"
[13] "    3990877 (20.67%) aligned >1 times"
[14] "65.02% overall alignment rate"                                                                                                                                                                                                                                                          

Preprocessing steps

The previous alignment step produced 8 BAM files. Two paired-end BAM files for each of the genome, ribosome, linear junctions, and scrambled junctions alignments. We are interested in junctional reads, those reads from the 8 files that align to a scrambled or linear junctions. Since the majority of the reads in the genome and ribosome BAM files are not junctional, we need to filter them out.

1. Get read IDs for genomic mate pairs

Since the next preprocessing steps are memory intensive, we can intially filter out reads in which both pair mates align to the genome or ribosome. We can also filter out unaligned reads and return lists of read IDs in which both pairs align to the indicies, or in which one pair mate aligns but the other doesn’t. To proceed we need to get the names of those reads that are found in both R1 and R2 genomic alignment files.

getReadIds <- function(bamfile) {
  scanBam(bamfile,
          param = ScanBamParam(flag = scanBamFlag(isUnmappedQuery = FALSE),
                               what = "qname"))[[1]][[1]]
  }
getCommonGenomicPairedReads <- function(genome_bamfiles, ribo_bamfiles) {

  genome_reads <- bplapply(genome_bamfiles,
                           getReadIds,
                           BPPARAM = MulticoreParam(workers = 2))
  ribosomal_reads <- bplapply(ribo_bamfiles,
                              getReadIds,
                              BPPARAM = MulticoreParam(workers = 2))

  genomic_reads_1 <- union(genome_reads[[1]], ribosomal_reads[[1]])
  genomic_reads_2 <- union(genome_reads[[2]], ribosomal_reads[[2]])
  common_reads <- intersect(genomic_reads_1, genomic_reads_2)
  r1_not_r2 <- setdiff(genomic_reads_1, genomic_reads_2)
  r2_not_r1 <- setdiff(genomic_reads_2, genomic_reads_1)

  return(
    list(
      common_reads = common_reads,
      r1_not_r2 = r1_not_r2,
      r1_not_r2 = r2_not_r1
    )
  )
}

Note that code above uses the parallel lapply function bplapply. If you system has less than 32G RAM, you can modify the code to run with one worker or use lapply instead. Here we retrieve the location of the aligned files and feed them into the function defined above.

genome_aligned <- outpaths(args_align_genome)
ribo_aligned <- outpaths(args_align_ribo)
reads_to_ignore <- getCommonGenomicPairedReads(genome_aligned, ribo_aligned)

2. Read alignment data and run initial filtering

Now we use these read IDs to create two filters. The first one filters out genomic reads in which both mates (R1 and R2) align. The second one is a junctional filter that keeps only junctional reads that overlap a particular junction by at least 8 nucleotides and that also satify the first filter.

genomic_reads_filter <-
  FilterRules(list(
    include = function(x) !(mcols(x)$qname %in% reads_to_ignore$common_reads)
  ))
MIDPOINT <- 150
JUNC_OVERLAP <- 8
junction_reads_filter <-
  FilterRules(list(
    include = function(x)
      (
        (MIDPOINT + JUNC_OVERLAP + 1 - width(x) <= start(x)) &
        (start(x) <= MIDPOINT - JUNC_OVERLAP + 1) &
        !(mcols(x)$qname %in% reads_to_ignore$common_reads)
      )
  ))

As part of the filtering step, we also update the alignment score for the junction reads such that the N-penalty is corrected.

countAmbigousBases <- function(align_obj) {
  sapply(gregexpr("N0|0N", mcols(align_obj)$MD),
         function(x)
           if (attributes(x)$match.length[1] == -1)
             return(0)
         else
           return(length(x))
  )
}
#MD and XN do agree
filterAlignments <-
  function(bamfile, param, filter, updateScore = TRUE) {
    x <- readGAlignments(file = bamfile, param = param)
    x <- subsetByFilter(x, filter)
    seqlevels(x) <- seqlevelsInUse(x)
    # correct for N-penalty in junction alignments and update alignment score
    if (updateScore) {
      mcols(x)$AS <- mcols(x)$AS + countAmbigousBases(x)
    }
    return(x)
  }
# get file paths
linear_junctions_aligned <- outpaths(args_align_juncts_reg)
scrambled_junctions_aligned <- outpaths(args_align_juncts_scram)
# run filters
param <-
  ScanBamParam(
    what = c("qname","mapq"), # get the read name and the mapping quality
    flag = scanBamFlag(isUnmappedQuery = FALSE),
    tag = c("AS", "XS", "XN", "MD") # get tags
  )
genomic_filtered <- bplapply(
  genome_aligned,
  filterAlignments,
  param = param,
  filter = genomic_reads_filter,
  updateScore = FALSE,
  BPPARAM = MulticoreParam(workers = 2)
)
linear_junctions_filtered <- bplapply(
  linear_junctions_aligned,
  filterAlignments,
  param = param,
  filter = junction_reads_filter,
  BPPARAM = MulticoreParam(workers = 2)
)
scrambled_junctions_filtered <- bplapply(
  scrambled_junctions_aligned,
  filterAlignments,
  param = param,
  filter = junction_reads_filter,
  BPPARAM = MulticoreParam(workers = 2)
)

Let’s take a look at first few rows of the GAlignment object for the linear junctions R1 aligned reads

linear_junctions_filtered[[1]][1:3,] %>% as.data.frame()

We can also examine the frequency of the cigar strings

linear_junctions_filtered[[1]] %>%
  cigar() %>%
  table %>%
  sort(decreasing = TRUE) %>%
  head() 

.
   60M    59M    58M    57M    56M    55M
397762 138477  45942  13278   3761   3041 

3. Convert GAlignment objects into dataframes

toDataFrame <- function(x) {
  x <- x %>%
    as.data.frame() %>%
    subset(select = -c(end, cigar, qwidth, njunc, XS, MD))
  colnames(x)[1:8] <-   c("junction", "str", "pos", "len",
                         "id", "MQ", "AS", "XN")
  x$junction <- as.character(x$junction)
  x$str <- as.character(x$str)
  return(x)
}
appendJunctionTags <- function(x) {
  df <-  x %>%
    seqnames() %>%
    as.character() %>%
    lapply(., function(y) unlist(strsplit(y, "[:]|[|]"))) %>%
    do.call(rbind, .) %>%
    data.frame(stringsAsFactors = FALSE)
  colnames(df) <- c("jChr", "jGene1", "jPos1", "jGene2",
                    "jPos2", "jType", "jStr")
  df$jPos1 <- as.integer(df$jPos1)
  df$jPos2 <- as.integer(df$jPos2)
  mcols(x) <- cbind(mcols(x), df)
  return(toDataFrame(x))
}
genomic_filtered_dt <- bplapply(
  genomic_filtered,
  toDataFrame,
  BPPARAM = MulticoreParam(workers = 2)
)
linear_junctions_filtered_dt <- bplapply(
  linear_junctions_filtered,
  appendJunctionTags,
  BPPARAM = MulticoreParam(workers = 2)
)
scrambled_junctions_filtered_dt <- bplapply(
  scrambled_junctions_filtered,
  appendJunctionTags,
  BPPARAM = MulticoreParam(workers = 2)
)

Let’s take a look at the dataframe for the linear junctions R1 aligned reads

linear_junctions_filtered_dt[[1]][1:3,]

4. Filter R1 junctional reads and pair them with R2 mates

A read (R1) is a linear junction read if it aligns to a linear junction and does not align to the genome or ribosome. In the code these reads are stored in linear_1. A read R1 is a scrambled junction read if it aligns to a scrambled junction and does not align to the genome or ribosome or to the linear junction. Scrambled junction reads are stored in scrambled_1. The read mate R2 can be a genomic, linear junction, or a scrambled junction read as determined by whehter the read name can be found in the R2 aligned reads. An R1 read can have a single or multiple R2 mates. Some have none at all. These are the unmapped reads.

createJunctionsReadsDF <- function(genomic,
                                   linear,
                                   scrambled,
                                   r1_not_r2,
                                   suffixes = c(".R1", ".R2g", ".R2l", ".R2s")) {

  linear_1 <- linear[[1]][!(linear[[1]]$id %in% r1_not_r2), ]
  scrambled_1 <- scrambled[[1]][!(scrambled[[1]]$id %in% c(r1_not_r2, linear_1$id)), ]

  genomic_mate <- genomic[[2]][(genomic[[2]]$id %in% c(linear_1$id,scrambled_1$id)), ]
  linear_mate <- linear[[2]][(linear[[2]]$id %in% c(linear_1$id,scrambled_1$id)), ]
  scrambled_mate <- scrambled[[2]][(scrambled[[2]]$id %in% c(linear_1$id, scrambled_1$id)), ]

  linear_genomic <- left_join(linear_1,
                              genomic_mate,
                              by = "id",
                              suffix = suffixes[c(1, 2)])

  linear_linear <- left_join(linear_1,
                             linear_mate ,
                             by = "id",
                             suffix = suffixes[c(1, 3)])

  linear_scrambled <- left_join(linear_1,
                                scrambled_mate,
                                by = "id",
                                suffix = suffixes[c(1, 4)])

  scrambled_genomic <- left_join(scrambled_1,
                                 genomic_mate,
                                 by = "id",
                                 suffix = suffixes[c(1, 2)])

  scrambled_linear <- left_join(scrambled_1,
                                linear_mate ,
                                by = "id",
                                suffix = suffixes[c(1, 3)])

  scrambled_scrambled <- left_join(scrambled_1,
                                   scrambled_mate,
                                   by = "id",
                                   suffix = suffixes[c(1, 4)])

  linear_df <- cbind(linear_linear[, 1:15],
                     linear_genomic[, 16:22],
                     linear_linear[, 16:29],
                     linear_scrambled[, 16:29])

  scrambled_df <- cbind(scrambled_linear[, 1:15],
                        scrambled_genomic[, 16:22],
                        scrambled_linear[, 16:29],
                        scrambled_scrambled[, 16:29])

  return(list(linear_df = linear_df, scrambled_df = scrambled_df))
}
R1R2MatesJunctionReads <- createJunctionsReadsDF(genomic_filtered_dt,
                                                 linear_junctions_filtered_dt,
                                                 scrambled_junctions_filtered_dt,
                                                 reads_to_ignore$r1_not_r2)

5. Classify Mate Reads

Here we classify mate reads based on alignment score. We also create a new variable ascore that averages the alignment scores of the two matched reads.

classsifyMateReads <- function(df, suffixes = c(".R1", ".R2g", ".R2l", ".R2s")) {

    scores_df <- df[, paste0("AS", suffixes)]
    nMates <- apply(scores_df[, 2:4] , 1, function(x) sum(!is.na(x)))
    max_score <- apply(scores_df[, 2:4] , 1, which.max)

    mate <- sapply(max_score, function(x)
      if (is.null(attributes(x)))
        {return(0)}
      else
        {return(as.integer(x))}
      )

    matescore <- scores_df[cbind(1:nrow(scores_df), (mate + 1))]
    ASmean <- 0.5 * (scores_df$AS.R1 + matescore)
    df <- cbind(df,
                ASmean = ASmean,
                nMates = nMates,
                mate = mate)
    return(df)
  }
linear_df <-
  R1R2MatesJunctionReads$linear_df %>%
  classsifyMateReads()
scrambled_df <-
  R1R2MatesJunctionReads$scrambled_df %>%
  classsifyMateReads()
cat( " number of linear junction reads",
     NROW(linear_df), "\n",
     "number of scrambled junction reads",
     NROW(scrambled_df), "\n",
  "total number of junctional reads:",
    NROW(linear_df) + NROW(scrambled_df), "\n",
    "number of junctions with aligned reads:",
    length(unique(linear_df$junction.R1)) +
      length(unique(scrambled_df$junction.R1)))

 number of linear junction reads 600658
 number of scrambled junction reads 1599
 total number of junctional reads: 602257
 number of junctions with aligned reads: 108949

Let’s take a look at the dataframe for the linear junctions R1 aligned reads

scrambled_df[1:3, ]  

Let’s also output the frequency of R1 reads that have 0, 1, 2, or 3 potential paired-end mate reads (R2)

linear_junction_mapcateg <-
  linear_df$nMates %>%
  table  %>%
  as.data.frame()
scrambled_junction_mapcateg <-
  scrambled_df$nMates %>%
  table %>%
  as.data.frame()
names_df <- c("Unmapped","OneMateMatch", "TwoMateMatchs", "ThreeMateMatchs")
rownames(linear_junction_mapcateg) <- names_df
rownames(scrambled_junction_mapcateg) <- names_df
linear_junction_mapcateg 

scrambled_junction_mapcateg 

6. Classify Junction Reads

We split the reads into reads that have a mapped mate R2 and those that do not. We intially assign all reads as uncategorized (i.e. class=0).

linear_reads_split <- split(linear_df, linear_df$nMates == 0)
scrambled_reads_split <- split(scrambled_df, scrambled_df$nMates == 0)
unmapped_reads <- rbind(linear_reads_split[[2]], scrambled_reads_split[[2]])[,1:15]
linear_reads_mapped <- cbind(linear_reads_split[[1]] , class = 0)
scrambled_reads_mapped <- cbind(scrambled_reads_split[[1]] , class = 0)

Let’s take a look at the unmapped reads

unmapped_reads[1:3,]

Now we proceed to classify scrambled junctional reads as circular or decoy and linear junctional reads as linear or anomalous.

classifyCircularPairedMates <- function(x) {

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- x$mate == 3
  cond4 <- x$junction.R1 == x$junction.R2s
  x$class[(cond1 & cond2 & cond3 & cond4)] <- 1 #(circ)

  return(x)
}
classifyCircularLinearMates <- function(x, BUFFER = 15, MIDPOINT = 150) {

  r1_start <- pmin(x$jPos1.R1, x$jPos2.R1)
  r1_end <- pmax(x$jPos1.R1, x$jPos2.R1)
  r2_start <- pmin(x$jPos1.R2l, x$jPos2.R2l)  - MIDPOINT +  x$pos.R2l
  r2_end <- pmin(x$jPos1.R2l, x$jPos2.R2l)  - MIDPOINT +  x$pos.R2l +  x$len.R2l - 1

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- r2_start >=  r1_start - BUFFER
  cond4 <- r2_start <=  r1_end + BUFFER
  cond5 <- r2_end >= r1_start - BUFFER
  cond6 <- r2_end <= r1_end + BUFFER
  cond7 <- x$mate == 2
  cond8 <- !(x$class %in% c(1))

  x$class[(cond1 & cond2 & cond3 & cond4 & cond5 & cond6 & cond7 & cond8)] <- 2 #(circ_lin)
  return(x)
}
classifyCircularGenomicMates <- function(x, BUFFER = 15, MIDPOINT = 150) {

  r1_start <- pmin(x$pos.R1, x$jPos2.R1)
  r1_end <- pmax(x$jPos1.R1, x$jPos2.R1)
  r2_start <- x$pos.R2g
  r2_end <- x$pos.R2g + x$len.R2g - 1

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- r2_start >=  r1_start - BUFFER
  cond4 <- r2_start <=  r1_end + BUFFER
  cond5 <- r2_end >= r1_start - BUFFER
  cond6 <- r2_end <= r1_end + BUFFER
  cond7 <- x$mate == 1
  cond8 <- !(x$class %in% c(1, 2))

  x$class[(cond1 & cond2 & cond3 & cond4 & cond5 & cond6 & cond7 & cond8 )] <- 3 #(circ_gen)
  return(x)
}
classifyCircularIgnore <- function(x, BUFFER = 15, MIDPOINT = 150) {

  r1_start <- pmin(x$jPos1.R1, x$jPos2.R1) - MIDPOINT + x$pos.R1
  r1_end <- r1_start + x$len.R1 - 1

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- r1_start >= pmin(x$jPos1.R2s, x$jPos2.R2s) - BUFFER
  cond4 <- r1_start <= pmax(x$jPos1.R2s, x$jPos2.R2s) + BUFFER
  cond5 <- r1_end >= pmin(x$jPos1.R2s, x$jPos2.R2s) - BUFFER
  cond6 <- r1_end <= pmax(x$jPos1.R2s, x$jPos2.R2s) + BUFFER
  cond7 <- x$mate == 3
  cond8 <- !(x$class %in% c(1, 2, 3))
  x$class[(cond1 & cond2 & cond3 & cond4 & cond5 & cond6 & cond7 & cond8)] <- 4
  return(x)
}
classifyCircularDecoy <- function(x) {
  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- !(x$class %in% c(1, 2, 3, 4))
  cond4 <- x$mate == 3
  x$class[(cond1 & cond2 & cond3) ] <- 5
  x$class[!(cond1 & cond2) ] <- 5

  return(x)
}
classifyLinearPairedMates <- function(x, BUFFER = 15, MIDPOINT = 150) {

  r1_start <- pmin(x$jPos1.R1, x$jPos2.R1) - MIDPOINT + x$pos.R1 # myCoord
  r2_start <- pmin(x$jPos1.R2l, x$jPos2.R2l)  - MIDPOINT +  x$pos.R2l #mateCoord

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3_p <- r2_start >= r1_start - BUFFER
  cond3_m <- r1_start >= r2_start - BUFFER
  cond3 <- ifelse(x$str.R1 == "+", cond3_p, cond3_m)
  cond4 <- x$mate == 2
  x$class[(cond1 & cond2 & cond3 & cond4)] <- 1
  return(x)

}
classifyLinearGenomicMates <- function(x, BUFFER = 15, MIDPOINT = 150) {

  r1_start <- pmin(x$jPos1.R1, x$jPos2.R1) - MIDPOINT + x$pos.R1
  r2_start <- x$pos.R2g

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3_p <- r2_start >= r1_start - BUFFER
  cond3_m <- r1_start >= r2_start - BUFFER
  cond3 <- ifelse(x$str.R1 == "+", cond3_p, cond3_m)
  cond4 <- x$mate == 1
  cond5 <- !(x$class %in% c(1))
  x$class[(cond1 & cond2 & cond3 & cond4 & cond5)] <- 2
  return(x)
}
classifyLinearIgnore <- function(x) {

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  cond3 <- x$mate == 3
  x$class[(cond1 & cond2 & cond3 )] <- 3
  return(x)
}
classifyLinearAnomaly <- function(x) {

  cond1 <- x$str.R1 != x[, c("str.R2g", "str.R2l", "str.R2s")][cbind(1:nrow(x), x$mate)]
  cond2 <- x$jChr.R1 == x[, c("junction.R2g", "jChr.R2l", "jChr.R2s")][cbind(1:nrow(x), x$mate)]
  x$class[!(cond1 & cond2)] <- 4
  x$class[x$class == 0] <- 4
  return(x)
}
linear_reads_mapped <-
  linear_reads_mapped %>%
  classifyLinearPairedMates() %>%
  classifyLinearGenomicMates() %>%
  classifyLinearIgnore()  %>%
  classifyLinearAnomaly()
scrambled_reads_mapped <-
  scrambled_reads_mapped %>%
  classifyCircularPairedMates() %>%
  classifyCircularLinearMates() %>%
  classifyCircularGenomicMates  %>%
  classifyCircularIgnore  %>%
  classifyCircularDecoy
linear_reads_mapped[]

The frequency of read categories are as follows:

linear_junction_reads <-
  linear_reads_mapped$class %>%
  tabulate  %>%
  as.data.frame()
scrambled_junction_reads <-
  scrambled_reads_mapped$class %>%
  tabulate %>%
  as.data.frame()
rownames(linear_junction_reads) <- c("Linear.l",
                                     "Linear.g",
                                     "Ignore",
                                     "Anomaly")
rownames(scrambled_junction_reads) <- c("Circular.s",
                                        "Circular.l",
                                        "Circular.g",
                                        "Ignore",
                                        "Decoy")
linear_junction_reads 

scrambled_junction_reads

cat( " number of linear reads:",
     sum(linear_junction_reads[1:2,]), "\n",
     "number of circular reads:",
     sum(scrambled_junction_reads[1:3,]), "\n",
  "total number of anomaly reads:",
    sum(linear_junction_reads[4,]), "\n",
    "number of  decoy reads:",
    sum(scrambled_junction_reads[5,]))

 number of linear reads: 540360
 number of circular reads: 1265
 total number of anomaly reads: 14474
 number of  decoy reads: 117

Reorder columns

refcols <- c("id","class", "mate")
linear_reads_mapped <- linear_reads_mapped[, c(refcols,
                                               setdiff(names(linear_reads_mapped),
                                                       refcols))]
scrambled_reads_mapped <- scrambled_reads_mapped[, c(refcols,
                                                     setdiff(names(scrambled_reads_mapped),
                                                             refcols))]
scrambled_reads_mapped

Save data

write.table(linear_reads_mapped,
            "linear_reads_mapped.txt",
            row.names = FALSE,
            sep = "\t")
write.table(scrambled_reads_mapped,
            "scrambled_reads_mapped.txt",
            row.names = FALSE,
            sep = "\t")
saveRDS(scrambled_reads_mapped, "scrambled_reads_mapped.rds")
saveRDS(linear_reads_mapped, "linear_reads_mapped.rds")
#scrambled_reads_mapped <- readRDS("scrambled_reads_mapped.rds")
#linear_reads_mapped <- readRDS("linear_reads_mapped.rds")

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] dplyr_0.5.0                systemPipeR_1.9.2          ShortRead_1.32.0           GenomicAlignments_1.10.0   SummarizedExperiment_1.4.0 Biobase_2.34.0
 [7] BiocParallel_1.8.1         Rsamtools_1.26.1           Biostrings_2.42.1          XVector_0.14.0             GenomicRanges_1.26.1       GenomeInfoDb_1.10.2
[13] IRanges_2.8.1              S4Vectors_0.12.1           BiocGenerics_0.20.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8            locfit_1.5-9.1         lattice_0.20-34        GO.db_3.4.0            assertthat_0.1         digest_0.6.10          R6_2.2.0
 [8] plyr_1.8.4             BatchJobs_1.6          backports_1.0.4        RSQLite_1.1-1          ggplot2_2.2.1          zlibbioc_1.20.0        GenomicFeatures_1.26.2
[15] lazyeval_0.2.0         annotate_1.52.1        Matrix_1.2-7.1         checkmate_1.8.2        GOstats_2.40.0         splines_3.3.2          stringr_1.1.0
[22] pheatmap_1.0.8         RCurl_1.95-4.8         biomaRt_2.30.0         munsell_0.4.3          sendmailR_1.2-1        rtracklayer_1.34.1     base64enc_0.1-3
[29] BBmisc_1.10            fail_1.3               tibble_1.2             edgeR_3.16.5           XML_3.98-1.5           AnnotationForge_1.16.0 bitops_1.0-6
[36] grid_3.3.2             RBGL_1.50.0            xtable_1.8-2           GSEABase_1.36.0        gtable_0.2.0           DBI_0.5-1              magrittr_1.5
[43] scales_0.4.1           graph_1.52.0           stringi_1.1.2          hwriter_1.3.2          genefilter_1.56.0      limma_3.30.7           latticeExtra_0.6-28
[50] brew_1.0-6             rjson_0.2.15           RColorBrewer_1.1-2     tools_3.3.2            Category_2.40.0        survival_2.40-1        AnnotationDbi_1.36.0
[57] colorspace_1.3-2       memoise_1.0.0          knitr_1.15.1          

References

Szabo, Linda, Robert Morey, Nathan J Palpant, Peter L Wang, Nastaran Afari, Chuan Jiang, Mana M Parast, Charles E Murry, Louise C Laurent, and Julia Salzman. 2015. “Statistically based splicing detection reveals neural enrichment and tissue-specific induction of circular RNA during human fetal development.” Genome Biology 16. Genome Biology: 126. doi:10.1186/s13059-015-0690-5.