forked from vinuesa/get_phylomarkers
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathadd_labels2tree.pl
executable file
·258 lines (210 loc) · 8.31 KB
/
add_labels2tree.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
#!/usr/bin/perl -w
# fastaheader2treetag_PhyMLConsenseTopols_V02.pl
# by Pablo Vinuesa, CCG-UNAM, Mexico. V01 : April 16th 2006
# V02: March 25th 2007
# V03: June 16th 2007 #Marfil's fasta header after filtering redundant sequences
# V04: March 8th 2008 # can now parse
# Used to edit the taxon labels on newick phylogenies generated by PhyML or CONSENSE
# converting numbers to full taxon labels. The latter are read from the
# fasta headers of the original *.faa files, or from the my_*.faa files
# in the formats of GABO, Neurogadgets and IH-PVs
# usage: run from within the directory containing the tree files,
# providing two args: 1.- an *.faa file and 2.- the extension of the topologies to edit phyml_tree.txt
# perl fastaheader2treetag_PhyMLTopols_V01.pl my_rpoA_NCTC11168_aln.faa phyml_tree.txt
use strict;
use File::Basename;
my $progname = basename($0);
my $VERSION = 1.4; # Sept 15th, 2013; added general case (/^>?(\d+)[\t\s+](.*$); again allows >ID [bichococcus ..] type labels to work!
# Jul 10th, 2013; added the general case: (/(\S+)[\t\s+](\S+)/)
# Aug. 4th, 2012;
# 9 Feb 2010 ahora lee también elsif(/^>?(\d+)\s+(.*)/)
# v 1.0 ahora lee if(/^>(\d+)\s+\[(.*?)\]/) y elsif(/^(\d+)\s+\[(.*?)\]/) ... es decir >ID [blah] y ID [blah]
die "\n\n to use $progname v. $VERSION you need to provide two args:
1) the name of an alignment file
2) the file extension name of a Clann supertree or consensus.ph file on which to edit the taxon tags
for example:
$0 rpoC-Rhizobiales_extended_aln.faa rpoC-Rhizobiales_extended_aln_nj.ph
or
$0 aln3.faa aln3_nj.ph\n\n" unless @ARGV ==2;
my($aln_file, $tree_file_ext) = (@ARGV);
my (@tree_files, $tree_file, $species, $strain,$tax_label, $seq_key, %seq_ID, $basename, $outfile, $ext, $gi);
open ALNFILE, $aln_file or die "$0 can't open aln_file $aln_file: $!\n";
print "the seq_key => seq_ID correspondences for aln $aln_file are:\n";
while(<ALNFILE>){
chomp;
if(/^>(\d+)\s+\[(.*?)\]/){
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}
elsif(/^(\d+)\s+\[(.*?)\]/)
{
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}
elsif(/^>?(\S+)[\t\s+](\S+)/)
{
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}
elsif(/^>?(\d+)[\t\s+](.*$)/)
{
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$strain =~ s/[\s\|;-]+/_/g;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}
elsif(/^>?(\d+)\s+\|\d+\|\[(.*)\]/)
{
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$strain =~ s/[\s\|;-]+/_/g; #print "# strain label is $strain\n"; exit;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}
elsif(/^>(\d+) \d+\|\[(.*)\]\|(.*)\|/)
{
#$GI_number = $1;
$seq_key = $1;
$species = $2;
$strain = $3;
$tax_label = $species . '_' . $strain;
$tax_label =~ s/[\s+|;|-]/_/g;
$seq_ID{$seq_key}=$tax_label;
print "$seq_key => $seq_ID{$seq_key}\n"; # format of our extract_CDSfromGBK_batch.pl
}elsif(/^>(\d+)\s+\S+\s+\|\[(.*?)\]/){ # Fasta header for get_homologues.pl
$seq_key=$1;
$strain = $2;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif(/^>(\d+) \d+_\[(\S+)\]/){
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif(/^>\d+_(\S+)/){
#$GI_number = $1;
$seq_key++;
$strain = $1;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif(/^>(\w+)\|/){
#$GI_number = $1;
$seq_key++;
$strain = $1;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif(/^>\d+\s+\S+\s+\[(\S+)\]/){ #GABO's fasta header
$seq_key++;
$strain = $1;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif(/^>(\d+) \d+ \S* (\d+)\|\[(\w+\s*\w*\s*\w*)\]/){ #Marfil's fasta header after filtering redundant sequences
$seq_key = $1; #print "saw seq_key $seq_key\n";
$gi = $2;
$strain = $3; #print "saw strain $strain\n";
$strain =~ s/\s+/_/g;
$strain = $strain . '_' . $gi;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}elsif (/^>(\w+) (.*)/){ # RDPII fasta header; most general fasta header ID w+ w+ w+ ...
#$GI_number = $1;
$seq_key = $1;
$strain = $2;
$strain =~ s/\s+/_/g;
$strain =~ s/,/_/g;
$strain =~ s/:/_/g;
$strain =~ s/".*"/_/g;
$strain =~ s/\(//g;
$strain =~ s/\)//g;
$seq_ID{$seq_key}=$strain;
print "$seq_key => $seq_ID{$seq_key}\n";
}else{
next;
}
}
print "\n";
@tree_files = <*$tree_file_ext>;
foreach $tree_file(@tree_files){
($basename, $ext) = (split(/\./, $tree_file))[0,1];
$outfile = $basename . "_ed." . "$ext";
open TREEFILE, $tree_file or die "$0 can't open aln_file $aln_file: $!\n";
while(<TREEFILE>){
foreach $seq_key (sort keys %seq_ID){
open OUT, ">$outfile" or die "$0 can't write to outfile $outfile: $!\n";
if ($_ =~ /\($seq_key:/){
#print "the match line is $_\n";
s/\($seq_key:/\($seq_ID{$seq_key}:/g;
#print "the substituted line is $_\n";
}elsif ($_ =~ /,$seq_key:/){
s/,$seq_key:/,$seq_ID{$seq_key}:/g;
}elsif ($_ =~ /\(($seq_key)(\__\d+):/){ #if we used Fasta2cluster, with tree labels such as 000000001__3
s/\($1$2:/\($seq_ID{$seq_key}$2:/g; #if we used Fasta2cluster, with tree labels such as 000000001__3
}elsif ($_ =~ /,($seq_key)(\__\d+):/){ #if we used Fasta2cluster, with tree labels such as 000000001__3
s/,$1$2:/,$seq_ID{$seq_key}$2:/g; #if we used Fasta2cluster, with tree labels such as 000000001__3
}
elsif ($_ =~ /\(($seq_key,)/){ #if we used Fasta2cluster, with tree labels such as 000000001__3
s/$1/$seq_ID{$seq_key},/g; #if we used Fasta2cluster, with tree labels such as 000000001__3
}
elsif ($_ =~ /(,$seq_key)\)/){
s/$1/,$seq_ID{$seq_key},/g;
}
elsif ($_ =~ /\)(,$seq_key)/) {
s/$1/,$seq_ID{$seq_key},/g;
}
elsif ($_ =~ /(,$seq_key)\)/) {
s/$1/,$seq_ID{$seq_key}/g;
}
}
print OUT;
print "the edited topology for treefile $tree_file is:\n $_\n";
}
close TREEFILE;
close OUT;
}
# User documentation
=head1 NAME
<fastaheader2treetag_PhyML_Consense_Topol_V03.pl> - <Used for editing of the taxon labels on a newick phylogeny generated by PhyML or CONSENSE converting seqID numbers to full taxon labels>
=head1 VERSION
This documentation refers to <fastaheader2treetag_PhyML_Consense_Topol_V03.pl> version 0.0.3
=head1 USAGE
usage: run from within the directory containing the tree files, providing two args:
1.- my_alignment_or_Fasta_File.faa
2.- the extension of the topologies to edit, like phyml_tree.txt or aln.ph
The latter should complement the basename of the aln files
Usage example:
path/to/fastaheader2treetag_PhyML_Consense_batch_V01.pl my_fasta_file.faa phyml_tree.txt
=head1 REQUIRED ARGUMENTS
1) name of alignment file
2) tree extension name, like phyml_tree.txt or aln.ph
=head1 DESCRIPTION
The program will process all trees in a directory, changing the numeric sequence identifiers to full taxon labels using the alignment provided
=head1 DIAGNOSTICS
Just standard diagnostics
=head1 DEPENDENCIES
No external dependiencies.
=head1 BUGS AND LIMITATIONS
There are no known bugs in this module.
Please report problems to <Maintainer name(s)> (<[email protected]>)
Patches are welcome.
=head1 AUTHOR
<Pablo Vinuesa> (<[email protected]>)
=head1 LICENCE AND COPYRIGHT
Copyright (c) <2007> <Pablo Vinuesa> (<[email protected]>). All rights reserved.
This program is free software; you can redistribute it and/or
modify it under the same terms as Perl itself. See L<perlartistic>.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.