使用find_circ识别环状RNA

欢迎关注”生信修炼手册”!

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语言开发，安装过程如下

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软件，代码如下

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软件提取没比对上的序列，代码如下

samtools view -hf 4 accepted_hits.bam | samtools view -Sb - > unmapped.bam

3. 从序列两端提取锚点序列

代码如下

unmapped2anchors.py unmapped.bam anchor.fq

4. 将锚点序列比对参考基因组

bowtie2 -p 16 \
--reorder  \
--mm \
--score-min=C,-15,0 \
-q -x human_bowtie2_index \
-U anchor.fq  \
-S align.sam

5. 预测circRNA

代码如下

cat align.sam | find_circ.py  -G hg19.fa -p hsa_ > splice_sites.bed

结果如下所示

-p参数指定的是第四列内容的前缀，建议指定为物种对应的三字母缩写，需要注意的是，在sites.bed中同时包含了环状RNA和线性RNA,环状RNA的名称用circ标识，线性RNA的名称用norm标识。

6. 结果筛选

根据以下规则对结果进行筛选

1. 根据关键词CIRCULAR筛选环状RNA
2. 去除线粒体上的环状RNA
3. 筛选unique junction reads数至少为2的环状RNA
4. 去除断裂点不明确的环状RNA
5. 过滤掉长度大于100kb的circRNA,这里的100kb为基因组长度，直接用环状RNA的头尾相减即可

代码如下

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·

—如果喜欢，快分享给你的朋友们吧—

扫描关注微信号，更多精彩内容等着你！