A pipeline for de novo circRNA identification
Switch branches/tags
Nothing to show
Clone or download
Latest commit 9651cda Feb 16, 2017
Permalink
Failed to load latest commit information.
CB_splice acfs v2.0 Aug 4, 2014
ACF_MAKE.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step1.pl v2.1 May 6, 2016
ACF_Step2.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step2_MuSeg.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step3.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step3_MuSeg.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step4.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step4_CBR.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step4_MEA.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_Step5.pl v2.1 May 18, 2016
ACF_Step5m.pl acfs v2.0 Aug 4, 2014
ACF_Step5m2.pl V2.2 Oct 11, 2016
ACF_Step6.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_fusion_circRNAs.pl v2.1 May 13, 2016
ACF_trans_splice_step1.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_trans_splice_step2.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
ACF_trans_splice_step3.pl v2.1 May 17, 2016
BASH_example.sh Changes to be committed: May 28, 2015
LICENSE acfs v2.0 Aug 4, 2014
README.md minor changes Feb 16, 2017
SPEC_example.txt acfs v2.0 Aug 4, 2014
Truseq_merge_unique_fa.pl minor changes Feb 16, 2017
Truseq_merge_unique_filelist.pl minor changes Feb 16, 2017
bwa-0.7.3a.tar Include bwa backage. Support for higher versions of bwa will be added… Mar 4, 2015
change_fastq_header.pl C Jan 27, 2016
convert_unmapped_SAM_to_fa_for_acfs.pl V2.2 Oct 11, 2016
get_bed12_from_refFlat.pl C Jan 27, 2016
get_circRNA_from_bed.pl V2.2.1 Can handel chromosome names with both Ensembl (1.fa,MT.fa) and… Dec 25, 2016
get_id_from_simulate_fusion_circRNA_gtf.pl v2.1 May 9, 2016
get_linear_id_from_psl.pl v2.1 May 13, 2016
get_pseudo_circle.pl v2.1 May 13, 2016
get_selected_fa_from_pool.pl v2.1 May 9, 2016
get_selected_from_pool_singleline.pl v2.1 May 9, 2016
get_separate_unmap_sam.pl V2.2 Oct 11, 2016
get_seq_from_agtf.pl V2.2 Oct 11, 2016
get_split_exon_border_biotype_genename.pl V2.2 Oct 11, 2016
reverse_complement.pl minor changes Feb 16, 2017
simulate_PE_reads_from_circRNA.pl v2.1 May 13, 2016
simulate_PE_reads_from_linear.pl v2.1 May 13, 2016
simulate_SE_reads_from_circRNA.pl v2.1 May 13, 2016
simulate_SE_reads_from_linear.pl v2.1 May 13, 2016
simulate_gtf_for_circRNA.pl v2.1 May 18, 2016
simulate_gtf_for_fusion_circRNA.pl v2.1 May 9, 2016
simulate_reads_for_fusion_circRNA.pl v2.1 May 13, 2016
test.tar.gz minor changes Feb 16, 2017

README.md

ACFS

Accurate CircRNA Finder Suite. Discovering circRNAs from RNA-Seq data.

Overview

CircRNAs are generated through splicing, or to be precise, back-splicing where the downstream splice donor attacks an upstream splice acceptor. Identifying the exact site of back-splice lies in the heart of circRNA discovery.

ACFS first examines and pinpoints the back-splice site from RNA-Seq alignment using a maximal entropy model. The expression of circRNAs is estimated from a second round of alignment to the inferred pseudo circular sequences.

No prior knowledge of gene annotation is needed for circRNA prediection, but reading the coordinates is far less interesting than reading gene names, so circRNAs are annotated using the gene annotation if provided. Also, given annotation, circRNA sequences can be better estimated (especially those contains short exons) and enhance the accuracy of circRNA expression quantification.

ACFS is designed for Single-end RNA-Seq reads. Seeing is believing, we would like to see read directly spanning the back-splice sites. Paired-end data is also supported, albeit with lower sensitivity (if neither of the read ends crosses the back-splice, but the read-pair does).

Change Log

  • Update on 2016-05-06 : Added extension for fining trans-splicing evidence, which can be used to identify fusion-circRNAs
  • Update on 2016-01-27 : Performance improvement and add scripts for simulation
  • Update on 2015-09-17 : Added a small real-world example
  • Update on 2015-08-20 : Added support for paired-end reads
  • Update on 2015-08-11 : Corrected the Tutorial section in README, thanks to Zol
  • Update on 2015-03-09 : Now ACF can include pre-defined circRNA annotations from a bed6 or bed12 file (and their authenticity will be checked, so please adjust (minJump, maxJump, minSplicingScore) accordingly ). This way, you can both predict novel circRNAs in your data and estimate the abundance of annotated circRNAs.
  • First release on 2013-11-01

Installation

Simply unpack the ACFS package.

Requirement

  • bwa-0.7.3a (included in the package, but you need to do "make")
  • perl
  • blat (not necessary, sometime it helps to rule out false positive fusion-circRNAs when there is a gene-duplication or gene-pseudogene dilemma)

Pipeline scheme

  1. Map all Tophat2-unmapped-reads to genome using BWA, seperate :
    • 1-part
    • 2-part-same-chromosome-same-strand (contain circRNAs)
    • 2-part-same-chromosome-diff-strand (possibly PolII backwalk)
    • >2-part-same-chromosome (contain circRNAs, if they are originated from short exons)
    • >=2-part-diff-chromosome (contain trans-splicing, or even fusion-circRNAs)
  2. Estimate splicing strength using Christoph Berge's method for "2-part-same-chromosome-same-strand". report:
    • forward-splice (not interesting)
    • back-splice and canonical splice-motif {[GT-AT],[GC-AG],[AT-AC]}, calculate strength using MaxEnt
    • back-splice and non-canonical splice-motif, find the maximal score using MaxEnt and the corresponding back-splice site
    • try to rescue circles using ">2-part-same-chromosome"
  3. Check if both splice sites of the cirRNAs are on known exon-borders if an annotation is provided; otherwise all are CBR. report:
    • both known exon-border (termed : MEA . M atch with E xisting A nnotation, which is somewhat more reliable)
    • at least one splice site sits on unknown exon-border (termed : CBR . C hris B erge R escued)
  4. Build gtf and pseudo-transcript for results from step3
  5. Map all unmapped-reads to pseudo-transcipts and estimate the abundance of circRNAs
  6. Generate bed track for visulization
  7. Optional search of trans-splicing and fusion-circRNA events

Before running ACFS, a few pre-process

  1. Map the RNA-Seq reads to genome and transcriptom, and extract the unmapped reads. This is recommended as those mapped reads will NOT span the back-splice sites, and therefore do NOT contribute to circRNA discovery.

  2. Change fasta/fastq header format to allow processing multiple samples in one run. This is IMPORTANT ! ACFS expects a special header format so that multiple samples can be processed in one run. Do change the default header such as ">HWUSI-EAS100R:6:73:941:1973" into ">Truseq_sample1_HWUSI-EAS100R:6:73:941:1973", where the "sample1" is the name of your choice describing the sample. Do the conversion as:

    perl change_fastq_header.pl SRR650317_1.fasta SRR650317_1.fa Truseq_SRR650317left
    perl change_fastq_header.pl SRR650317_2.fasta SRR650317_2.fa Truseq_SRR650317right
    

    Make sure there is No underline within the sample name. e.g. ">Truseq_ctrl.1_Default_header" and ">Truseq_AGO2KO.1_Default_header" are OK; ">Truseq_ctrl_1_Default_header" is BAD because ACF will register "ctrl" as the sample name instead of "ctrl_1".

  3. Merge sequences from multiple fasta/fastq files into one fasta file, which saves time for mapping.

    perl Truseq_merge_unique_fa.pl UNMAP SRR650317_1.fa SRR650317_2.fa
    

    Alternatively, if there are many files to merge, generate a file (named filelist for example) contains the full-path of each file

    perl Truseq_merge_unique_filelist.pl <UNMAP> <filelist>
    

    UNMAP is the collasped fasta file which will be processed by ACFS, and UNMAP_expr contains the readcount per sequence in all the samples. If you change the name of UNMAP to SomethingElse, then the readcount file will be automatically named as SomethingElse_expr

    However, one can bypass the previous and this step to run ACFS sample by sample. This way, no fasta header reformatting and reads collapsing is needed. For each sample, set the value of UNMAP to the name of fasta/fastq in the SPEC_example.txt file, and set the value of UNMAP_expr to "no".

  4. Build BWA index, using verion 0.73a (currently not support for other versions as the output format changes between versions)

    /bin/bwa073a/bwa index /data/iGenome/human/Ensembl/GRCh37/Sequence/BWAIndex/genome.fa
    
  5. Prepare for annotation (recommended) Download the gtf file from iGenome package or download ensembl gtf here : ftp://ftp.ensembl.org/pub/current_gtf/ Then run:

    perl get_split_exon_border_biotype_genename.pl </data/iGenome/human/Ensembl/GRCh37/Annotation/Genes/genes.gtf> </data/iGenome/human/Ensembl/GRCh37/Annotation/Genes/Homo_sapiens.GRCh37.71_split_exon.gtf>
    

    The first argument is the input gtf file, the second argument is the output

  6. Strandedness is assumed as from the Truseq Stranded RNA-Seq, so the reads are reverse-complementary to mRNAs.

    • If the reads are actually sense (the same 5'->3' direction as mRNA), please reverse-complement all reads.
    • If the reads are actually stransless or pair-ended, please run in parallel original reads and reverse-complemented reads.
    • No pair-end information is used, as the exact junction-site must be supported by a single read. Seeing is believing.
  7. For Paired-end data, it is highly recommended to align the reads to genome+transcriptome first (e.g. using Tophat2), and extract the unmapped read (read pairs) using the following script (the fasta header line is also modified)

    perl convert_unmapped_SAM_to_fa_for_acfs.pl <output_file_name> <unmapped.sam> sample_id
    

Parameters

There are nine mandatory parameters to run ACFS in a basic mode. Searching for fusion-circRNAs is disabled by default. Please modify the config file "SPEC_example.txt" according to your specific organism, experimental design and sequence specs. The config file "SPEC_example.txt" is a two-column tab-delimited file : <name>\t<value>

Mandatory paramters:

Parameter value Note
BWA_folder /home/bin/bwa037a/ path of the folder of bwa
BWA_genome_Index /data/.../BWAIndex/genome.fa full path to the index files
BWA_genome_folder /data/.../Chromosomes/ full path to the folder containing separeted chromosome files
ACF_folder /home/bin/ACFS/ path of the folder of ACFS
CBR_folder /home/bin/ACFS/CB_splice/ path of the folder for MaxEnt
Agtf /data/.../Homo_sapiens.GRCh37.71_split_exon.gtf processed annotation, see the previous section point-4
UNMAP UNMAP the collapsed fasta file
UNMAP_expr UNMAP_expr the expression of the collasped reads
Seq_len 150 length of sequencing reads

Optional parameters, the values in below are set as default:

Parameter value Note
Thread 16 number of threads used in bwa
BWA_seed_length 16 bwa seed length
BWA_min_score 20 bwa min score to trigger report. For shorter reads, e.g. 50nt, set to 10 or lower could report more circRNAs at risk of higher FDR
minJump 100 the minimum distance of a back-splice. The smaller, the more likely you can find circles from short exons
maxJump 2500000 the maximum distance of a back-splice. The larger, the more likely you can find circles from long genes. The longest human gene is CNTNAP2 which spans 2.3M bp
minSplicingScore 10 the minimum score for the sum of splicing strength at both splice site, 10 corresponds to 95% of all human/mouse splice site pairs. One could also set it to a lower value, e.g. zero, and do a post-filtering after running acfs
minSampleCnt 1 the minimum number of samples that detect any given circle
minReadCnt 1 the minimum number of reads (from all samples) that detect any given circle
minMappingQuality 20 the minimum mapping quality of any given sequence
minSpanJunc 6 the minimum number of bases reach beyond the back-splice-site. The larger the less likely of false-positive
Coverage 0.9 the minimum percentage of any given read is aligned. The larger the more conserved the results are
ErrorRate 0.05 the maximum error rate for re-alignment. The smaller the better
Strandness - the strand information of sequencing, must be one of the {+, -, no}
pre_defined_circle_bed no pre-defined circle annotation in bed6 or bed12 format (to increase sensitivity for lowly expressed circRNAs please include bed12 files of annotated one, e.g. merge the bed files from GSE61991)
Search_trans_splicing no set to "yes" to seach for trans-splicing reads
blat_search yes use blat to discard false positives results from gene duplication, turn off by "no"
blat_path blat full path to the executable, such as "/usr/bin/blat/blat", it is ignored if the blat_search option is set to no
trans_splicing_coverage 0.9 see Coverage
trans_splicing_minMappingQuality 0 see minMappingQuality
trans_splicing_minSplicingScore 10 see minSplicingScore
trans_splicing_maxSpan 2500000 the maximum distance between the junctions on the same gene for fusion circRNAs

run ACFS

  1. make Pipeline Bash file
    perl ACF_MAKE.pl <SPEC_example.txt> <BASH_example.sh>
    
  2. find circles
    bash <BASH_example.sh>
    

Results

  1. circRNAs are stored in two files, which can be visualzed using UCSC Genome Browser:

    • circle_candidates_MEA.bed12
    • circle_candidates_CBR.bed12
      The higher the value in 5th column, the more "likely" that circle is true.
      The name in 4th column can be seperated by "|" into four segments: circle-ID, sum-of-Splicing-Strength, number-of-supporting-samples, number-of-supporting-reads.
      The sequence for the "middle exons" in almost every circle is hypothetically filled in using the annotations provided, true sequences could be determined by integrating mRNA-Seq data and validated by inward and outward PCRs.
  2. The circRNA expression is stored here:

    • circle_candidates_expr (merged file of <circle_candidates_MEA.expr> and <circle_candidates_CBR.expr>) The second column denote the name of the gene from which this circRNA is derived.
  3. Fusion-circRNAs, if enabled:

    • fusion_circRNAs
      The Junctional sequences are stored in "unmap.trans.splicing.tsloci.fa".

Tutorial

  1. We provide a toy example in the test folder after unpacking test.tar.gz. One can run the pipeline in few minutes.

  2. pre-processing for sequencing reads

    1A. for single-end RNA-Seq Raw data only. Note this sample is from a stranded RNA-Seq experiment from mouse, therefore mouse annotation should be used. (Feasible but not a good idea in practice, since you don't want to map all reads again. Use unmapped reads instead.)

    wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX852%2FSRX852583/SRR1772422/SRR1772422.sra  
    fastq-dump.2 --fasta 0 --split-files SRR1772422.sra  
    perl change_fastq_header.pl SRR1772422.fasta SRR1772422.fa Truseq_HippoSyn  
    perl Truseq_merge_unique_fa.pl UNMAP SRR1772422.fa
    

    1B. for paired-end RNA-Seq Raw data only. Note this sample is from an unstranded RNA-Seq experiment from human, therefore human annotation should be used. (Feasible but not a good idea in practice, since you don't want to map all reads again. Use unmapped reads instead.) As the RNA-Seq was un-stranded, we need to make reverse complement for all reads. For paired-end stranded RNA-Seq data, we need to make reverse complement for Read-2 if the Read-1 is antisense to mRNA (dUTP, NSR, NNSR protocol).

    wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX218%2FSRX218203/SRR650317/SRR650317.sra  
    fastq-dump.2 --fasta 0 --split-files SRR650317.sra
    perl reverse_complement.pl SRR650317_1.fasta
    perl reverse_complement.pl SRR650317_2.fasta
    perl change_fastq_header.pl SRR650317_1.fasta SRR650317_1.fa Truseq_SRR650317left  
    perl change_fastq_header.pl SRR650317_2.fasta SRR650317_2.fa Truseq_SRR650317right  
    perl change_fastq_header.pl SRR650317_1.fasta.rc SRR650317_1.fa.rc Truseq_SRR650317left  
    perl change_fastq_header.pl SRR650317_2.fasta.rc SRR650317_2.fa.rc Truseq_SRR650317right
    perl Truseq_merge_unique_fa.pl UNMAP SRR650317_1.fa SRR650317_1.fa.rc SRR650317_2.fa SRR650317_2.fa.rc
    

    1C. for single-end and paired-end RNA-Seq unmapped reads.
    If you have already aligned raw RNA-Seq reads to some reference (using say Bowtie/BWA/STAR/Tophat2), obtain a SAM file unmapped.sam containing all the unmapped reads. For example if you used Tophat2, simply convert the file unmapped.bam to unmapped.sam and then run the following command:

    perl convert_unmapped_SAM_to_fa_for_acfs.pl <USER_converted.fa> <unmapped.sam> sample_id
    

    where sample_id is the new sample ID that is added to the fasta/fastq header according to ACFS style. It could be "Ctrl12" or "TreatA.rep2", your choice as long as there is NO underscore as "_" in it.
    USER_converted.fa is the name of the output file. Then collapse the reads using the following command:

    perl Truseq_merge_unique_fa.pl <UNMAP> <USER_converted.fa>
    

    Or if you have more than one samples:

    perl Truseq_merge_unique_fa.pl <UNMAP> <USER_converted_1.fa> <USER_converted_2.fa> <USER_converted_3.fa>  ...
    
  3. prepare for annotation, if you haven't done so

    perl get_split_exon_border_biotype_genename.pl /data/iGenome/human/Ensembl/GRCh37/Annotation/Genes/genes.gtf /data/iGenome/human/Ensembl/GRCh37/Annotation/Genes/Homo_sapiens.GRCh37.71_split_exon.gtf
    
  4. modify the config file SPEC_example.txt accordingly.

  5. generate ACFS pipeline

    perl ACF_MAKE.pl <SPEC_example.txt> <BASH_example.sh>
    
  6. run ACFS

    nohup bash <BASH_example.sh> &
    
  7. take a look at the results when finished :

    • circle_candidates_MEA.bed12 # circRNAs back-splice at annotated boundarys of exon(s)
    • circle_candidates_CBR.bed12 # circRNAs back-splice at un-annotated boundarys of exon(s)
      And the expression(readcounts) table for circRNAs, with circRNAs in rows and samples in columns
    • circle_candidates.expr
      And the potential fusion circRNAs
    • fusion_circRNAs

A few useful scripts for simulation

To see the usage, simply run the perl scripts with no arguments.

  1. simulating SE reads from linear transcripts

    simulate_SE_reads_from_linear.pl
    
  2. simulating PE reads from linear transcripts

    simulate_PE_reads_from_linear.pl
    
  3. simulating circRNAs

    simulate_gtf_for_circRNA.pl
    get_split_exon_border_biotype_genename.pl
    get_seq_from_agtf.pl
    
  4. simulating SE reads from circRNAs

    simulate_SE_reads_from_circRNA.pl
    
  5. simulating PE reads from circRNAs

    simulate_PE_reads_from_circRNA.pl
    
  6. simulating fusion-circRNAs

    simulate_gtf_for_fusion_circRNA.pl
    get_split_exon_border_biotype_genename.pl
    get_seq_from_agtf.pl
    get_id_from_simulate_fusion_circRNA_gtf.pl
    
  7. simulating SE reads from fusion-circRNAs (or PE reads where the insert-size is the read length)

    simulate_reads_for_fusion_circRNA.pl
    

Contact

This pipeline is developed and is maintained by Arthur Xintian You: arthur.yxt@gmail.com. I will do my best to respond in a timely manner.

Cite

Nat Neurosci 2015 Apr;18(4):603-10 [PMID: 25714049]