forked from JCVenterInstitute/PanACEA
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmake_panacea.pl
More file actions
executable file
·11834 lines (9658 loc) · 424 KB
/
Copy pathmake_panacea.pl
File metadata and controls
executable file
·11834 lines (9658 loc) · 424 KB
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
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env perl
#Copyright (C) 2016-2021 The J. Craig Venter Institute (JCVI). All rights reserved
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#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. See the
#GNU General Public License for more details.
#You should have received a copy of the GNU General Public License
#along with this program. If not, see <http://www.gnu.org/licenses/>
=pod
=head1 NAME
make_panacea.pl
=head1 SYNOPSIS
make_panacea.pl [options] -i PanACEA.InputFile.txt
Options:
-i --input Input PanACEA Flat File. Required.
-o --output Output Directory for all the files. Default is current directory.
-n --name Header name of HTML and SVG files. Default is \"PanHTML\"
-d --dir_name Root directory of the HTML for use in the browser. Default is
current directory.
-f --function File containing information about additional gene cluster annotation.
-a --fasta Either directory containing the cluster multiple alignment files
or a single fasta file. Multiple inputs could be added in a comma separated
list
-t --tree File containing the phylogeny of the genomes in newick format.
(Requires BioPerl to use)
-m --metafiles File containing group metadata for the genomes. Only used in
connection with tree.
-g --graphic File containing the graphics configure file to change the
size of the whole image and parts of the image
=head1 DESCRIPTION
make_panacea.pl generates an interconnected set of html, svg, javascript and json files
from a PanACEA flat file generated from a possible suite of programs, including Pangenome and
Roary, that allows users to explore multiple related prokayrotic genomes through their
pan-chromosomes, including core and flexible regions. While the basic version only shows
the location and relative position of the core and flexible regions, the user can also
add information on the gene cluster annotation and alignment as well as the genomic
annotation and phylogeny to obtain a more complete view of the relationship between the
genomes, their function and any regions with their associated genes that might drive
this relationship. The end result is a central html file that can be opened in many web
browsers.
=head2 Notes
=head3 'function' annotation config file format (-f flag)
PanACEA can use many types of annotations , but requires each to be inputed in modular format
in a single configure file with the format shown below. Each annotation block begins with
START and ends with END, and includes:
* START ID a string used on the table tab
* mapfile: a tab-delineated file with gene or centroid ID in the first column
and a comma-separated term list in the second
* ontology: an obo type onotology file containing definitions for the terms in
the mapfile
* name: a tab-delineated file with gene or centroid ID in the first column
and a function in the second, OR a string that is assigned to each of the genes in the mapfile
* color: a tab-delineated file with the function in the first column and a color
in the second, OR a rgb color string that is assigned to each of the gene
with the functionq
More information can be found in the manual
=head3 'graphic' configure file (-g flag)
The graphical annotation config file allows the user to change the dimensions and colors
of the image. The file consists of unique lines with the variable name followed by a space
and then followed by the value for one of 6 different variables. Variable names not listed
below will be ignored. Repeated variables will be given the last value. A list of the
variables along with a description follows:
* BACKGROUND color of the background of the outer ring/core regions on the main
view
* fGR color of the fGRs in the outer ring/core regions on the main view
* BORDER_SIZE size of the border between the outer ring and the edge of the
circle in pixels
* CIRCLE_SIZE size of the radius of the outer edge of the outer ring in pixels
* TEXT_SIZE size of the text in pixels
* GENE_SIZE size of the gene in pixels
=head1 AUTHOR
Thomas Clarke (tclarke@jcvi.org)
2017 J. Craig Venter Institute
=cut
#Detailed decriptions of the sub-ruetines are found at the end of the script
use strict;
use Getopt::Long;
use Cwd;
use warnings;
use Pod::Usage;
#load information from the files
my $GENOME_NUM = 0;
#array function
my @funcInfo;
#This an image of the disk
#Used on a couple of images
my $disk_image_orig = disk_image("default");
#This is the image of the disk appropiate for javascript
my $disk_image_js = disk_image("js");
my $core_cnt = 0; #Running count of the number of core regions to make the figures
my $CORE_PERCT = 0.75; #The core percentage- used to get a minimal y-axis in the core region graphs
my %jso; #hash used to make the tables to be made into JSON format for javascript to read
my %head; #keeps track of the table headers for each of the tables
my %col_width; # the widths
my %gene_num; #Matches the gene cluster name with the corresponding array number in the JSON matrix
my %seqs; #keeps the sequences associated with genes
my %gene_func; #gene function
my %clusters; #stores names and number of genomes in cluster
my %fgi_grps; #type = circular/linear; bound = the start and stop core genes; parent; sz = size; hit
my %fgi_member; #the fGI that a cluster belongs to
my %cores; #cores and their type
my %list; #used to make fGI gene array: contains bounding genes, and the different fGRs
my $seq_len = 0; #The length of genome
my %core_list; #list of core genes
my %tick; #gets information for ticks around core
my %core_size; #gets the core size for the number
my %fgi_gene_list; #hash of arrays
my %geneLenInfo; # hash containing the gene length mean, min, max and standard deviation for each cluster;
my %jsHead; #json version of the table headers
my %go_name; #names for the go terms (and any other added funnction) and the genes associated with it
my $tree_level = 0; #Level at which trees are initial shown with 1 = leaves. O = tree not loaded
my %meta_types; #Differnt genome metadata variables
my $tree_size = 500; #Size of the tree image square
my %func_list; #list of the functions associated with the genes
my $defaultFunct; #The default funtional value
#### User Inputed variables on Command line###
my $mult_align_dir = ""; #The directory containing the cluster multiple alignments
my $dir = cwd(); #Directory to write the output
my $out_dir = cwd(); #Default value for the output head directory
my $output_id = "PanACEA"; #Default value for -n
my $func_file; #Location of the configure file with the
my $file_in1; #Location of the PanACEA flat file
my $verb = 1; #Verbosity level: not currently used
my $graphic_config_file; #Configure file with graphical output information
my $tree_file; #File with newick tree
my $meta_data_file; #file(s) with group information with the newick trees
#Graphical values
# default values:
my $PI = 3.14159265358979323846;
my $max_radius = 400; #Radius of the chromosome image
my $border = 120; #Size of the border between chromosome and the edge of the screen
my $layer_height = $max_radius / 12; #
my $gene_height = 24; #height of the gene icons in the image
my $tick_height = $gene_height / 3; # tick size on the edge of the image
my $text_height = 6; #Height of the
my $thumb_dif = 0.15; #The proportional size of the thumbnail sketches
my $DPI = 3; #gives the DPI for outputed PNG images
my $SVGHEIGHT = 2 * ( $border + $max_radius ); #Height of the total image
my $SVGWIDTH = 2 * ( $border + $max_radius ); #Width of the total image
my $thumb_size_h = $SVGHEIGHT * $thumb_dif; # size of the thumbnails
my $thumb_size_w = $SVGWIDTH * $thumb_dif; # size of thumbnails
my $chng_w = $gene_height + $thumb_size_w; #
my $disp = "inline-block"; #The
my $bw = ( $max_radius * 2 - 2.25 * $border ); #body width
my $bodyTop = "30px"; #Sets the location of the table body
my $bodyHeight = ( $max_radius / 4 - 30 ) . "px"; #height
my $bodyWidth = "100%"; #width of the table body
my $cr_start = 0; #starting number of chromosomes
my $cr_cnt; #number of chromosomes
#Default name for files generated by the program
my $out_name = "OUT"; #
my $filenameSVG = $out_name . ".svg"; #written svgs
my $filenameTXT = $out_name . ".txt"; #written txt files
my $filenamePNG = $out_name . ".png"; #written png files
my %term2type;
my %type2term;
my $cur_st; #current start value when drawing the pan chromosomes: this allows the horseshoe linear chromosomes to be drawn
my $default_table = "ButtonIDRegion"; #The starting table selection Javascript ID
my $default_table_type = "Region"; #The starting table selection
my %ont; #Connect term IDs to their OBO definitions
my %rgb; #Connects terms IDS to their color
my $help;
#Getting all the user inputed usage
GetOptions(
"i|input:s" => \$file_in1,
"t|tree:s" => \$tree_file,
"m|metadata:s" => \$meta_data_file,
"a|fasta:s" => \$mult_align_dir,
"o|output:s" => \$out_dir,
"d|dir:s" => \$dir,
"f|function:s" => \$func_file,
"n|name:s" => \$output_id,
"g|graphic:s" => \$graphic_config_file,
"help|h|?" => sub { pod2usage( { -verbose => 2, -exitval => 1 } ) }
) || die("Cannot recognize usage. Please use --help to see options");
if ( !$file_in1 ) {
pod2usage(2);
exit();
}
#Initializes & runs the
main();
#This is the main sub reutine that performs all the task
sub main {
if ($verb) {
warn "Welcome to PanNav Creation Software\n\n";
warn "The PanGenome File located in $file_in1 will be used..\n\n\n";
}
#Makes the output directories
make_output_dirs();
# This is the hash of the header of the tables
#The default gene and region tabs
$head{Region} = [ "ID", "Type", "# of Genes", "Associated Terms" ];
$head{Gene} = [ "ID", "Name", "# of Genomes", "Region", "Associated Terms" ];
#Setting the column widths of the file
$col_width{Region} = [ .20, .20, .10, .50 ];
$col_width{Gene} = [ .10, .20, .10, .20, .40 ];
foreach my $a ( keys(%head) ) {
for ( my $i = 0 ; $i < scalar( @{ $head{$a} } ) ; $i++ ) {
$jsHead{$a}->{ $head{$a}->[$i] } = 1;
}
}
#Checks to see if a tree is included, and if so loads it in
#If the tree is loaded in, then tree_level is set to the defaut value (4): when it is not zero this means
#that a tree is included
if ($tree_file) {
warn "Using $tree_file to incorporate tree in visualization\n\n";
make_tree_json($tree_file);
}
#This is the default function if and only if there is a function file to will assign the rest of the
$defaultFunct = "Hypothetical";
if ( !$func_file ) { $defaultFunct = "NA"; } #This is the case where no functions have been assigned. All will be NA
$func_list{$defaultFunct} = 1; #This is the default assigment
#file list;
if ( -e $file_in1 ) {
if ($verb) { warn "Reading in PanGenome information...\n"; }
get_gene_info_new();
if ($verb) { warn "Finished\n\n\n"; }
} else {
die("Cannot find $file_in1 for reading in files. Please double check...\n\n");
}
$rgb{$defaultFunct} = "dddddd"; #Default Function color is grey
$rgb{"Background"} = "e0e0e0"; #Background color is the
$rgb{"fGR"} = "606060"; #Default fGR color is red: should I change this? Yes to dark gray
#Reading in functional information and annotation
if ($func_file) {
if ($verb) {
warn "Reading in gene function information...\n";
}
read_function_file();
#
foreach my $funcStr (@funcInfo) {
load_functions($funcStr);
}
if ($verb) {
warn "Finished ...\n\n\n";
}
}
#If the user wants to change some of the graphic size and color values
if ($graphic_config_file) {
read_config_file();
}
if ($verb) {
warn "Starting the main image..\n";
}
#Write the javascript files used by the scripts: main (used by home page), fgi (fgi/core pages), and gene page
fgi_javascript();
make_main_javascript();
write_geneFile_javascript();
#making the main figure: the sub-script also uses this secondary and tertiary pages
make_main_figure();
if ($verb) {
warn "Finished all\n";
}
}
#All the files used in the image
sub make_main_figure() {
#This is the main html file. Already the directory has been created
open( SVG_MAIN, ">", $out_dir . "/$output_id/main.html" ) or die();
#%out is the main variable to store all the SVG information for the main pages images with the chromosome number
#as the hash index.
#Since SVGs doen't have z scores, we use staggered hash format for each chromosome with 100 being the highest, ie first drawn
#
#Number 100 is the maximum level allowed
#Number 100 is the back ground
#Number 2 is for the preview windows. It is also in a seperate hash key "preview"
#Number 1 is for the central ring + region
#Only levels 100 and 1 are used in the thumbnail
#Other keys include:
#key "load" shows the loading screen (currenly rotating exapanding and diminishing JCVI square logo)
#key "script" loads the javascript
#There are two scripts: first which loads the default variables and second which calls the script
#Key "beg" is the opening for the HTML: this includes the javascript, legend, save
#Key "end" is the closing for the HTML
#Key "load" is the animation for the loading screen
#Key "script" is the script
#Key "table" is the table
#Key "thumb" contains the thumbnail sketches
my %out;
#The ratio of the thumbnail to the main images
my $assem_core_num = scalar( keys(%cores) );
my %thumb_list;
#Making the html images. Ordered similarly as %out -->
my %html;
$html{beg} = "<!DOCTYPE html><html lang=\"en\">";
#getting the javascript functions
$html{beg} .= "<head> <meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\"/>";
#Initializeing the image by drawing
%out = &svg_init_image( 2 * ( $border + $max_radius ), 2 * ( $border + $max_radius ), \%out );
#
$out{beg} = sprintf(
"<div id=\"mainDiv\" style=\"visibility: hidden; left:%fpx; top:0px; position:absolute;\" onmousemove=\"mouseMove()\">",
1 * $gene_height )
. $out{beg};
#This is the loading screen described above. The blocks are the squares and twirl is the central screen
$out{load} = make_loading_screen();
#Setting the loading of the initial javascript variables
my $background = $rgb{"Background"}; #setting the background color for the arc
my $core_num = 1; #starting the number of the core genes
#This is making the save disk icons for SVG and PNGs
$out{beg} .= sprintf(
"<path d=\"M%f %f L%f %f A%f %f 0 0 1 %f %f Z\" fill=\"lightgray\"/>\n",
( $border + $max_radius ),
( $border + $max_radius ),
$border * 2,
( $border + $max_radius ),
( $max_radius - $border ),
( $max_radius - $border ),
2 * $max_radius,
$border + $max_radius
);
#Making the Disk image by: (1) shrinking it with $trans; (2) adding the action; and (3) making a new variable containing the
#This gives the location(translate) and appropriate size for the main page save SVG image
my $trans = sprintf(
"translate(%f,%f) scale(%f, %f)",
( $border + $max_radius - 20 ),
( $border + $max_radius - $border * 1.90 ),
0.045, 0.045
);
my $action = "onclick=\"saveSVG(\'svg\')\"";
my $disk_image = $disk_image_orig;
$disk_image =~ s/TRANS/$trans/g;
$disk_image =~ s/ACTION/$action/g;
$disk_image =~ s/TEXT/SVG/g;
$disk_image =~ s/FS/250/g;
#Adding disk image to the SVG
$out{beg} .= $disk_image;
#Same as above but for PNG
$trans = sprintf(
"translate(%f,%f) scale(%f, %f)",
( $border + $max_radius + 20 ),
( $border + $max_radius - $border * 1.90 ),
0.045, 0.045
);
$action = "onclick=\"saveSVG(\'png\')\"";
$disk_image = $disk_image_orig;
$disk_image =~ s/TRANS/$trans/g;
$disk_image =~ s/ACTION/$action/g;
$disk_image =~ s/TEXT/PNG/g;
$disk_image =~ s/FS/250/g;
$out{beg} .= $disk_image;
#Legend Title
$out{beg} .= "<g id=\"Legend1\" visiblity=\"visible\">\n";
$out{beg} .= sprintf(
"<text id=\"LegendText\" x=\"%f\" y=\"%f\" fill=\"black\" font-size=\"40\" text-anchor=\"middle\" numChrom=\"$core_cnt\" onclick=\"makeThumbnails(evt)\">%s</text>\n",
( $border + $max_radius ),
( $border + $max_radius - $border * 1.95 ), "Legend"
);
#Legend: The different functions
my @curr_y; #Keeps track of current y position of the function in the table
my @curr_x; #Keeps track of current x position of the function in the table
my $num_leg_cols = 3; #Setting the number of columns in the legend
#Making the starting x and y position for each column
my $x1 = ( ( $max_radius - $border ) * sin( deg2rad(310) ) ) + ( $border + $max_radius );
my $y1 = ( ( $max_radius - $border ) * cos( deg2rad(310) ) * -1 ) + ( $border + $max_radius );
for ( my $i = 0 ; $i < $num_leg_cols ; $i++ ) {
$curr_x[$i] = $x1 + ( $i * $border * 1.15 );
$curr_y[$i] = $y1;
}
my $i = 0;
if ($verb) {
warn "Found ", scalar( keys(%func_list) ), " different functions\n\n";
}
#func_list: this cycles through by column, and once it
foreach my $a ( keys(%func_list) ) {
$x1 = $curr_x[$i];
$y1 = $curr_y[$i];
#in case of long names, this function (1) switchs words for equivilant symbols (ie and -> &) and making multi-line
my $id_name = split_id( $a, $x1 + $gene_height * 0.6, 22 );
#Getting the number of lines
my $n = -1;
while ( $id_name =~ /\<\/tspan\>/g ) { $n++; }
#Adding the name, the color rectangle. All should be clickable to turn it on/off
$out{beg} .= sprintf( "<g onclick=\"runType(\'GeneType\',evt,\'\')\" href=\"%s\" fillcol=\"%s\">", $a, "#" . $rgb{$a} );
#The rectangle around: thought that I might add a box if it's clicked. Currently turrned off (instead all others fade when one is clicked)
$out{beg} .= sprintf(
"<rect href=\"%s\" id=\"%s\" class=\"LegRect\" x = \"%f\" y=\"%f\" width=\"%f\" height=\"%f\" fill=\"none\" visiblility=\"hidden\" fillcol=\"%s\"/>",
$a, $a . "rect",
$x1, $y1, $gene_height * 0.75 + $border * 1.15, $gene_height, "#" . $rgb{$a}
);
#The colored rectange before
$out{beg} .= sprintf(
"<rect href=\"%s\" id=\"%s\" class=\"LegCol\" x=\"%f\" y=\"%f\" width=\"%f\" height=\"%f\" fill=\"%s\" fillcol=\"%s\"/>\n",
$a, $a . "col", $x1,
$y1 + $gene_height * ($n) / 5,
$gene_height * 0.5,
$gene_height / 3.0,
"#" . $rgb{$a}, "#" . $rgb{$a}
);
#The text printed out. Can be multiple lines thanks to col_names
$out{beg} .= sprintf(
"<text href=\"%s\" id=\"%s\" class=\"LegType\" x=\"%f\" y=\"%f\" width=\"%f\" height=\"auto\" font-size=\"10\" dominant-baseline=\"hanging\" text-anchor=\"start\" fill=\"black\" fillcol=\"%s\>%s</text>\n",
$a,
$a . "Text",
$x1 + $gene_height * 0.6,
$y1, $border * 1.15, "#" . $rgb{$a}, $id_name
);
$out{beg} .= sprintf("</g>");
$y1 += $gene_height * ( ( $n + 1 ) / 2.5 + 0.15 ); #cycling through colunms
$curr_y[$i] = $y1;
$i++;
if ( $i == 3 ) { $i = 0; }
}
$out{beg} .= "</g>";
####################################################################################
# Cycling through the pan-chromosome by:
# 1. Looking at each core gene and:
# 1.1 If it is a core gene,
# 1.1.1 add information to cor_reg hash and goto 1
# 1.2 If it it is a puedo-core gene (i.e. flexible region)
# 1.2.1 If core genes in the cor_reg hash:
# 1.2.1.1 Make core region preview pane and page from cor_reg
# 1.2.1.2 Draw core region clickable segment on the chromosomes
# 1.2.1.2 Reset cor_reg hash
# 1.2.2 Find all genes in associated with fGRs and build fgi_reg hash
# 1.2.3 Make fGR region preview pane and page from fgi_reg
# 2. If core genes in cor_reg hash
# 2.2.1 Make core region preview pane and page from cor_reg
# 2.2.2 Draw core region clickable segment on the chromosomes
####################################################################################
my %cor_reg; #Hash with a list of core region information passed to the svg draw page
#Keys: st -> starting bp
# end -> ending bp
# coords -> array with gene
my %reg2gene; #Hash that maps the region to a gene
my $core_type = ""; # contains the annotations for each core gene: used for table
my $core_oth = ""; # contains any additional annotation for each core gene; ibid
my $cur_gene_list = ""; #Contains a semi-colon sperated list of gene IDs for each region
my $old_st = 0; #Used to keep track of the position of the core region start
my $even = 0; # used to iterate the position of the core regions to shuffle the verticle postion of the fGRs
for my $cur_chr ( keys(%cores) ) {
$out{script_chr} .= sprintf( "\n\nvar %s = \'{", "fastaJSONPlasmid" . $cur_chr );
if ($verb) {
warn "Starting the image for Core $cur_chr..\n";
}
my $c3 = 0; #Number of chromosomes in the thumbnails
if ( $thumb_list{Assembly_core} ) { $c3 = scalar( @{ $thumb_list{Assembly_core} } ); }
$thumb_list{Assembly_core}->[$c3] = $cur_chr;
if ( $cores{$cur_chr}->{type} eq "cycle" ) {
#Make a complete circle
$cur_st = 0;
#Setting the size of the chromosome. This is to make sure that the ticks/genes are correctly spaced
$seq_len = $cores{$cur_chr}->{sz};
#Drawing the ticks around the chromosome
#Trying to get ticks that have have intervals of a factor of 10 and number between 10 to 100
my $tick_space = 10**( int( log($seq_len) / log(10) ) - 1 );
$out{$cur_chr} = &place_svg_ticks( $tick_space, $seq_len, $out{$cur_chr} );
#Drawing the background of a full circle
$out{$cur_chr} = &svg_draw_arc_full(
100, &convert_layer(1),
&convert_layer(1) + $gene_height,
"#" . $background,
0, 0, $out{$cur_chr}
);
}
if ( $cores{$cur_chr}->{type} eq "chain" ) {
#Change to incomplete circle- 300 degree in a horseshoe shape
$cur_st = $cores{$cur_chr}->{sz} * 7.5 / 300;
#This is to make sure that the ticks/genes are in the appropiate place
#The functions are designed to work on a full
$seq_len = $cores{$cur_chr}->{sz} * 315 / 300;
#svg_draw_arc_seg-> change starting position on arc but not genome.
$out{$cur_chr} = &svg_draw_arc_seg(
100, 0, $cores{$cur_chr}->{sz},
&convert_layer(1),
&convert_layer(1) + $gene_height,
"#" . $background,
"", 0, "", "", $cur_chr, $out{$cur_chr}
);
#Trying to get ticks that have have intervals of a factor of 10 and number between 10 to 100
my $tick_space = 10**( int( log($seq_len) / log(10) ) - 1 );
#A modified verion of the above
$out{$cur_chr} = &place_svg_ticks_arc( $tick_space, $cores{$cur_chr}->{sz}, $out{$cur_chr} );
}
my $terms;
my %term_cnt;
my $out_cnt = 0;
$cr_start = 0;
#Counting through the core genes: can also include psuedo-core genes, that is fGRs
foreach my $core ( @{ $core_list{$cur_chr} } ) {
#Each corde gene includes the chromosome ($chr), the ID, the starting BP, the ending BP, the long gene name, the type (core/fgi), and the length in BP
my ( $chr, $ID, $st, $end, $def, $type, $len ) = @$core;
#Setting the start value for the next group
$old_st = $end + 1;
if ( $type eq "fGI" ) {
#Getting the heights of the Regions on the pan-chromosome
my $h1 = &convert_layer(1);
my $h2 = $h1 + ( $gene_height * 0.75 );
#This gets the fGR staggering variable
if ($even) { $h1 += $gene_height * .25; $h2 += $gene_height * .25; $even = 0; }
else { $even = 1; }
#Checking to see if there is a core-region (ie area with core-genes)
if (%cor_reg) {
#The location of the Core Region Detail Page
my $new = "core/Core_Region.$core_cnt.html";
my $c = 0; #The row of the Region in the Region Table Array
if ( $jso{Region} ) { $c = scalar( @{ $jso{Region} } ); }
#The values to be added to the Region Table
my @tmp = (
"Core_Region.$core_cnt", "Core Region", scalar( @{ $cor_reg{cords} } ),
$core_type, "CORE" . $core_num
);
for ( my $i = 0 ; $i < scalar( @{ $head{Region} } ) ; $i++ ) {
$jso{Region}->[$c]->{ $head{Region}->[$i] } = $tmp[$i];
}
$reg2gene{ "CORE" . $core_num } = $cur_gene_list;
#Adding information to json array for the Region to be used in the javascript but not shown in the tba
$jso{Region}->[$c]->{type_ref} = $core_type;
$jso{Region}->[$c]->{oth_ref} = $core_oth;
$jso{Region}->[$c]->{core_clust} = $cur_chr;
$jso{Region}->[$c]->{href} = "CORE" . $core_num;
#Counting the number of regions as so to name them
$core_cnt++;
#Drawing the core region preview pane
#$tmp is a string holding the SVG and then the HTML values for the pane and then the page
my $tmp;
$tmp = &make_core_region_svg( $new, "core" . $core_cnt, \%cor_reg, \%ont, $tmp );
#Adding the preview pane to the SVG $out variable in the preview key at level 2
$out{$chr}->{preview}->{2} .= $tmp;
#Making the core region detail page
$tmp = &make_core_region_page( "core" . $core_cnt, \%cor_reg, $tmp );
$tmp =~ s/hidden/visible/g; #Probably can change this
open( TMP_OUT, ">", $out_dir . "/$output_id/" . $new );
print TMP_OUT $tmp;
close(TMP_OUT);
#Drawing the clickable area on the pan-chromosome
$out{$chr} = &svg_draw_arc_seg(
1, $cr_start, $st - 1, &convert_layer(1), &convert_layer(1) + $gene_height, "#" . $background, $new,
sprintf( "CORE REGION $cr_start to %d\n%d Core Genes", ( $st - 1 ), $cr_cnt ), "CORE" . $core_num,
"Region", $chr, $out{$chr}
);
#Counting the number of core regions
$core_num++;
}
#Resetting the core_region values
$cur_gene_list = "";
$core_type = "";
$core_oth = "";
#String with the file id of the
my $new = $out_dir . "/$output_id/FGI.$core_cnt.html";
#The FGI_region hash used to make the preview pane and the detail page
#Keys: st -> starting bp
# end -> ending bp
# number -> count of the different gene presence and absences paths observed in the fGR
# order -> array with the different gene paths sorted by the number of genomes with that order
# cnt -> hash with key of fgi ID with number of genomes containing the path
# The ID for most of these is the path ID
# st_id -> hash with key of fgi ID with the bounding core gene ID
# st_sz -> hash with key of fgi ID with the bounding core gene size
# st_dir -> hash with key of fgi ID with the bounding core gene direction
# st_col -> hash with key of fgi ID with the bounding core gene color
# end_id -> hash with key of fgi ID with the bounding core gene ID
# end_sz -> hash with key of fgi ID with the bounding core gene size
# end_dir -> hash with key of fgi ID with the bounding core gene direction
# end_col -> hash with key of fgi ID with the bounding core gene color
# gen -> string with semi-colon separated list of genomes with
# coords -> array with gene
my %fgi_reg;
my %fgi_ids; #Keeps track of the gene IDs in each gene path using a simple 1 if its in the path
$fgi_reg{st} = min( $st, $end ); #Ignoring the direction
$fgi_reg{end} = max( $st, $end );
$fgi_reg{number} = scalar( keys( %{ $list{$ID} } ) ); #Getting the counts of the genomes
#@sort is the based on the number of genes in the path this is subsequently used to sort the genes on the x-axis
#@sort2 is based on the number of genomes containing this path. This is the order used on the pane/page y-axis
my @sort =
sort { scalar( @{ $list{$ID}->{$b}->{arr} } ) <=> scalar( @{ $list{$ID}->{$a}->{arr} } ) }
keys( %{ $list{$ID} } );
my @sort2 = sort { $list{$ID}->{$b}->{cnt} <=> $list{$ID}->{$a}->{cnt} } keys( %{ $list{$ID} } );
$fgi_reg{order} = \@sort2;
#Going through each of the paths
foreach my $a (@sort2) {
$fgi_reg{cnt}->{$a} = ( $list{$ID}->{$a}->{cnt} + 0 );
#Adding information about the starting bounding core gene if it exists for this gene path
#See list above for the variables being collected
if ( $list{$ID}->{$a}->{st} ) {
$fgi_reg{st_id}->{$a} = $list{$ID}->{$a}->{st};
$fgi_reg{st_sz}->{$a} = abs( $core_size{ $list{$ID}->{$a}->{st} } );
$fgi_reg{st_dir}->{$a} = 1;
if ( $core_size{ $list{$ID}->{$a}->{st} } < 0 ) {
$fgi_reg{st_dir}->{$a} = -1;
}
if ( $gene_func{ $list{$ID}->{$a}->{st} } ) {
$fgi_reg{st_col}->{$a} = $rgb{ $gene_func{ $list{$ID}->{$a}->{st} } };
} else {
#Setting empty/blank color values to the default function
$fgi_reg{st_col}->{$a} = $rgb{$defaultFunct};
}
} else {
#The blank start core gene/break values
$fgi_reg{st_id}->{$a} = 0;
$fgi_reg{st_sz}->{$a} = 0;
$fgi_reg{st_dir}->{$a} = 0;
$fgi_reg{st_col}->{$a} = "ffffff";
}
if ( $list{$ID}->{$a}->{end} ) {
#Adding information about the starting bounding core gene if it exists for this gene path
#See list above for the variables being collected
$fgi_reg{end_id}->{$a} = $list{$ID}->{$a}->{end};
$fgi_reg{end_sz}->{$a} = abs( $core_size{ $list{$ID}->{$a}->{end} } );
$fgi_reg{end_dir}->{$a} = 1;
if ( $gene_func{ $list{$ID}->{$a}->{end} } ) {
$fgi_reg{end_col}->{$a} = $rgb{ $gene_func{ $list{$ID}->{$a}->{end} } };
} else {
$fgi_reg{end_col}->{$a} = $rgb{$defaultFunct};
}
} else {
#No bounding end core gene
$fgi_reg{end_id}->{$a} = 0;
$fgi_reg{end_sz}->{$a} = 0;
$fgi_reg{end_dir}->{$a} = 0;
$fgi_reg{end_col}->{$a} = "ffffff";
}
$fgi_reg{end_type}->{$a} = $gene_func{ $list{$ID}->{$a}->{end} };
$fgi_reg{st_type}->{$a} = $gene_func{ $list{$ID}->{$a}->{st} };
$fgi_reg{gen}->{$a} = $list{$ID}->{$a}->{gen};
}
my %mat;
#mat is hash to store the array with the genes so it can be ordered and displayed
#keys are arr -> array with the genes in order; loc -> hash where the gene id points to its order
#For the first gene in the fgi
for ( my $j = 0 ; $j < scalar( @{ $list{$ID}->{ $sort[0] }->{arr} } ) ; $j++ ) {
$mat{arr}->[$j] = $list{$ID}->{ $sort[0] }->{arr}->[$j];
my $a = $list{$ID}->{ $sort[0] }->{arr}->[$j];
if ($a) {
if ( !$clusters{$a} ) {
die("Cannot find a cluster for $a. Ending PanACEA...\n\n");
}
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{href} = "CL_$a";
#Adding the region information to the table. This could allow for clicking on a gene to highlight it's location in the region
if (! $jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{$head{Gene}->[3]})
{
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{ $head{Gene}->[3] } = $ID;
}
else
{
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{ $head{Gene}->[3] } .= ",". $ID;
}
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{core_clust} = $cur_chr;
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{detail} =
"CL_" . $list{$ID}->{ $sort[0] }->{arr}->[$j];
$cur_gene_list .= "CL_" . $list{$ID}->{ $sort[0] }->{arr}->[$j] . ";";
$mat{loc}->{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } = $j + 1;
#Adding values to the functional annotation of the Region, both the type ref and the other ref
while (
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{type_ref} =~ /([^\;]+)/g ) {
my $in1 = $1;
if ( $core_type !~ /$in1/ ) { $core_type .= $in1 . ";"; }
}
while (
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{oth_ref} =~ /([^\;]+)/g ) {
my $in1 = $1 . ";";
if ( $ont{$in1} ) {
$term_cnt{$in1}++;
$terms->{$in1} = $ont{$in1};
}
if ( !$term2type{$in1} ) {
my $type_ref =
$jso{Gene}->[ $gene_num{ $list{$ID}->{ $sort[0] }->{arr}->[$j] } ]->{type_ref};
$term2type{$in1} = $type_ref;
if ( !$type2term{$type_ref} ) {
$type2term{$type_ref}->[0] = $in1;
} else {
push @{ $type2term{$type_ref} }, $in1;
}
}
if ( $core_oth !~ /$in1/ ) { $core_oth .= $in1; }
}
#Flags the gene as being in the fgi
$fgi_ids{ $list{$ID}->{ $sort[0] }->{arr}->[$j] }->{ $sort[0] } = 1;
}
}
#For all subsequent genes in the fGI
for ( my $i = 1 ; $i < scalar(@sort) ; $i++ ) {
#Prev is the order value, which is used to insert a previously unseen gene into the order
my $prev = 0;
#For each of the paths with the Gene
for ( my $j = 0 ; $j < scalar( @{ $list{$ID}->{ $sort[$i] }->{arr} } ) ; $j++ ) {
my $a = $list{$ID}->{ $sort[$i] }->{arr}->[$j];
if ($a) {
$fgi_ids{$a}->{ $sort[$i] } = 1;
#If there is an order position for the gene, use it
if ( $mat{loc}->{$a} ) {
$prev = $mat{loc}->{$a};
} else {
#Otherwise, cycle through all the previously seen genes to get a count
for ( my $k = scalar( @{ $mat{arr} } ) ; $k > ($prev) ; $k-- ) {
$mat{arr}->[$k] = $mat{arr}->[ $k - 1 ];
$mat{loc}->{ $mat{arr}->[$k] }++;
}
$mat{arr}->[$prev] = $a;
$mat{loc}->{$a} = $prev + 1;
$prev++;
}
}
}
}
#If the matrix array exists
if ( $mat{arr} ) {
#Go through the array of genes in the fgi
for ( my $i = 0 ; $i < scalar( @{ $mat{arr} } ) ; $i++ ) {
my $a = $mat{arr}->[$i];
if ($a) {
$fgi_grps{ $fgi_member{$a} }->{hit} = $ID;
if ( !$clusters{$a} ) { die("$ID $i $a"); }
#Put the gene information into the coordinate array
push(
@{ $fgi_reg{cords} },
[
$i, $clusters{$a}->{num_of_members},
$clusters{$a}->{protein_name}, $fgi_member{$a},
$fgi_ids{$a}, "CL_$a",
$rgb{ $jso{Gene}->[ $gene_num{$a} ]->{type_ref} }, $clusters{$a}->{start},
$clusters{$a}->{end}, $jso{Gene}->[ $gene_num{$a} ]
]
);
}
}
my $c = 0;
#If there are any genes in the fGR...
if ( scalar( @{ $fgi_reg{cords} } ) > 0 ) {
#if there are- make the pane and the page
my $fig;
#drawing the fgi page and writing it to a html file
$fig = &make_fgi_page_svg( "fgi" . $ID, \%fgi_reg, \%ont, $fig );
$new = "fgi/$ID.html";
open( OUT, ">", $out_dir . "/$output_id/" . $new );
print OUT "$fig\n";
close(OUT);
#drawing the preview pane
$out{$chr}->{preview}->{2} .= &make_fgi_region_svg( $new, \%fgi_reg, \%ont, $out{$chr}->{2} );
}
#Counting the number of regions to keep in the table
if ( $jso{Region} ) { $c = scalar( @{ $jso{Region} } ); }
#Putting the region information into a table row for the table
my @tmp = ( "FGR.$core_cnt", "fGR", scalar( @{ $fgi_reg{cords} } ), $core_type, $ID );
for ( my $i = 0 ; $i < scalar( @{ $head{Region} } ) ; $i++ ) {
$jso{Region}->[$c]->{ $head{Region}->[$i] } = $tmp[$i];
}
#More table information for put in the table jso
$jso{Region}->[$c]->{href} = $ID;
$jso{Region}->[$c]->{core_clust} = $cur_chr;
$jso{Region}->[$c]->{oth_ref} = $core_oth;