最近在看circRNA的文章，检测工具从2013年到现在已经有了几款，但是每个软件的检测结果都差异很大。

在biostar上，有人在最近几天讨论了2个国人做的circRNA检测工具，其中的问题还是很大的。

现在就我使用过的circRNA detection tools 做一个简要的介绍。

1. find_circ

这款工具是在2013年随着Circular RNAs are a large class of animal RNAs with regulatory potency这篇文章发布的，从而吹起了一阵研究circRNA的风潮～。

现有的所有circRNA工具的第一步都是mapping，将RNA reads mapping到基因组上，但是普通线性RNA是可以很好的mapping回去的，对于circRNA来说，成环的剪接位点不能直接mapping回基因组上。 首先，就是用比对mapping的方法筛选出这样的序列，也就是把mapping不上的序列都体取出来。 然后，将这些mapping不上的reads取两头20bp，变成短序列，重新mapping到基因组上。

这里有个参数20bp的reads，为什么要用20bp的？是否可以用更长的，或者更短的？

我们来做个简单计算。多长的一条序列可以唯一比对到基因组的某个地方？ 首先，基因组的大小为3G，也就是三十亿字节。\(4^{15}=1073741824\)，约为十亿字节。 那么任意取出一条15bp的read，它的碱基排列情况是4^15的其中一种，任意取出两条15bp的reads 它们完全相同的概率为\(\frac{1}{4^{15}}*\frac{1}{4^{15}}\)。 从基因组上的任意取出15bp长的一段序列，有三十亿种可能性，取两次， 取出同一条15bp长的序列的可能性为\(\frac{1}{3*10^{9}}*\frac{1}{3*10^{9}}\)， 这种可能性(概率)小于任意从15bp序列的全部碱基排列情况中抽出2次相同碱基序列的概率，所以说15bp的read可以不唯一比对到基因组的某个地方。 同理计算16bp reads \(4^{16}=4294967296\)，所以任意抽出一条16bp的序列可以唯一比对到基因组上。（可能算得不对，如果有人发现错误，请立即指正，谢谢）

这样看来20bp的序列，可以唯一比对到基因组上，也可以根据测序reads的长度适当放大缩小（我觉得缩小还是没必要的，可以放大，取25bp，30bp之类的）。

接下来，将短序列mapping到基因组上后，他们开发了一个方法，检测这些短序列是否是circRNA的短序列。

需要检测的条件如下：

  1. GU/AG 在剪接位点的两侧出现
  2. 可以检测到清晰的breakpoint
  3. 只支持2个mismatch
  4. breakpoint不能再短序列（anchor）以里2 nucleotides 之外的地方出现，也就是最多2nt （这条可能理解有偏差）
  5. 至少有两条reads支持这个junction
  6. mapping正确的一个短序列的位置要比它mapping到其他位置的分值高35分以上。

通过这些条件，就筛选出了文章认为正确的circRNA。

另外在这篇文章里重点介绍了circRNA的一个功能，富集microRNA，某些circRNA上有很多microRNA的seed。

2. CIRCexplorer

这款工具是在2014年随着Complementary Sequence-Mediated Exon Circularization的发布而问世，这篇文章我非常欣赏，我认为它是去年做RNA领域最值得读的一篇文章。

CIRCexplorer这个工具巧妙的用fusion gene这个思路去检测circRNA。 首先，过滤出Tophat无法mapping上的reads，再把这些reads用Tophat-Fusion mapping到基因组上。那些用Tophat-Fusion mapping到基因组上的非线性候选位置的reads就是潜在的back-splied juction reads。 接下来，这些reads会在有基因注释的帮助下，确定更加精确的donor, acceptor位置。 最后，对circular RNA进行注释。

这篇文章里介绍了在intron区域的反向重复Alu序列是引起circRNA形成的一个原因（即“内含子配对驱动环化”（intron-pair-driven circularization）模型，除此之外还有一种模型是“套索驱动环化”（lariat-driven circularization）模型），这些序列在跨越外显子的内含子区间上可以形成互补序列，所以在转录过程中容易形成茎环结构。环的部分就是外显子，然后在剪接酶的参与下，被剪接成了环形。

另外，circRNA也存在很多的可变环形结构，这也是由于在基因组中广泛存在的内含子上的Alu序列造成的，不同的Alu序列互补，形成含有不同exon的环形结构。

强烈建议研究circRNA的工作者好好学习一下这篇文章。附薛宇老师对这篇文章的评价，薛宇老师的博客是我看科学网唯一的动力。

3. CIRI

这款工具是中科院北京生命科学研究院的赵老师组的工作，今年发表在genome biology上， CIRI: an efficient and unbiased algorithm for de novo circular RNA identification，是一篇方法工具类的文章。

文章中总结了目前从RNA-seq数据中识别circRNA所遇到的问题：

  1. circRNA比其他RNA在细胞中的比例低，一般RNA-seq的实验步骤中，不包括circRNA富集的步骤（例如RNase R treatment，我觉得他们要表到的意思是在RNA-seq中能探测到的circRNA比用RNase R处理过的数据中少），所以在RNA-seq中的circRNA比较低，假阳性也高。
  2. 目前的注释文件都是根据线性RNA进行的，所以也不适合circRNA的识别，另外非模式动物的注释信息更少，就更不用谈在那些物种里找circRNA了。
  3. 由于RNA测序数据的reads长度差别大，也对检测circRNA工作带来了不便
  4. 套索结构和融合基因同circRNA的reads结果类似，不好区分。

文章中总结了从2012年到2014年出现的几种检测circRNA的方法：

  • Salzman Cell-type specific features of circular RNA expression 用了一个依赖注释信息的方法来检测circRNA，通过搜索已知注释的外显子边界来查找，并在最近的工作中更新了方法，加入了false discovery rate 控制比对的质量分数。这种方法存在上述描述的第二个问题。
  • Memczak Circular RNAs are a large class of animal RNAs with regulatory potency （也就是第一部分介绍的那个软件）用了GT-AG信号来找寻splicing 位点，也有其他工作用类似的方法筛选micorRNA-sponge 的候选circRNA。这种方法会找不到“长外显子1-短外显子-长外显子2”形成的环形结构，这种结构中一条测序Read上会有三个部分，第一部分序列属于长外显子1，第二部分序列属于短外显子，第三部分序列属于长外显子2。 Memczak的方法只是把一条序列切成两部分，这种算法会把“长外显子1-短外显子-长外显子2”丢掉，或者识别成“长外显子1-长外显子2”。
  • Jeck Circular RNAs are abundant, conserved, and associated with ALU repeats 采用了比较的方法，比较没有经过RNase处理和经过RNase处理的序列的结果，用来确定潜在的circRNA，消除假阳性。这种方法在富集circRNA阶段会有系统误差。

说了那么多其他人的工作，赵老师组的工作是采用sam格式中的CIGAR值进行分析的，从sam文件中扫描PCC信号（paired chiastic clipping signals）。 CIGAR值在junction read的特征是xS/HyM或者xMyS/H，其中x,y代表碱基数目，M是mapping上的，S是soft clipping，H是hard clipping。

对于单外显子成环，或者“长外显子1-短外显子-长外显子2”形成的环形结构，CIGAR值应该是xS/HyMzS/H以及(x+y)S/HzM或者xM(y+z)S/H，CIRI软件可以很好的将这两种情况分开。

对于paired-end reads，CIRI算法会考虑一对reads，其中一条可以mapping到circRNA上，另一条也需要mapping到circRNA的区间内。

对于splicing 信号(GT,AG) CIRI也会考虑其他弱splicing 信息例如（AT-AC），算法会从GTF/GFF文件中抽取外显子边界位置，并用已知的边界来过滤假阳性。

CIRI方法消耗的内存比较大，我跑了个12G的sam文件，内存用了20G。

从mapping工具中找backsplice junction

除了上述直接找circRNA的工具外，还可以用mapping工具直接找backsplice的junction，然后在人工判断这些junction，例如segemehl。

#Create index:
./segemehl.x -d hg19.fa -x hg18.idx
#Map you reads:
./segemehl.x -d hg19.fa -i chr1.idx -q reads.fastq -S | samtools view -bS - | samtools sort -o - deleteme | samtools view -h - > mapped.sam
#Call different splice junctions:
./testrealign.x -d hg19.fa -q mapped.sam -n

最后一步可以生成两个文件，一个是splicesites.bed，另一个是transrealigned.bed，从splicesites.bed文件中可以得到backsplice junction，在文件中的指示标志是”C”。

从我的测试样本中可以找到4085个通过检验的backsplice junction。

方法比较

从一开始的biostar链接中就可以看到不同方法的差异非常大，我测试了一个实际数据，从第一个软件有126163行结果，筛选后剩下491行，第二个软件只有327行，第三个软件有687行。

我认为不同的mapping策略从一开始就会导致搜索circRNA的空间不一样。

下游软件

画图

https://github.com/dieterich-lab/CircTest

中文介绍

环形RNA分子揭秘——张 杨 1 ,王海滨 1,2 ,陈玲玲 1 *

（未完待续）