bedtools miRanda——妈妈再也不用担心断网了我不能预测circRNA相结合的miR...

问题的提出

解决这个问题的思路其实很简单。大家都知道核酸分子间相互结合的理论基础就是碱基互补配对，并且已经有了很好的基于序列信息预测结合位点的生物信息学软件，如miRanda。该软件最初是为miRNA靶向mRNA设计的，需要miRNA序列和mRNA UTR序列作为输入。在这个基础上，我们只要把mRNA相关序列换成环状RNA的序列就解决问题了。即最终的问题转换成了如何获取环状RNA全长序列。

首先要知道环状RNA在基因组上的位置

这里推荐使用CIRCexplorer2流程完成环状RNA的鉴定工作（只针对RNA-seq分析）。因为CIRCexplorer2给出了相当丰富的环状RNA位置信息，如下图：

mark

由于大多数环状RNA来源于可变剪切的外显子区域，我们在获得其首尾相接的位点后还需要考虑去除中间可能包含的内含子序列。利用以上信息，我们只要摘取外显子序列就可以避免引入内含子区域序列。

下一步就是在CIRCexplorer2的输出基础上制作bed格式文件，主要目的是满足bedtools getfasta的输入要求。假如有一批18个样本的CIRCexplorer2输出结果，每一个样本放置在自己的文件夹下面，如下图：

mark

mark

红框标记是C01号样本的CIRCexplorer2鉴定结果输出文件，具体文件格式和内容请参考CIRCexplorer2网站。

我们要做的就是从18个这样的文件中提取全部环状RNA的位置区间信息，并以bed文件格式输出。实现方式用python,代码如下：

import os
import glob
import sys

CIRCexplorer2_path = sys.argv[1]
file_name = sys.argv[2]
f = open(file_name, 'w+')

circ_list = []
for i in next(os.walk(CIRCexplorer2_path))[1]:
    for file in glob.glob(os.path.join(CIRCexplorer2_path, i, '*known.txt')):
        with open(file, 'r') as filereader:
            for line in filereader:
                l = line.strip('\n').split()
                if len(l) <>10: continue
                if int(l[12]) <>2: continue
                circ_id = '%s_%s_%s'%(l[0], l[1], l[2])
                if circ_id in circ_list:
                    continue
                else:
                    circ_list.append(circ_id)
                    f.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' % (l[0], l[1], l[2], circ_id, l[4], l[5], l[6], l[7], l[8], l[9], l[10], l[11]))


bedtools getfasta获取环状RNA全长序列

执行完上一步中最后的python脚本，得到如下形式的示例文件：

mark

随后执行bedtools getfasta代码，从参考基因组中直接获取环状RNA全长序列。

bedtools getfasta -fi /data/yaoyh/public_data/hg38.fa -bed my_project.txt -split -s -name | fold -w 70 > my_circRNA.fa

miRanda压轴

有了环状RNA全长序列后，直接就能运行miRanda了（前提是你已经下载好了miRNA序列），代码很简单：

miranda mature.hsa.fa my_circRNA.fa -en 20 -strict -out my_predication.txt


最后的表演来自对miRanda输出文件的处理。miRanda的输出结果想必多多少少会得到大家的吐槽，有的没的全都展示。没关系，我们同样用几行简单的python代码就可以提取出想要的结果，并且控制输出最少的结合位点个数，代码如下：

import sys

input_file = sys.argv[1]
output_file = sys.argv[2]
num_site = sys.argv[3]

f = open(output_file, 'w+')
f.write('Seq1\tSeq2\tTot Score\tTot Energy\tMax Score\tMax Energy\tStrand\tLen1\tLen2\tPositions\n')

with open(input_file, 'r') as filereader:
    for line in filereader:
        if line.startswith('>>'):
            line = line.strip('>>')
            l = line.split('\t')
            binding_sites = l[9].strip().split()
            if len(binding_sites) > int(num_site):
                f.write(line)

学会了以上操作，就可以在本地部署自己的预测流程，不再依赖网络，并且可以得到未知环状RNA的靶点信息。