circRNA分析软件CIRCexplorer2之初体验

日常瞎掰

  时间过得有点快啊,不知不觉躺平了这么多天,内心还是有些好不情愿的内疚。不过,话说回来,说到躺平让我想起来了不知道是谁说的一段非常有意思的话,与大家一起分享:“现在内卷的这么严重,当你被别人打倒的时候,心里不要愤怒,而是应该心存感激,因为当你倒下的时候,你会发现原来躺平真的很舒服!”。
  虽然躺了这么久,但咱也没有闲着啊,这不,让我刷到了一些别人挣钱的外门邪道,下面让我们一起来鉴别一下这到底有多毒:
(1)今天下大雪,刚才出门见一大爷摔倒了,我过去问道:大爷,我一月工资不到2000块钱,能扶您起来么?大爷:小伙子,你走吧,我再等一会儿。我:好勒!天气虽冷,大爷的话却是暖暖的,满满的都是正能量……
(2)走在路上看到大爷躺在地上,我赶紧上去要扶,大爷看了看我说:后生你别动,看你也是个打工的,你走吧我再等等。我被感动了赶紧说:大爷,有个开路虎的停好车就快过来了。大爷也激动了:你这后生还算实在,你干脆别走了给我做个证人,完事儿也买个车上班开。
(3)走在路上看到个大爷躺在地上,我赶紧上陪大爷躺下,大爷看了看我说:后生你别跟我抢,我儿子要买房娶媳妇。我说:我也要买房娶媳妇啊,但我不啃老。大爷说:后生有志气,这地给你了,我换个地儿!

走起

  看到人家大爷和小伙都那么玩命地挣钱,那咱也不能闲着是不!下面咱们就来说点正事。这段时间收到一个circRNA的数据分析,虽然以前也没有分析过这个数据类型,但咱可以学习啊。于是乎,先到网上一搜,还别说关于circRNA的帖子还真不少,不管对错先看了再说,至少让咱对circRNA有个初步的认识了,关于circRNA内容大致如下:
  环状 RNA (circRNA)是一类由 mRNA 前体经反向可变剪切(back-splicing)而来、共价闭合且保守的单链转录本,通过 miRNA 海绵功能、干扰可变剪切、结合蛋白等方式调控来源基因及线性 mRNA 的表达。
  主流的circRNA形成模型主要有以下5种:
  1、反向可变剪切环化;
  2、内含子驱动反向互补序列形成环化;
  3、单个基因内不同内含子序列配对环化;
  4、受到RNA结合蛋白调控的外显子环化;
  5、依赖于RBPs环化模式。

  circRNA的生物功能:
  1、调控转录及剪切过程;
  2、编码蛋白;
  3、circRNA驱动的假基因;
  4、miRNA海绵及蛋白互作;
  5、引起免疫响应;
  6、生物标记物。
  circRNA由特殊可变剪切产生,主要来源于外显子,少部分来源于内含子,表达水平具有种属、组织、时间特异性,具有一定序列保守性。虽然circRNA的表达量很低,但其在生物过程的功能却不容忽视。

鉴定软件

  目前,已经有不少公开发表的circRNA鉴定程序,这些软件根据原理以及需要是否注释信息可分成两大类,如不需要注释信息的软件有CIRIquant、CIRI2、find_circ2、CIRCexplorer2、MapSplice2等,需要注释信息的有KNIFE、NCLscan等,当然有注释信息的软件会有更高一点的准确性。不过,都是基于算法来鉴定circRNA肯定会有假阳性,想要提高准确性最好结合两种或多种不同软件预测的结果来有效提高 circRNA 鉴定。关于这些软件也有专业的测评,感兴趣的可以自己去读一读<<Computational Strategies for Exploring Circular RNAs>>,有助于选择合适的软件。

各软件预测率

CIRCexplorer2

  下面我们来介绍一下CIRCexplorer2的使用方法,为什么介绍这个软件呢?原因很简单,因为咱关注的circRNA存在于人家的数据库里,人家数据库里的circRNA鉴定使用的就是自己开发的CIRCexplorer2,所以咱也想尝试一下用这款软件能不能在咱们的数据中也鉴定出想要的circRNA。

可以看出这款软件功能还是很厉害的,共分为5个功能模块:
1、Align用于将序列比对到参考基因组上,比对软件可以基于TopHat2、STAR 、MapSplice、BWA、segemehl。PS: 如果需要做Assemble、Denovo最好使用TopHat2来比对;
2、Parse用于从比对结果中识别junction reads;
3、Annotate用于预测环状RNA;
4、Assemble用于组装环状RNA的转录本序列;
5、Denovo根据序列组装结果,识别新的环状RNA和分析环状RNA上的可变剪切事件。

下面咱们来看看使用该软件来鉴定circRNA的流程:

# Align
tophat2 -p 5 --fusion-search --keep-fasta-order --bowtie1 --no-coverage-search  -o tophat2 bwtidx read1 read2

# parse
CIRCexplorer2 parse --pe -t TopHat-Fusion -b sample_bsj.bed tophat2/accepted_hits.bam >CIRCexplorer2_parse.log

# annotate
CIRCexplorer2 annotate -r gene_anno.txt -g ref_fasta -b sample_bsj.bed -o sample_circ.txt >CIRCexplorer2_annotate.log

# assemble
CIRCexplorer2 assemble -r gene_anno.txt -m tophat -o assemble >CIRCexplorer2_assemble.log

# denovo
CIRCexplorer2 denovo -r gene_anno.txt -g ref_fasta  -b sample_bsj.bed --abs abs --as as -m tophat -n pAplus_tophat -o denovo >CIRCexplorer2_denovo.log

gene_anno.txt: 基因注释文件,可用软件自带脚本获取,如fetch_ucsc.py hg38 ref out
其实,人家也贴心地把上面的模块做了一个封装,比如可以一行代码从fastq做注释或者鉴定cicrRNA:

# 注释流程
fast_circ.py annotate -r gene_anno.txt -g ref_fasta -G gtf -p 5 -o outdir  -f fastq
# 鉴定流程
fast_circ.py denovo -r gene_anno.txt -g ref_fasta -G gtf  -p 5 -o outdir  -f fastq

踩坑实录

  在使用该软件的时候遇到了一点问题,记录一下,以便遇到同样问题的人。咱当时走的是上面分步骤的流程,使用tophat2做比对加了--bowtie1参数就会调用bowtie1,服务器bowtie1默认版本为version 1.2.1.1,然后遇到了下面的错误:

[sam_read1] missing header? Abort!
        [FAILED]
Error running bowtie:
Error while flushing and closing output
terminate called after throwing an instance of 'int'

  网上一搜,遇到同类问题的人还真不少,一开始以为是存储空间不够导致的,后来提升了一下存储空间还是同样出错。那就只能朝其他的原因努力了,最终发现原因主要是tophat2调用bowtie1做比对时生成的temp.samheader.sam文件里没有@SQ头部信息。有人就给出了解决方法,在tophat2 源码里找到get_index_sam_header这个函数(tophat2是由python写的软件),然后添加如下语句:

if params.bowtie2:
            bowtie_header_cmd += ["-x"]

然后我就按照这个方法去操作了,打开一看,傻眼了,咱软件的源码里有这句代码了:

def get_index_sam_header(params, idx_prefix, name = ""):
    noSkip = currentStage >= resumeStage
    try:
        temp_sam_header_filename = tmp_dir + "temp.samheader.sam"
        temp_sam_header_file = None
        if noSkip:
          temp_sam_header_file = open(temp_sam_header_filename, "w")

        bowtie_header_cmd = [bowtie_path]

        read_params = params.read_params
        if not params.bowtie2:
            bowtie_header_cmd += ["--sam"]

        if read_params.color:
            bowtie_header_cmd.append('-C')

        if params.bowtie2:
            bowtie_header_cmd += ["-x"]

  心里顿时不知所措,原因找到了,但不知道怎么处理,此时好像是趴在玻璃窗上的苍蝇,前途一片光明,就是找不到出路!最后,实在没招了只能换一个bowtie1版本(bowtie version 1.1.2)试一试了,真的是不试不知道,一试吓一跳。试完以后,当时的心情有点稍微的复杂,有点像“做错了事的孩子找到了弥补的方法,内疚与喜悦交加,愤怒与激动并存”。如果一定要找一个词来形容当时的心情,那只能是言简意赅却特别富有感染力和共情性的词语——卧槽

结束语

  好吧,咱承认今天废话有点多了,内容有点水可以实锤了。但习惯还是要保持一致,按照惯例,为了大家的方便,最后附上官网的链接和参考文献。

CIRCexplorer2: https://circexplorer2.readthedocs.io/en/latest
<<Computational Strategies for Exploring Circular RNAs>>

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,470评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,393评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,577评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,176评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,189评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,155评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,041评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,903评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,319评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,539评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,703评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,417评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,013评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,664评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,818评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,711评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,601评论 2 353

推荐阅读更多精彩内容