欢迎关注”生信修炼手册”!
circBase是环状RNA数据库的开山鼻祖，其中的环状RNA都是采用find_circ这个软件预测得到的。该软件的原理如下
和参考基因组比对完之后，首先剔除和基因组完全比对的reads,只保留没比对上的reads。这部分reads直接比是比对不上基因组的，因为其来自不同的外显子区域，直接比对的话不允许这么大片段的缺失。
那么如何区分剪切的spliced read 和 来自环状RNA的junction read呢，从上面的示意图我们可以直接看出，spliced read 的两部分比对在基因组上的前后位置和转录本中的位置保持一致，而来自circRNA的junction read 其比对的位置是相反的。
具体操作的时候，首先从junction read的5’端和3’端取一部分序列，分别叫做5’ anchor 和 3’ anchor, 如果两个序列比对的位置是相反的，这条reads 就是一个可能的junction read, 然后将anchor read一直延伸，直到连接处为止，如果到连接处为止序列都能够完全匹配，再看连接点处的剪切模式是否符合AG-GT的剪切模式，如果以上条件都满足，就认定这是一个circRNA。
该软件采用python语言开发，安装过程如下
代码语言：javascript复制wget https://github.com/marvin-jens/find_circ/archive/v1.2.tar.gz
tar xzvf v1.2.tar.gz需要注意的是，该软件是基于python2的语法开发的，依赖pysam和numpy两个模块。具体的运行流程可以分成以下几步
1. 比对参考基因组
官方的pipeline使用的是bowtie2软件，代码如下
代码语言：javascript复制bowtie2 -p16
--very-sensitive
--score-min=C,-15,0
--mm
-x hg19 -q
-1 R1.fastq.gz -2 R2.fastq.gz
2> bowtie2.log
| samtools view -hbuS -
| samtools sort - accepted_hits最终生成了一个排序之后的bam文件，其实这一步选择其他的比对软件，比如hisat也是可以的，只需要产生bam文件就可以了。
2. 提取没比上参考基因组的序列
采用samtools软件提取没比对上的序列，代码如下
代码语言：javascript复制samtools view -hf 4 accepted_hits.bam | samtools view -Sb - > unmapped.bam3. 从序列两端提取锚点序列
代码如下
代码语言：javascript复制unmapped2anchors.py unmapped.bam anchor.fq4. 将锚点序列比对参考基因组
代码语言：javascript复制bowtie2 -p 16
--reorder
--mm
--score-min=C,-15,0
-q -x human_bowtie2_index
-U anchor.fq
-S align.sam5. 预测circRNA
代码如下
代码语言：javascript复制cat align.sam | find_circ.py -G hg19.fa -p hsa_ > splice_sites.bed结果如下所示
-p参数指定的是第四列内容的前缀，建议指定为物种对应的三字母缩写，需要注意的是，在sites.bed中同时包含了环状RNA和线性RNA,环状RNA的名称用circ标识，线性RNA的名称用norm标识。
6. 结果筛选
根据以下规则对结果进行筛选
- 根据关键词CIRCULAR筛选环状RNA
- 去除线粒体上的环状RNA
- 筛选unique junction reads数至少为2的环状RNA
- 去除断裂点不明确的环状RNA
- 过滤掉长度大于100kb的circRNA,这里的100kb为基因组长度，直接用环状RNA的头尾相减即可
代码如下
代码语言：javascript复制grep CIRCULAR splice_sites.bed |
grep -v chrM |
awk '$5>=2' |
grep UNAMBIGUOUS_BP |
grep ANCHOR_UNIQUE |
./maxlength.py 100000
> circ_candidates.bed该软件只需要基因组的fasta序列就可以识别环状RNA，适用于绝大多数物种的环状RNA分析。
·end·
—如果喜欢，快分享给你的朋友们吧—
扫描关注微信号，更多精彩内容等着你！
