快速提取序列的Perl脚本

是这么回事的,整理一下一些有用的提问。感谢大家的提问。

网友sun:

我从miRbase数据库下载到了,matrue.fa,即里面有测出来的物种的成熟miRNA,里面有human 我想用perl把human的提出来另外存到一个文件.这是其中人的一条的
>hsa-miR-96 MIMAT0000095 Homo sapiens miR-96
UUUGGCACUAGCACAUUUUUGCU
怎么把所有的都提出来.
ps:正则表达式不是这个$s=~/^>hsa/期待你的回信.

有接触过miRbase数据库的,都知道miRNA的命名Human是用hsa开头,Mouse是用mmu开头,Rat是用rno开头。

解释一下

忘了大伙对Fasta格式可能有不明白的,我来稍微解释一下上面的话。比如有一个文件:matrue.fa。里面存放的内容的格式是这样子的:

>hsa-miR-96
UUUGGCACUAGCACAUUUUUGCU
>mmu-miR-96
AAAGGCACUAGCACUCGAUCGA
>rno-miR-85
UAGUCUGCAUGUUGCACAGUGUA

上面只是举一个为例子,数据都是批量的。想要取出以hsa命名开头的数据,取出来的结果是:

>hsa-miR-96
UUUGGCACUAGCACAUUUUUGCU
……

 

这样子就有了下面的Perl脚本:

提取Fasta序列的Perl脚本

#!/usr/bin/perl
$/ = '>';
print '>';

while(<>){
    print if(/hsa-/);
}

这样就实现了从Fasta格式的文件中提取序列的目的,如果你再稍微改改,类似GenBank格式的也是可以提取的。还有好多,只要有唯一分隔符的都可以,不一一列出。

提醒一下:

$/ = '>';

强大之处是用了这句,自行慢慢体会。

27 回复
  1. gahoo says:

    这句是什么意思?
    $/ = ‘>’;

    柳城 回复:

    你可以自己测试一下啊。你找一个fasta格式的序列。

    #!/usr/bin/perl
    while(<>){
    print;
    }

    第二种是:

    #!/usr/bin/perl
    $/='>';
    while(<>){
    print;
    }

    试了你就知道了呗~~ [ok]

  2. niche says:

    非常感谢您的回复,PAML是算核苷酸替换率的,http://www.bioperl.org/wiki/HOWTO:PAML
    非常喜欢你的博客,很多内容都很有用。

  3. motianfeng says:

    想请问一下各位大虾,我现在有一个达6G的txt数据,也是用’>’分开的,采用上面的方法提取每条染色体下面的数据,结果生成了几千个只有’>’的文件,不知是咋回事?另外,如果想要达到我的目的,有没有什么不太耗内存的方法?多谢啦先!~~

    柳城 回复:

    具体问题具体分析嘛。可到yunbio.com提问哦

  4. tempo8 says:

    类似于http://liucheng.name/1205/
    #Usage:perl fasta2.pl in_fasta out_file
    in_fasta,out_file写在哪行?如何读文件以及如何输出结果?

  5. tempo8 says:

    楼主还有一个问题:
    如何批量定点提取序列,即CDS-Genome作BLAST比对之后,要提取ATG上游2000bp的序列!
    如:
    正链/正链(plus/plus)
    >subject1
    Query 1 ATCG
    Subject 5186 tagc
    这里有个特征是Query 1下一行对应的Subject后的数字(5186)减2000,即subject1中提取3186到5186

    此外,如何提取CDS下游2000bp的序列?
    若(plus/plus),在下一个”>”到上一个”>”之间在所有Subject这一行取最大数,直接加2000
    若(plus/minus)在Query这一行取最大数,再看Subject这一行对应的数值,然后减2000

    望不吝赐教!谢谢

    柳城 回复:

    请到yunbio.com提问。。

  6. Mana Kimzey says:

    Hi there, I found your website via Google while looking for a related topic, your site came up, it looks great. I have bookmarked it in my google bookmarks.

  7. 尹志远 says:

    柳大神,你好!拜读你的网站已久,可恕在下鲁钝,还是看不太懂,你的perl脚本貌似是直接把序列输出到屏幕,怎样用perl脚本按指定ID搜索,然后把这条序列输出为一个单独的fasta文件?谢谢啊!由于我刚学perl,还望能说得详细一点。

    柳城 回复:

    可到yunbio.com提问

评论已关闭。