1. find_circ鉴定circRNA的原理
find_circ 的基本原理: find_circ根据Bowtie2比对结果，从没有比对到参考序列的 reads 的两端各提取 20nt 的 anchor 序列，将每一对 anchor 序列再次与参考序列比对。如果 anchor 序列的 5' 端比对到参考序列（起始与终止位点分别记为 A3，A4），anchor 序列的 3' 端比对到此位点的上游（起始与终止位点分别记为 A1，A2），并且在参考序列的 A2 到 A3 之间存在剪接位点（GT-AG），则将此 read 作为候选 circRNA。最后将 read count 大于等于 2 的候选 circRNA 作为鉴定的 circRNA。
2. find_circ的安装
find_circ需要运行在装有python 2.7的64位系统上，同时需要安装numpy和pysam这两个python模块。其运行需要借助bowtie2和samtools来完成基因组mapping的过程。
1 wget https://github.com/marvin-jens/find_circ/archive/v1.2.tar.gz
2 tar -xzvf v1.2.tar.gz
3. 参考基因组的下载
通过fetch_ucsc.py下载ucsc最新版本的参考基因组
1 fetch_ucsc.py hg19/hg38/mm9/mm10 ref/kg/ens/fa out
4. bowtie2建立参考基因组索引
1 bowtie2_build hg38.fa hg38
5. 基于RNA-Seq的基因组比对（pair-end模式）
bowtie2参数介绍
-p 使用多线程
--very-sensitive 允许多重比对，报告出最好的一个
--score-min=C,-15,0 设置比对分数函数
--mm 设置I/O模式。
samtools view参数介绍
-h 文件包含header line
-b 输出bam格式
-u 输出非压缩的bam格式
–S 忽略版本兼容
1 bowtie2 -p 16 --very-sensitive --score-min=C,-15,0 --mm -x /path/to/bowtie2_index -q -1 reads1.fq -2 reads2.fq | samtools view -hbuS - | samtools sort - -o output.bam
6. 挑出没有比对上的序列，各取两头20bp短序列（anchor)
1 samtools view -hf 4 output.bam | samtools view -Sb - > unmapped.bam
2 /path/to/unmapped2anchors.py unmapped.bam | gzip > anchors.fq.gz
7.根据anchor比对基因组情况寻找潜在的circRNA
find_circ.py参数介绍
--prefix参数指定的是spliced_sites.bed文件中第四列内容的前缀，建议指定为物种对应的三字母缩写，需要注意的是，在spliced_sites_bed中同时包含了环状RNA和线性RNA,环状RNA的名称用circ标识，线性RNA的名称用norm标识，这里设置为--prefix=hsa_
--name参数会在生成的spliced_sites.bed文件中指定tissues列的名字
--reads参数会生成包含spliced reads的fa文件
--stats参数会生成包含数值统计信息的txt文件
1 bowtie2 -p 16 --reorder --mm --score-min=C,-15,0 -q -x /path/to/bowtie2_index -U anchors.fq.gz | /path/to/find_circ.py --genome=/path/to/hg38.fa --prefix=hsa_ --name=my_test_sample --stats=<run folder>/stats.txt --reads=<run folder>/splice_reads.fa > <run folder>/spliced_sites.bed
根据以下规则对结果进行筛选
1.根据关键词CIRCULAR筛选环状RNA
2.去除线粒体上的环状RNA
3.筛选unique junction reads数至少为2的环状RNA
4.去除断裂点不明确的环状RNA
5.过滤掉长度大于100kb的circRNA,这里的100kb为基因组长度，直接用环状RNA的头尾相减即可
1 grep CIRCULAR spliced_sites.bed | grep -v chrM | gawk '$5>=2' | grep UNAMBIGUOUS_BP | grep ANCHOR_UNIQUE | /path/to/maxlength.py 100000 > find_circ.candidates.bed
8.分析多个样本
如果有多个样本，需要分别用find_circ.py运行，然后将各自的结果合并
1 /path/to/merge_bed.py sample1.bed sample2.bed [...] >combined.bed
9. 输出的spliced_sites_bed文件格式
输出的spliced_sites.bed文件前六列为标准的BED文件格式，剩余的12列关于junction的一些信息
| column | name | description |
|---|---|---|
| 1 | chrom | chromosome/contig name |
| 2 | start | left splice site (zero-based) |
| 3 | end | right splice site (zero-based).(Always: end > start. 5' 3' depends on strand) |
| 4 | name | (provisional) running number/name assigned to junction |
| 5 | n_reads | number of reads supporting the junction (BED 'score') |
| 6 | strand | genomic strand (+ or -) |
| 7 | n_uniq | number of distinct read sequences supporting the junction |
| 8 | uniq_bridges | number of reads with both anchors aligning uniquely |
| 9 | best_qual_left | alignment score margin of the best anchor alignment supporting the left splice junction (max=2 * anchor_length) |
| 10 | best_qual_right | same for the right splice site |
| 11 | tissues | comma-separated, alphabetically sorted list of supporting the left splice junction (max=2 * anchor_length) |
| 12 | tiss_counts | comma-separated list of corresponding read-counts |
| 13 | edits | number of mismatches in the anchor extension process |
| 14 | anchor_overlap | number of nucleotides the breakpoint resides within one anchor |
| 15 | breakpoints | number of alternative ways to break the read with flanking GT/AG |
| 16 | signal | flanking dinucleotide splice signal (normally GT/AG) |
| 17 | strandmatch | 'MATCH', 'MISMATCH' or 'NA' for non-stranded analysis |
| 18 | category | list of keywords describing the junction. Useful for quick grep filtering |
参考资料
