Question

生物信息学perl程序，求助。

rps-0,B0393.1
Exon12CAATCTTAACTTCCAAATGCAGCAATACGTCTACAAGAGACGTTTTGACGGACCAAATATCATCAACGTCAAGAAGACCTGGGAGAAGCTTCTTCTCGCT
Exon34TCAAGGAGCCAAGACTTCTCGTCATCTCTGACCCACGTATCGATCATCAGGCTGTCACTGAGGCTTCCTACGTCGGAGTTCCAGTCATCTCCTTCGTCAA
Intron1gttagtttaaaaaaaattcaattcatttttcttgtgaattatttaataaaattataacaattatcattacag
Intron2gtaagacttttctaaatagttgttttttctaattatttaaaatttag
Intron3gtatatatttaatttataaaaaatgtaaatatgtatttgaattgtag
rps-2
Exon12GGACGTGGAGGACGTGGTGGCCGCGGAGGACGTGGTGGACGCGGAGGAAGAGCCGGACGCGGAGGAGAGAAGGAGACCGAATGGACTCCAGTCACCAAGC
Intron1gtagttattttattcttaaaccttcagattttgaatatttaaatttcag
我想让基因rps-0,B0393.1中的Exon12和Intron1比对，Exon34和Intron2比对；同理基因rps-2也是Exon12和Intron1比对。方法调用BLAST ，用Bioperl 来解析。现在我只能处理2条序列间的比对，一批量就出问题。这是我的一部分代码怎么修改可批量处理。
use strict;
use Bio::SearchIO;
my $seq1="exon12.txt"; #两个输入文件
my $seq2="intron1.txt";
my $tempoutput="bl2seq.txt"; #输出文件
my $bl2seq="C:/blast/bin/bl2seq";
my $system_check=system("$bl2seq -i $seq1 -j $seq2 -p blastn -F F -o $tempoutput -g T");
my ($result,$hit,$hsp);
my $in = new Bio::SearchIO(-format => 'blast',
-file => "$tempoutput",
-report_type => 'tblastn');
while( $result = $in->next_result ) {
while( $hit = $result->next_hit ) {
while( $hsp = $hit->next_hsp ) {
if( $hsp->length('total') > 1500 ) {
if ( $hsp->percent_identity >= 80 ) {
print "Hit= ", $hit->name,
",Length=", $hsp->length('total'),
",Percent_id=", $hsp->percent_identity, "\n";
}
}
}
}
}

Answer 1

　　老实说我也不是很明白，我之前的一段程序也是：只做1个的时候，就可以。一旦大批量的文件IO了，就挂掉了……而且我里面也是要用system来打开一个外部线程的。所以看到你这个我觉得真心虎躯一震，菊花一紧啊……

 

1. 不过既然你的“my $system_check=system("$bl2seq -i $seq1 -j $seq2 -p blastn -F F -o $tempoutput -g T");”这句都已经取得返回值了，为什么不检查一下$system_check的值是不是0这个正常返回值呢？

2. 话说NCBI是不是不提倡用bl2seq了？还是我的记忆牙碜了？

3. 看着这段程序的样子，貌似你不是让基因rps-0,B0393.1中的Exon12和Intron1比对，Exon34和Intron2比对；同理基因rps-2也是Exon12和Intron1比对吧？而是基因rps-0,B0393.1中的Exon12和它自己的Intron们加上基因rps-2的intron1比对吧，然后其他的外显子也是一样的……你这里是不是得调整成只让它的文件1里面只有一个exon的而且文件2里面只有对应的intron的时候才行啊？

    

    

   我觉得这个要是调整起来就大发了……你加油吧……

追问

这位老兄雄起呀，分析真是透彻。

追答

没有雄起……你做好了么？要是还有需要讨论的，咱们可以再聊？744684934

Answer 2

代码没问题，把这段代码写成函数，将$seq1, $seq2, $output设置成参数，用一个循环往里面传参数就行了

追问

大侠，可否给出代码看看。