-
Notifications
You must be signed in to change notification settings - Fork 1
/
6.CIP_CALP.pl
37 lines (31 loc) · 1.32 KB
/
6.CIP_CALP.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#####filter blastp output for ancestral karyotype analysis by CIP & CALP parameters###
#####Both CIP CALP should >=0.5###
##e.g.##
##makeblastdb -in species2.pep -dbtype prot -out species2.pep.db##
##blastp -query species1.pep -db species2.pep.db -num_threads 12 -evalue 1E-5 -outfmt 0 -out species1-species2.blastout###
##perl cip_calp.pl species1-species2.blastout>species1-species2.blastout.filtered###
use Bio::SearchIO;
my $blastout = shift;
my $filter = Bio::SearchIO->new(-file => $blastout,
-format => 'blast');
print "Query\tlength_Q\tnumber_of_hits\tHit\tlength_H\te-value\tCIP\tCALP\n";
while( my $result = $filter->next_result ) {
while(my $hit=$result->next_hit){
my $al=0;
my @id=();
while(my $hsp=$hit->next_hsp){
$al=$al+($hsp->hsp_length);
my $ai=($hsp->num_identical)/($hsp->hsp_length);
push @id,$ai;
}
my $cip=0;
foreach (@id){
$cip=$cip+($_/$al*100);
}
my $calp=$al/($result->query_length);
if($cip ge 0.5 and $calp ge 0.5)
{
print $result->query_name,"\t",$result->query_length,"\t",$result->num_hits,"\t",$hit->name,"\t",$hit->hit_length,"\t",$hit->significance,"\t",$cip,"\t",$calp,"\n";
}
}
}