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

use Bio::SeqIO;
use Bio::SearchIO;
my $in1=Bio::SeqIO->new(-file=>"exon12.txt",-format=>"fasta");
my $seq_num_1=1;
my $out="out";
open OUT,'>>',"bl2seq.txt";
my $bl2seq="C:/blast/bin/bl2seq";
while(my $seq1=$in1->next_seq())
{
open TMP1,'>',"tmp1";
print TMP1 ">",$seq1->display_id,"\n",$seq1->seq,"\n";
my $seq_num_2=1;
my $in2=Bio::SeqIO->new(-file=>"intron1.txt",-format=>"fasta");
while(my $seq2=$in2->next_seq())
{
if($seq_num_2==$seq_num_1)
{
open TMP2,'>',"tmp2";
print TMP2 ">",$seq2->display_id,"\n",$seq2->seq,"\n";
system("$bl2seq -i tmp1 -j tmp2 -p blastn -F F -o $out -g T");
system("cat $out >>blast");
}
$seq_num_2++;
}
$seq_num_1++;
close TMP1;
close TMP2;
}
my ($result,$hit,$hsp);
my $in = new Bio::SearchIO(-format => 'blast', -file => "blast", -report_type => 'blastn');
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 OUT "Hit= ",$hit->name,",Length=",$hsp->length('total'),",Percent_id=", $hsp->percent_identity,"\n";
}
}
}
}
}

追问

呵呵，谢谢你。
1 不过不能运行出批量运行结果，也是只能2条比对。
2 还有代码里的cut（system("cat $out >>blast")我这里怎么不能识别;
3 给我你的邮箱，我给你传个附件你帮我运行下，我的邮箱zhangqiang829@163.com

追答

不好意思 我在Linux下用的 都是shell的命令～～

追问

那在windos系统有相应的命令吗？怎么替换，谢谢。

Answer 2

真心不懂生物....不知道什么操作才是比较两序列

Answer 3

把他拆分成两个文件，然后blast不就行了？

追问

在线blast我能弄了，可是基因比较多，批量不会呀。

追答

不是有个多序列比对嘛我记得，你找找好像叫什么clusterx