annotate PGAP-1.2.1/inparanoid.pl @ 0:83e62a1aeeeb draft

Uploaded
author dereeper
date Thu, 24 Jun 2021 13:51:52 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1 #! /usr/bin/perl
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
2 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
3 # InParanoid version 4.1
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
4 # Copyright (C) Erik Sonnhammer, Kristoffer Forslund, Isabella Pekkari,
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
5 # Ann-Charlotte Berglund, Maido Remm, 2007
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
6 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
7 # This program is provided under the terms of a personal license to the recipient and may only
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
8 # be used for the recipient's own research at an academic insititution.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
9 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
10 # Distribution of the results of this program must be discussed with the authors.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
11 # For using this program in a company or for commercial purposes, a commercial license is required.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
12 # Contact Erik.Sonnhammer@sbc.su.se in both cases
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
13 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
14 # Make sure that Perl XML libraries are installed!
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
15 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
16 # NOTE: This script requires blastall (NCBI BLAST) version 2.2.16 or higher, that supports
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
17 # compositional score matrix adjustment (-C2 flag).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
18
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
19 my $usage =" Usage: inparanoid.pl <FASTAFILE with sequences of species A> <FASTAFILE with sequences of species B> [FASTAFILE with sequences of species C]
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
20
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
21 ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
22
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
23 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
24 # The program calculates orthologs between 2 datasets of proteins
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
25 # called A and B. Both datasets should be in multi-fasta file
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
26 # - Additionally, it tries to assign paralogous sequences (in-paralogs) to each
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
27 # thus forming paralogous clusters.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
28 # - Confidence of in-paralogs is calculated in relative scale between 0-100%.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
29 # This confidence value is dependent on how far is given sequence from the
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
30 # seed ortholog of given group
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
31 # - Confidence of groups can be calculated with bootstrapping. This is related
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
32 # to score difference between best hit and second best hit.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
33 # - Optionally it can use a species C as outgroup.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
34
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
35 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
36 # You may need to run the following command manually to increase your
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
37 # default datasize limit: 'limit datasize 500000 kb'
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
38
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
39 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
40 # Set following variables: #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
41 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
42
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
43 # What do you want the program to do? #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
44 $run_blast = 1; # Set to 1 if you don't have the 4 BLAST output files #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
45 # Requires 'blastall', 'formatdb' (NCBI BLAST2) #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
46 # and parser 'blast_parser.pl' #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
47 $blast_two_passes = 1; # Set to 1 to run 2-pass strategy #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
48 # (strongly recommended, but slower) #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
49 $run_inparanoid = 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
50 $use_bootstrap = 1;# Use bootstrapping to estimate the confidence of orthologs#
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
51 # Needs additional programs 'seqstat.jar' and 'blast2faa.pl'
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
52 $use_outgroup = 0; # Use proteins from the third genome as an outgroup #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
53 # Reject best-best hit if outgroup sequence is MORE #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
54 # similar to one of the sequences #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
55 # (by more than $outgroup_cutoff bits) #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
56
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
57 # Define location of files and programs:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
58 #$blastall = "blastall -VT"; #Remove -VT for blast version 2.2.12 or earlier
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
59 $blastall = "$ARGV[0] -a $ARGV[1]"; #Add -aN to use N processors
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
60 $formatdb = "$ARGV[2]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
61 $seqstat = "seqstat.jar";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
62 $blastParser = "blast_parser.pl";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
63
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
64 #$matrix = "BLOSUM62"; # Reasonable default for comparison of eukaryotes.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
65 $matrix = "BLOSUM45"; #(for prokaryotes),
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
66 #$matrix = "BLOSUM80"; #(orthologs within metazoa),
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
67 #$matrix = "PAM70";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
68 #$matrix = "PAM30";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
69
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
70 # Output options: #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
71 $output = 0; # table_stats-format output #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
72 $table = 0; # Print tab-delimited table of orthologs to file "table.txt" #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
73 # Each orthologous group with all inparalogs is on one line #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
74 $mysql_table = 1; # Print out sql tables for the web server #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
75 # Each inparalog is on separate line #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
76 $html = 0; # HTML-format output #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
77
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
78 # Algorithm parameters:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
79 # Default values should work without problems.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
80 # MAKE SURE, however, that the score cutoff here matches what you used for BLAST!
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
81 $score_cutoff = $ARGV[3]; # In bits. Any match below this is ignored #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
82 $outgroup_cutoff = 50; # In bits. Outgroup sequence hit must be this many bits#
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
83 # stronger to reject best-best hit between A and B #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
84 $conf_cutoff = 0.05; # Include in-paralogs with this confidence or better #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
85 $group_overlap_cutoff = 0.5; # Merge groups if ortholog in one group has more #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
86 # than this confidence in other group #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
87 $grey_zone = 0; # This many bits signifies the difference between 2 scores #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
88 $show_times = 0; # Show times spent for execution of each part of the program #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
89 # (This does not work properly) #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
90 $debug = 0; # Print debugging messages or not. Levels 0,1,2 and 4 exist #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
91
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
92 my $seq_overlap_cutoff = $ARGV[4]; # Match area should cover at least this much of longer sequence. Match area is defined as area from start of
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
93 # first segment to end of last segment, i.e segments 1-10 and 90-100 gives a match length of 100.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
94 my $segment_coverage_cutoff = $ARGV[5]; # Actually matching segments must cover this much of longer sequence.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
95 # For example, segments 1-10 and 90-100 gives a total length of 20.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
96
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
97 splice(@ARGV,0,6);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
98
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
99 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
100 # No changes should be required below this line #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
101 ###############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
102 $ENV{CLASSPATH} = "./$seqstat" if ($use_bootstrap);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
103
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
104 if (!@ARGV){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
105 print STDERR $usage;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
106 exit 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
107 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
108
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
109 if ((@ARGV < 2) and ($run_inparanoid)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
110 print STDERR "\n When \$run_inparanoid=1, at least two distinct FASTA files have to be specified.\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
111
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
112 print STDERR $usage;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
113 exit 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
114 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
115
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
116 if ((!$run_blast) and (!$run_inparanoid)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
117 print STDERR "run_blast or run_inparanoid has to be set!\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
118 exit 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
119 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
120
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
121 # Input files:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
122 $fasta_seq_fileA = "$ARGV[0]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
123 $fasta_seq_fileB = "$ARGV[1]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
124 $fasta_seq_fileC = "$ARGV[2]" if ($use_outgroup); # This is outgroup file
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
125
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
126 my $blast_outputAB = $fasta_seq_fileA . "-" . $fasta_seq_fileB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
127 my $blast_outputBA = $fasta_seq_fileB . "-" . $fasta_seq_fileA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
128 my $blast_outputAA = $fasta_seq_fileA . "-" . $fasta_seq_fileA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
129 my $blast_outputBB = $fasta_seq_fileB . "-" . $fasta_seq_fileB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
130
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
131 if ($use_outgroup){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
132 $blast_outputAC = $fasta_seq_fileA . "-" . $fasta_seq_fileC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
133 $blast_outputBC = $fasta_seq_fileB . "-" . $fasta_seq_fileC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
134 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
135 my %idA; # Name -> ID combinations for species 1
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
136 my %idB; # Name -> ID combinations for species 2
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
137 my @nameA; # ID -> Name combinations for species 1
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
138 my @nameB; # ID -> Name combinations for species 2
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
139 my @nameC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
140 my %scoreAB; # Hashes with pairwise BLAST scores (in bits)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
141 my %scoreBA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
142 my %scoreAA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
143 my %scoreBB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
144 my @hitnAB; # 1-D arrays that keep the number of pairwise hits
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
145 my @hitnBA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
146 my @hitnAA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
147 my @hitnBB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
148 my @hitAB; # 2-D arrays that keep the actual matching IDs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
149 my @hitBA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
150 my @hitAA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
151 my @hitBB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
152 my @besthitAB; # IDs of best hits in other species (may contain more than one ID)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
153 my @besthitBA; # IDs of best hits in other species (may contain more than one ID)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
154 my @bestscoreAB; # best match A -> B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
155 my @bestscoreBA; # best match B -> A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
156 my @ortoA; # IDs of ortholog candidates from species A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
157 my @ortoB; # IDs of ortholog candidates from species B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
158 my @ortoS; # Scores between ortoA and ortoB pairs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
159 my @paralogsA; # List of paralog IDs in given cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
160 my @paralogsB; # List of paralog IDs in given cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
161 my @confPA; # Confidence values for A paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
162 my @confPB; # Confidence values for B paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
163 my @confA; # Confidence values for orthologous groups
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
164 my @confB; # Confidence values for orthologous groups
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
165 my $prev_time = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
166
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
167 $outputfile = "Output." . $ARGV[0] . "-" . $ARGV[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
168 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
169 open OUTPUT, ">$outputfile" or warn "Could not write to OUTPUT file $filename\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
170 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
171
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
172 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
173 # Assign ID numbers for species A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
174 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
175 open A, "$fasta_seq_fileA" or die "File A with sequences in FASTA format is missing
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
176 Usage $0 <FASTAFILE with sequences of species A> <FASTAFILE with sequences of species B> <FASTAFILE with sequences of species C>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
177 $id = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
178 while (<A>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
179 if(/^\>/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
180 ++$id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
181 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
182 s/\>//;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
183 @tmp = split(/\s+/);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
184 #$name = substr($tmp[0],0,25);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
185 $name = $tmp[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
186 $idA{$name} = int($id);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
187 $nameA[$id] = $name;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
188 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
189 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
190 close A;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
191 $A = $id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
192 print "$A sequences in file $fasta_seq_fileA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
193
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
194 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
195 print OUTPUT "$A sequences in file $fasta_seq_fileA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
196 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
197
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
198 if (@ARGV >= 2) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
199 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
200 # Assign ID numbers for species B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
201 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
202 open B, "$fasta_seq_fileB" or die "File B with sequences in FASTA format is missing
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
203 Usage $0 <FASTAFILE with sequences of species A> <FASTAFILE with sequences of species B> <FASTAFILE with sequences of species C>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
204 $id = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
205 while (<B>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
206 if(/^\>/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
207 ++$id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
208 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
209 s/\>//;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
210 @tmp = split(/\s+/);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
211 #$name = substr($tmp[0],0,25);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
212 $name = $tmp[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
213 $idB{$name} = int($id);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
214 $nameB[$id] = $name;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
215 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
216 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
217 $B = $id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
218 print "$B sequences in file $fasta_seq_fileB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
219 close B;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
220
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
221 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
222 print OUTPUT "$B sequences in file $fasta_seq_fileB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
223 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
224 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
225 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
226 # Assign ID numbers for species C (outgroup)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
227 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
228 if ($use_outgroup){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
229 open C, "$fasta_seq_fileC" or die "File C with sequences in FASTA format is missing
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
230 Usage $0 <FASTAFILE with sequences of species A> <FASTAFILE with sequences of species B> <FASTAFILE with sequences of species C>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
231 $id = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
232 while (<C>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
233 if(/^\>/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
234 ++$id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
235 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
236 s/\>//;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
237 @tmp = split(/\s+/);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
238 #$name = substr($tmp[0],0,25);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
239 $name = $tmp[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
240 $idC{$name} = int($id);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
241 $nameC[$id] = $name;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
242 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
243 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
244 $C = $id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
245 print "$C sequences in file $fasta_seq_fileC\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
246 close C;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
247 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
248 print OUTPUT "$C sequences in file $fasta_seq_fileC\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
249 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
250 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
251 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
252 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
253 printf ("Indexing sequences took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
254 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
255 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
256
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
257 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
258 # Run BLAST if not done already
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
259 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
260 if ($run_blast){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
261 print "Trying to run BLAST now - this may take several hours ... or days in worst case!\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
262 print STDERR "Formatting BLAST databases\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
263 system ("$formatdb -i $fasta_seq_fileA");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
264 system ("$formatdb -i $fasta_seq_fileB") if (@ARGV >= 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
265 system ("$formatdb -i $fasta_seq_fileC") if ($use_outgroup);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
266 print STDERR "Done formatting\nStarting BLAST searches...\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
267
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
268 # Run blast only if the files do not already exist is not default.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
269 # NOTE: you should have done this beforehand, because you probably
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
270 # want two-pass blasting anyway which is not implemented here
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
271 # this is also not adapted to use specific compositional adjustment settings
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
272 # and might not use the proper blast parser...
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
273
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
274 do_blast ($fasta_seq_fileA, $fasta_seq_fileA, $A, $A, $blast_outputAA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
275
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
276 if (@ARGV >= 2) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
277 do_blast ($fasta_seq_fileA, $fasta_seq_fileB, $B, $B, $blast_outputAB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
278 do_blast ($fasta_seq_fileB, $fasta_seq_fileA, $A, $A, $blast_outputBA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
279 do_blast ($fasta_seq_fileB, $fasta_seq_fileB, $B, $B, $blast_outputBB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
280 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
281
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
282 if ($use_outgroup){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
283
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
284 do_blast ($fasta_seq_fileA, $fasta_seq_fileC, $A, $C, $blast_outputAC);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
285 do_blast ($fasta_seq_fileB, $fasta_seq_fileC, $B, $C, $blast_outputBC);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
286 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
287
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
288 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
289 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
290 printf ("BLAST searches took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
291 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
292 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
293 print STDERR "Done BLAST searches. ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
294
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
295 } else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
296 print STDERR "Skipping blast! \n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
297 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
298
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
299 if ($run_inparanoid){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
300 print STDERR "Starting ortholog detection...\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
301 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
302 # Read in best hits from blast output file AB
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
303 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
304 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
305 open AB, "$blast_outputAB" or die "Blast output file A->B is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
306 $old_idQ = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
307 while (<AB>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
308 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
309 @Fld = split(/\s+/); # Get query, match and score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
310
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
311 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
312 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
313 print STDERR "AB ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
314 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
315 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
316 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
317
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
318 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
319 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
320 $idQ = $idA{$q}; # ID of query sequence
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
321 $idM = $idB{$m}; # ID of match sequence
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
322 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
323
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
324 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
325
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
326 # Score must be equal to or above cut-off
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
327 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
328
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
329 if(!$count || $q ne $oldq){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
330 print "Match $m, score $score, ID for $q is missing\n" if ($debug == 2 and !(exists($idA{$q})));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
331 $hitnAB[$idA{$oldq}] = $hit if($count); # Record number of hits for previous query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
332 $hit = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
333 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
334 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
335 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
336 ++$hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
337 $hitAB[$idQ][$hit] = int($idM);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
338 # printf ("hitAB[%d][%d] = %d\n",$idQ,$hit,$idM);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
339 $scoreAB{"$idQ:$idM"} = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
340 $scoreBA{"$idM:$idQ"} = $score_cutoff; # Initialize mutual hit score - sometimes this is below score_cutoff
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
341 $old_idQ = $idQ;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
342 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
343 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
344 $hitnAB[$idQ] = $hit; # For the last query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
345 #printf ("hitnAB[1] = %d\n",$hitnAB[1]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
346 #printf ("hitnAB[%d] = %d\n",$idQ,$hit);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
347 close AB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
348 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
349 print OUTPUT "$count sequences $fasta_seq_fileA have homologs in dataset $fasta_seq_fileB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
350 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
351 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
352 # Read in best hits from blast output file BA
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
353 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
354 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
355 open BA, "$blast_outputBA" or die "Blast output file B->A is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
356 $old_idQ = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
357 while (<BA>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
358 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
359 @Fld = split(/\s+/); # Get query, match and score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
360
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
361 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
362 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
363 print STDERR "BA ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
364 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
365 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
366 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
367
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
368 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
369 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
370 $idQ = $idB{$q};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
371 $idM = $idA{$m};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
372 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
373
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
374 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
375
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
376 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
377
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
378 if(!$count || $q ne $oldq){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
379 print "ID for $q is missing\n" if ($debug == 2 and (!exists($idB{$q})));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
380 $hitnBA[$idB{$oldq}] = $hit if($count); # Record number of hits for previous query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
381 $hit = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
382 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
383 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
384 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
385 ++$hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
386 $hitBA[$idQ][$hit] = int($idM);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
387 # printf ("hitBA[%d][%d] = %d\n",$idQ,$hit,$idM);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
388 $scoreBA{"$idQ:$idM"} = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
389 $scoreAB{"$idM:$idQ"} = $score_cutoff if ($scoreAB{"$idM:$idQ"} < $score_cutoff); # Initialize missing scores
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
390 $old_idQ = $idQ;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
391 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
392 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
393 $hitnBA[$idQ] = $hit; # For the last query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
394 #printf ("hitnBA[%d] = %d\n",$idQ,$hit);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
395 close BA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
396 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
397 print OUTPUT "$count sequences $fasta_seq_fileB have homologs in dataset $fasta_seq_fileA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
398 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
399 ##################### Equalize AB scores and BA scores ##########################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
400
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
401 ###################################################################################################################################### Modification by Isabella 1
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
402
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
403 # I removed the time consuming all vs all search and equalize scores for all pairs where there was a hit
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
404
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
405 foreach my $key (keys %scoreAB) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
406
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
407 my ($a, $b) = split(':', $key);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
408 my $key2 = $b . ':' . $a;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
409
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
410 # If debugg mod is 5 and the scores A-B and B-A are unequal
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
411 # the names of the two sequences and their scores are printed
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
412 if ($scoreAB{$key} != $scoreBA{$key2}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
413 printf ("%-20s\t%-20s\t%d\t%d\n",$nameA[$a], $nameB[$b], $scoreAB{$key}, $scoreBA{$key2}) if ($debug == 5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
414 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
415
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
416 # Set score AB and score BA to the mean of scores AB and BA.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
417 # The final score is saved as an integer so .5 needs to be added to avoid rounding errors
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
418 $scoreAB{$key} = $scoreBA{$key2} = int(($scoreAB{$key} + $scoreBA{$key2})/2.0 +.5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
419 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
420
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
421 # For all ids for sequences from organism A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
422 #for $a(1..$A){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
423 #For all ids for sequences from organism B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
424 #for $b(1..$B){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
425
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
426 # No need to equalize score if there was no match between sequence with id $a from species A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
427 # and sequence with id $b from species B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
428 # next if (!$scoreAB{"$a:$b"});
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
429
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
430 # If debugg mod is 5 and the scores A-B and B-A are unequal
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
431 # the names of the two sequences and their scores are printed
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
432 # if ($scoreAB{"$a:$b"} != $scoreBA{"$b:$a"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
433 # printf ("%-20s\t%-20s\t%d\t%d\n",$nameA[$a], $nameB[$b], $scoreAB{"$a:$b"}, $scoreBA{"$b:$a"}) if ($debug == 5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
434 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
435
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
436 # Set score AB and score BA to the mean of scores AB and BA.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
437 # The final score is saved as an integer so .5 needs to be added to avoid rounding errors
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
438 # $scoreAB{"$a:$b"} = $scoreBA{"$b:$a"} = int(($scoreAB{"$a:$b"} + $scoreBA{"$b:$a"})/2.0 +.5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
439
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
440 # printf ("scoreAB{%d: %d} = %d\n", $a, $b, $scoreAB{"$a:$b"});
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
441 # printf ("scoreBA{%d: %d} = %d\n", $b, $a, $scoreBA{"$a:$b"});
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
442 #}
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
443 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
444
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
445 ####################################################################################################################################### End modification by Isabella 1
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
446
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
447 ##################### Re-sort hits, besthits and bestscore #######################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
448 for $idA(1..$A){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
449 # print "Loop index = $idA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
450 # printf ("hitnAB[%d] = %d\n",$idA, $hitnAB[$idA]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
451 next if (!($hitnAB[$idA]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
452 for $hit (1..($hitnAB[$idA]-1)){ # Sort hits by score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
453 while($scoreAB{"$idA:$hitAB[$idA][$hit]"} < $scoreAB{"$idA:$hitAB[$idA][$hit+1]"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
454 $tmp = $hitAB[$idA][$hit];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
455 $hitAB[$idA][$hit] = $hitAB[$idA][$hit+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
456 $hitAB[$idA][$hit+1] = $tmp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
457 --$hit if ($hit > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
458 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
459 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
460 $bestscore = $bestscoreAB[$idA] = $scoreAB{"$idA:$hitAB[$idA][1]"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
461 $besthitAB[$idA] = $hitAB[$idA][1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
462 for $hit (2..$hitnAB[$idA]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
463 if ($bestscore - $scoreAB{"$idA:$hitAB[$idA][$hit]"} <= $grey_zone){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
464 $besthitAB[$idA] .= " $hitAB[$idA][$hit]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
465 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
466 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
467 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
468 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
469 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
470 undef $is_besthitAB[$idA]; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
471 grep (vec($is_besthitAB[$idA],$_,1) = 1, split(/ /,$besthitAB[$idA]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
472 # printf ("besthitAB[%d] = hitAB[%d][%d] = %d\n",$idA,$idA,$hit,$besthitAB[$idA]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
473
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
474 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
475
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
476 for $idB(1..$B){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
477 # print "Loop index = $idB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
478 next if (!($hitnBA[$idB]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
479 for $hit (1..($hitnBA[$idB]-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
480 # Sort hits by score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
481 while($scoreBA{"$idB:$hitBA[$idB][$hit]"} < $scoreBA{"$idB:$hitBA[$idB][$hit+1]"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
482 $tmp = $hitBA[$idB][$hit];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
483 $hitBA[$idB][$hit] = $hitBA[$idB][$hit+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
484 $hitBA[$idB][$hit+1] = $tmp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
485 --$hit if ($hit > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
486 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
487 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
488 $bestscore = $bestscoreBA[$idB] = $scoreBA{"$idB:$hitBA[$idB][1]"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
489 $besthitBA[$idB] = $hitBA[$idB][1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
490 for $hit (2..$hitnBA[$idB]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
491 if ($bestscore - $scoreBA{"$idB:$hitBA[$idB][$hit]"} <= $grey_zone){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
492 $besthitBA[$idB] .= " $hitBA[$idB][$hit]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
493 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
494 else {last;}
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
495 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
496 undef $is_besthitBA[$idB]; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
497 grep (vec($is_besthitBA[$idB],$_,1) = 1, split(/ /,$besthitBA[$idB]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
498 # printf ("besthitBA[%d] = %d\n",$idA,$besthitAB[$idA]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
499 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
500
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
501 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
502 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
503 printf ("Reading and sorting homologs took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
504 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
505 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
506
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
507 ######################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
508 # Now find orthologs:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
509 ######################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
510 $o = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
511
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
512 for $i(1..$A){ # For each ID in file A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
513 if (defined $besthitAB[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
514 @besthits = split(/ /,$besthitAB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
515 for $hit (@besthits){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
516 if (vec($is_besthitBA[$hit],$i,1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
517 ++$o;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
518 $ortoA[$o] = $i;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
519 $ortoB[$o] = $hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
520 $ortoS[$o] = $scoreAB{"$i:$hit"}; # Should be equal both ways
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
521 # --$o if ($ortoS[$o] == $score_cutoff); # Ignore orthologs that are exactly at score_cutoff
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
522 print "Accept! " if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
523 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
524 else {print " " if ($debug == 2);}
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
525 printf ("%-20s\t%d\t%-20s\t", $nameA[$i], $bestscoreAB[$i], $nameB[$hit]) if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
526 print "$bestscoreBA[$hit]\t$besthitBA[$hit]\n" if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
527 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
528 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
529 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
530 print "$o ortholog candidates detected\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
531 #####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
532 # Sort orthologs by ID and then by score:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
533 #####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
534
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
535 ####################################################################################################### Modification by Isabella 2
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
536
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
537 # Removed time consuiming bubble sort. Created an index array and sort that according to id and score.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
538 # The put all clusters on the right place.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
539
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
540 # Create an array used to store the position each element shall have in the final array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
541 # The elements are initialized with the position numbers
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
542 my @position_index_array = (1..$o);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
543
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
544 # Sort the position list according to id
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
545 my @id_sorted_position_list = sort { ($ortoA[$a]+$ortoB[$a]) <=> ($ortoA[$b] + $ortoB[$b]) } @position_index_array;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
546
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
547 # Sort the list according to score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
548 my @score_id_sorted_position_list = sort { $ortoS[$b] <=> $ortoS[$a] } @id_sorted_position_list;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
549
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
550 # Create new arrays for the sorted information
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
551 my @new_ortoA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
552 my @new_ortoB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
553 my @new_orthoS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
554
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
555 # Add the information to the new arrays in the orer specifeid by the index array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
556 for (my $index_in_list = 0; $index_in_list < scalar @score_id_sorted_position_list; $index_in_list++) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
557
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
558
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
559 my $old_index = $score_id_sorted_position_list[$index_in_list];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
560 $new_ortoA[$index_in_list + 1] = $ortoA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
561 $new_ortoB[$index_in_list + 1] = $ortoB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
562 $new_ortoS[$index_in_list + 1] = $ortoS[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
563 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
564
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
565 @ortoA = @new_ortoA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
566 @ortoB = @new_ortoB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
567 @ortoS = @new_ortoS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
568
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
569 # Use bubblesort to sort ortholog pairs by id
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
570 # for $i(1..($o-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
571 # while(($ortoA[$i]+$ortoB[$i]) > ($ortoA[$i+1] + $ortoB[$i+1])){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
572 # $tempA = $ortoA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
573 # $tempB = $ortoB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
574 # $tempS = $ortoS[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
575 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
576 # $ortoA[$i] = $ortoA[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
577 # $ortoB[$i] = $ortoB[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
578 # $ortoS[$i] = $ortoS[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
579 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
580 # $ortoA[$i+1] = $tempA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
581 # $ortoB[$i+1] = $tempB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
582 # $ortoS[$i+1] = $tempS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
583 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
584 # --$i if ($i > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
585 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
586 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
587 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
588 # # Use bubblesort to sort ortholog pairs by score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
589 # for $i(1..($o-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
590 # while($ortoS[$i] < $ortoS[$i+1]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
591 # # Swap places:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
592 # $tempA = $ortoA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
593 # $tempB = $ortoB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
594 # $tempS = $ortoS[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
595 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
596 # $ortoA[$i] = $ortoA[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
597 # $ortoB[$i] = $ortoB[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
598 # $ortoS[$i] = $ortoS[$i+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
599 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
600 # $ortoA[$i+1] = $tempA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
601 # $ortoB[$i+1] = $tempB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
602 # $ortoS[$i+1] = $tempS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
603 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
604 # --$i if ($i > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
605 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
606 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
607
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
608 ###################################################################################################### End modification bt Isabella 2
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
609
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
610 @all_ortologsA = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
611 @all_ortologsB = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
612 for $i(1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
613 push(@all_ortologsA,$ortoA[$i]); # List of all orthologs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
614 push(@all_ortologsB,$ortoB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
615 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
616 undef $is_ortologA; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
617 undef $is_ortologB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
618 grep (vec($is_ortologA,$_,1) = 1, @all_ortologsA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
619 grep (vec($is_ortologB,$_,1) = 1, @all_ortologsB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
620
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
621 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
622 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
623 printf ("Finding and sorting orthologs took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
624 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
625 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
626 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
627 # Read in best hits from blast output file AC
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
628 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
629 if ($use_outgroup){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
630 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
631 open AC, "$blast_outputAC" or die "Blast output file A->C is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
632 while (<AC>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
633 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
634 @Fld = split(/\s+/); # Get query, match and score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
635
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
636 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
637 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
638 print STDERR "AC ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
639 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
640 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
641 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
642
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
643 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
644 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
645 $idQ = $idA{$q};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
646 $idM = $idC{$m};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
647 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
648 next unless (vec($is_ortologA,$idQ,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
649
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
650 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
651
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
652 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
653
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
654 next if ($count and ($q eq $oldq));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
655 # Only comes here if this is the best hit:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
656 $besthitAC[$idQ] = $idM;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
657 $bestscoreAC[$idQ] = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
658 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
659 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
660 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
661 close AC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
662 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
663 # Read in best hits from blast output file BC
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
664 #################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
665 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
666 open BC, "$blast_outputBC" or die "Blast output file B->C is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
667 while (<BC>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
668 chomp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
669 @Fld = split(/\s+/); # Get query, match and score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
670
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
671 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
672 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
673 print STDERR "BC ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
674 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
675 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
676 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
677
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
678 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
679 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
680 $idQ = $idB{$q};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
681 $idM = $idC{$m};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
682 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
683 next unless (vec($is_ortologB,$idQ,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
684
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
685 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
686
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
687 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
688
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
689 next if ($count and ($q eq $oldq));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
690 # Only comes here if this is the best hit:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
691 $besthitBC[$idQ] = $idM;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
692 $bestscoreBC[$idQ] = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
693 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
694 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
695 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
696 close BC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
697 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
698 # Detect rooting problems
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
699 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
700 $rejected = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
701 @del = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
702 $file = "rejected_sequences." . $fasta_seq_fileC;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
703 open OUTGR, ">$file";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
704 for $i (1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
705 $diff1 = $diff2 = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
706 $idA = $ortoA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
707 $idB = $ortoB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
708 $diff1 = $bestscoreAC[$idA] - $ortoS[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
709 $diff2 = $bestscoreBC[$idB] - $ortoS[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
710 if ($diff1 > $outgroup_cutoff){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
711 print OUTGR "Ortholog pair $i ($nameA[$idA]-$nameB[$idB]).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
712 $nameA[$idA] from $fasta_seq_fileA is closer to $nameC[$besthitAC[$idA]] than to $nameB[$idB]\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
713 print OUTGR " $ortoS[$i] < $bestscoreAC[$idA] by $diff1\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
714 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
715 if ($diff2 > $outgroup_cutoff){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
716 print OUTGR "Ortholog pair $i ($nameA[$idA]-$nameB[$idB]).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
717 $nameB[$idB] from $fasta_seq_fileB is closer to $nameC[$besthitBC[$idB]] than to $nameA[$idA]\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
718 print OUTGR " $ortoS[$i] < $bestscoreBC[$idB] by $diff2\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
719 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
720 if (($diff1 > $outgroup_cutoff) or ($diff2 > $outgroup_cutoff)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
721 ++$rejected;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
722 $del[$i] = 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
723 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
724 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
725 print "Number of rejected groups: $rejected (outgroup sequence was closer by more than $outgroup_cutoff bits)\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
726 close OUTGR;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
727 } # End of $use_outgroup
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
728 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
729 # Read inside scores from AA
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
730 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
731 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
732 $max_hit = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
733 open AA, "$blast_outputAA" or die "Blast output file A->A is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
734 while (<AA>) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
735 chomp; # strip newline
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
736
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
737 @Fld = split(/\s+/); # Get query and match names
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
738
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
739 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
740 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
741 print STDERR "AA ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
742 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
743 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
744 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
745
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
746 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
747 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
748 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
749 next unless (vec($is_ortologA,$idA{$q},1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
750
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
751 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
752
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
753 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
754
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
755 if(!$count || $q ne $oldq){ # New query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
756 $max_hit = $hit if ($hit > $max_hit);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
757 $hit = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
758 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
759 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
760 ++$hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
761 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
762 $scoreAA{"$idA{$q}:$idA{$m}"} = int($score + 0.5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
763 $hitAA[$idA{$q}][$hit] = int($idA{$m});
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
764 $hitnAA[$idA{$q}] = $hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
765 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
766 close AA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
767 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
768 print OUTPUT "$count $fasta_seq_fileA-$fasta_seq_fileA matches\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
769 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
770 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
771 # Read inside scores from BB
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
772 ################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
773 $count = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
774 open BB, "$blast_outputBB" or die "Blast output file B->B is missing\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
775 while (<BB>) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
776 chomp; # strip newline
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
777
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
778 @Fld = split(/\s+/); # Get query and match names
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
779
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
780 if( scalar @Fld < 9 ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
781 if($Fld[0]=~/done/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
782 print STDERR "BB ok\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
783 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
784 next;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
785 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
786
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
787 $q = $Fld[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
788 $m = $Fld[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
789 $score = $Fld[2];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
790 next unless (vec($is_ortologB,$idB{$q},1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
791
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
792 next if (!overlap_test(@Fld));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
793
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
794 next if ($score < $score_cutoff);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
795
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
796 if(!$count || $q ne $oldq){ # New query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
797 $max_hit = $hit if ($hit > $max_hit);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
798 $oldq = $q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
799 $hit = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
800 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
801 ++$count;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
802 ++$hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
803 $scoreBB{"$idB{$q}:$idB{$m}"} = int($score + 0.5);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
804 $hitBB[$idB{$q}][$hit] = int($idB{$m});
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
805 $hitnBB[$idB{$q}] = $hit;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
806 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
807 close BB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
808 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
809 print OUTPUT "$count $fasta_seq_fileB-$fasta_seq_fileB matches\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
810 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
811 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
812 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
813 printf ("Reading paralogous hits took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
814 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
815 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
816 print "Maximum number of hits per sequence was $max_hit\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
817 #####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
818 # Find paralogs:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
819 #####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
820 for $i(1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
821 $merge[$i] = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
822 next if($del[$i]); # If outgroup species was closer to one of the seed orthologs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
823 $idA = $ortoA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
824 $idB = $ortoB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
825 local @membersA = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
826 local @membersB = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
827
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
828 undef $is_paralogA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
829 undef $is_paralogB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
830
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
831 print "$i: Ortholog pair $nameA[$idA] and $nameB[$idB]. $hitnAA[$idA] hits for A and $hitnBB[$idB] hits for B\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
832 # Check if current ortholog is already clustered:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
833 for $j(1..($i-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
834 # Overlap type 1: Both orthologs already clustered here -> merge
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
835 if ((vec($is_paralogA[$j],$idA,1)) and (vec($is_paralogB[$j],$idB,1))){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
836 $merge[$i] = $j;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
837 print "Merge CASE 1: group $i ($nameB[$idB]-$nameA[$idA]) and $j ($nameB[$ortoB[$j]]-$nameA[$ortoA[$j]])\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
838 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
839 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
840 # Overlap type 2: 2 competing ortholog pairs -> merge
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
841 elsif (($ortoS[$j] - $ortoS[$i] <= $grey_zone)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
842 and (($ortoA[$j] == $ortoA[$i]) or ($ortoB[$j] == $ortoB[$i]))
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
843 # and ($paralogsA[$j])
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
844 ){ # The last condition is false if the previous cluster has been already deleted
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
845 $merge[$i] = $j;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
846 print "Merge CASE 2: group $i ($nameA[$ortoA[$i]]-$nameB[$ortoB[$i]]) and $j ($nameA[$ortoA[$j]]-$nameB[$ortoB[$j]])\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
847 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
848 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
849 # Overlap type 3: DELETE One of the orthologs belongs to some much stronger cluster -> delete
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
850 elsif (((vec($is_paralogA[$j],$idA,1)) or (vec($is_paralogB[$j],$idB,1))) and
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
851 ($ortoS[$j] - $ortoS[$i] > $score_cutoff)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
852 print "Delete CASE 3: Cluster $i -> $j, score $ortoS[$i] -> $ortoS[$j], ($nameA[$ortoA[$j]]-$nameB[$ortoB[$j]])\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
853 $merge[$i]= -1; # Means - do not add sequences to this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
854 $paralogsA[$i] = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
855 $paralogsB[$i] = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
856 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
857 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
858 # Overlap type 4: One of the orthologs is close to the center of other cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
859 elsif (((vec($is_paralogA[$j],$idA,1)) and ($confPA[$idA] > $group_overlap_cutoff)) or
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
860 ((vec($is_paralogB[$j],$idB,1)) and ($confPB[$idB] > $group_overlap_cutoff))){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
861 print "Merge CASE 4: Cluster $i -> $j, score $ortoS[$i] -> $ortoS[$j], ($nameA[$ortoA[$j]]-$nameB[$ortoB[$j]])\n" if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
862 $merge[$i] = $j;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
863 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
864 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
865 # Overlap type 5:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
866 # All clusters that were overlapping, but not catched by previous "if" statements will be DIVIDED!
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
867 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
868 next if ($merge[$i] < 0); # This cluster should be deleted
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
869 ##### Check for paralogs in A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
870 $N = $hitnAA[$idA];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
871 for $j(1..$N){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
872 $hitID = $hitAA[$idA][$j]; # hit of idA
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
873 # print "Working with $nameA[$hitID]\n" if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
874 # Decide whether this hit is inside the paralog circle:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
875 if ( ($idA == $hitID) or ($scoreAA{"$idA:$hitID"} >= $bestscoreAB[$idA]) and
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
876 ($scoreAA{"$idA:$hitID"} >= $bestscoreAB[$hitID])){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
877 if ($debug == 2){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
878 print " Paralog candidates: ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
879 printf ("%-20s: %-20s", $nameA[$idA], $nameA[$hitID]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
880 print "\t$scoreAA{\"$idA:$hitID\"} : $bestscoreAB[$idA] : $bestscoreAB[$hitID]\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
881 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
882 $paralogs = 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
883 if ($scoreAA{"$idA:$idA"} == $ortoS[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
884 if ($scoreAA{"$idA:$hitID"} == $scoreAA{"$idA:$idA"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
885 $conf_here = 1.0; # In the center
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
886 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
887 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
888 $conf_here = 0.0; # On the border
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
889 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
890 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
891 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
892 $conf_here = ($scoreAA{"$idA:$hitID"} - $ortoS[$i]) /
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
893 ($scoreAA{"$idA:$idA"} - $ortoS[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
894 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
895 # Check if this paralog candidate is already clustered in other clusters
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
896 for $k(1..($i-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
897 if (vec($is_paralogA[$k],$hitID,1)){ # Yes, found in cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
898 if($debug == 2){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
899 print " $nameA[$hitID] is already in cluster $k, together with:";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
900 print " $nameA[$ortoA[$k]] and $nameB[$ortoB[$k]] ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
901 print "($scoreAA{\"$ortoA[$k]:$hitID\"})";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
902 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
903 if (($confPA[$hitID] >= $conf_here) and
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
904 ($j != 1)){ # The seed ortholog CAN NOT remain there
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
905 print " and remains there.\n" if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
906 $paralogs = 0; # No action
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
907 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
908 else { # Ortholog of THIS cluster is closer than ortholog of competing cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
909 print " and should be moved here!\n" if ($debug == 2); # Remove from other cluster, add to this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
910 @membersAK = split(/ /, $paralogsA[$k]); # This array contains IDs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
911 $paralogsA[$k] = "";# Remove all paralogs from cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
912 @tmp = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
913 for $m(@membersAK){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
914 push(@tmp,$m) if ($m != $hitID); # Put other members back
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
915 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
916 $paralogsA[$k] = join(' ',@tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
917 undef $is_paralogA[$k]; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
918 grep (vec($is_paralogA[$k],$_,1) = 1, @tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
919 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
920 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
921 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
922 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
923 next if (! $paralogs); # Skip that paralog - it is already in cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
924 push (@membersA,$hitID); # Add this hit to paralogs of A
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
925 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
926 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
927 # Calculate confidence values now:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
928 @tmp = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
929 for $idP (@membersA){ # For each paralog calculate conf value
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
930 if($scoreAA{"$idA:$idA"} == $ortoS[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
931 if ($scoreAA{"$idA:$idP"} == $scoreAA{"$idA:$idA"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
932 $confPA[$idP] = 1.00;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
933 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
934 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
935 $confPA[$idP] = 0.00;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
936 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
937 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
938 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
939 $confPA[$idP] = ($scoreAA{"$idA:$idP"} - $ortoS[$i]) /
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
940 ($scoreAA{"$idA:$idA"} - $ortoS[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
941 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
942 push (@tmp, $idP) if ($confPA[$idP] >= $conf_cutoff); # If one wishes to use only significant paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
943 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
944 @membersA = @tmp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
945 ########### Merge if necessary:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
946 if ($merge[$i] > 0){ # Merge existing cluster with overlapping cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
947 @tmp = split(/ /,$paralogsA[$merge[$i]]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
948 for $m (@membersA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
949 push (@tmp, $m) unless (vec($is_paralogA[$merge[$i]],$m,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
950 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
951 $paralogsA[$merge[$i]] = join(' ',@tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
952 undef $is_paralogA[$merge[$i]];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
953 grep (vec($is_paralogA[$merge[$i]],$_,1) = 1, @tmp); # Refresh index of paralog array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
954 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
955 ######### Typical new cluster:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
956 else{ # Create a new cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
957 $paralogsA[$i] = join(' ',@membersA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
958 undef $is_paralogA; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
959 grep (vec($is_paralogA[$i],$_,1) = 1, @membersA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
960 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
961 ##### The same procedure for species B:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
962 $N = $hitnBB[$idB];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
963 for $j(1..$N){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
964 $hitID = $hitBB[$idB][$j];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
965 # print "Working with $nameB[$hitID]\n" if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
966 if ( ($idB == $hitID) or ($scoreBB{"$idB:$hitID"} >= $bestscoreBA[$idB]) and
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
967 ($scoreBB{"$idB:$hitID"} >= $bestscoreBA[$hitID])){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
968 if ($debug == 2){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
969 print " Paralog candidates: ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
970 printf ("%-20s: %-20s", $nameB[$idB], $nameB[$hitID]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
971 print "\t$scoreBB{\"$idB:$hitID\"} : ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
972 print "$bestscoreBA[$idB] : $bestscoreBA[$hitID]\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
973 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
974 $paralogs = 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
975 if ($scoreBB{"$idB:$idB"} == $ortoS[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
976 if ($scoreBB{"$idB:$hitID"} == $scoreBB{"$idB:$idB"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
977 $conf_here = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
978 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
979 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
980 $conf_here = 0.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
981 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
982 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
983 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
984 $conf_here = ($scoreBB{"$idB:$hitID"} - $ortoS[$i]) /
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
985 ($scoreBB{"$idB:$idB"} - $ortoS[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
986 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
987
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
988 # Check if this paralog candidate is already clustered in other clusters
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
989 for $k(1..($i-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
990 if (vec($is_paralogB[$k],$hitID,1)){ # Yes, found in cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
991 if($debug == 2){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
992 print " $nameB[$hitID] is already in cluster $k, together with:";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
993 print " $nameB[$ortoB[$k]] and $nameA[$ortoA[$k]] ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
994 print "($scoreBB{\"$ortoB[$k]:$hitID\"})";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
995 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
996 if (($confPB[$hitID] >= $conf_here) and
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
997 ($j != 1)){ # The seed ortholog CAN NOT remain there
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
998 print " and remains there.\n" if ($debug == 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
999 $paralogs = 0; # No action
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1000 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1001 else { # Ortholog of THIS cluster is closer than ortholog of competing cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1002 print " and should be moved here!\n" if ($debug == 2); # Remove from other cluster, add to this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1003 @membersBK = split(/ /, $paralogsB[$k]); # This array contains names, not IDs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1004 $paralogsB[$k] = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1005 @tmp = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1006 for $m(@membersBK){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1007 push(@tmp,$m) if ($m != $hitID); # Put other members back
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1008 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1009 $paralogsB[$k] = join(' ',@tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1010 undef $is_paralogB[$k]; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1011 grep (vec($is_paralogB[$k],$_,1) = 1, @tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1012 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1013 last; # Don't search in other clusters
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1014 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1015 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1016 next if (! $paralogs); # Skip that paralog - it is already in cluster $k
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1017 push (@membersB,$hitID);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1018 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1019 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1020 # Calculate confidence values now:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1021 @tmp = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1022 for $idP (@membersB){ # For each paralog calculate conf value
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1023 if($scoreBB{"$idB:$idB"} == $ortoS[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1024 if ($scoreBB{"$idB:$idP"} == $scoreBB{"$idB:$idB"}){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1025 $confPB[$idP] = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1026 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1027 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1028 $confPB[$idP] = 0.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1029 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1030 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1031 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1032 $confPB[$idP] = ($scoreBB{"$idB:$idP"} - $ortoS[$i]) /
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1033 ($scoreBB{"$idB:$idB"} - $ortoS[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1034 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1035 push (@tmp, $idP) if ($confPB[$idP] >= $conf_cutoff); # If one wishes to use only significant paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1036 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1037 @membersB = @tmp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1038 ########### Merge if necessary:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1039 if ($merge[$i] > 0){ # Merge existing cluster with overlapping cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1040 @tmp = split(/ /,$paralogsB[$merge[$i]]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1041 for $m (@membersB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1042 push (@tmp, $m) unless (vec($is_paralogB[$merge[$i]],$m,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1043 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1044 $paralogsB[$merge[$i]] = join(' ',@tmp);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1045 undef $is_paralogB[$merge[$i]];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1046 grep (vec($is_paralogB[$merge[$i]],$_,1) = 1, @tmp); # Refresh index of paralog array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1047 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1048 ######### Typical new cluster:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1049 else{ # Create a new cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1050 $paralogsB[$i] = join(' ',@membersB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1051 undef $is_paralogB; # Create index that we can check later
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1052 grep (vec($is_paralogB[$i],$_,1) = 1, @membersB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1053 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1054 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1055 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1056 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1057 printf ("Finding in-paralogs took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1058 $prev_time = $user_time;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1059 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1060 #####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1061 &clean_up(1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1062 ####################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1063 # Find group for orphans. If cluster contains only one member, find where it should go:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1064 for $i (1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1065 @membersA = split(/ /, $paralogsA[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1066 @membersB = split(/ /, $paralogsB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1067 $na = @membersA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1068 $nb = @membersB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1069 if (($na == 0) and $nb){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1070 print "Warning: empty A cluster $i\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1071 for $m (@membersB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1072 $bestscore = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1073 $bestgroup = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1074 $bestmatch = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1075 for $j (1..$o) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1076 next if ($i == $j); # Really need to check against all 100% members of the group.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1077 @membersBJ = split(/ /, $paralogsB[$j]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1078 for $k (@membersBJ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1079 next if ($confPB[$k] != 1); # For all 100% in-paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1080 $score = $scoreBB{"$m:$k"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1081 if ($score > $bestscore){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1082 $bestscore = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1083 $bestgroup = $j;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1084 $bestmatch = $k;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1085 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1086 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1087 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1088 print "Orphan $nameB[$m] goes to group $bestgroup with $nameB[$bestmatch]\n" ;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1089 @members = split(/ /, $paralogsB[$bestgroup]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1090 push (@members, $m);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1091 $paralogsB[$bestgroup] = join(' ',@members);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1092 $paralogsB[$i] = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1093 undef $is_paralogB[$bestgroup];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1094 undef $is_paralogB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1095 grep (vec($is_paralogB[$bestgroup],$_,1) = 1, @members); # Refresh index of paralog array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1096 # grep (vec($is_paralogB[$i],$_,1) = 1, ());
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1097 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1098 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1099 if ($na and ($nb == 0)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1100 print "Warning: empty B cluster $i\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1101 for $m (@membersA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1102 $bestscore = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1103 $bestgroup = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1104 $bestmatch = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1105 for $j (1..$o) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1106 next if ($i == $j);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1107 @membersAJ = split(/ /, $paralogsA[$j]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1108 for $k (@membersAJ){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1109 next if ($confPA[$k] != 1); # For all 100% in-paralogs
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1110 $score = $scoreAA{"$m:$k"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1111 if ($score > $bestscore){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1112 $bestscore = $score;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1113 $bestgroup = $j;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1114 $bestmatch = $k;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1115 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1116 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1117 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1118 print "Orphan $nameA[$m] goes to group $bestgroup with $nameA[$bestmatch]\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1119 @members = split(/ /, $paralogsA[$bestgroup]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1120 push (@members, $m);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1121 $paralogsA[$bestgroup] = join(' ',@members);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1122 $paralogsA[$i] = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1123 undef $is_paralogA[$bestgroup];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1124 undef $is_paralogA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1125 grep (vec($is_paralogA[$bestgroup],$_,1) = 1, @members); # Refresh index of paralog array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1126 # grep (vec($is_paralogA[$i],$_,1) = 1, ());
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1127 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1128 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1129 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1130
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1131 &clean_up(1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1132 ###################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1133 $htmlfile = "orthologs." . $ARGV[0] . "-" . $ARGV[1] . ".html";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1134 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1135 open HTML, ">$htmlfile" or warn "Could not write to HTML file $filename\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1136 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1137
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1138
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1139 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1140 print OUTPUT "\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1141 print OUTPUT "$o groups of orthologs\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1142 print OUTPUT "$totalA in-paralogs from $fasta_seq_fileA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1143 print OUTPUT "$totalB in-paralogs from $fasta_seq_fileB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1144 print OUTPUT "Grey zone $grey_zone bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1145 print OUTPUT "Score cutoff $score_cutoff bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1146 print OUTPUT "In-paralogs with confidence less than $conf_cutoff not shown\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1147 print OUTPUT "Sequence overlap cutoff $seq_overlap_cutoff\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1148 print OUTPUT "Group merging cutoff $group_overlap_cutoff\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1149 print OUTPUT "Scoring matrix $matrix\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1150 print OUTPUT "\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1151 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1152 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1153 print HTML "<pre>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1154 print HTML "\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1155 print HTML "$o groups of orthologs\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1156 print HTML "$totalA in-paralogs from $fasta_seq_fileA\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1157 print HTML "$totalB in-paralogs from $fasta_seq_fileB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1158 print HTML "Grey zone $grey_zone bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1159 print HTML "Score cutoff $score_cutoff bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1160 print HTML "In-paralogs with confidence less than $conf_cutoff not shown\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1161 print HTML "Sequence overlap cutoff $seq_overlap_cutoff\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1162 print HTML "Group merging cutoff $group_overlap_cutoff\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1163 print HTML "Scoring matrix $matrix\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1164 print HTML "\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1165 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1166 # ##############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1167 # Check for alternative orthologs, sort paralogs by confidence and print results
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1168 # ##############################################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1169 if ($use_bootstrap and $debug){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1170 open FF, ">BS_vs_bits" or warn "Could not write to file BS_vs_bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1171 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1172 for $i(1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1173 @membersA = split(/ /, $paralogsA[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1174 @membersB = split(/ /, $paralogsB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1175 $message = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1176 $htmlmessage = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1177
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1178 $idB = $ortoB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1179 $nB = $hitnBA[$idB];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1180 for $idA(@membersA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1181 next if ($confPA[$idA] != 1.0);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1182 $nA = $hitnAB[$idA];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1183 $confA[$i] = $ortoS[$i]; # default
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1184 $bsA[$idA] = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1185 ##############
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1186 for $j(1..$nB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1187 $idH = $hitBA[$idB][$j];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1188 ################ Some checks for alternative orthologs:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1189 # 1. Don't consider sequences that are already in this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1190 next if (vec($is_paralogA[$i],$idH,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1191 next if ($confPA[$idH] > 0); # If $conf_cutoff > 0 idH might be incide circle, but not paralog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1192
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1193 # 2. Check if candidate for alternative ortholog is already clustered in stronger clusters
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1194 $in_other_cluster = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1195 for $k(1..($i-1)){ # Check if current ortholog is already clustered
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1196 if (vec($is_paralogA[$k],$idH,1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1197 $in_other_cluster = $k;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1198 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1199 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1200 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1201 # next if ($in_other_cluster); # This hit is clustered in cluster $k. It cannot be alternative ortholog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1202
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1203 # 3. The best hit of candidate ortholog should be ortoA or at least to belong into this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1204 @besthits = split (/ /,$besthitAB[$idH]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1205 $this_family = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1206 for $bh (@besthits){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1207 $this_family = 1 if ($idB == $bh);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1208 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1209 # next unless ($this_family); # There was an alternative BA match but it's best match did not belong here
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1210 ################# Done with checks - if sequence passed, then it could be an alternative ortholog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1211 $confA[$i] = $ortoS[$i] - $scoreBA{"$idB:$idH"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1212 if ($use_bootstrap){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1213 if ($confA[$i] < $ortoS[$i]){ # Danger zone - check for bootstrap
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1214 $bsA[$idA] = &bootstrap($fasta_seq_fileB,$idB,$idA,$idH);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1215 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1216 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1217 $bsA[$idA] = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1218 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1219 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1220 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1221 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1222 $message .= sprintf("Bootstrap support for %s as seed ortholog is %d%%.", $nameA[$idA], 100*$bsA[$idA]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1223 $message .= sprintf(" Alternative seed ortholog is %s (%d bits away from this cluster)", $nameA[$idH], $confA[$i]) if ($bsA[$idA] < 0.75);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1224 $message .= sprintf("\n");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1225 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1226 if ($bsA[$idA] < 0.75){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1227 $htmlmessage .= sprintf("<font color=\"red\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1228 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1229 elsif ($bsA[$idA] < 0.95){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1230 $htmlmessage .= sprintf("<font color=\"\#FFCC00\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1231 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1232 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1233 $htmlmessage .= sprintf("<font color=\"green\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1234 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1235 $htmlmessage .= sprintf("Bootstrap support for %s as seed ortholog is %d%%.\n", $nameA[$idA], 100*$bsA[$idA]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1236 $htmlmessage .= sprintf("Alternative seed ortholog is %s (%d bits away from this cluster)\n", $nameA[$idH], $confA[$i]) if ($bsA[$idA] < 0.75);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1237 $htmlmessage .= sprintf("</font>");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1238 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1239 printf (FF "%s\t%d\t%d\n", $nameA[$idA], $confA[$i], 100*$bsA[$idA]) if ($use_bootstrap and $debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1240 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1241 ########
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1242 $idA = $ortoA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1243 $nA = $hitnAB[$idA];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1244 for $idB(@membersB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1245 next if ($confPB[$idB] != 1.0);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1246 $nB = $hitnBA[$idB];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1247 $confB[$i] = $ortoS[$i]; # default
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1248 $bsB[$idB] = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1249
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1250 for $j(1..$nA){ # For all AB hits of given ortholog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1251 $idH = $hitAB[$idA][$j];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1252 # ############### Some checks for alternative orthologs:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1253 # 1. Don't consider sequences that are already in this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1254 next if (vec($is_paralogB[$i],$idH,1));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1255 next if ($confPB[$idH] > 0); # If $conf_cutoff > 0 idH might be incide circle, but not paralog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1256
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1257 # 2. Check if candidate for alternative ortholog is already clustered in stronger clusters
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1258 $in_other_cluster = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1259 for $k(1..($i-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1260 if (vec($is_paralogB[$k],$idH,1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1261 $in_other_cluster = $k;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1262 last; # out from this cycle
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1263 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1264 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1265 # next if ($in_other_cluster); # This hit is clustered in cluster $k. It cannot be alternative ortholog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1266
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1267 # 3. The best hit of candidate ortholog should be ortoA
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1268 @besthits = split (/ /,$besthitBA[$idH]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1269 $this_family = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1270 for $bh (@besthits){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1271 $this_family = 1 if ($idA == $bh);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1272 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1273 # next unless ($this_family); # There was an alternative BA match but it's best match did not belong here
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1274 # ################ Done with checks - if sequence passed, then it could be an alternative ortholog
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1275 $confB[$i] = $ortoS[$i] - $scoreAB{"$idA:$idH"};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1276 if ($use_bootstrap){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1277 if ($confB[$i] < $ortoS[$i]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1278 $bsB[$idB] = &bootstrap($fasta_seq_fileA,$idA,$idB,$idH);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1279 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1280 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1281 $bsB[$idB] = 1.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1282 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1283 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1284 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1285 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1286 $message .= sprintf("Bootstrap support for %s as seed ortholog is %d%%.", $nameB[$idB], 100*$bsB[$idB]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1287 $message .= sprintf(" Alternative seed ortholog is %s (%d bits away from this cluster)", $nameB[$idH],$confB[$i]) if ($bsB[$idB] < 0.75);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1288 $message .= sprintf("\n");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1289 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1290 if ($bsB[$idB] < 0.75){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1291 $htmlmessage .= sprintf("<font color=\"red\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1292 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1293 elsif ($bsB[$idB] < 0.95){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1294 $htmlmessage .= sprintf("<font color=\"\#FFCC00\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1295 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1296 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1297 $htmlmessage .= sprintf("<font color=\"green\">");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1298 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1299 $htmlmessage .= sprintf("Bootstrap support for %s as seed ortholog is %d%%.\n", $nameB[$idB], 100*$bsB[$idB]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1300 $htmlmessage .= sprintf("Alternative seed ortholog is %s (%d bits away from this cluster)\n", $nameB[$idH],$confB[$i]) if ($bsB[$idB] < 0.75);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1301 $htmlmessage .= sprintf("</font>");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1302 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1303 printf (FF "%s\t%d\t%d\n", $nameB[$idB], $confB[$i], 100*$bsB[$idB]) if ($use_bootstrap and $debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1304 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1305 close FF;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1306 ########### Print header ###############
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1307 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1308 print OUTPUT "___________________________________________________________________________________\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1309 print OUTPUT "Group of orthologs #" . $i .". Best score $ortoS[$i] bits\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1310 print OUTPUT "Score difference with first non-orthologous sequence - ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1311 printf (OUTPUT "%s:%d %s:%d\n", $fasta_seq_fileA,$confA[$i],$fasta_seq_fileB,$confB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1312 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1313
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1314 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1315 print HTML "</pre>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1316 print HTML "<hr WIDTH=\"100%\">";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1317 print HTML "<h3>";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1318 print HTML "Group of orthologs #" . $i .". Best score $ortoS[$i] bits<br>\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1319 print HTML "Score difference with first non-orthologous sequence - ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1320 printf (HTML "%s:%d %s:%d</h3><pre>\n", $fasta_seq_fileA,$confA[$i],$fasta_seq_fileB,$confB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1321 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1322 ########### Sort and print members of A ############
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1323 $nA = @membersA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1324 $nB = @membersB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1325 $nMAX = ($nA > $nB) ? $nA : $nB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1326 # Sort membersA inside the cluster by confidence:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1327 for $m (0..($nA-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1328 while($confPA[$membersA[$m]] < $confPA[$membersA[$m+1]]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1329 $temp = $membersA[$m];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1330 $membersA[$m] = $membersA[$m+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1331 $membersA[$m+1] = $temp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1332 --$m if ($m > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1333 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1334 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1335 $paralogsA[$i] = join(' ',@membersA); # Put them back together
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1336 # Sort membersB inside the cluster by confidence:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1337 for $m (0..($nB-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1338 while($confPB[$membersB[$m]] < $confPB[$membersB[$m+1]]){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1339 $temp = $membersB[$m];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1340 $membersB[$m] = $membersB[$m+1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1341 $membersB[$m+1] = $temp;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1342 --$m if ($m > 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1343 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1344 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1345 $paralogsB[$i] = join(' ',@membersB); # Put them back together
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1346 # Print to text file and to HTML file
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1347 for $m (0..($nMAX-1)){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1348 if ($m < $nA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1349 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1350 printf (OUTPUT "%-20s\t%.2f%%\t\t", $nameA[$membersA[$m]], (100*$confPA[$membersA[$m]]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1351 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1352 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1353 print HTML "<B>" if ($confPA[$membersA[$m]] == 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1354 printf (HTML "%-20s\t%.2f%%\t\t", $nameA[$membersA[$m]], (100*$confPA[$membersA[$m]]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1355 print HTML "</B>" if ($confPA[$membersA[$m]] == 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1356 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1357 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1358 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1359 printf (OUTPUT "%-20s\t%-7s\t\t", " ", " ");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1360 printf (HTML "%-20s\t%-7s\t\t", " ", " ") if ($html);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1361 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1362 if ($m < $nB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1363 if ($output){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1364 printf (OUTPUT "%-20s\t%.2f%%\n", $nameB[$membersB[$m]], (100*$confPB[$membersB[$m]]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1365 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1366 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1367 print HTML "<B>" if ($confPB[$membersB[$m]] == 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1368 printf (HTML "%-20s\t%.2f%%", $nameB[$membersB[$m]], (100*$confPB[$membersB[$m]]));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1369 print HTML "</B>" if ($confPB[$membersB[$m]] == 1);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1370 print HTML "\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1371 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1372 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1373 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1374 printf (OUTPUT "%-20s\t%-7s\n", " ", " ") if($output);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1375 print HTML "\n" if ($html);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1376 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1377 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1378 print OUTPUT $message if ($use_bootstrap and $output);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1379 print HTML "$htmlmessage" if ($use_bootstrap and $html);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1380 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1381 if ($output) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1382 close OUTPUT;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1383 print "Output saved to file $outputfile\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1384 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1385 if ($html){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1386 close HTML;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1387 print "HTML output saved to $htmlfile\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1388 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1389
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1390 if ($table){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1391 $filename = "table." . $ARGV[0] . "-" . $ARGV[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1392 open F, ">$filename" or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1393 print F "OrtoID\tScore\tOrtoA\tOrtoB\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1394 for $i(1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1395 print F "$i\t$ortoS[$i]\t";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1396 @members = split(/ /, $paralogsA[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1397 for $m (@members){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1398 $m =~ s/://g;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1399 printf (F "%s %.3f ", $nameA[$m], $confPA[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1400 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1401 print F "\t";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1402 @members = split(/ /, $paralogsB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1403 for $m (@members){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1404 $m =~ s/://g;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1405 printf (F "%s %.3f ", $nameB[$m], $confPB[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1406 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1407 print F "\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1408 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1409 close F;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1410 print "Table output saved to $filename\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1411 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1412 if ($mysql_table){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1413 $filename2 = "sqltable." . $ARGV[0] . "-" . $ARGV[1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1414 open F2, ">$filename2" or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1415 for $i(1..$o){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1416 @membersA = split(/ /, $paralogsA[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1417 for $m (@membersA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1418 # $m =~ s/://g;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1419 if ($use_bootstrap && $bsA[$m]) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1420 printf (F2 "%d\t%d\t%s\t%.3f\t%s\t%d%\n", $i, $ortoS[$i], $ARGV[0], $confPA[$m], $nameA[$m], 100*$bsA[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1421 } else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1422 printf (F2 "%d\t%d\t%s\t%.3f\t%s\n", $i, $ortoS[$i], $ARGV[0], $confPA[$m], $nameA[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1423 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1424 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1425 @membersB = split(/ /, $paralogsB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1426 for $m (@membersB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1427 # $m =~ s/://g;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1428 if ($use_bootstrap && $bsB[$m]) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1429 printf (F2 "%d\t%d\t%s\t%.3f\t%s\t%d%\n", $i, $ortoS[$i], $ARGV[1], $confPB[$m], $nameB[$m], 100*$bsB[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1430 }else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1431 printf (F2 "%d\t%d\t%s\t%.3f\t%s\n", $i, $ortoS[$i], $ARGV[1], $confPB[$m], $nameB[$m]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1432 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1433 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1434 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1435 close F2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1436 print "mysql output saved to $filename2\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1437 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1438 if ($show_times){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1439 ($user_time,,,) = times;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1440 printf ("Finding bootstrap values and printing took %.2f seconds\n", ($user_time - $prev_time));
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1441 printf ("The overall execution time: %.2f seconds\n", $user_time);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1442 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1443 if ($run_blast) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1444 unlink "formatdb.log";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1445 unlink "$fasta_seq_fileA.phr", "$fasta_seq_fileA.pin", "$fasta_seq_fileA.psq";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1446 unlink "$fasta_seq_fileB.phr", "$fasta_seq_fileB.pin", "$fasta_seq_fileB.psq" if (@ARGV >= 2);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1447 unlink "$fasta_seq_fileC.phr", "$fasta_seq_fileC.pin", "$fasta_seq_fileC.psq" if ($use_outgroup);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1448 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1449 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1450
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1451 ##############################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1452 # Functions:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1453 ##############################################################
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1454 sub clean_up { # Sort members within cluster and clusters by size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1455 ############################################################################################### Modification by Isabella 3
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1456
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1457 # Sort on index arrays with perl's built in sort instead of using bubble sort.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1458
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1459 $var = shift;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1460 $totalA = $totalB = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1461 # First pass: count members within each cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1462 foreach $i (1..$o) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1463 @membersA = split(/ /, $paralogsA[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1464 $clusnA[$i] = @membersA; # Number of members in this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1465 $totalA += $clusnA[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1466 $paralogsA[$i] = join(' ',@membersA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1467
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1468 @membersB = split(/ /, $paralogsB[$i]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1469 $clusnB[$i] = @membersB; # Number of members in this cluster
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1470 $totalB += $clusnB[$i];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1471 $paralogsB[$i] = join(' ',@membersB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1472
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1473 $clusn[$i] = $clusnB[$i] + $clusnA[$i]; # Number of members in given group
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1474 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1475
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1476 # Create an array used to store the position each element shall have in the final array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1477 # The elements are initialized with the position numbers
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1478 my @position_index_array = (1..$o);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1479
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1480 # Sort the position list according to cluster size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1481 my @cluster_sorted_position_list = sort { $clusn[$b] <=> $clusn[$a]} @position_index_array;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1482
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1483 # Create new arrays for the sorted information
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1484 my @new_paralogsA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1485 my @new_paralogsB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1486 my @new_is_paralogA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1487 my @new_is_paralogB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1488 my @new_clusn;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1489 my @new_ortoS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1490 my @new_ortoA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1491 my @new_ortoB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1492
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1493
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1494 # Add the information to the new arrays in the orer specifeid by the index array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1495 for (my $index_in_list = 0; $index_in_list < scalar @cluster_sorted_position_list; $index_in_list++) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1496
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1497 my $old_index = $cluster_sorted_position_list[$index_in_list];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1498
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1499 if (!$clusn[$old_index]) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1500 $o = (scalar @new_ortoS) - 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1501 last;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1502 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1503
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1504 $new_paralogsA[$index_in_list + 1] = $paralogsA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1505 $new_paralogsB[$index_in_list + 1] = $paralogsB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1506 $new_is_paralogA[$index_in_list + 1] = $is_paralogA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1507 $new_is_paralogB[$index_in_list + 1] = $is_paralogB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1508 $new_clusn[$index_in_list + 1] = $clusn[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1509 $new_ortoA[$index_in_list + 1] = $ortoA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1510 $new_ortoB[$index_in_list + 1] = $ortoB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1511 $new_ortoS[$index_in_list + 1] = $ortoS[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1512 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1513
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1514 @paralogsA = @new_paralogsA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1515 @paralogsB = @new_paralogsB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1516 @is_paralogA = @new_is_paralogA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1517 @is_paralogB = @new_is_paralogB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1518 @clusn = @new_clusn;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1519 @ortoS = @new_ortoS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1520 @ortoA = @new_ortoA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1521 @ortoB = @new_ortoB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1522
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1523 # Create an array used to store the position each element shall have in the final array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1524 # The elements are initialized with the position numbers
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1525 @position_index_array = (1..$o);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1526
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1527 # Sort the position list according to score
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1528 @score_sorted_position_list = sort { $ortoS[$b] <=> $ortoS[$a] } @position_index_array;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1529
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1530 # Create new arrays for the sorted information
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1531 my @new_paralogsA2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1532 my @new_paralogsB2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1533 my @new_is_paralogA2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1534 my @new_is_paralogB2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1535 my @new_clusn2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1536 my @new_ortoS2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1537 my @new_ortoA2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1538 my @new_ortoB2 = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1539
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1540 # Add the information to the new arrays in the orer specifeid by the index array
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1541 for (my $index_in_list = 0; $index_in_list < scalar @score_sorted_position_list; $index_in_list++) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1542
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1543 my $old_index = $score_sorted_position_list[$index_in_list];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1544 $new_paralogsA2[$index_in_list + 1] = $paralogsA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1545 $new_paralogsB2[$index_in_list + 1] = $paralogsB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1546 $new_is_paralogA2[$index_in_list + 1] = $is_paralogA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1547 $new_is_paralogB2[$index_in_list + 1] = $is_paralogB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1548 $new_clusn2[$index_in_list + 1] = $clusn[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1549 $new_ortoA2[$index_in_list + 1] = $ortoA[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1550 $new_ortoB2[$index_in_list + 1] = $ortoB[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1551 $new_ortoS2[$index_in_list + 1] = $ortoS[$old_index];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1552 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1553
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1554 @paralogsA = @new_paralogsA2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1555 @paralogsB = @new_paralogsB2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1556 @is_paralogA = @new_is_paralogA2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1557 @is_paralogB = @new_is_paralogB2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1558 @clusn = @new_clusn2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1559 @ortoS = @new_ortoS2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1560 @ortoA = @new_ortoA2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1561 @ortoB = @new_ortoB2;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1562
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1563 #################################################################################### End modification by Isabella 3
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1564
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1565 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1566 sub bootstrap{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1567 my $species = shift;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1568 my $seq_id1 = shift; # Query ID from $species
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1569 my $seq_id2 = shift; # Best hit ID from other species
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1570 my $seq_id3 = shift; # Second best hit
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1571 # Retrieve sequence 1 from $species and sequence 2 from opposite species
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1572 my $significance = 0.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1573
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1574 if ($species eq $fasta_seq_fileA){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1575 $file1 = $fasta_seq_fileA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1576 $file2 = $fasta_seq_fileB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1577 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1578 elsif ($species eq $fasta_seq_fileB){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1579 $file1 = $fasta_seq_fileB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1580 $file2 = $fasta_seq_fileA;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1581 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1582 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1583 print "Bootstrap values for ortholog groups are not calculated\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1584 return 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1585 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1586
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1587 open A, $file1 or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1588 $id = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1589 $print_this_seq = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1590 $seq1 = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1591 $seq2 = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1592
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1593 $query_file = $seq_id1 . ".faq";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1594 open Q, ">$query_file" or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1595
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1596 while (<A>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1597 if(/^\>/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1598 ++$id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1599 $print_this_seq = ($id == $seq_id1)?1:0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1600 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1601 print Q if ($print_this_seq);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1602 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1603 close A;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1604 close Q;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1605 ###
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1606 open B, $file2 or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1607 $db_file = $seq_id2 . ".fas";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1608 open DB, ">$db_file" or die;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1609 $id = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1610 $print_this_seq = 0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1611
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1612 while (<B>){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1613 if(/^\>/){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1614 ++$id;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1615 $print_this_seq = (($id == $seq_id2) or ($id == $seq_id3))?1:0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1616 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1617 print DB if ($print_this_seq);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1618 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1619 close B;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1620 close DB;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1621
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1622 system "$formatdb -i $db_file";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1623 # Use soft masking in 1-pass mode for simplicity.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1624 system "$blastall -F\"m S\" -i $query_file -z 5000000 -d $db_file -p blastp -M $matrix -m7 | ./$blastParser 0 -a > $seq_id2.faa";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1625
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1626 # Note: Changed score cutoff 50 to 0 for blast2faa.pl (060402).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1627 # Reason: after a cluster merger a score can be less than the cutoff (50)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1628 # which will remove the sequence in blast2faa.pl. The bootstrapping will
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1629 # then fail.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1630 # AGAIN, updaye
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1631
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1632 if (-s("$seq_id2.faa")){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1633
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1634 system("java -jar $seqstat -m $matrix -n 1000 -i $seq_id2.faa > $seq_id2.bs"); # Can handle U, u
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1635
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1636 if (-s("$seq_id2.bs")){
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1637 open BS, "$seq_id2.bs" or die "pac failed\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1638 $_ = <BS>;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1639 ($dummy1,$dummy2,$dummy3,$dummy4,$significance) = split(/\s+/);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1640 close BS;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1641 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1642 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1643 print STDERR "pac failed\n"; # if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1644 $significance = -0.01;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1645 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1646 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1647 else{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1648 print STDERR "blast2faa for $query_file / $db_file failed\n"; # if ($debug);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1649 $significance = 0.0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1650 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1651
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1652 unlink "$seq_id2.fas", "$seq_id2.faa", "$seq_id2.bs", "$seq_id1.faq";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1653 unlink "formatdb.log", "$seq_id2.fas.psq", "$seq_id2.fas.pin", "$seq_id2.fas.phr";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1654
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1655 return $significance;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1656 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1657
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1658 sub overlap_test{
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1659 my @Fld = @_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1660
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1661 # Filter out fragmentary hits by:
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1662 # Ignore hit if aggregate matching area covers less than $seq_overlap_cutoff of sequence.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1663 # Ignore hit if local matching segments cover less than $segment_coverage_cutoff of sequence.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1664 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1665 # $Fld[3] and $Fld[4] are query and subject lengths.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1666 # $Fld[5] and $Fld[6] are lengths of the aggregate matching region on query and subject. (From start of first matching segment to end of last matching segment).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1667 # $Fld[7] and $Fld[8] are local matching length on query and subject (Sum of all segments length's on query).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1668
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1669 $retval = 1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1670 # if ($Fld[3] >= $Fld[4]) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1671 if ($Fld[5] < ($seq_overlap_cutoff * $Fld[3])) {$retval = 0};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1672 if ($Fld[7] < ($segment_coverage_cutoff * $Fld[3])) {$retval = 0};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1673 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1674
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1675 # if ($Fld[4] >= $Fld[3]) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1676 if ($Fld[6] < ($seq_overlap_cutoff * $Fld[4])) {$retval = 0};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1677 if ($Fld[8] < ($segment_coverage_cutoff * $Fld[4])) {$retval = 0};
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1678 # }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1679
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1680 # print "$Fld[3] $Fld[5] $Fld[7]; $Fld[4] $Fld[6] $Fld[8]; retval=$retval\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1681
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1682 return $retval;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1683 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1684
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1685 sub do_blast
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1686 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1687 my @parameter=@_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1688 my $resultfile=@parameter[@parameter-1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1689 my $go_to_blast=1;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1690 my $resultfilesize;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1691 if (-e $resultfile)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1692 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1693 $resultfilesize= -s "$resultfile";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1694 if ($resultfilesize >10240)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1695 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1696 $go_to_blast=0;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1697 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1698 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1699 if ($go_to_blast)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1700 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1701 if ($blast_two_passes)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1702 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1703 do_blast_2pass(@parameter);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1704 } else
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1705 {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1706 do_blast_1pass(@parameter);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1707 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1708 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1709 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1710
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1711 sub do_blast_1pass {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1712 my @Fld = @_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1713
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1714 # $Fld [0] is query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1715 # $Fld [1] is database
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1716 # $Fld [2] is query size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1717 # $Fld [3] is database size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1718 # $Fld [4] is output name
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1719
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1720 # Use soft masking (low complexity masking by SEG in search phase, not in alignment phase).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1721 system ("$blastall -F\"m S\" -i $Fld[0] -d $Fld[1] -p blastp -v $Fld[3] -b $Fld[3] -M $matrix -z 5000000 -m7 | ./$blastParser $score_cutoff > $Fld[4]");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1722 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1723
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1724 sub do_blast_2pass {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1725
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1726 my @Fld = @_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1727
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1728 # $Fld [0] is query
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1729 # $Fld [1] is database
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1730 # $Fld [2] is query size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1731 # $Fld [3] is database size
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1732 # $Fld [4] is output name
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1733
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1734 # assume the script has already formatted the database
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1735 # we will now do 2-pass approach
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1736
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1737 # load sequences
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1738
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1739 %sequencesA = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1740 %sequencesB = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1741
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1742 open (FHA, $Fld [0]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1743 while (<FHA>) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1744
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1745 $aLine = $_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1746 chomp ($aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1747
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1748 $seq = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1749
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1750 if ($aLine =~ />/) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1751 @words = split (/\s/, $aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1752 $seqID = $words[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1753 $sequencesA {$seqID} = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1754 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1755 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1756 $sequencesA {$seqID} = $sequencesA {$seqID}.$aLine;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1757 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1758 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1759 close (FHA);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1760
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1761 open (FHB, $Fld [1]);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1762 while (<FHB>) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1763 $aLine = $_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1764 chomp ($aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1765
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1766 $seq = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1767
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1768 if ($aLine =~ />/) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1769 @words = split (/\s/, $aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1770 $seqID = $words[0];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1771 $sequencesB {$seqID} = "";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1772 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1773 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1774 $sequencesB {$seqID} = $sequencesB {$seqID}.$aLine;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1775 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1776 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1777 close (FHB);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1778
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1779 # Do first pass with compositional adjustment on and soft masking.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1780 # This efficiently removes low complexity matches but truncates alignments,
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1781 # making a second pass necessary.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1782 print STDERR "\nStarting first BLAST pass for $Fld[0] - $Fld[1] on ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1783 system("date");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1784 open FHR, "$blastall -C3 -F\"m S\" -i $Fld[0] -d $Fld[1] -p blastp -v $Fld[3] -b $Fld[3] -M $matrix -z 5000000 -m7 | ./$blastParser $score_cutoff|";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1785
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1786 %theHits = ();
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1787 while (<FHR>) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1788 $aLine = $_;
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1789 chomp ($aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1790 @words = split (/\s+/, $aLine);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1791
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1792 if (exists ($theHits {$words [0]})) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1793 $theHits {$words [0]} = $theHits {$words [0]}." ".$words [1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1794 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1795 else {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1796 $theHits {$words [0]} = $words [1];
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1797 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1798
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1799 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1800 close (FHR);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1801
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1802 $tmpdir = "."; # May be slightly (5%) faster using the RAM disk "/dev/shm".
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1803 $tmpi = "$tmpdir/tmpi";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1804 $tmpd = "$tmpdir/tmpd";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1805
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1806 # Do second pass with compositional adjustment off to get full-length alignments.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1807 print STDERR "\nStarting second BLAST pass for $Fld[0] - $Fld[1] on ";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1808 system("date");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1809 unlink "$Fld[4]";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1810 foreach $aQuery (keys % theHits) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1811
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1812 # Create single-query file
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1813 open (FHT, ">$tmpi");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1814 print FHT ">$aQuery\n".$sequencesA {">$aQuery"}."\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1815 close (FHT);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1816
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1817 # Create mini-database of hit sequences
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1818 open (FHT, ">$tmpd");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1819 foreach $aHit (split (/\s/, $theHits {$aQuery})) {
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1820 print FHT ">$aHit\n".$sequencesB {">$aHit"}."\n";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1821 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1822 close (FHT);
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1823
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1824 # Run Blast and add to output
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1825 system ("$formatdb -i $tmpd");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1826 system ("$blastall -C0 -FF -i $tmpi -d $tmpd -p blastp -v $Fld[3] -b $Fld[3] -M $matrix -z 5000000 -m7 | ./$blastParser $score_cutoff >> $Fld[4]");
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1827 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1828
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1829 unlink "$tmpi", "$tmpd", "formatdb.log", "$tmpd.phr", "$tmpd.pin", "$tmpd.psq";
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1830 }
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1831
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1832
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1833 # Date Modification
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1834 # -------- ---------------------------------------------------
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1835 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1836 # 2006-04-02 [1.36] - Changed score cutoff 50 to 0 for blast2faa.pl.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1837 # Reason: after a cluster merger a score can be less than the cutoff (50)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1838 # which will remove the sequence in blast2faa.pl. The bootstrapping will
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1839 # then fail.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1840 # - Fixed bug with index variable in bootstrap routine.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1841 #
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1842 # 2006-06-01 [2.0] - Fixed bug in blast_parser.pl: fields 7 and 8 were swapped,
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1843 # it was supposed to print match_area before HSP_length.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1844 # - Fixed bug in blastall call: -v param was wrong for the A-B
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1845 # and B-A comparisons.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1846 # -
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1847 # - Changed "cluster" to "group" consistently in output.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1848 # - Changed "main ortholog" to "seed ortholog" in output.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1849 # - Replace U -> X before running seqstat.jar, otherwise it crashes.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1850 # 2006-08-04 [2.0] - In bootstrap subroutine, replace U with X, otherwise seqstat
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1851 # will crash as this is not in the matrix (should fix this in seqstat)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1852 # 2006-08-04 [2.1] - Changed to writing default output to file.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1853 # - Added options to run blast only.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1854 # - Fixed some file closing bugs.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1855 # 2007-12-14 [3.0] - Sped up sorting routines (by Isabella).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1856 # - New XML-based blast_parser.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1857 # - New seqstat.jar to handle u and U.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1858 # - Modified overlap criterion for rejecting matches. Now it agrees with the paper.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1859 # 2009-04-01 [4.0] - Further modification of overlap criteria (require that they are met for both query and subject).
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1860 # - Changed bit score cutoff to 40, which is suitable for compositionally adjusted BLAST.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1861 # - Added in 2-pass algorithm.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1862 # 2009-06-11 [4.0] - Moved blasting out to subroutine.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1863 # - Changed blasting in bootstrap subroutine to use unconditional score matrix adjustment and SEG hard masking,
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1864 # to be the same as first step of 2-pass blast.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1865 # 2009-06-17 [4.0] - Compensated a Blast "bug" that sometimes gives a self-match lower score than a non-identical match.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1866 # This can happen with score matrix adjustment and can lead to missed orthologs.
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1867 # 2009-08-18 [4.0] - Consolidated Blast filtering parameters for 2-pass (-C3 -F\"m S\"; -C0 -FF)
83e62a1aeeeb Uploaded
dereeper
parents:
diff changeset
1868 # 2009-10-09 [4.1] - Fixed bug that caused failure if Fasta header lines had more than one word.