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.gztar xzvf v1.2.tar.gz
需要注意的是,该软件是基于python2的语法开发的,依赖pysam
和numpy
两个模块。具体的运行流程可以分成以下几步
官方的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文件就可以了。
采用samtools软件提取没比对上的序列,代码如下
samtools view -hf 4 accepted_hits.bam | samtools view -Sb - > unmapped.bam
代码如下
unmapped2anchors.py unmapped.bam anchor.fq
bowtie2 -p 16 --reorder --mm --score-min=C,-15,0 -q -x human_bowtie2_index -U anchor.fq -S align.sam
代码如下
cat align.sam | find_circ.py -G hg19.fa -p hsa_ > splice_sites.bed
结果如下所示
-p
参数指定的是第四列内容的前缀,建议指定为物种对应的三字母缩写,需要注意的是,在sites.bed
中同时包含了环状RNA和线性RNA,环状RNA的名称用circ
标识,线性RNA的名称用norm
标识。
根据以下规则对结果进行筛选
根据关键词CIRCULAR筛选环状RNA
去除线粒体上的环状RNA
筛选unique junction reads数至少为2的环状RNA
去除断裂点不明确的环状RNA
过滤掉长度大于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·
联系客服