||
适用情形:需要从包含多个FASTA序列的文件中,根据序列名提取其中部分FASTA序列。一般情况,这个序列文件中包含的序列数量比较多,例如水稻基因组全部预测基因的序列,每个序列为一个fasta文件。主要为两步:第一步,建立索引文件。第二步,提取需要的序列。在这两步间还有一个准备步骤。
第1步的脚本文件:
####从一个或多个FASTA文件中建立索引
use Bio::Index::Fasta;
use strict;
my $Index_File_Name ='IRGSPv5_gene.fasta';###索引文件名,该文件在程序运行的当前目录中
my $inx = Bio::Index::Fasta->new(-filename => $Index_File_Name,
-write_flag => 1);
$inx->make_index('/home/bioysy/rice_genome_seq/IRGSPv5_gene.fasta');###建立索引文件的原始文件。
第2步的脚本文件:
###利用索引文件,按序列名提取序列
use Bio::Index::Fasta;
use Bio::SeqIO;
use Bio::Seq;
open(GENID,'/home/bioysy/bioperl_study/Bio_Index_Fasta /get_genesid');####包含基因名的文件,我建的这个文件中基因名为:Os05t0149400-01,
############如果用Os05g0149400去搜索的话是没有结果的,用Os05t0149400也没结果。
###GENID为文件句柄
@genename=<GENID>; ##将文件内容储存在数组@genename中。
my $Index_File_Name ='IRGSPv5_gene.fasta'; ##已经建立的索引文件IRGSPv5_gene.fasta。
my $inx = Bio::Index::Fasta->new(-filename=>$Index_File_Name);
my $out = Bio::SeqIO->new(-format => 'Fasta',
-fh =>\*STDOUT); ###结果输出到屏幕上,存到文件中,我用的是LINUX系统的重定向,如>gene.fasta
foreach my $id(@genename) {
chomp $id;
my $seq = $inx->fetch($id); # Returns Bio::Seq object
$out->write_seq($seq);
}
准备步骤:得准备下第2步里面的名为get_genesid的文件。
当然,要运行上面的脚本得用bioperl.
以前鼓捣过这个东西,失败了。这次算是成功了吧。所以,写篇短文,以念之。当有许多序列需要处理的时候,就会感到现成的软件不够用了,所以得写点小程序才能应付自如。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 04:29
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社