Mercurial > repos > dereeper > pgap

comparison PGAP-1.2.1/Blast_Filter.pl @ 0:83e62a1aeeeb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Thu, 24 Jun 2021 13:51:52 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:83e62a1aeeeb
1 #!/usr/bin/perl
2 use strict;
3 use warnings;
4
5
6 my $BLASTResult=$ARGV[0];
7 my $SEQFile=$ARGV[1];
8 my $coverage=$ARGV[2];
9 my $identity=$ARGV[3]*100;#percents
10 my $score=$ARGV[4];
11 my %hash;
12 my @row;
13 my $line;
14
15 #my $coverage=0.5;
16 #my $identity=50;#percents
17 #my $score=50;
18
19 &ReadSeqLength("$SEQFile",\%hash);
20
21 open(F,$BLASTResult);
22 while ($line=<F>)
23 {
24 if ($line!~/\#/)
25 {
26 @row=split(/\t/,$line);
27 if ($row[2]>=$identity and $row[11]>=$score)
28 {
29 if ((($row[7]-$row[6]+1)/$hash{$row[0]})>=$coverage)
30 {
31 if ((($row[9]-$row[8]+1)/$hash{$row[1]})>=$coverage)
32 {
33 print "$row[0]\t$row[1]\t$row[11]";
34 }
35 }
36 }
37 }
38 }
39 close(F);
40
41 sub ReadSeqLength()
42 {
43 use Bio::SeqIO;
44 (my $file,my $hash)=@_;
45 my $seq;
46 my $in=Bio::SeqIO->new(-file=>"$file",-format=>"fasta");
47 while ($seq=$in->next_seq())
48 {
49 $$hash{$seq->id}=length($seq->seq());
50 }
51 }