/* Last edited: May 6 10:53 2003 (edgrif) */ /* belvu.c - Pretty colour viewing of multiple alignments ------------------------------------------------------------- | File: belvu.c | | Author: Erik.Sonnhammer@cgb.ki.se | | Copyright (C) E Sonnhammer | ------------------------------------------------------------- Date Modification -------- --------------------------------------------------- 94-08-06 Created 95-01-23 Removed multiple windows - call xfetch instead. Made sure all slx & mul formats are read properly Built-in ruler 95-04-25 Changed internal storage system into an array of ALN structs. This allowed direct selex support and any size of alignment. 95-04-28 Added MSF support, highlighting of whole line and new conserved residue colouring. 95-05-01 Changed Conservation Colouring to work by colorMap matrix. 95-05-08 Added 'Draw Wrapped Alignment' 95-08-25 [1.4] Added "Make non-redundant" and "Remove highlighted sequence". 95-08-29 [1.4] Added "Remove columns". 95-09-21 [1.4] Added "Remove columns that are all-gaps" 95-10-10 [2.0] Completely rewrote layout and made own scrollbars. 95-10-23 [2.1] Re-implemented "average similarity" coloring. This has the effect over %ID that columns with many different types of similar residues get higher scores, while columns with some predominant residue type will get lower scores if the other residues are very different. There is a fudge-factor of 2 for Blosum62, which is totally ad-hoc. In the old implementation it was set to 4 (lowest Blosum62 diagonal) but that is too weak. With a factor of 2, the highest values becomes 550% (W-W)... Whether this is better than %ID God knows. 95-10-25 [2.2] Added optional column with scores. [2.2] Added 'Ask to save' if modified. 95-10-30 [2.2] Added -t Title [2.2] Added Centering with middle button. NOTE: graphBoxShift was too slow. [2.2] Replaced lastboxAbs with highlight() - lights up all boxes of protein. [2.2] Added checkmarks in menu and fixed Toby's colours. 95-11-06 [2.2] Fixed bug that lost Highlight_alnp when re-sorting alignment. 95-11-10 [2.2] Window sized to fit small alignments. 95-12-11 [2.2] Added "Add sequences with matching segments". 96-01-25 [2.2] Fixed -m crashes by actually writing countInserts() (had forgotten). 96-03-08 [2.2] Fixed bug that anticipated '*' to mean 'any insertion'. 96-03-21 [2.2] Enabled efetching when double clicking in alignment. Fixed unwanted scrolling when clipping to the left. Added 'Remove outliers'. Same directory-buffer for all I/O operations. 96-04-01 [2.3] Adapts to any name length up to 50 residues. Added "Only colour residues above %id threshold". Added "Save current colour scheme to file". 96-05-20 [2.3] "Use gray shades (for printing)" with reverse video. 96-08-23 [2.3] Improved WWW efetch. 96-08-28 [2.4] Made conservation colours and cutoffs changeable. 96-11-19 [2.4] Added "Read labels of highlighted sequence and spread them". 97-03-12 [2.4] Added "Ignore gaps in conservation levels". 97-03-18 [2.4] Removed dirName and fileName from filqueryopen to avoid crashes. 97-06-01 [2.4] Changed setConservColors() for color_by_similarity to not count self-pair residues. This led to ghost conservation for e.g. single C's and H's at small samples (<~3). 97-07-29 [2.5] Added tree-construction, tree-sorting and tree-drawing. Changed toolbar to more microsoft'ish look and rephrased various options. Fixed bug that would screw up color schemes when changing coloring mode with color codes window up. Tried to make it blocking but got strange crashes. The graph libarary does not like to redraw other windows from a blocking window - gets confused about which one to block. Solved by temporarily changing to harmless colour menu, which prevents changing mode. Added 'Plot conservation profile'. Added conservation values to '-c Print Conservation table'. 97-10-06 [2.5] Added '-' as legal gap character (in addition to '.'). Added toggle "Automatically remove all-gap columns after sequence removals" Added pickable 'markup' sequences, for consensus or markup lines. 97-12-18 [2.6] Added hiding/unhiding of picked line. 97-12-30 [2.7] Added MSF writing and multiple-choice of output formats on command line. Added "Save As..." tool and changed save routines to use the global variables saveCoordsOn, saveSeparator, and saveFormat for the format details. Added parsing of aligned Fasta. Discontinued ProDom parsing support in favor of Aligned Fasta format. 97-12-30 [2.7] Fixed bug so negative scores can be read in. Added command-line "make non-redundant". 99-03-02 [2.7] Added "Remove partial sequences". Fixed bug in rmOutliers that skipped one seq after a deletion. 99-03-02 [2.7] Added "Remove gappy sequences" and "Remove gappy columns". 99-06-05 [2.8] Added support for Stockholm format GR and GC markup lines. (In the belvu code, Stockholm format is still called MUL). 99-06-23 [2.8] Added "Sort by similarity/identity to highlighted sequence". 99-08-05 [2.8] Added "Print tree in New Hampshire format". 99-08-08 [2.8] Code changes in sorting routintes so that #=GR markup sticks to the sequence line it belongs to (markups are removed before sorting by separateMarkupLines() and put back in the right place afterwards by reInsertMarkupLines(). Added commandline option to sort by similarity to 1st sequence. 99-10-26 [2.8] Changed tree calculation and drawing to handle multiple trees simultaneously that can still be printed. Moved tree line width to tree option (treesettings) control window. 00-01-24 [2.8] Added Collapse functionality to Wrap-around display (residues between []) (Should change this to a markup-line, much easier) Removed right-hand coordinate in wrap-around display. 00-01-25 [2.8] Added "select gap character". 00-03-03 [2.8] Fixed bug in parseMul so that seqnames can contain '-'. 00-03-06 [2.8] Changed code for foreground coloring to sync normal and wrapped displays. 00-03-21 [2.8] Fixed bug in Tree-mode selection from command line (string didn't change). 00-04-27 [2.9] Center on highlighted after sorting. 00-05-02 [2.9] Show Highlighted in tree. 00-05-15 [2.9] Color tree by organism. User-choosable colors (from color setup file). Default tree color scheme with decent colors. Reroot tree by clicking. "Find orthologs option in tree - marks up nodes with different left->org and right->org. 00-05-20 [2.9] Print out list of found orthologs. 00-05-23 [2.9] Added back dirName and fileName to filqueryopen - what was the problem? Also added default file name = the original file name. Reroot tree on "best balanced" node. 00-08-04 [2.9] Added output tree only in batch mode. Fixed BUG: "sort by tree" -> tree highlighting finds wrong sequence in alignment. 00-08-10 [2.9] Added batch-mode generation of bootstrapped trees. Added trees without sequence coordinates. 00-08-16 [2.9] Fixed bugs in treeFindBalance and treeSize3way that chose the wrong root. - treeSize3way had bugs in the branchlength selection when going up parents. - treeFindBalance had bugs that ameliorated bugs in treeSize3way. - treeFindBalance now chooses "best subtree balance" tree when more than one rooting give perfect balance. 00-08-16 [2.9] Fixed bug in bootstrapping (the original wasn't kept separately). 00-10-19 [2.9] Added "core name" (without coordinates) parsing of #=GS OS lines. 00-10-19 [2.9] Fixed bug in 'find orthologs': if the node was originally the root, the organism spreading was incorrect. Fix: spread organisms to higher nodes after the tree is calculated. Added average conservation to "Plot conservation" window. 00-10-25 [2.9] Added "Show Organisms" window. Changed mknonredundant to remove only inclusive matches. 00-11-13 [2.9] Fixed bug in #=GS OS parsing (messed up the original). 00-12-11 [2.9] Added batch mode "Remove Gappy Columns". 01-09-11 [2.10] Fixed bugs in sequence coordinate handling. 01-11-28 [2.11] Added Sonnhammer-Storm distance correction. 01-12-04 [2.12] Fixed a bug in parsing alignments without coordinates (from 2.10). Added treeScale to scale tree according to mode (1.0 normal, 0.3 distcorr). 02-03-22 [2.13] Added reading of user-choosable label for organism info (e.g. OC). 02-04-14 [2.14] Fixed small but segfaulting bug in stripCoordTokens. 02-05-24 [2.15] Added "Show organism" in the tree. 02-07-08 [2.16] Converted to GTK and ACEIN. Pending: Internal bootstrapping. Rationale for group bootstrapping: For each sequence, store array of sequences in merging order in bootstrap tree. For each node in original tree, check if group of sequences corresponds to the first n sequences in bootstraptree array. Rationale for exact bootstrapping: Traverse bootstrap tree from leaf to root: oritree nodes = {same, different, 0} if issequence(parent->other_child) { if (parent->other_child == boottree(parent->other_child)) parent = same else parent = different if isseen(parent->other_child) { if (parent->other_child == boottree(parent->other_child)) parent = same else parent = diferent if same, continue to parent gap penalties: add in scoring function. Clean up the parsing code with str2aln calls. Bug reported by Jerome Gracy: Linux version can crash in readMul. read in other trees with bootstraps use confidence cutoff in find orthologs . Use similiarity -> distance in tree calc. . Include gaps in similarity calculation. make alignment collapsing easier to use, see above. Keyboard edit of one sequence. How to find the right place ? !! Undo edits. Would like to abort parsing if two sequence lines accidentally have same name - How? - this is used in selex format... Koonin ideas: Consensus pattern lines (Default 0% and 100% lines, maybe 90% too). Add/remove any nr of lines cutoffs 0 - 100%. Tool to define groups of residues and corresponding symbol. These are shown on pattern lines (priority to smaller groups). Low priority: "Add sequences with matching segments" for more than one sequence. Will probably be very deleterious to the alignment. Ability to choose both foreground and background color - makes it twice as slow - see "Black/white for printing. */ /* Description of color_by_similarity algorithm in setConservColors(): ( Corresponds to summing all pairwise scores) for each residue i { for each residue j { if (i == j) score(i) += (count(i)-1)*count(j)*matrix(i,j) else score(i) += count(i)*count(j)*matrix(i,j) } } if (ignore gaps) n = nresidues(pos) else n = nsequences if (n == 1) id = 0.0 else id = score/(n*(n-1)) } */ #include #include #include #include #include #include #include #include #include /*#include / * Needed for RAND_MAX but clashes with other stuff */ #include #include #define MAXLINE 512 #define MAXLENGTH 10000 #define MAXNAMESIZE 256 #define boxColor WHITE /* LIGHTGRAY */ #define DBL_CLK_DELAY 2 /* seconds */ #define HSCRWID 1.0 /* Width of Horizontal scroll bar */ #define SCRBACKCOLOR PALEGRAY #define NOCOLOR -1 #define NN NOCOLOR #define BG WHITE enum { MUL, RAW }; /* Input formats for IN_FORMAT */ enum { dummy, GF, GC, GS, GR }; /* Markup types in Stockholm format */ enum { UNCORR, KIMURA, STORMSONN}; /* Distance correction for tree building */ /* Residue colors */ static int color[] = { NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN, NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN,NN }; /* Residue colors */ static int markupColor[] = { BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG, BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG,BG }; static int oldColor[256]; /* BLOSUM62 930809 A R N D C Q E G H I L K M F P S T W Y V B Z X \* */ static int BLOSUM62[24][24] = { 4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0, -4, -1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1, -4, -2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1, -4, -2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1, -4, 0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2, -4, -1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1, -4, -1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1, -4, 0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1, -4, -2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1, -4, -1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1, -4, -1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1, -4, -1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1, -4, -1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1, -4, -2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1, -4, -1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2, -4, 1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0, -4, 0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0, -4, -3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2, -4, -2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1, -4, 0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1, -4, -2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1, -4, -1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1, -4, 0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -1 }; #define NA 0 static int a2b[] /* ASCII-to-binary translation table */ = { NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA, 1,NA, 5, 4, 7,14, 8, 9,10,NA,12,11,13, 3,NA, 15, 6, 2,16,17,NA,20,18,NA,19,NA,NA,NA,NA,NA,NA, NA, 1,NA, 5, 4, 7,14, 8, 9,10,NA,12,11,13, 3,NA, 15, 6, 2,16,17,NA,20,18,NA,19,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA }; static int a2b_sean[] = { NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA, 1,NA, 2, 3, 4, 5, 6, 7, 8,NA, 9,10,11,12,NA, 13,14,15,16,17,NA,18,19,NA,20,NA,NA,NA,NA,NA,NA, NA, 1,NA, 2, 3, 4, 5, 6, 7, 8,NA, 9,10,11,12,NA, 13,14,15,16,17,NA,18,19,NA,20,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA }; typedef struct alnStruct { char name[MAXNAMESIZE+1]; int start; int end; int len; /* Length of this sequence - Do not use! Use maxLen instead ! */ char *seq; int nr; /* Ordinal number in array - must be sorted on this */ char fetch[MAXNAMESIZE+11]; float score; int color; /* Background color of name */ int markup; /* Markup line */ int hide; /* Hide this line */ int nocolor; /* Exclude from coloring */ char *organism; } ALN; typedef struct segStruct { int start; int end; int qstart; int qend; struct segStruct *next; } SEG; typedef struct treenode { /* KEEP IN SYNC WITH TREECPY !!!! */ float dist; /* Absolute distance position */ float branchlen; /* Length of branch to higher node */ float boot; struct treenode *left; struct treenode *right; struct treenode *parent; char *name; char *organism; ALN *aln; int box; int color; /* KEEP IN SYNC WITH TREECPY !!!! */ } treeNode; typedef struct treestruct { treeNode *head; int lastNodeBox; /* Last box used by a node - name boxes come next */ int currentPickedBox; } treeStruct; static char b2a[] = ".ARNDCQEGHILKMFPSTWYV" ; static char *colorNamesOld[] = { "WHITE", "BLACK", "LIGHTGRAY", "DARKGRAY", "RED", "GREEN", "BLUE", "YELLOW", "CYAN", "MAGENTA", "LIGHTRED", "LIGHTGREEN", "LIGHTBLUE", "DARKRED", "DARKGREEN", "DARKBLUE"}; static char *colorNames[NUM_TRUECOLORS] = { "WHITE", "BLACK", "LIGHTGRAY", "DARKGRAY", "RED", "GREEN", "BLUE", "YELLOW", "CYAN", "MAGENTA", "LIGHTRED", "LIGHTGREEN", "LIGHTBLUE", "DARKRED", "DARKGREEN", "DARKBLUE", "PALERED", "PALEGREEN", "PALEBLUE", "PALEYELLOW", "PALECYAN", "PALEMAGENTA", "BROWN", "ORANGE", "PALEORANGE", "PURPLE", "VIOLET", "PALEVIOLET", "GRAY", "PALEGRAY", "CERISE", "MIDBLUE" }; static int treeColors[16] = { RED, BLUE, DARKGREEN, ORANGE, MAGENTA, BROWN, PURPLE, CYAN, VIOLET, MIDBLUE, CERISE, LIGHTBLUE, DARKRED, GREEN, DARKBLUE, GRAY }; static void externalCommand(); /* Global variables *******************/ enum { COLORBYRESIDUE, COLORSCHEMESTANDARD, COLORSCHEMEGIBSON, COLORSCHEMECYS, COLORSCHEMEEMPTY, COLORSIM, COLORID, COLORIDSIM }; /* Modes - used for checkmarks */ enum { NOLL, SORT_ALPHA, SORT_ORGANISM, SORT_TREE, SORT_SCORE, SORT_SIM, SORT_ID }; /* Initial sorting modes */ enum {UPGMA, NJ}; /* Tree building methods */ static int **conservCount=0, /* Matrix of conservation values - 21 x maxLen */ **colorMap=0, /* Matrix of conservation colours - 21 x maxLen */ *conservResidues=0, /* Array of number of residues present in each column */ colorScheme = COLORSIM, /* Current colour scheme mode (used for checkmarks) */ color_by_conserv = 1, /* 1 => by id or similarity; 0 => by residue */ id_blosum=1, /* Paint similar residues too */ color_by_similarity = 1, /* 0 => by id */ maxLen=0, /* Columns in alignment */ nseq=0, /* Rows in alignment */ verbose=0, debug=0, xmosaic = 0, ip, /* Int pointer used in array operations */ init_sort=0, init_rmPartial=0, rmPickLeftOn, pickedCol = 0, /* Last picked column */ colorRectangles=1, maxNameLen, /* Max string length of any sequence name */ maxStartLen=0, /* Max string length of any sequence start */ maxEndLen=0, /* Max string length of any sequence end */ maxScoreLen=0, /* Max string length of any sequence score */ AlignXstart=0, /* x Start position of alignment in alignment window */ AlignYstart=0, /* y Start position of alignment in alignment window */ AlignXend, /* x End position of alignment in alignment window */ AlignYend, /* y End position of alignment in alignment window */ Alignwid, /* Witdh of shown alignment */ Alignheig, /* Height of shown alignment */ statsBox, HscrollSliderBox, VscrollSliderBox, scrLeftBox, scrRightBox, scrUpBox, scrDownBox, scrLeftBigBox, scrRightBigBox, scrUpBigBox, scrDownBigBox, colorButtonBox, sortButtonBox, editButtonBox, maxBox, midBox, lowBox, gridOn = 0, IN_FORMAT = MUL, saved = 1, displayScores = 0, Vcrosshair, Highlight_matches, colorByResIdOn = 0, /* colour by residue type above identity cutoff */ matchFooter = 0, maxfgColor = BLACK, midfgColor = BLACK, lowfgColor = BLACK, maxbgColor = CYAN, midbgColor = MIDBLUE, lowbgColor = LIGHTGRAY, oldmaxfgColor, oldmidfgColor, oldlowfgColor, oldmaxbgColor, oldmidbgColor, oldlowbgColor, printColorsOn = 0, wrapLinelen, ignoreGapsOn = 0, treeOrderNr, /* Tree node number */ conservationWindow=1, rmEmptyColumnsOn=1, saveCoordsOn=1, treeSampleOn=0, treeMethod = NJ, maxTreeWidth = 0, treeColorsOn = 1, treeShowBranchlen = 0, treeShowOrganism = 1, treeCoordsOn = 1, treebootstraps = 0, treebootstrapsDisplay = 0, stripCoordTokensOn = 1; static float Xpos, /* x screen position of alignment window */ Aspect, /* Aspect ratio of coordinate system (width/height) */ VSCRWID, /* Width of Vertical scroll bar - depends on Aspect */ HscrollStart, /* Start of legal Horizontal scrollbar slider area */ HscrollEnd, /* End of -''- */ VscrollStart, /* as above */ VscrollEnd, /* as above */ HsliderLen, /* Length of Horizontal slider */ VsliderLen, /* Length of Vertical slider */ SliderOffset, colorByResIdCutoff = 20.0, /* Colour by residue + id cutoff */ lowIdCutoff = 0.4, /* %id cutoff for lowest colour */ midIdCutoff = 0.6, /* %id cutoff for medium colour */ maxIdCutoff = 0.8, /* %id cutoff for maximum colour */ lowSimCutoff = 0.5, /* %id cutoff for lowest colour */ midSimCutoff = 1.5, /* %id cutoff for medium colour */ maxSimCutoff = 3.0, /* %id cutoff for maximum colour */ oldMax, oldMid, oldLow, screenWidth, screenHeight, fontwidth, fontheight, tree_y, treeLinewidth = 0.3, conservationLinewidth = 0.2, *conservation, /* The max conservation in each column [0..maxLen] */ conservationRange, conservationXScale=1, conservationYScale=2, oldlinew, treeBestBalance, treeBestBalance_subtrees, treeDistCorr = UNCORR, treeScale = 1.0; static char belvuVersion[] = "2.16", line[MAXLENGTH+1], /* General purpose array */ Title[256] = "", /* Window title */ *ruler=0, /* Residue ruler on top */ stats[256], /* Status bar string */ dirName[DIR_BUFFER_SIZE+1], /* Default directory for file browser */ fileName[FIL_BUFFER_SIZE+1], maxText[11], /* Cutoff text arrays for boxes */ midText[11], lowText[11], saveFormat[50], treeMethodString[50] = "NJ", treeDistString[50] = "UNCORR", saveSeparator='/', gapChar='.', *OrganismLabel = "OS"; static Graph belvuGraph, showColorGraph, treeGraph, conservationGraph, saveGraph, annGraph, treeGUIGraph, gapCharGraph, organismsGraph; static Stack stack ; static Stack AnnStack ; static Array Align, MarkupAlign, organismArr; static ALN aln, /* General purpose ALN */ *alnp, /* General purpose ALNP */ *Highlight_alnp=0; static treeStruct *treestruct; /* General purpose treeStruct */ static treeNode *treeHead, /* Global current treehead */ *treeBestBalancedNode; ACEIN ace_in; static int stripCoordTokens(char *cp); static void readColorCodesMenu(void); static void readColorCodes(FILE *file, int *colorarr); static void showColorCodesRedraw(void); static void saveColorCodes(void); static void setColorCodes(void); static void belvuRedraw(void); static void colorSchemeStandard(void); static void colorByResidue(void); static void colorSchemeGibson(void); static void colorSchemeCys(void); static void colorSchemeEmpty(void); static void colorId(void); static void colorIdSim(void); static void colorSim(void); static void colorRect(void); static void xmosaicFetch(void); static void saveAlignment(void); static void saveMul(void); static void saveFasta(void); static void saveMSF(void); static void saveRedraw(void); static void saveMSFSelect(void); static void saveMulSelect(void); static void saveFastaAlnSelect(void); static void saveFastaSelect(void); static void treeUPGMAselect(void); static void treeNJselect(void); static void treeUNCORRselect(void); static void treeKIMURAselect(void); static void treeSTORMSONNselect(void); static void wrapAlignmentWindow(void); static void mkNonRedundant(float cutoff); static void mkNonRedundantPrompt(void); static void rmOutliers(void); static void rmPartial(void); static void rmGappySeqs(float cutoff); static void rmGappySeqsPrompt(void); static void readLabels(void); static void selectGaps(void); static void listIdentity(void); static float identity(char *s1, char *s2); static void rmPicked(void); static void rmPickLeft(void); static void Help(void); static void rmColumnPrompt(void); static void rmColumnCutoff(void); static void rmColumnLeft(void); static void rmColumnRight(void); static void rmEmptyColumns(float cutoff); static void rmEmptyColumnsInteract(void); static void rmScore(void); static void alphaSort(void); static void simSort(void); static void idSort(void); static void organismSort(void); static void scoreSortBatch(void); static void scoreSort(void); static void treeSettings(void); static void treeDisplay(void); static void treeSortBatch(void); static void treeSort(void); static void printScore(void); static void belvuDestroy(void); static void readMatch(FILE *fil); static void readMatchPrompt(void); static int findCommand (char *command, char **retp); static void colorByResId(void); static void printColors (void); static void ignoreGaps (void); static void finishShowColorOK(void); static void finishShowColorCancel(void); static void initResidueColors(void); static void initMarkupColors(void); static int isAlign(char c); static int isGap(char c); static void treePrintNH_init(void); static void treeFindOrthologs(void); static void treeSetLinewidth(char *cp); static void treeSetScale(char *cp); static void conservationPlot(void); static void conservationSetWindow(void); static void conservationSetLinewidth(void); static void conservationSetScale(void); static void rmEmptyColumnsToggle(void); static void markup(void); static void hide(void); static void unhide(void); static int score(char *s1, char *s2); static void arrayOrder(void); static void arrayOrder10(void); static void highlightScoreSort(char mode); static int bg2fgColor (int bkcol); static void alncpy(ALN *dest, ALN *src); static int alignFind( Array arr, ALN *obj, int *ip); static void centerHighlighted(void); static void str2aln(char *src, ALN *alnp); static treeNode *treeReroot(treeNode *node); static void treeDraw(treeNode *treeHead); static treeNode *treecpy(treeNode *node); static void treeTraverse(treeNode *node, void (*func)()); static void treeBalanceByWeight(treeNode *node1, treeNode *node2, float *llen, float *rlen); static float treeSize3way(treeNode *node, treeNode *fromnode); static void showOrganisms(void); static void parseMulLine(char *line, ALN *aln); #define xmosaicFetchStr "Use efetch via WWW" static MENU belvuMenu; static MENUOPT mainMenu[] = { belvuDestroy, "Quit", Help, "Help", wrapAlignmentWindow, "Wrap-around alignment window (pretty print)", graphPrint, "Print this window", menuSpacer, "", treeDisplay, "Make tree from current alignment", treeSettings, "Tree options", menuSpacer, "", conservationPlot, "Plot conservation profile", menuSpacer, "", saveAlignment, "Save alignment", saveRedraw, "Save alignment as...", printScore, "Output score and coords of line", /* readMatchPrompt, "Add sequence with matching segments", */ /* Too difficult to use interactively - rely solely on matchFooters instead */ menuSpacer, "", xmosaicFetch, xmosaicFetchStr, listIdentity, "Compare all vs. all, list identity", graphCleanUp, "Clean up windows", 0, 0 }; static MENUOPT treeMenu[] = { graphDestroy, "Quit", graphPrint, "Print", treePrintNH_init, "Save tree in New Hampshire format", treeFindOrthologs, "Find putative orthologs", showOrganisms, "Show current organisms", 0, 0 }; static MENUOPT conservationMenu[] = { graphDestroy, "Quit", graphPrint, "Print", conservationSetWindow, "Set window size", conservationSetScale, "Set plot scale factors", conservationSetLinewidth, "Set line width of tree", 0, 0 }; static MENUOPT wrapMenu[] = { graphDestroy, "Quit", graphPrint, "Print", wrapAlignmentWindow, "New Wrap-around alignment window", 0, 0 }; #define colorByResidueStr "Colour scheme: By residue, Last palette" #define colorSchemeStandardStr "Colour scheme: By residue, Erik's" #define colorSchemeGibsonStr "Colour scheme: By residue, Toby's" #define colorSchemeCysStr "Colour scheme: By residue, Cys/Gly/Pro" #define colorSchemeEmptyStr "Colour scheme: By residue, Clean slate" #define colorSimStr "Colour scheme: Average similarity by Blosum62" #define colorIdStr "Colour scheme: Percent identity only" #define colorIdSimStr "Colour scheme: Percent identity + Blosum62" #define colorRectStr "Display colors (faster without)" #define printColorsStr "Use gray shades (for printing)" #define ignoreGapsStr "Ignore gaps in conservation calculation" #define thresholdStr "Only colour residues above %id threshold" static MENU colorMenu; static MENUOPT colorMENU[] = { colorSchemeStandard, colorSchemeStandardStr, colorSchemeGibson, colorSchemeGibsonStr, colorSchemeCys, colorSchemeCysStr, colorSchemeEmpty, colorSchemeEmptyStr, colorByResidue, colorByResidueStr, colorByResId, thresholdStr, /* setColorCodes, "Set colours manually",*/ saveColorCodes, "Save current colour scheme", readColorCodesMenu, "Read colour scheme from file", menuSpacer, "", colorSim, colorSimStr, colorId, colorIdStr, colorIdSim, colorIdSimStr, ignoreGaps, ignoreGapsStr, printColors, printColorsStr, menuSpacer, "", markup, "Exclude highlighted from calculations on/off", colorRect, colorRectStr, showColorCodesRedraw, "Open window to edit colour scheme", 0, 0 }; static MENU colorEditingMenu; static MENUOPT colorEditingMENU[] = { showColorCodesRedraw, "Close Colour Codes window first", 0, 0 }; #define rmEmptyColumnsStr "Automatically remove all-gap columns after sequence removals" static MENU editMenu; static MENUOPT editMENU[] = { rmPicked, "Remove highlighted line", rmPickLeft, "Remove many sequences", rmGappySeqsPrompt, "Remove gappy sequences", rmPartial, "Remove partial sequences", mkNonRedundantPrompt,"Make non-redundant", rmOutliers, "Remove outliers", rmScore, "Remove sequences below given score", menuSpacer, "", rmColumnPrompt, "Remove columns", rmColumnLeft, "<- Remove columns to the left of cursor (inclusive)", rmColumnRight, "Remove columns to the right of cursor (inclusive) -> ", rmColumnCutoff, "Remove columns according to conservation", rmEmptyColumnsInteract, "Remove gappy columns", rmEmptyColumnsToggle, rmEmptyColumnsStr, menuSpacer, "", menuSpacer, "", readLabels, "Read labels of highlighted sequence and spread them", selectGaps, "Select gap character", hide, "Hide highlighted line", unhide, "Unhide all hidden lines", 0, 0 }; static MENU showColorMenu; static MENUOPT showColorMENU[] = { finishShowColorCancel, "CANCEL", finishShowColorOK, "OK", 0, 0 }; static MENU sortMenu; static MENUOPT sortMENU[] = { scoreSort, "Sort by score", alphaSort, "Sort alphabetically", organismSort, "Sort by swissprot organism", treeSort, "Sort by tree order", simSort, "Sort by similarity to highlighted sequence", idSort, "Sort by identity to highlighted sequence", 0, 0 }; #define MSFStr "MSF" #define MulStr "Stockholm (Pfam/HMMER)" #define FastaAlnStr "Aligned Fasta" #define FastaStr "Unaligned Fasta" static MENU saveMenu; static MENUOPT saveMENU[] = { saveMSFSelect, MSFStr, saveMulSelect, MulStr, saveFastaAlnSelect, FastaAlnStr, saveFastaSelect, FastaStr, 0, 0 }; #define UPGMAstr "UPGMA" #define NJstr "NJ" static MENU treeGUIMenu; static MENUOPT treeGUIMENU[] = { treeUPGMAselect, UPGMAstr, treeNJselect, NJstr, 0, 0 }; #define UNCORRstr "UNCORR" #define KIMURAstr "KIMURA" #define STORMSONNstr "STORM & SONNAMMER" static MENU treeDistMenu; static MENUOPT treeDistMENU[] = { treeUNCORRselect, UNCORRstr, /* treeKIMURAselect, KIMURAstr,*/ treeSTORMSONNselect, STORMSONNstr, 0, 0 }; static void fatal(char *format, ...) { va_list ap; printf("\nFATAL ERROR: "); va_start(ap, format); vprintf(format, ap); va_end(ap); printf("\n"); exit(-1); } static void menuCheck(MENU menu, int mode, int thismode, char *str) { if (mode == thismode) menuSetFlags(menuItem(menu, str), MENUFLAG_TOGGLE_STATE); else menuUnsetFlags(menuItem(menu, str), MENUFLAG_TOGGLE_STATE); } static void setMenuCheckmarks(void) { menuCheck(colorMenu, colorScheme, COLORBYRESIDUE, colorByResidueStr); menuCheck(colorMenu, colorScheme, COLORSCHEMESTANDARD, colorSchemeStandardStr); menuCheck(colorMenu, colorScheme, COLORSCHEMEGIBSON, colorSchemeGibsonStr); menuCheck(colorMenu, colorScheme, COLORSCHEMECYS, colorSchemeCysStr); menuCheck(colorMenu, colorScheme, COLORSCHEMEEMPTY, colorSchemeEmptyStr); menuCheck(colorMenu, colorScheme, COLORSIM, colorSimStr); menuCheck(colorMenu, colorScheme, COLORID, colorIdStr); menuCheck(colorMenu, colorScheme, COLORIDSIM, colorIdSimStr); menuCheck(colorMenu, 1, colorRectangles, colorRectStr); menuCheck(colorMenu, 1, printColorsOn, printColorsStr); menuCheck(colorMenu, 1, ignoreGapsOn, ignoreGapsStr); menuCheck(editMenu, 1, rmEmptyColumnsOn, rmEmptyColumnsStr); if (!graphExists(showColorGraph)) graphNewBoxMenu(colorButtonBox, colorMenu); else graphNewBoxMenu(colorButtonBox, colorEditingMenu); graphNewBoxMenu(sortButtonBox, sortMenu); graphNewBoxMenu(editButtonBox, editMenu); menuCheck(belvuMenu, 1, xmosaic, xmosaicFetchStr); graphNewMenu(belvuMenu); } static void Help(void) { graphMessage (messprintf("\ Belvu - a multiple alignment viewer.\n\ \n\ Version %s\n\ Copyright (C) Erik Sonnhammer, 1995\n\ \n\ \n\ ALIGNMENT:\n\ %d sequences, %d residues wide.\n\ \n\ LEFT MOUSE BUTTON:\n\ Click on residue to show coordinate.\n\ Double click on sequence to fetch.\n\ Scrollbars: drag and page.\n\ Toolbar buttons:\n\ File: Display this text.\n\ Edit: Remove many sequences.\n\ Colour: Edit colour scheme window.\n\ Sort: Sort alphabetically.\n\ \n\ MIDDLE MOUSE BUTTON:\n\ In alignment: centre on crosshair.\n\ Scrollbars: drag and jump.\n\ \n\ RIGHT MOUSE BUTTON:\n\ Toolbar menus (file menu elsewhere).\n\ \n\ KEYBOARD:\n\ Arrow keys: Move screenfulls.\n\ Home: Go to Top.\n\ End: Go to Bottom.\n\ Insert: Go to Start of line.\n\ Delete: Go to End of line.\n\ For more details, see http://www.sanger.ac.uk/~esr/Belvu.html\n\ ", belvuVersion, nseq, maxLen)); } static void fetchAln(ALN *alnp) { if (xmosaic) { static char *browser=0; if (!browser) { printf("Looking for WWW browsers ...\n"); if (!findCommand("netscape", &browser) && !findCommand("Netscape", &browser) && !findCommand("Mosaic", &browser) && !findCommand("mosaic", &browser) && !findCommand("xmosaic", &browser)) { messout("Couldn't find any WWW browser. Looked for " "netscape, Netscape, Mosaic, xmosaic & mosaic"); return; } } printf("Using WWW browser %s\n", browser); fflush(stdout); system(messprintf("%s http://www.sanger.ac.uk/cgi-bin/seq-query?%s&", browser, alnp->fetch)); /* OLD: system(messprintf("xfetch '%s' &", alnp->fetch));*/ } else externalCommand(messprintf("efetch '%s' &", alnp->fetch)); } static void scrRight() { AlignXstart++; belvuRedraw(); } static void scrLeft() { AlignXstart--; belvuRedraw(); } static void scrUp() { AlignYstart--; belvuRedraw(); } static void scrDown() { AlignYstart++; belvuRedraw(); } static void scrRightBig() { AlignXstart += Alignwid; belvuRedraw(); } static void scrLeftBig() { AlignXstart -= Alignwid; belvuRedraw(); } static void scrUpBig() { AlignYstart -= Alignheig; belvuRedraw(); } static void scrDownBig() { AlignYstart += Alignheig; belvuRedraw(); } static void noAction(void) {} static void unregister(void) { graphRegister(LEFT_UP, noAction); graphRegister(MIDDLE_UP, noAction); graphRegister(LEFT_DRAG, noAction); graphRegister(MIDDLE_DRAG, noAction); } static void HscrollUp(double x, double y) { float X = x - SliderOffset; AlignXstart = (X - HscrollStart)/(HscrollEnd - HscrollStart) * maxLen; unregister(); belvuRedraw(); } static void HscrollDrag(double x, double y) { float X = x - SliderOffset; if (X < HscrollStart) X = HscrollStart; if (X + HsliderLen > HscrollEnd) X = HscrollEnd - HsliderLen; graphBoxShift(HscrollSliderBox, X, VscrollEnd+HSCRWID); graphRegister(LEFT_UP, HscrollUp); } static void VscrollUp(double x, double y) { float Y = y - SliderOffset; AlignYstart = (Y - VscrollStart)/(VscrollEnd - VscrollStart) * nseq; unregister(); belvuRedraw(); } static void VscrollDrag(double x, double y) { float Y = y - SliderOffset; if (Y < VscrollStart) Y = VscrollStart; if (Y + VsliderLen > VscrollEnd) Y = VscrollEnd - VsliderLen; graphBoxShift(VscrollSliderBox, HscrollEnd+VSCRWID, Y); graphRegister(LEFT_UP, VscrollUp); } static int x2col(double *x) { /* Round off to half units */ *x = floor(*x)+0.5; if (*x < HscrollStart - VSCRWID+1) *x = HscrollStart - VSCRWID+1; else if (*x > HscrollStart - VSCRWID + Alignwid) *x = HscrollStart - VSCRWID + Alignwid; return (int) (*x - (HscrollStart - VSCRWID + 0.5) + AlignXstart + 1); } static void updateStatus(ALN *alnp, int col) { int i, gaps=0, star; if (!alnp) return; strcpy(stats, " "); if (col) strcat(stats, messprintf("Column %d: ", col)); strcat(stats, messprintf("%s/%d-%d", alnp->name, alnp->start, alnp->end)); if (col) { strcat(stats, messprintf(" %c = ", alnp->seq[col-1])); for (star = i = 0; i < col; i++) { if (isGap(alnp->seq[i])) gaps++; else if (alnp->seq[i] == '*') star = 1; } if (star) { strcat(stats, "(unknown position due to insertion)"); } else { strcat(stats, messprintf("%d", col-1 + alnp->start - gaps)); } } strcat(stats, messprintf(" (%d match", Highlight_matches)); if (Highlight_matches > 1) strcat(stats, "es"); strcat(stats, ")"); graphBoxDraw(statsBox, BLACK, LIGHTBLUE); } static void xorVcrosshair(int mode, double x) { static double xold; if (mode == 1) graphXorLine(xold, VscrollStart-HSCRWID, xold, VscrollEnd+HSCRWID); graphXorLine(x, VscrollStart-HSCRWID, x, VscrollEnd+HSCRWID); xold = x; } static void middleDrag(double x, double y) { pickedCol = x2col(&x); if (pickedCol == Vcrosshair) return; updateStatus(Highlight_alnp, pickedCol); Vcrosshair = pickedCol; xorVcrosshair(1, x); } static void middleUp(double x, double y) { pickedCol = x2col(&x); AlignXstart = pickedCol - Alignwid/2; unregister(); belvuRedraw(); } static void leftDown(double x, double y) { float oldlinew = graphLinewidth(0.05); graphRectangle(0.5, floor(y), wrapLinelen+0.5, floor(y)+1); graphLinewidth(oldlinew); graphRedraw(); } static void middleDown(double x, double y) { if (x > HscrollStart && x < HscrollEnd && y > VscrollEnd+HSCRWID && y < VscrollEnd+2*HSCRWID) { graphRegister(MIDDLE_DRAG, HscrollDrag); graphRegister(MIDDLE_UP, HscrollUp); SliderOffset = HsliderLen/2; HscrollDrag(x, y); } else if (y > VscrollStart && y < VscrollEnd && x > HscrollEnd+VSCRWID && x < HscrollEnd+2*VSCRWID) { graphRegister(MIDDLE_DRAG, VscrollDrag); graphRegister(MIDDLE_UP, VscrollUp); SliderOffset = VsliderLen/2; VscrollDrag(x, y); } else if (x > HscrollStart-VSCRWID && x < HscrollEnd+VSCRWID && y > VscrollStart-HSCRWID && y < VscrollEnd+HSCRWID) { /* In alignment */ graphRegister(MIDDLE_DRAG, middleDrag); graphRegister(MIDDLE_UP, middleUp); pickedCol = x2col(&x); updateStatus(Highlight_alnp, pickedCol); Vcrosshair = pickedCol; xorVcrosshair(0, x); } } static int nrorder(ALN *x, ALN *y) { if (x->nr < y->nr) return -1; else if (x->nr > y->nr) return 1; else return 0; } static int scoreorder(ALN *x, ALN *y) { if (x->score < y->score) return -1; else if (x->score > y->score) return 1; else return 0; } static int scoreorderRev(ALN *x, ALN *y) { if (x->score > y->score) return -1; else if (x->score < y->score) return 1; else return 0; } static int organism_order(ALN *x, ALN *y) { return strcmp(x->organism, y->organism); } static int alphaorder(ALN *x, ALN *y) { int retval=0; if (!(retval = strcmp(x->name, y->name))) { if (x->start == y->start) { if (x->end < y->end) retval = -1; else if (x->end > y->end) retval = 1; } else if (x->start < y->start) retval = -1; else if (x->start > y->start) retval = 1; } /* printf("Comparing %10s %4d %4d with %10s %4d %4d = %d\n", x->name, x->start, x->end, y->name, y->start, y->end, retval); */ return retval; } /* Separate markuplines to another array before resorting */ static void separateMarkupLines(void) { int i, count=0; MarkupAlign = arrayReCreate(MarkupAlign, 100, ALN); if (Highlight_alnp) alncpy(&aln, Highlight_alnp); arrayOrder(); for (i = 0; i+count < nseq; ) { alnp = arrp(Align, i, ALN); if (alnp->markup) { arrayInsert(MarkupAlign, alnp, (void*)alphaorder); arrayRemove(Align, alnp, (void*)nrorder); count++; /* printf ("Moving line %d, nseq=%d\n", i, nseq); */ } else i++; } nseq -= count; arrayOrder(); if (Highlight_alnp) { if (!alignFind(Align, &aln, &ip)) Highlight_alnp = 0; else Highlight_alnp = arrp(Align, ip, ALN); } } /* Reinsert markuplines after mother sequence or if orphan at bottom */ static void reInsertMarkupLines(void) { int i, j; char tmpname[MAXNAMESIZE+1], *cp; for (i = arrayMax(MarkupAlign)-1; i >=0 ; i--) { alnp = arrp(MarkupAlign, i, ALN); strcpy(tmpname, alnp->name); if (cp = strchr(tmpname, ' ')) *cp = 0; arrayOrder10(); for (j = 0; j < arrayMax(Align); j++) if (!strcmp(tmpname, arrp(Align, j, ALN)->name)) break; alnp->nr = (j+1)*10+5; if (0) { printf ("Current alignment:\n"); { int k; ALN *alnk; for (k = 0; k < nseq; k++) { alnk = arrp(Align, k, ALN); printf("%s %d \n", alnk->name, alnk->nr); } } printf ("inserting %s at %d \n\n", alnp->name, alnp->nr); } arrayInsert(Align, alnp, (void*)nrorder); nseq++; } arrayOrder(); } static int organismorder(ALN *x, ALN *y) { int retval=0; char *p1 = strchr(x->name, '_'), *p2 = strchr(y->name, '_'); if (!p1 && !p2) return alphaorder(x, y); if (!p1) return 1; if (!p2) return -1; if (!(retval = strcmp(p1, p2))) return alphaorder(x, y); return retval; } static void showOrganisms(void) { int i; if (!graphActivate(organismsGraph)) { organismsGraph = graphCreate (TEXT_FULL_SCROLL, "Organisms", 0, 0, 0.3, 0.5); graphTextFormat(FIXED_WIDTH); /* graphRegister(PICK, savePick); */ } graphPop(); graphClear(); graphTextBounds(50, arrayMax(organismArr)+ 2); for (i = 0; i < arrayMax(organismArr); i++) { graphColor(arrp(organismArr, i, ALN)->color); graphText(arrp(organismArr, i, ALN)->organism, 1, 1+i); } graphRedraw(); } /* Highlight the box of Highlight_alnp */ static void highlight(int ON) { int i, box; if (!Highlight_alnp) return; box = Highlight_alnp->nr - AlignYstart; if (box > 0 && box <= Alignheig) { /* The alignment */ if (ON) graphBoxDraw(box, WHITE, BLACK); else graphBoxDraw(box, BLACK, WHITE); } /* All names * / for (i = Highlight_matches = 0; i < Alignheig; i++) if (!strcmp(arrp(Align, AlignYstart+i, ALN)->name, Highlight_alnp->name)) { if (ON) { graphBoxDraw(Alignheig+i+1, WHITE, BLACK); Highlight_matches++; } else graphBoxDraw(Alignheig+i+1, BLACK, WHITE); } */ /* All names - also count all matches */ for (i = Highlight_matches = 0; i < nseq; i++) if (!strcmp(arrp(Align, i, ALN)->name, Highlight_alnp->name)) { Highlight_matches++; if (i >= AlignYstart && i < AlignYstart+Alignheig) { if (ON) graphBoxDraw(i-AlignYstart+1+Alignheig, WHITE, BLACK); else graphBoxDraw(i-AlignYstart+1+Alignheig, BLACK, WHITE); } } } /* General purpose routine to convert a string to ALN struct. Note: only fields Name, Start, End are filled! */ static void str2aln(char *src, ALN *alnp) { char *tmp = messalloc(strlen(src)+1), *cp; strcpy(tmp, src); stripCoordTokens(tmp); if (sscanf(tmp, "%s%d%d", alnp->name, &alnp->start, &alnp->end) != 3) { messout("Name to field conversion failed for %s (%s)", src, tmp); return; } messfree(tmp); } /* General purpose routine to convert an ALN struct to a string. */ static char *aln2str(ALN *alnp) { char *cp = messalloc(strlen(alnp->name)+maxStartLen+maxEndLen+10); sprintf(cp, "%s/%d-%d\n", alnp->name, alnp->start, alnp->end); return cp; } /* Action to picking a box in a tree window */ static void treeboxPick (int box, double x, double y, int modifier_unused) { treeNode *treenodePicked; if (!box) return; if (!graphAssFind(assVoid(1), &treestruct)) { messout("Can not find treeStruct"); return; }; if (!graphAssFind(assVoid(100+box), &treenodePicked)) { messout("Picked box %d not associated", box); return; } /* Invert box */ graphBoxDraw(treestruct->currentPickedBox, BLACK, WHITE); graphBoxDraw(box, WHITE, BLACK); treestruct->currentPickedBox = box; if (box > treestruct->lastNodeBox) { /* sequence box */ /* Highlight it in the alignment */ Highlight_alnp = treenodePicked->aln; centerHighlighted(); belvuRedraw(); } else { /* Tree node box - reroot */ if (treenodePicked->parent->parent) { /* Do not reroot with the same root */ treestruct->head = treecpy(treestruct->head); treeDraw(treeReroot(treenodePicked)); } } } /* Click on name -> fetch * Click on sequence -> update status bar * * NOTE: the first Alignheig boxes are the alignment ones, then come Alignheig name ones. */ static void boxPick (int box, double x, double y, int modifier_unused) { static int j, lastbox=0; static Graph ga, lastga, statGraph; static time_t lasttime; pickedCol = (int)x+1+AlignXstart; ga = graphActive(); /* printf("Picked box %d\n", box); */ if (box > Alignheig*2) { /* Click outside alignment - check for scrollbar controls */ if (box == scrLeftBox) scrLeft(); else if (box == scrRightBox) scrRight(); else if (box == scrUpBox) scrUp(); else if (box == scrDownBox) scrDown(); else if (box == scrLeftBigBox) scrLeftBig(); else if (box == scrRightBigBox) scrRightBig(); else if (box == scrUpBigBox) scrUpBig(); else if (box == scrDownBigBox) scrDownBig(); else if (box == HscrollSliderBox) { graphRegister(LEFT_DRAG, HscrollDrag); SliderOffset = x; } else if (box == VscrollSliderBox) { graphRegister(LEFT_DRAG, VscrollDrag); SliderOffset = y; } return; } /* Alignment clicked */ /* Turn last box off */ if (lastga && lastbox) { graphActivate(lastga); } highlight(0); /* Turn all selections off */ if (!box) { Highlight_alnp = 0; return; } aln.nr = (box > Alignheig ? box-Alignheig : box); aln.nr += AlignYstart; if (!arrayFind(Align, &aln, &ip, (void*)nrorder)) { messout("boxPick: Cannot find row %d in alignment array", aln.nr); return; } Highlight_alnp = alnp = arrp(Align, ip, ALN); /* Double click */ if (box == lastbox && (time(0) - lasttime < DBL_CLK_DELAY)) { lasttime = time(0) - DBL_CLK_DELAY; /* To avoid triple clicks */ /* 'Remove many sequences' mode */ if (rmPickLeftOn) { rmPicked(); return; } /* if (box > Alignheig) */ fetchAln(alnp); } else lasttime = time(0); /* Single click */ /* Turn this box on */ graphActivate(ga); highlight(1); lastbox = box; lastga = ga; if (box <= Alignheig) { /* In alignment - update status bar */ updateStatus(alnp, pickedCol); } else updateStatus(alnp, 0); } static void keyboard(int key, int unused) { switch (key) { /* case UP_KEY: scrUp(); break; case DOWN_KEY: scrDown(); break; case LEFT_KEY: scrLeft(); break; case RIGHT_KEY: scrRight(); break; */ case PAGE_UP_KEY: scrUpBig(); break; /* Doesn't work */ case PAGE_DOWN_KEY: scrDownBig(); break; /* Doesn't work */ case UP_KEY: scrUpBig(); break; case DOWN_KEY: scrDownBig(); break; case LEFT_KEY: scrLeftBig(); break; case RIGHT_KEY: scrRightBig(); break; case INSERT_KEY: AlignXstart = 0; belvuRedraw(); break; case DELETE_KEY: AlignXstart = maxLen; belvuRedraw(); break; case HOME_KEY: AlignYstart = 0; belvuRedraw(); break; case END_KEY: AlignYstart = nseq; belvuRedraw(); break; } } static void initConservMtx(void) { int i; conservCount = (int **)messalloc(21*sizeof(int *)); colorMap = (int **)messalloc(21*sizeof(int *)); for (i = 0; i < 21; i++) { conservCount[i] = (int *)messalloc(maxLen*sizeof(int)); colorMap[i] = (int *)messalloc(maxLen*sizeof(int)); } conservResidues = (int *)messalloc(maxLen*sizeof(int)); conservation = (float *)messalloc(maxLen*sizeof(float)); } /* Calculate conservation in each column */ static int countResidueFreqs(void) { int i, j, max, consensus, nseqeff=0; float id; if (!conservCount) initConservMtx(); /* Must reset since this routine may be called many times */ for (i = 0; i < maxLen; i++) { for (j = 0; j < 21; j++) { conservCount[j][i] = 0; colorMap[j][i] = 0; } conservResidues[i] = 0; } for (j = 0; j < nseq; j++) { alnp = arrp(Align, j, ALN); if (alnp->markup) continue; nseqeff++; for (i = 0; i < maxLen; i++) { conservCount[a2b[alnp->seq[i]]][i]++; if (isalpha(alnp->seq[i]) || alnp->seq[i] == '*') conservResidues[i]++; } } return nseqeff; } /* Return 1 if c1 has priority over c2, 0 otherwise */ static int colorPriority(int c1, int c2) { if (c2 == WHITE) return 1; if (c2 == maxbgColor) return 0; if (c2 == lowbgColor) { if (c1 == lowbgColor) return 0; else return 1; } if (c2 == midbgColor) { if (c1 == maxbgColor) return 1; else return 0; } fatal("Unknown colour %s", colorNames[c2]); } static void setConservColors() { int i, j, k, l, colornr, simCount, n, nseqeff; float id, maxid; if (!conservCount) initConservMtx(); nseqeff = countResidueFreqs(); for (i = 0; i < maxLen; i++) for (k = 1; k < 21; k++) colorMap[k][i] = WHITE; for (i = 0; i < maxLen; i++) { maxid = -100.0; for (k = 1; k < 21; k++) { if (color_by_similarity) { /* Convert counts to similarity counts */ simCount = 0; for (j = 1; j < 21; j++) { if (j == k) { if (1) simCount += (conservCount[j][i]-1)*conservCount[k][i]*BLOSUM62[j-1][k-1]; else /* Alternative, less good way */ simCount += (int)floor(conservCount[j][i]/2.0)* (int)ceil(conservCount[k][i]/2.0)* BLOSUM62[j-1][k-1]; } else simCount += conservCount[j][i]*conservCount[k][i]*BLOSUM62[j-1][k-1]; } if (ignoreGapsOn) n = conservResidues[i]; else n = nseqeff; if (n < 2) id = 0.0; else { if (1) id = (float)simCount/(n*(n-1)); else /* Alternative, less good way */ id = (float)simCount/(n/2.0 * n/2.0); } /* printf("%d, %c: simCount= %d, id= %.2f\n", i, b2a[k], simCount, id); */ if (id > lowSimCutoff) { if (id > maxSimCutoff) colornr = maxbgColor; else if (id > midSimCutoff) colornr = midbgColor; else colornr = lowbgColor; if (colorPriority(colornr, colorMap[k][i])) colorMap[k][i] = colornr; /* Colour all similar residues too */ for (l = 1; l < 21; l++) { if (BLOSUM62[k-1][l-1] > 0 && colorPriority(colornr, colorMap[l][i])) { /*printf("%d: %c -> %c\n", i, b2a[k], b2a[l]);*/ colorMap[l][i] = colornr; } } } } else { if (ignoreGapsOn && conservResidues[i] != 1) id = (float)conservCount[k][i]/conservResidues[i]; else id = (float)conservCount[k][i]/nseqeff; if (colorByResIdOn) { if (id*100.0 > colorByResIdCutoff) colorMap[k][i] = color[b2a[k]]; else colorMap[k][i] = WHITE; } else if (id > lowIdCutoff) { if (id > maxIdCutoff) colornr = maxbgColor; else if (id > midIdCutoff) colornr = midbgColor; else colornr = lowbgColor; colorMap[k][i] = colornr; if (id_blosum) { /* Colour all similar residues too */ for (l = 1; l < 21; l++) { if (BLOSUM62[k-1][l-1] > 0 && colorPriority(colornr, colorMap[l][i])) { /*printf("%d: %c -> %c\n", i, b2a[k], b2a[l]);*/ colorMap[l][i] = colornr; } } } } } if (id > maxid) maxid = id; } conservation[i] = maxid; } } static int findResidueBGcolor(ALN* alnp, int i) { if (alnp->nocolor) return boxColor; else if (alnp->markup) return markupColor[alnp->seq[i]]; else if (color_by_conserv || colorByResIdOn) return colorMap[a2b[alnp->seq[i]]][i]; else return color[alnp->seq[i]]; } static void colorCons(void) { setConservColors(); menuSetFlags(menuItem(colorMenu, thresholdStr), MENUFLAG_DISABLED); menuUnsetFlags(menuItem(colorMenu, printColorsStr), MENUFLAG_DISABLED); menuUnsetFlags(menuItem(colorMenu, ignoreGapsStr), MENUFLAG_DISABLED); colorByResIdOn = 0; belvuRedraw(); } static void colorId(void) { colorScheme = COLORID; color_by_conserv = 1; color_by_similarity = 0; id_blosum = 0; colorCons(); } static void colorIdSim(void) { colorScheme = COLORIDSIM; color_by_conserv = 1; color_by_similarity = 0; id_blosum = 1; colorCons(); } static void colorSim(void) { colorScheme = COLORSIM; color_by_conserv = 1; color_by_similarity = 1; colorCons(); } static void graphScrollBar(float x1, float y1, float x2, float y2) { int oldColor = graphColor(LIGHTBLUE); graphFillRectangle(x1, y1, x2, y2); graphColor(BLACK); graphRectangle(x1, y1, x2, y2); graphColor(oldColor); } static void graphTriangle(float x1, float y1, float x2, float y2, float x3, float y3) { Array temp = arrayCreate(6, float) ; int oldColor = graphColor(LIGHTBLUE); array(temp, 0, float) = x1; array(temp, 1, float) = y1; array(temp, 2, float) = x2; array(temp, 3, float) = y2; array(temp, 4, float) = x3; array(temp, 5, float) = y3; graphPolygon(temp); arrayDestroy(temp); graphColor(BLACK); graphLine(x1, y1, x2, y2); graphLine(x2, y2, x3, y3); graphLine(x3, y3, x1, y1); graphColor(oldColor); } static void belvuDestroy() { if (graphActive() == belvuGraph && !saved && messQuery ("Alignment was modified - save ?")) { saveMul(); } graphDestroy(); } static void centerHighlighted(void) { if (Highlight_alnp) AlignYstart = Highlight_alnp->nr - Alignheig/2; } static void belvuRedraw(void) { char seq[MAXLINE+1], *word1end, *c, pads[10], ichar[10], db[5], widBylen[99], ch[2] = " "; int ac, i, j, k, l, box, nx, ny, seqlen, statsStart, bkcol; float Ypos=3, Ysave, x1, x2, y1, y2; if (!graphActivate(belvuGraph)) return; /* Batch usage */ graphClear(); graphTextFormat(FIXED_WIDTH); graphFitBounds(&nx, &ny); /* Development grid lines */ if (gridOn) { graphColor(ORANGE); for (i = 0; i <= ny; i++) graphLine(0, i, nx, i); for (i=0; i <= nx; i++) graphLine(i, 0, i, ny); graphColor(BLACK); } Xpos = maxNameLen+1 + 1; if (1 /* coordsOn */ ) Xpos += maxStartLen+1 + maxEndLen+1; if (displayScores) Xpos += maxScoreLen+1; Alignwid = nx-0.5 - VSCRWID - Xpos; if (Alignwid > maxLen) { Alignwid = maxLen; AlignXstart = 0; } if (AlignXstart < 0) AlignXstart = 0; if (AlignXstart + Alignwid > maxLen) AlignXstart = maxLen - Alignwid; seqlen = (Alignwid > MAXLINE ? MAXLINE : Alignwid); /* RULER *******/ { if (!ruler) { ruler = messalloc(maxLen+1); memset(ruler, '-', maxLen); for (i = 10 ; i <= maxLen; i += 10) { sprintf(ichar, "%d", i); strncpy(ruler + i - strlen(ichar), ichar, strlen(ichar)); } } strncpy(seq, ruler+AlignXstart, seqlen); seq[seqlen] = 0; graphText(seq, Xpos, Ypos); Ypos += 2; } /* Matching sequences ********************/ { /* No, bad idea - have to rewrite all boxPick functions, and horizontal scrolling... */ /* graphLine(0.5, Ypos-0.5, nx-0.5, Ypos-0.5); graphText("Query goes here", 1, Ypos); graphTextFormat (PLAIN_FORMAT); Ypos += 2; */ } Alignheig = ny - 0.5 - HSCRWID - Ypos - 0.5; if (Alignheig > nseq) { Alignheig = nseq; AlignYstart = 0; } if (AlignYstart < 0) AlignYstart = 0; if (AlignYstart + Alignheig > nseq) AlignYstart = nseq - Alignheig; /* Alignment ***********/ Ysave = Ypos; for (j = AlignYstart; j < AlignYstart+Alignheig && j < nseq; j++) { alnp = arrp(Align, j, ALN); if (alnp->hide) { /* Need dummy box to get box-counting right. Ugly but cheap */ graphBoxStart(); graphBoxEnd(); continue; } graphBoxStart(); if (colorRectangles) { for (i = AlignXstart; i < AlignXstart+Alignwid && i < maxLen; i++) { if (!isGap(alnp->seq[i]) && alnp->seq[i] != ' ') { bkcol = findResidueBGcolor(alnp, i); graphColor(bkcol); graphFillRectangle(Xpos+i-AlignXstart, Ypos, Xpos+i-AlignXstart+1, Ypos+1); } else bkcol = NOCOLOR; /* Foreground color */ if (color_by_conserv) { graphColor(bg2fgColor(bkcol)); *ch = alnp->seq[i]; graphText(ch, Xpos+i-AlignXstart, Ypos); } } } graphColor(BLACK); if (!color_by_conserv) { strncpy(seq, alnp->seq+AlignXstart, seqlen); seq[seqlen] = 0; graphText(seq, Xpos, Ypos); } graphBoxEnd(); Ypos++; } graphColor(BLACK); /* Names, coordinates and scores *******************************/ Ypos = Ysave; for (j = AlignYstart; j < AlignYstart+Alignheig && j < nseq; j++) { alnp = arrp(Align, j, ALN); if (alnp->hide) { /* Need dummy box to get box-counting right. Ugly but cheap */ graphBoxStart(); graphBoxEnd(); continue; } box = graphBoxStart(); graphText(alnp->name, 1, Ypos); graphBoxEnd(); graphBoxDraw(box, BLACK, alnp->color); if (alnp->markup != GC /* && coordsOn */ ) { graphText(messprintf("%*d", maxStartLen, alnp->start), maxNameLen+2, Ypos); graphText(messprintf("%*d", maxEndLen, alnp->end), maxNameLen+maxStartLen+3, Ypos); } if (displayScores && !alnp->markup) { graphText(messprintf("%*.1f", maxScoreLen, alnp->score), maxNameLen+maxStartLen+maxEndLen+4, Ypos); } /* Fetch command */ /* Sort out database here since short PIR ones will be seen as * Swissprot otherwise - this overrides Efetch' default * parsing. I do this since normally belvu has complete seq * names - no need to guess them. * / if (strchr(alnp->name, ':')) *db = 0; else if (strchr(alnp->name, '_')) strcpy(db, "SW:"); else if (strchr(alnp->name, '.')) strcpy(db, "WP:"); else strcpy(db, "PIR:"); sprintf(alnp->fetch, "%s%s", db, alnp->name); */ sprintf(alnp->fetch, "%s", alnp->name); Ypos++; } /* Scrollbars ************/ HscrollStart = Xpos - 0.5 + VSCRWID; HscrollEnd = nx - 0.5 - 2*VSCRWID; VscrollStart = Ysave - 0.5 + HSCRWID; VscrollEnd = ny - 0.5 - 2*HSCRWID; x1 = HscrollStart + ((float)AlignXstart/maxLen)*(HscrollEnd-HscrollStart); x2 = HscrollStart + ((float)(AlignXstart+Alignwid)/maxLen)*(HscrollEnd-HscrollStart); HsliderLen = x2-x1; y1 = VscrollStart + ((float)AlignYstart/nseq)*(VscrollEnd-VscrollStart); y2 = VscrollStart + ((float)(AlignYstart+Alignheig)/nseq)*(VscrollEnd-VscrollStart); VsliderLen = y2-y1; /* Horizontal scrollbar */ /* arrows */ scrLeftBox = graphBoxStart(); graphTriangle(HscrollStart-VSCRWID, ny-0.5-0.5*HSCRWID, HscrollStart, ny-0.5-HSCRWID, HscrollStart, ny-0.5); graphBoxEnd(); graphBoxDraw(scrLeftBox, BLACK, SCRBACKCOLOR); scrRightBox = graphBoxStart(); graphTriangle(HscrollEnd+VSCRWID, ny-0.5-0.5*HSCRWID, HscrollEnd, ny-0.5-HSCRWID, HscrollEnd, ny-0.5); graphBoxEnd(); graphBoxDraw(scrRightBox, BLACK, SCRBACKCOLOR); /* big-scroll boxes */ scrLeftBigBox = graphBoxStart(); graphRectangle(HscrollStart, ny-0.5-HSCRWID, x1, ny-0.5); graphBoxEnd(); graphBoxDraw(scrLeftBigBox, BLACK, SCRBACKCOLOR); scrRightBigBox = graphBoxStart(); graphRectangle(x2, ny-0.5-HSCRWID, HscrollEnd, ny-0.5); graphBoxEnd(); graphBoxDraw(scrRightBigBox, BLACK, SCRBACKCOLOR); /* slider */ HscrollSliderBox = graphBoxStart(); graphScrollBar(x1, ny-0.5-HSCRWID, x2, ny-0.5); graphBoxEnd(); /* Vertical scrollbar */ /* arrows */ scrUpBox = graphBoxStart(); graphTriangle(nx-0.5-0.5*VSCRWID, VscrollStart-HSCRWID, nx-0.5-VSCRWID, VscrollStart, nx-0.5, VscrollStart); graphBoxEnd(); graphBoxDraw(scrUpBox, BLACK, SCRBACKCOLOR); scrDownBox = graphBoxStart(); graphTriangle(nx-0.5-0.5*VSCRWID, VscrollEnd+HSCRWID, nx-0.5-VSCRWID, VscrollEnd, nx-0.5, VscrollEnd); graphBoxEnd(); graphBoxDraw(scrDownBox, BLACK, SCRBACKCOLOR); /* big-scroll boxes */ scrUpBigBox = graphBoxStart(); graphRectangle(nx-0.5-VSCRWID, VscrollStart, nx-0.5, y1); graphBoxEnd(); graphBoxDraw(scrUpBigBox, BLACK, SCRBACKCOLOR); scrDownBigBox = graphBoxStart(); graphRectangle(nx-0.5-VSCRWID, y2, nx-0.5, VscrollEnd); graphBoxEnd(); graphBoxDraw(scrDownBigBox, BLACK, SCRBACKCOLOR); /* slider */ VscrollSliderBox = graphBoxStart(); graphScrollBar(nx-0.5-VSCRWID, y1, nx-0.5, y2); graphBoxEnd(); /* patch up frame of scroll boxes */ graphRectangle(HscrollStart-VSCRWID, ny-0.5-HSCRWID, HscrollEnd+VSCRWID, ny-0.5); graphRectangle(nx-0.5-VSCRWID, VscrollStart-HSCRWID, nx-0.5, VscrollEnd+HSCRWID); /* Buttons *********/ statsStart = 1; graphButton("File", Help, statsStart, 0.9); /* Uses main menu, haha */ statsStart += 6; editButtonBox = graphButton("Edit", rmPickLeft, statsStart, 0.9); statsStart += 6; colorButtonBox = graphButton("Colour", showColorCodesRedraw, statsStart, 0.9); statsStart += 8; sortButtonBox = graphButton("Sort", alphaSort, statsStart, 0.9); statsStart += 6; /* Status bar ************/ sprintf(widBylen, "(%dx%d)", nseq, maxLen); graphText(widBylen, 1, 3); /* statsStart += strlen(widBylen)+1; */ graphText("Picked:", statsStart, 1); statsStart += 8; statsBox = graphBoxStart(); graphTextPtr(stats, statsStart, 1, nx-statsStart-1); graphBoxEnd(); graphBoxDraw(statsBox, BLACK, LIGHTBLUE); /* Frame lines *************/ graphLine(0.5, 0.5, nx-0.5, 0.5); graphLine(0.5, 0.5, 0.5, ny-0.5); graphLine(nx-0.5, 0.5, nx-0.5, ny-0.5); graphLine(0.5, ny-0.5, nx-0.5, ny-0.5); graphLine(0.5, 2.5, nx-0.5, 2.5); graphLine(0.5, Ysave-0.5, nx-0.5, Ysave-0.5); Xpos = 0.5 + maxNameLen + 1; graphLine(Xpos, Ysave-0.5, Xpos, ny-0.5); Xpos += maxStartLen + 1; graphLine(Xpos, Ysave-0.5, Xpos, ny-0.5); Xpos += maxEndLen + 1; graphLine(Xpos, Ysave-0.5, Xpos, ny-0.5); if (displayScores) { Xpos += maxScoreLen + 1; graphLine(Xpos, Ysave-0.5, Xpos, ny-0.5); } /* Highlighted sequence **********************/ highlight(1); setMenuCheckmarks(); graphRedraw() ; } static int isGap(char c) { if (c == '.' || c == '-' || c == '[' || c == ']' /* Collapse-control chars */ ) return 1; return 0; } static int isAlign(char c) { if (isalpha(c) || isGap(c) || c == '*') return 1; else return 0; } static int getMatchStates(void) { int i, j, n, retval=0; for (j = 0; j < maxLen; j++) { n = 0; for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (isalpha(alnp->seq[j])) n++; } if (n > nseq/2) retval++; } return retval; } /* Output a.a. probabilities in HMMER format. Disabled counting of match-state residues only (script HMM-random does that already) Disabled pseudocounts - not sure what the best way is. */ static void outputProbs(FILE *fil) { /* From swissprot 33: Ala (A) 7.54 Gln (Q) 4.02 Leu (L) 9.31 Ser (S) 7.19 Arg (R) 5.15 Glu (E) 6.31 Lys (K) 5.94 Thr (T) 5.76 Asn (N) 4.54 Gly (G) 6.86 Met (M) 2.36 Trp (W) 1.26 Asp (D) 5.29 His (H) 2.23 Phe (F) 4.06 Tyr (Y) 3.21 Cys (C) 1.70 Ile (I) 5.72 Pro (P) 4.91 Val (V) 6.52 */ float f[21] = { 0.0, 0.0754, 0.0170, 0.0529, 0.0631, 0.0406, 0.0686, 0.0223, 0.0572, 0.0594, 0.0931, 0.0236, 0.0454, 0.0491, 0.0402, 0.0515, 0.0719, 0.0576, 0.0652, 0.0126, 0.0321}; float f_sean[21] = { 0.0, .08713, .03347, .04687, .04953, .03977, .08861, .03362, .03689, .08048, .08536, .01475, .04043, .05068, .03826, .04090, .06958, .05854, .06472, .01049, .02992 }; float p; int i, j, c[21], n=0, nmat; char *cp; for (i = 1; i <= 20; i++) c[i] = 0; for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); cp = alnp->seq; for (; *cp; cp++) { if (a2b_sean[*cp]) { c[a2b_sean[*cp]]++; n++; } } } if (!n) fatal("No residues found"); if (0) { nmat = getMatchStates(); /* Approximation, HMM may differ slightly */ if (!nmat) fatal("No match state columns found"); printf("Amino\n"); for (i = 1; i <= 20; i++) { /* One way: p = ((float)c[i]/nmat + 20*f_sean[i]) / (float)(n/nmat+20);*/ /* Other way: */ p = ((float)c[i] + 20*f_sean[i]) / (float)(n+20); if (p < 0.000001) p = 0.000001; /* Must not be zero */ printf("%f ", p); /* printf(" %d %f %d \n", c[i], f[i], n);*/ } } else { printf("Amino\n"); for (i = 1; i <= 20; i++) { p = (float)c[i] / (float)n; if (p < 0.000001) p = 0.000001; /* Must not be zero */ printf("%f ", p); } } printf("\n"); fclose(fil); fflush(fil); saved = 1; } /* Format the name/start-end string For convenience, used by writeMul */ char *writeMulName(ALN *aln) { static char name[MAXNAMESIZE+50]; char *cp, *namep, GRname[MAXNAMESIZE*2+2], GRfeat[MAXNAMESIZE+1]; if (aln->markup == GC) return messprintf("#=GC %s", aln->name); /* NOTE: GR lines have the feature concatenated in aln->name - must separate */ if (aln->markup == GR) { namep = GRname; strncpy(namep, aln->name, 50); cp = strchr(namep, ' '); strncpy(GRfeat, cp+1, 50); *cp = 0; } else namep = aln->name; if (!saveCoordsOn) { strcpy(name, messprintf("%s", namep)); } else { strcpy(name, messprintf("%s%c%d-%d", namep, saveSeparator, aln->start, aln->end)); } if (aln->markup == GR) return messprintf("#=GR %s %s", name, GRfeat); else return name; } static void writeMul(FILE *fil) { static char *cp; int i, w, W; /* Write Annotation */ if (!stackEmpty(AnnStack)) { stackCursor(AnnStack, 0); while (cp = stackNextText(AnnStack)) { fprintf(fil, "%s\n", cp); } } /* Find max width of name column */ for (i = w = 0; i < nseq; i++) if ( (W = strlen(writeMulName(arrp(Align, i, ALN)))) > w) w = W; /* Write alignment */ for (i = 0; i < nseq; i++) fprintf(fil, "%-*s %s\n", w, writeMulName(arrp(Align, i, ALN)), arrp(Align, i, ALN)->seq); fprintf(fil, "//\n"); fclose(fil); fflush(fil); saved = 1; } static void saveMul(void) { FILE *fil; if (!(fil = filqueryopen(dirName, fileName, "","w", "Save as Stockholm file:" ))) return; writeMul(fil); } static void writeFasta(FILE *pipe) { int i, n; char *cp; for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (alnp->markup) continue; if (saveCoordsOn) fprintf(pipe, ">%s%c%d-%d\n", alnp->name, saveSeparator, alnp->start, alnp->end); else fprintf(pipe, ">%s\n", alnp->name); for (n=0, cp = alnp->seq; *cp; cp++) { if (!strcmp(saveFormat, FastaAlnStr)) { fputc(*cp, pipe); n++; } else { if (!isGap(*cp)) { fputc(*cp, pipe); n++; } } if (n && !(n % 80) ) { fputc('\n', pipe); n = 0; } } if (n) fputc('\n', pipe); } fclose(pipe); fflush(pipe); } static void saveFasta(void) { FILE *fil; if (!(fil = filqueryopen(dirName, fileName, "","w", "Save as unaligned Fasta file:"))) return; writeFasta(fil); saved = 1; } static void saveAlignment(void) { if (!strcmp(saveFormat, MSFStr)) saveMSF(); else if (!strcmp(saveFormat, FastaAlnStr)) saveFasta(); else if (!strcmp(saveFormat, FastaStr)) saveFasta(); else saveMul(); } static void graphButtonCheck(char* text, void (*func)(void), float x, float *y, int On) { graphButton(messprintf("%s %s", On ? "*" : " ", text), func, x, *y); *y += 2; } static void saveMSFSelect(void) { strcpy(saveFormat, MSFStr); saveRedraw(); } static void saveMulSelect(void){ strcpy(saveFormat, MulStr); saveRedraw(); } static void saveFastaAlnSelect(void){ strcpy(saveFormat, FastaAlnStr); saveRedraw(); } static void saveFastaSelect(void){ /* unaligned */ strcpy(saveFormat, FastaStr); saveRedraw(); } static void treeUPGMAselect(void){ strcpy(treeMethodString, UPGMAstr); treeMethod = UPGMA; treeSettings(); } static void treeNJselect(void){ strcpy(treeMethodString, NJstr); treeMethod = NJ; treeSettings(); } static void treeUNCORRselect(void){ strcpy(treeDistString, UNCORRstr); treeDistCorr = UNCORR; treeScale = 1.0; treeSettings(); } static void treeKIMURAselect(void){ strcpy(treeDistString, KIMURAstr); treeDistCorr = KIMURA; treeSettings(); } static void treeSTORMSONNselect(void){ strcpy(treeDistString, STORMSONNstr); treeDistCorr = STORMSONN; treeScale = 0.3; treeSettings(); } static void saveCoordsToggle(void) { saveCoordsOn = !saveCoordsOn; saveRedraw(); } static void saveSeparatorToggle(void) { if (saveSeparator == '/') saveSeparator = '='; else saveSeparator = '/'; saveRedraw(); } static void saveOK(void) { saveAlignment(); graphDestroy(); graphActivate(belvuGraph); } static void saveRedraw(void) { float x = 1.0, y = 1.0; if (!graphActivate(saveGraph)) { saveGraph = graphCreate (TEXT_FIT, "Save As", 0, 0, 0.5, 0.25); graphTextFormat(FIXED_WIDTH); /* graphRegister(PICK, savePick); */ } graphPop(); graphClear(); graphText("Format:", x, y); graphNewBoxMenu(graphButton(saveFormat, saveRedraw, x+9, y), saveMenu); y += 2.0; graphButtonCheck("Save with coordinates", saveCoordsToggle, x, &y, saveCoordsOn); if (saveCoordsOn) { graphText("Separator char between name and coordinates:", x, y); graphButton(messprintf(" %c ", saveSeparator), saveSeparatorToggle, x+45, y); y += 1.0; graphText("(Use = for GCG)", x, y); y += 1.0; } y += 2; graphButton("OK", saveOK, x, y); graphButton("Cancel", graphDestroy, x+4, y); graphRedraw(); } static void treeShowBranchlenToggle(void) { treeShowBranchlen = !treeShowBranchlen; treeSettings(); } static void treeShowOrganismToggle(void) { treeShowOrganism = !treeShowOrganism; treeSettings(); } static void treeSettings(void) { float x = 1.0, y = 1.0; static char linewidthstr[256], treeScalestr[256]; if (!graphActivate(treeGUIGraph)) { treeGUIGraph = graphCreate (TEXT_FIT, "tree GUI", 0, 0, 0.5, 0.35); graphTextFormat(FIXED_WIDTH); /* graphRegister(PICK, savePick); */ } graphPop(); graphClear(); graphText("Select tree building method:", x, y); graphNewBoxMenu(graphButton(treeMethodString, treeSettings, x+30, y), treeGUIMenu); y += 2.0; graphText("Select distance correction method:", x, y); graphNewBoxMenu(graphButton(treeDistString, treeSettings, x+35, y), treeDistMenu); y += 2.0; sprintf(treeScalestr, "%.2f", treeScale); graphText("Tree Scale:", x, y); graphTextEntry (treeScalestr, 5, x+12, y, (void*)(treeSetScale)); y += 2.0; sprintf(linewidthstr, "%.2f", treeLinewidth); graphText("Line width:", x, y); graphTextEntry (linewidthstr, 5, x+12, y, (void*)(treeSetLinewidth)); y += 2.0; graphButtonCheck("Display branch lenghts", treeShowBranchlenToggle, x, &y, treeShowBranchlen); graphButtonCheck("Display organism", treeShowOrganismToggle, x, &y, treeShowOrganism); /*graphButtonCheck("Sample trees in intervals", treeSampleToggle, x, &y, treeBootstrapOn); if (treeBootstrapOn) { graphText("Number of bootstrapsSeparator char between name and coordinates:", x, y); graphButton(messprintf(" %c ", saveSeparator), saveSeparatorToggle, x+45, y); y += 1.0; graphText("(Use = for GCG)", x, y); y += 1.0; }*/ y += 2; graphButton("Make tree", treeDisplay, x, y); y += 2; graphButton("Quit", graphDestroy, x, y); graphRedraw(); } static void alncpy(ALN *dest, ALN *src) { strncpy(dest->name, src->name, MAXNAMESIZE); dest->start = src->start; dest->end = src->end; dest->len = src->len; /* dest->seq = messalloc(strlen(src->seq)); strcpy(dest->seq, src->seq); */ dest->seq = src->seq; dest->nr = src->nr; strncpy(dest->fetch, src->fetch, MAXNAMESIZE+10); dest->score = src->score; dest->color = src->color; } static void resetALN(ALN *alnp) { *alnp->name = *alnp->fetch = 0; alnp->start = alnp->end = alnp->len = alnp->nr = 0; alnp->seq = 0; alnp->score = 0.0; alnp->color = WHITE; alnp->markup = 0; } /* Set the ->nr field in Align array in ascending order */ static void arrayOrder(void) { int i; for (i = 0; i < nseq; i++) arrp(Align, i, ALN)->nr = i+1; } /* Set the ->nr field in Align array in ascending order with increments of 10 Necessary if one wants to insert items. */ static void arrayOrder10(void) { int i; for (i = 0; i < nseq; i++) arrp(Align, i, ALN)->nr = (i+1)*10; } static void init_sort_do(void) { arraySort(Align, (void*)nrorder); if (init_sort) { switch(init_sort) { case SORT_ALPHA : arraySort(Align, (void*)alphaorder); break; case SORT_ORGANISM : arraySort(Align, (void*)organismorder); break; case SORT_SCORE : scoreSortBatch(); break; case SORT_SIM : Highlight_alnp = arrp(Align, 0, ALN); highlightScoreSort('P'); break; case SORT_ID : Highlight_alnp = arrp(Align, 0, ALN); highlightScoreSort('I'); break; case SORT_TREE : treeSortBatch(); break; } } } /* Convert plot coordinates to screen coordinates */ static float xTran(float x) { return(5 + x*conservationXScale); } static float yTran(float y) { return(3 + conservationRange - y*conservationYScale); } static void conservationSetWindow(void) { ACEIN cons_in; if (!(cons_in = messPrompt("Window size:", messprintf("%d", conservationWindow), "iz", 0))) return; aceInInt(cons_in, &conservationWindow); aceInDestroy(cons_in); conservationPlot(); } static void conservationSetLinewidth(void) { ACEIN cons_in; if (!(cons_in = messPrompt("Line width:", messprintf("%.2f", conservationLinewidth), "fz", 0))) return; aceInFloat(cons_in, &conservationLinewidth); aceInDestroy(cons_in); conservationPlot(); } static void conservationSetScale(void) { if (!(ace_in = messPrompt("X and Y scales:", messprintf("%.1f %.1f", conservationXScale, conservationYScale), "ffz", 0))) return; aceInFloat(ace_in, &conservationXScale); aceInFloat(ace_in, &conservationYScale); aceInDestroy(ace_in); conservationPlot(); } static void conservationPlot() { int i; float oldlinew, scalex, scaley, mincons, maxcons, avgcons, range, XscaleBase, YscaleBase, *smooth, sum, f; /* init window */ if (!graphActivate(conservationGraph)) { conservationGraph = graphCreate (TEXT_FULL_SCROLL, "Belvu conservation profile", 0, 0, 1, 0.4); graphTextFormat(FIXED_WIDTH); } oldlinew = graphLinewidth(conservationLinewidth); graphClear(); graphMenu(conservationMenu); /* smooth conservation by applying a window */ smooth = messalloc(maxLen*sizeof(float)); sum = 0.0; for (i = 0; i < conservationWindow; i++) sum += conservation[i]; smooth[conservationWindow/2] = sum/conservationWindow; for ( ; i < maxLen; i++) { sum -= conservation[i-conservationWindow]; sum += conservation[i]; smooth[i-conservationWindow/2] = sum/conservationWindow; } /* Find max and min and avg conservation */ maxcons = -1; mincons = 10000; avgcons = 0; for (i = 0; i < maxLen; i++) { if (smooth[i] > maxcons) maxcons = smooth[i]; if (smooth[i] < mincons) mincons = smooth[i]; avgcons += conservation[i]; } avgcons /= maxLen*1.0; conservationRange = (maxcons - mincons)*conservationYScale; graphTextBounds(maxLen*conservationXScale+30, conservationRange*conservationYScale+10); graphText(messprintf("Window = %d", conservationWindow), 10, 1); /* Draw x scale */ YscaleBase = -1; graphLine(xTran(0), yTran(YscaleBase), xTran(maxLen), yTran(YscaleBase)); for (i=0; i < maxLen; i += 10) { graphLine(xTran(i), yTran(YscaleBase), xTran(i), yTran(YscaleBase)+0.5); if (i+5 < maxLen-1) graphLine(xTran(i+5), yTran(YscaleBase), xTran(i+5), yTran(YscaleBase)+0.25); graphText(messprintf("%d", i), xTran(i)-0.5, yTran(YscaleBase)+1); } /* Draw y scale */ XscaleBase = -1; graphLine(xTran(XscaleBase), yTran(0), xTran(XscaleBase), yTran(maxcons)); for (f=mincons; f < maxcons; f += 1) { graphLine(xTran(XscaleBase), yTran(f), xTran(XscaleBase-0.5), yTran(f)); if (f+0.5 < maxcons) graphLine(xTran(XscaleBase), yTran(f+0.5), xTran(XscaleBase-0.25), yTran(f+0.5)); graphText(messprintf("%.0f", f), xTran(XscaleBase)-3, yTran(f)-0.5); } /* Draw average line */ graphColor(RED); graphLine(xTran(0), yTran(avgcons), xTran(maxLen), yTran(avgcons)); graphText("Average conservation", xTran(maxLen), yTran(avgcons)); graphColor(BLACK); /* Plot conservation */ for (i=1; i 15) color = treeColors[15];*/ arrp(organismArr, i, ALN)->color = color; } } /* Convert swissprot name suffixes to organisms */ static void suffix2organism(void) { int i; char *cp; for (i=0; imarkup && (cp = strchr(alnp->name, '_'))) { char *suffix = messalloc(strlen(cp)+1); strcpy(suffix, cp+1); /* Add organism to table of organisms. This is necessary to make all sequences of the same organism point to the same place and to make a non-redundant list of present organisms */ alnp->organism = suffix; arrayInsert(organismArr, alnp, (void*)organism_order); /* Store pointer to unique organism in ALN struct */ arrayFind(organismArr, alnp, &ip, (void*)organism_order); alnp->organism = arrp(organismArr, ip, ALN)->organism; } } } static void organismSort(void) { if (Highlight_alnp) alncpy(&aln, Highlight_alnp); arraySort(Align, (void*)organismorder); if (Highlight_alnp) { if (!alignFind(Align, &aln, &ip)) { messout("Cannot find back highlighted seq after sort - probably a bug"); Highlight_alnp = 0; } else Highlight_alnp = arrp(Align, ip, ALN); } arrayOrder(); centerHighlighted(); belvuRedraw(); } static void showAnnotation(void) { int maxwidth=0, nlines=0, i=0; char *cp; if (stackEmpty(AnnStack)) return; stackCursor(AnnStack, 0); while (cp = stackNextText(AnnStack)) { nlines++; if (strlen(cp) > maxwidth) maxwidth = strlen(cp); } if (!graphActivate(annGraph)) { annGraph = graphCreate (TEXT_FULL_SCROLL, "Annotations", 0, 0, (maxwidth+1)/fontwidth*screenWidth, (nlines+1)/fontheight*screenHeight); graphTextFormat(FIXED_WIDTH); } graphClear(); graphTextBounds(maxwidth+1, nlines+1); stackCursor(AnnStack, 0); while (cp = stackNextText(AnnStack)) { graphText(cp, 0, i++); } graphRedraw(); } /* Traverse tree to calculate max distance, needed to set scale */ static void setTreeScale(treeNode *node) { return; } /* The actual tree drawing routine. Note: must be in sync with treeDrawNodeBox, which draws clickable boxes first. */ static float treeDrawNode(treeStruct *tree, treeNode *node, float x) { float y, yl, yr; int leftcolor; if (!node) return 0.0; yl = treeDrawNode(tree, node->left, x + node->branchlen*treeScale); leftcolor = graphColor(BLACK); yr = treeDrawNode(tree, node->right, x + node->branchlen*treeScale); if (yl) { /* internal node */ y = (yl + yr) / 2.0; /* connect children */ graphLine(x + node->branchlen*treeScale, yr, x + node->branchlen*treeScale, y); graphColor(leftcolor); graphLine(x + node->branchlen*treeScale, yl, x + node->branchlen*treeScale, y); if (node->left->organism != node->right->organism) graphColor(BLACK); } else { /* Sequence name */ int box, i; y = tree_y++; if (treeColorsOn && node->organism) { aln.organism = node->organism; if (arrayFind(organismArr, &aln, &ip, (void*)organism_order)) { graphColor(arrp(organismArr, ip, ALN)->color); } } else graphColor(BLACK); /* Make clickable box for sequence */ box = graphBoxStart(); graphText(node->name, x + node->branchlen*treeScale + 1, y-0.5); graphBoxEnd(); graphAssociate(assVoid(100+box), node); if (Highlight_alnp && node->aln == Highlight_alnp) { graphBoxDraw(box, WHITE, BLACK); tree->currentPickedBox = box; } else if (node->aln) { graphBoxDraw(box, BLACK, node->aln->color); } if (treeShowOrganism && node->organism) graphText(node->organism, x + node->branchlen*treeScale + 2 + strlen(node->name), y-0.5); { int pos = x + node->branchlen*treeScale + strlen(node->name); if (pos > maxTreeWidth) maxTreeWidth = pos; } } /* Horizontal branches */ graphLine(x + node->branchlen*treeScale, y, x, y); if (treeShowBranchlen && node->branchlen) graphText(messprintf("%.1f", node->branchlen*treeScale), x+.5, y); return y; } /* Draws clickable boxes for horizontal lines. The routine must be in sync with treeDrawNode, but can not be integrated since a box may accidentally overwrite some text. Therefore all boxes must be drawn before any text or lines. */ static float treeDrawNodeBox(treeStruct *tree, treeNode *node, float x) { float y, yl, yr; int box; if (!node) return 0.0; yl = treeDrawNodeBox(tree, node->left, x + node->branchlen*treeScale); yr = treeDrawNodeBox(tree, node->right, x + node->branchlen*treeScale); if (yl) { y = (yl + yr) / 2.0; } else { y = tree_y++; } /* Make box around horizontal lines */ box = graphBoxStart(); graphLine(x + node->branchlen*treeScale, y, x, y); oldlinew = graphLinewidth(0.0); graphColor(BG); graphRectangle(x + node->branchlen*treeScale, y-0.5, x, y+0.5); graphColor(BLACK); graphLinewidth(oldlinew); graphBoxEnd(); graphAssociate(assVoid(100+box), node); node->box = box; tree->lastNodeBox = box; return y; } static void treeDraw(treeNode *treeHead) { int i; float oldlinew; treeStruct *treestruct = messalloc(sizeof(treeStruct)); treeGraph = graphCreate (TEXT_FULL_SCROLL, messprintf("Belvu %s of %s", treeMethod == NJ ? "Neighbor-joining tree (unrooted)" : "UPGMA tree", Title), 0, 0, (treeMethod == UPGMA ? 130 : 110)/fontwidth*screenWidth, (nseq+7)/fontheight*screenHeight); graphTextFormat(FIXED_WIDTH); graphRegister(PICK, treeboxPick); treestruct->head = treeHead; graphAssociate(assVoid(1), treestruct); oldlinew = graphLinewidth(treeLinewidth); graphClear(); graphMenu(treeMenu); /* setTreeScale(treestruct->head); */ maxTreeWidth = 0; tree_y = 1; treeDrawNodeBox(treestruct, treestruct->head, 1.0); tree_y = 1; treeDrawNode(treestruct, treestruct->head, 1.0); graphColor(BLACK); graphTextBounds(maxTreeWidth+2 + (treeShowOrganism ? 25 : 0), nseq+6); /* Draw scale */ if (treeMethod == UPGMA) { tree_y += 2; graphLine(1*treeScale, tree_y, 101*treeScale, tree_y); for (i=1; i<=101; i += 10) { graphLine(i*treeScale, tree_y, i*treeScale, tree_y+0.5); if (i < 101) graphLine((i+5)*treeScale, tree_y, (i+5)*treeScale, tree_y+0.5); graphText(messprintf("%d", i-1), (i-0.5)*treeScale, tree_y+1); } } else { tree_y += 3; graphText("0.1", 5*treeScale, tree_y-1); graphLine(1*treeScale, tree_y, 11*treeScale, tree_y); graphLine(1*treeScale, tree_y-0.5, 1*treeScale, tree_y+0.5); graphLine(11*treeScale, tree_y-0.5, 11*treeScale, tree_y+0.5); } graphLinewidth(oldlinew); if (treeMethod == NJ) { float bal, lweight, rweight; treeNode *tree = treestruct->head; lweight = treeSize3way(tree->left, tree); rweight = treeSize3way(tree->right, tree); graphText((debug ? messprintf("Tree balance = %.1f (%.1f-%.1f)", fabsf(lweight - rweight), lweight, rweight) : messprintf("Tree balance = %.1f", fabsf(lweight - rweight))), 14, tree_y-0.5); } graphRedraw(); } static void treeSetLinewidth(char *cp) { treeLinewidth = atof(cp); treeSettings(); } static void treeSetScale(char *cp) { treeScale = atof(cp); treeSettings(); } /* Note: this treecpy does not reallocate any strings, assuming they will always exist and never change. Serious problems await if this is not true */ static treeNode *treecpy(treeNode *node) { treeNode *newnode; if (!node) return 0; newnode = messalloc(sizeof(treeNode)); newnode->dist = node->dist; newnode->branchlen = node->branchlen; newnode->boot = node->boot; newnode->name = node->name; newnode->organism = node->organism; newnode->aln = node->aln; newnode->box = node->box; newnode->color = node->color; newnode->left = treecpy(node->left); newnode->right = treecpy(node->right); return newnode; } static void fillParents(treeNode *parent, treeNode *node) { if (!node) return; node->parent = parent; fillParents(node, node->left); fillParents(node, node->right); } static char *fillOrganism(treeNode *node) { char *leftorganism, *rightorganism; if (node->name) return node->organism; leftorganism = fillOrganism(node->left); rightorganism = fillOrganism(node->right); /* printf("\nFill: (left=%s, right=%s): ", leftorganism, rightorganism); */ node->organism = (leftorganism == rightorganism ? leftorganism : 0); return node->organism; } static treeNode *treeParent2leaf(treeNode *newparent, treeNode *curr) { if (!curr->parent) { /* i.e. the old root */ if (curr->left == newparent) { newparent->branchlen += curr->right->branchlen; /* Add second part of this vertex */ return curr->right; } else { newparent->branchlen += curr->left->branchlen; return curr->left; } } else { if (curr->left == newparent) { /* Change the link to the new parent to the old parent */ curr->left = treeParent2leaf(curr, curr->parent); curr->left->branchlen = curr->branchlen; } else { curr->right = treeParent2leaf(curr, curr->parent); curr->right->branchlen = curr->branchlen; } } return curr; } /* Rerooting works roughly this way: - A new node is created with one child being the node chosen as new root and the other child the previous parent of it. The sum of left and right branch lengths of the new root should equal the branch length of the chosen node. - All previous parent nodes are visited and converted so that: 1. The previous parent node becomes a child. 2. The new parent is the node that calls. 3. The branchlength of the previous parent node is assigned to its parent. - When the previous root is reached, it is deleted and the other child becomes the child of the calling node. Note that treeReroot destroys the old tree, since it reuses the nodes. Use treecpy if you still need it for something later. */ static treeNode *treeReroot(treeNode *node) { treeNode *newroot = messalloc(sizeof(treeNode)); newroot->left = node; newroot->right = treeParent2leaf(node, node->parent); fillParents(newroot, newroot->left); fillParents(newroot, newroot->right); newroot->left->branchlen = newroot->right->branchlen = node->branchlen / 2.0; treeBalanceByWeight(newroot->left, newroot->right, &newroot->left->branchlen, &newroot->right->branchlen); return newroot; } static void treeTraverse(treeNode *node, void (*func)()) { if (!node) return; treeTraverse(node->left, func); func(node); treeTraverse(node->right, func); } /* Sum branchlengths, allow going up parents If going up parents, arg1 = parent ; arg2 = child. If you're not interested in going up parents, simply call me with the same node in both arguments. */ static float treeSize3way(treeNode *node, treeNode *fromnode) { int root = 0; treeNode *left, *right; float len; if (!node) return 0.0; /* Get the correct branch length */ if (node->left == fromnode || node->right == fromnode) /* Going up the tree */ len = fromnode->branchlen; else len = node->branchlen; if (node->left == fromnode) { left = node->parent; if (!left) root = 1; } else { left = node->left; } if (node->right == fromnode) { right = node->parent; if (!right) root = 1; } else { right = node->right; } if (root) { float retval; /* go across tree root */ if (left) /* Coming from right */ retval = treeSize3way(left, left); else /* Coming from left */ retval = treeSize3way(right, right); if (debug) printf("Returning (%.1f + %.1f) = %.1f\n", fromnode->branchlen, retval, fromnode->branchlen + retval); return fromnode->branchlen + retval; } else { float l = treeSize3way(left, node), r = treeSize3way(right, node), retval = (l + r)/2.0 + len; if (debug) printf("Returning (%.1f + %.1f)/2 + %.1f = %.1f\n", l, r, len, retval); return retval; } } /* Calculate the difference between left and right trees if the tree were to be rerooted at this node. What is calculated (bal) is the difference between 'left' and 'right' subtrees minus the length of the branch itself. If this difference is negative, a perfectly balanced tree can be made. For imperfectly balanced trees we want to root at the branch that gives the best balance. However, perfectly balanced trees are all 'perfect' so here we chose the branch with most equal subtrees. Actually it is not "left" and "right" but "down" and "up" subtrees. */ static void treeCalcBalance(treeNode *node) { float bal, lweight, rweight; if (node == treeBestBalancedNode) return; if (debug) printf("Left/Downstream weight\n"); lweight = treeSize3way(node, node); if (debug) printf("Right/Upstream weight\n"); rweight = treeSize3way(node->parent, node); bal = fabsf(lweight - rweight) - node->branchlen; if (debug) printf("Node=%s (branchlen = %.1f). Weights = %.1f %.1f. Bal = %.1f\n", node->name, node->branchlen, lweight, rweight, bal); if (bal < treeBestBalance) { /* better balance */ if (treeBestBalance > 0.0 || /* If previous tree was not perfectly balanced, or If previous tree was perfectly balanced - choose root with best subtree balance */ fabsf(lweight - rweight) < treeBestBalance_subtrees) { if (debug) printf(" %s has better balance %.1f < %.1f\n", node->name, bal, treeBestBalance); treeBestBalancedNode = node; treeBestBalance = bal; treeBestBalance_subtrees = fabsf(lweight - rweight); } } } /* Find the node which has most equal balance, return new tree with this as root. */ static treeNode *treeFindBalance(treeNode *tree) { float bal, lweight, rweight; treeNode *node = tree; lweight = treeSize3way(tree->left, tree->left); rweight = treeSize3way(tree->right, tree->right); treeBestBalancedNode = tree; treeBestBalance = fabsf(lweight - rweight); treeBestBalance_subtrees = fabsf((lweight - tree->left->branchlen) - (rweight - tree->right->branchlen)); if (debug) printf("Initial weights = %.1f %.1f. Bal = %.1f\n", lweight, rweight, treeBestBalance); treeTraverse(tree, treeCalcBalance); if (treeBestBalancedNode == tree) return tree; else return treeReroot(treeBestBalancedNode); } /* Balance the two sides of a tree branch by the weights of the two subtrees. The branch has two sides: llen and rlen. Rationale: The correct center point balances the treesizes on both sides. Method: Calculate half the unbalance, add it to the appropriate side. No real theory for this, but it seems to work in easy cases */ static void treeBalanceByWeight(treeNode *lnode, treeNode *rnode, float *llen, float *rlen) { float adhocRatio = 0.95; float halfbal, branchlen = *rlen+*llen; float lweight = treeSize3way(lnode, lnode) /*- lnode->branchlen*/; float rweight = treeSize3way(rnode, rnode) /*- rnode->branchlen*/; halfbal = fabsf(lweight-rweight) / 2.0; if (halfbal < *llen && halfbal < *rlen) { if (lweight < rweight) { *llen += halfbal; *rlen -= halfbal; } else { *llen -= halfbal; *rlen += halfbal; } } else { /* The difference is larger than the final branch - give nearly all weight to the shorter one (cosmetic hack) */ if (lweight < rweight) { *llen = branchlen*adhocRatio; *rlen = branchlen*(1.0 - adhocRatio); } else { *rlen = branchlen*adhocRatio; *llen = branchlen*(1.0 - adhocRatio); } } } /* Print tree in New Hampshire format */ static void treePrintNH(treeStruct *tree, treeNode *node, FILE *file) { if (!node) return; if (node->left && node->right) { fprintf(file, "(\n"); treePrintNH(tree, node->left, file); fprintf(file, ",\n"); treePrintNH(tree, node->right, file); fprintf(file, ")\n"); if (node != tree->head) /* Not exactly sure why this is necessary, but njplot crashes otherwise */ fprintf(file, "%.3f", node->boot); } else { fprintf(file, "%s", node->name); } if (node != tree->head) /* Not exactly sure why this is necessary, but njplot crashes otherwise */ fprintf(file, ":%.3f", node->branchlen/100.0); } /* Save tree in New Hampshire format */ static void treePrintNH_init(void) { FILE *file; treeStruct *treestruct; if (!(file = filqueryopen(dirName, fileName, "","w", "Write New Hampshire format tree to file:"))) return; if (!graphAssFind(assVoid(1), &treestruct)) { messout("Could not find tree for this graph"); return; } treePrintNH(treestruct, treestruct->head, file); fprintf(file, ";\n"); fclose(file); } static void treePrintNode(treeNode *node) { if (node->name) printf("%s ", node->name); } static void treeFindOrthologsRecur(treeStruct *treestruct, treeNode *node, char *org1, char *org2) { if (!node || !node->left || !node->right) return; if (debug) { printf("\n 1 (%s, seq=%s): ", node->left->organism, node->left->name); printf("\n 2 (%s, seq=%s)\n: ", node->right->organism, node->right->name); } /* The open-minded way */ if (node->left->organism && node->right->organism && node->left->organism != node->right->organism) { graphBoxDraw(node->box, WHITE, BLACK); printf("\nSpecies 1 (%s): ", node->left->organism); treeTraverse(node->left, treePrintNode); printf("\nSpecies 2 (%s): ", node->right->organism); treeTraverse(node->right, treePrintNode); printf("\n"); } /* The narrow-minded way * / if ((node->left->organism == org1 && node->right->organism == org2) || (node->left->organism == org2 && node->right->organism == org1)) { printf("Found orthologs"); }*/ else { treeFindOrthologsRecur(treestruct, node->left, org1, org2); treeFindOrthologsRecur(treestruct, node->right, org1, org2); } } static void treeFindOrthologs(void) { treeStruct *treestruct; char *org1, *org2; if (!graphAssFind(assVoid(1), &treestruct)) { messout("Could not find tree for this graph"); return; } /* Dialog to define org1 and org2 */ treeFindOrthologsRecur(treestruct, treestruct->head, org1, org2); } static void treeOrder(treeNode *node) { if (!node) return; treeOrder(node->left); if (node->aln) node->aln->nr = treeOrderNr++; treeOrder(node->right); } static void printMtx(float **mtx) { int i, j; printf ("\n"); for (i = 0; i < nseq; i++) { for (j = 0; j < nseq; j++) printf("%6.2f ", mtx[i][j]); printf ("\n"); } } static float treeKimura(float i) { return i; } /* Correct an observed distance to a 'true' distance. od = observed distance in percent */ static float treeSTORMSONN(float od) { float cd; /* Corrected distance */ if (od > 91.6) return 1000.0; /* Otherwise log of negative value below */ od /= 100; cd= -log(1 -0.95844*od -0.69957*od*od +2.4955*od*od*od -4.6353*od*od*od*od +2.8076*od*od*od*od*od); return cd*100; } static treeNode *treeMake(void) { int i, j, iter, maxi, maxj, count=0; ALN *alni, *alnj; float id, totid=0, maxid, pmaxid, **pairmtx, **Dmtx, **mtx, *src, *trg, *avgdist, /* vector r in Durbin et al */ llen, rlen; treeNode **node, /* Array of (primary) nodes. Value=0 => stale column */ *newnode; static STORE_HANDLE treeHandle = 0; /* Allocate memory */ if (treeHandle) handleDestroy(treeHandle); treeHandle = handleCreate(); pairmtx = handleAlloc(0, treeHandle, nseq*sizeof(float *)); Dmtx = handleAlloc(0, treeHandle, nseq*sizeof(float *)); node = handleAlloc(0, treeHandle, nseq*sizeof(treeNode *)); avgdist = handleAlloc(0, treeHandle, nseq*sizeof(float)); for (i = 0; i < nseq; i++) { pairmtx[i] = handleAlloc(0, treeHandle, nseq*sizeof(float)); Dmtx[i] = handleAlloc(0, treeHandle, nseq*sizeof(float)); node[i] = messalloc(sizeof(treeNode)); node[i]->name = messalloc(strlen(arrp(Align, i, ALN)->name)+50); if (!treeCoordsOn) { sprintf(node[i]->name, "%s", arrp(Align, i, ALN)->name); } else { sprintf(node[i]->name, "%s/%d-%d", arrp(Align, i, ALN)->name, arrp(Align, i, ALN)->start, arrp(Align, i, ALN)->end); } node[i]->aln = arrp(Align, i, ALN); node[i]->organism = arrp(Align, i, ALN)->organism; } /* Calculate pairwise distance matrix */ arrayOrder(); for (i = 0; i < nseq-1; i++) { alni = arrp(Align, i, ALN); for (j = i+1; j < nseq; j++) { alnj = arrp(Align, j, ALN); pairmtx[i][j] = 100.0 - identity(alni->seq, alnj->seq); if (treeDistCorr == KIMURA) pairmtx[i][j] = treeKimura(pairmtx[i][j]); else if (treeDistCorr == STORMSONN) pairmtx[i][j] = treeSTORMSONN(pairmtx[i][j]); } } /* Construct tree */ for (iter = 0; iter < nseq-1; iter++) { if (treeMethod == NJ) { int count; /* Calculate vector r (avgdist) */ for (i = 0; i < nseq; i++) { if (!node[i]) continue; avgdist[i] = 0.0; for (count=0, j = 0; j < nseq; j++) { if (!node[j]) continue; avgdist[i] += (i < j ? pairmtx[i][j] : pairmtx[j][i]); count++; } if (count == 2) /* Hack, to accommodate last node */ avgdist[i] = 1; else avgdist[i] /= 1.0*(count - 2); } /* Calculate corrected matrix Dij */ if (1 /* gjm */) { for (i = 0; i < nseq-1; i++) { if (!node[i]) continue; for (j = i+1; j < nseq; j++) { if (!node[j]) continue; Dmtx[i][j] = pairmtx[i][j] - (avgdist[i] + avgdist[j]); } } } else { /* Felsenstein */ float Q = 0.0; for (i = 0; i < nseq-1; i++) { if (!node[i]) continue; for (j = i+1; j < nseq; j++) { if (!node[j]) continue; Q += pairmtx[i][j]; } } for (i = 0; i < nseq-1; i++) { if (!node[i]) continue; for (j = i+1; j < nseq; j++) { if (!node[j]) continue; Dmtx[i][j] = (pairmtx[i][j] + (2.0*Q)/(count-(count == 2 ? 1 : 2)) - avgdist[i] - avgdist[j]) / 2.0; } } } mtx = Dmtx; } else mtx = pairmtx; if (debug) { printf("Node status, Avg distance:\n"); for (i = 0; i < nseq; i++) printf("%6d ", (node[i] ? 1 : 0)); printf("\n"); for (i = 0; i < nseq; i++) printf("%6.2f ", avgdist[i]); printf("\n\nPairdistances, corrected pairdistances:"); printMtx(pairmtx); printMtx(Dmtx); printf("\n"); } /* Find smallest distance pair in pairmtx */ maxi = -1; maxj = -1; maxid = pmaxid = 1000000; for (i = 0; i < nseq-1; i++) { if (!node[i]) continue; for (j = i+1; j < nseq; j++) { if (!node[j]) continue; /* printf("%d %d - %d, dist= %f\n", iter, i, j, mtx[i][j]);*/ if (mtx[i][j] < maxid) { maxid = mtx[i][j]; pmaxid = pairmtx[i][j]; maxi = i; maxj = j; } else if (treeMethod == NJ && Dmtx[i][j] == maxid && pairmtx[i][j] < pmaxid) { /* To resolve ties - important for tree look! */ maxi = i; maxj = j; } } } maxid = pairmtx[maxi][maxj]; /* Don't want to point to Dmtx in NJ */ /* Merge rows & columns of maxi and maxj into maxi Recalculate distances to other nodes */ for (i = 0; i < nseq; i++) { if (!node[i]) continue; if (i < maxi) trg = &pairmtx[i][maxi]; else if (i > maxi) trg = &pairmtx[maxi][i]; else continue; if (i < maxj) src = &pairmtx[i][maxj]; else if (i > maxj) src = &pairmtx[maxj][i]; else continue; if (treeMethod == UPGMA) *trg = (*trg + *src) / 2.0; else *trg = (*trg + *src - maxid) / 2.0; } /* Create node for maxi and maxj */ newnode = messalloc(sizeof(treeNode)); if (treeMethod == UPGMA) { /* subtract lower branch lengths from absolute distance Horribly ugly, only to be able to share code UPGMA and NJ */ treeNode *tmpnode = node[maxi]->left; for (llen = maxid; tmpnode; tmpnode = tmpnode->left) llen -= tmpnode->branchlen; tmpnode = node[maxj]->right; for (rlen = maxid; tmpnode; tmpnode = tmpnode->right) rlen -= tmpnode->branchlen; } else { llen = (maxid + avgdist[maxi] - avgdist[maxj]) / 2.0; rlen = maxid - llen; if (iter == nseq-2) { /* Root node */ /* Not necessary anymore, the tree is re-balanced at the end which calls this too treeBalanceByWeight(node[maxi], node[maxj], &llen, &rlen);*/ /* Put entire length of root branch in one leg so the rebalancing will work properly (otherwise it is hard to take this branch into account */ rlen += llen; llen = 0; } } if (debug) { printf("Iter %d: Merging %d and %d, dist= %f\n", iter, maxi+1, maxj+1, mtx[maxi][maxj]); printf("maxid= %f llen= %f rlen= %f\n", maxid, llen, rlen); printf("avgdist[left]= %f avgdist[right]= %f\n\n", avgdist[maxi], avgdist[maxj]); } newnode->left = node[maxi]; newnode->left->branchlen = llen; newnode->right = node[maxj]; newnode->right->branchlen = rlen; newnode->organism = (node[maxi]->organism == node[maxj]->organism ? node[maxi]->organism : 0); node[maxi] = newnode; node[maxj] = 0; } fillParents(newnode, newnode->left); fillParents(newnode, newnode->right); if (debug) treeDraw(newnode); if (treeMethod == UPGMA) newnode->branchlen = 100-maxid; if (treeMethod == NJ) newnode = treeFindBalance(newnode); fillOrganism(newnode); return newnode; } static void treeDisplay(void) { separateMarkupLines(); treeHead = treeMake(); treeDraw(treeHead); reInsertMarkupLines(); } static void treeSortBatch(void) { separateMarkupLines(); treeHead = treeMake(); treeOrderNr = 1; treeOrder(treeHead); /* Set nr field according to tree order */ arraySort(Align, (void*)nrorder); reInsertMarkupLines(); } /* Find back the ALN pointers lost in treeSortBatch */ static void treeFindAln(treeNode *node) { if (!node->name) return; str2aln(node->name, &aln); if (!alignFind(Align, &aln, &ip)) { messout("Cannot find back seq after sort - probably a bug"); } else node->aln = arrp(Align, ip, ALN); } static void treeSort(void) { ALN aln; if (Highlight_alnp) alncpy(&aln, Highlight_alnp); treeSortBatch(); /* After treeSortBatch the treeNodes point to the wrong ALN, because arraySort only copies contents and not entire structs. This can be fixed by either a linear string search or by remaking the tree... */ treeTraverse(treeHead, treeFindAln); if (Highlight_alnp) { if (!alignFind(Align, &aln, &ip)) { messout("Cannot find back highlighted seq after sort - probably a bug"); Highlight_alnp = 0; } else Highlight_alnp = arrp(Align, ip, ALN); centerHighlighted(); } treeDraw(treeHead); belvuRedraw(); } static void columnCopy(Array Aligndest, int dest, Array Alignsrc, int src) { int i; for (i = 0; i < nseq; i++) arrp(Aligndest, i, ALN)->seq[dest] = arrp(Alignsrc, i, ALN)->seq[src]; } static void treeBootstrap(void) { Array Aligntmp; int i, iter, src; treeStruct *treestruct = messalloc(sizeof(treeStruct)); separateMarkupLines(); Aligntmp = arrayCopy(Align); for (i = 0; i < nseq; i++) { arrp(Aligntmp, i, ALN)->seq = messalloc(strlen(arrp(Align, i, ALN)->seq)+1); strcpy(arrp(Aligntmp, i, ALN)->seq, arrp(Align, i, ALN)->seq); } #if defined(__CYGWIN__) || defined(DARWIN) srand(time(0)); #else srand48(time(0)); #endif for (iter = 0; iter < treebootstraps; iter++) { /* Make bootstrapped Alignment */ for (i = 0; i < maxLen; i++) { #if defined(__CYGWIN__) || defined(DARWIN) src = (int)((((float)rand())/RAND_MAX)*maxLen); #else src = (int)(drand48()*maxLen); #endif printf("%d %d\n", i, src); columnCopy(Align, i, Aligntmp, src); } /* printf("&\n"); */ treestruct->head = treeMake(); if (treebootstrapsDisplay) treeDraw(treestruct->head); else { printf("\n"); treePrintNH(treestruct, treestruct->head, stdout); printf(";\n"); } } /* Restore alignment */ for (i = 0; i < nseq; i++) { strcpy(arrp(Align, i, ALN)->seq, arrp(Aligntmp, i, ALN)->seq); } } static void highlightScoreSort(char mode) { int i, len; if (!Highlight_alnp) { messout("Please highlight a sequence first"); return; } if (Highlight_alnp->markup) { messout("Please do not highlight a markup line"); return; } alncpy(&aln, Highlight_alnp); separateMarkupLines(); /* if (displayScores) { if (!(graphQuery("This will erase the current scores, do you want to continue?"))) return; }*/ displayScores = 1; /* Calculate score relative to highlighted sequence */ for (i = 0; i < nseq; i++) { if (mode == 'P') arrp(Align, i, ALN)->score = score(aln.seq, arrp(Align, i, ALN)->seq); else if (mode == 'I') arrp(Align, i, ALN)->score = identity(aln.seq, arrp(Align, i, ALN)->seq); if ((len=strlen(messprintf("%.1f", arrp(Align, i, ALN)->score))) > maxScoreLen) maxScoreLen = len; } arraySort(Align, (void*)scoreorderRev); reInsertMarkupLines(); if (!alignFind(Align, &aln, &i)) { messout("Cannot find back highlighted seq after sort - probably a bug"); Highlight_alnp = 0; } else Highlight_alnp = arrp(Align, i, ALN); arrayOrder(); AlignYstart = 0; belvuRedraw(); } static void simSort(void) { highlightScoreSort('P'); } static void idSort(void) { highlightScoreSort('I'); } static void scoreSortBatch(void) { if (!displayScores) messcrash("No scores available"); arraySort(Align, (void*)scoreorder); arrayOrder(); } static void scoreSort(void) { if (Highlight_alnp) alncpy(&aln, Highlight_alnp); arraySort(Align, (void*)scoreorder); if (Highlight_alnp) { if (!alignFind(Align, &aln, &ip)) { messout("Cannot find back highlighted seq after sort - probably a bug"); Highlight_alnp = 0; } else Highlight_alnp = arrp(Align, ip, ALN); } arrayOrder(); AlignYstart = 0; belvuRedraw(); } static void printScore(void) { if (!Highlight_alnp) { messout("Select a line first"); return; } printf("%.1f %s/%d-%d\n", Highlight_alnp->score, Highlight_alnp->name, Highlight_alnp->start, Highlight_alnp->end); fflush(stdout); } int GCGchecksum(char *seq) { int check = 0, count = 0, i; for (i = 0; i < maxLen; i++) { count++; check += count * toupper((int) seq[i]); if (count == 57) count = 0; } return (check % 10000); } int GCGgrandchecksum(void) { int i, grand_checksum=0; for(i=0; i < nseq; i++) grand_checksum += GCGchecksum(arrp(Align, i, ALN)->seq); return (grand_checksum % 10000); } int writeMSF(FILE *pipe) /* c = separator between name and coordinates */ { int i, j, maxfullnamelen = maxNameLen+1+maxStartLen+1+maxEndLen, paragraph=0, alnstart, alnend, alnlen, linelen=50, blocklen=10; if (saveCoordsOn) maxfullnamelen = maxNameLen+1+maxStartLen+1+maxEndLen; else maxfullnamelen = maxNameLen; /* Title */ fprintf(pipe, "PileUp\n\n %s MSF: %d Type: X Check: %d ..\n\n", Title, maxLen, GCGgrandchecksum()); /* Names and checksums */ for(i=0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (alnp->markup) continue; fprintf(pipe, " Name: %-*s Len: %5d Check: %5d Weight: %.4f\n", maxfullnamelen, (saveCoordsOn ? messprintf("%s%c%d-%d", alnp->name, saveSeparator, alnp->start, alnp->end) : messprintf("%s", alnp->name)), maxLen, GCGchecksum(alnp->seq), 1.0); } fprintf(pipe, "\n//\n\n"); /* Alignment */ while (paragraph*linelen < maxLen) { for (i = 0; i < nseq; i++) { alnstart = paragraph*linelen; alnlen = ( (paragraph+1)*linelen < maxLen ? linelen : maxLen-alnstart ); alnend = alnstart + alnlen; alnp = arrp(Align, i, ALN); if (alnp->markup) continue; fprintf(pipe, "%-*s ", maxfullnamelen, (saveCoordsOn ? messprintf("%s%c%d-%d", alnp->name, saveSeparator, alnp->start, alnp->end) : messprintf("%s", alnp->name))); for (j = alnstart; j < alnend; ) { fprintf(pipe, "%c", alnp->seq[j]); j++; if ( !((j-alnstart)%blocklen) ) fprintf(pipe, " "); } fprintf(pipe, "\n"); } fprintf(pipe, "\n"); paragraph++; } fclose(pipe); fflush(pipe); } static void saveMSF(void) { FILE *fil; if (!(fil = filqueryopen(dirName, fileName, "","w", "Save as MSF (/) file:"))) return; writeMSF(fil); saved = 1; } static void readMSF(FILE *pipe) { char junk[1001], seq[1001], *cp, *cq; int len, i; fprintf(stderr, "\nDetected MSF format\n"); /* Read sequence names */ while (!feof (pipe)) { if (!fgets (line, MAXLENGTH, pipe)) break; if (!strncmp(line, "//", 2)) break; else if (strstr(line, "Name:") && (cp = strstr(line, "Len:")) && strstr(line, "Check:")) { int len; sscanf(cp+4, "%d", &len); if (maxLen < len) maxLen = len; cp = strstr(line, "Name:")+6; while (*cp && *cp == ' ') cp++; parseMulLine(cp, &aln); aln.len = len; aln.seq = messalloc(aln.len+1); aln.nr = ++nseq; if (!arrayInsert(Align, &aln, (void*)alphaorder)) messout("Failed to arrayInsert %s", aln.name); } } /* Read sequence alignment */ while (!feof (pipe)) { if (!fgets (line, MAXLENGTH, pipe)) break; cp = line; cq = seq; while (*cp && *cp == ' ') cp++; while (*cp && *cp != ' ') cp++; /* Spin to sequence */ while (*cp && cq-seq < 1000) { if (isAlign(*cp)) *cq++ = *cp; cp++; } *cq = 0; if (*seq) { cp = line; while (*cp && *cp == ' ') cp++; parseMulLine(cp, &aln); if (arrayFind(Align, &aln, &i, (void*)alphaorder)) { alnp = arrp(Align, i, ALN); strcat(alnp->seq, seq); } else { messout("Cannot find back %s %d %d seq=%s.\n", aln.name, aln.start, aln.end, seq); } } } strcpy(saveFormat, MSFStr); } static void readFastaAlnFinalise(char *seq) { aln.len = strlen(seq); aln.seq = messalloc(aln.len+1); strcpy(aln.seq, seq); if (maxLen) { if (aln.len != maxLen) fatal("Differing sequence lengths: %s %s", maxLen, aln.len); } else maxLen = aln.len; if (arrayFind(Align, &aln, &ip, (void*)alphaorder)) fatal ("Sequence name occurs more than once: %s%c%d-%d", aln.name, saveSeparator, aln.start, aln.end); else { aln.nr = ++nseq; arrayInsert(Align, &aln, (void*)alphaorder); } } /* readFastaAln() * * Read aligned fasta format, assumed to be: * * >5H1A_HUMAN/53-400 more words * GNACVVAAIAL...........ERSLQ.....NVANYLIG..S.LAVTDL * MVSVLV..LPMAAL.........YQVL..NKWTL......GQVT.CDL.. * >5H1A_RAT more words * GNACVVAAIAL...........ERSLQ.....NVANYLIG..S.LAVTDL * MVSVLV..LPMAAL.........YQVL..NKWTL......GQVT.CDL.. * ... */ static void readFastaAln(FILE *pipe) { char *cp, *cq, *seq, *seqp; int filesize=0, started=0; stack = stackReCreate (stack, 100); while (!feof (pipe)) { if (!fgets (line, MAXLENGTH, pipe)) break; if (cp = strchr(line, '\n')) *cp = 0; pushText(stack, line); filesize += strlen(line); } seq = messalloc(filesize); while (cp = stackNextText(stack)) { if (*cp == '>') { if (started) readFastaAlnFinalise(seq); parseMulLine(cp+1, &aln); seqp = seq; started = 1; } else { strcpy(seqp, cp); seqp += strlen(seqp); } } if (started) readFastaAlnFinalise(seq); messfree(seq); strcpy(saveFormat, FastaAlnStr); } /* Convenience routine for converting "name/start-end" to "name start end". Used by parsing routines. Return 1 if both tokens were found, otherwise 0. */ static int stripCoordTokens(char *cp) { char *end = cp; while (*end && !isspace(*end)) end++; if ((cp = strchr(cp, saveSeparator)) && cp < end) { *cp = ' '; if ((cp = strchr(cp, '-')) && cp < end) { *cp = ' '; IN_FORMAT = MUL; return 1; } } IN_FORMAT = RAW; return 0; } /* Parse name, start and end of a Mul format line Convenience routine, part of readMul and other parsers */ static void parseMulLine(char *line, ALN *aln) { char line2[MAXLENGTH+1], *cp=line2, *cq, *space, GRfeat[MAXNAMESIZE+1]; strcpy(cp, line); resetALN(aln); if (!strncmp(cp, "#=GC", 4)) { aln->markup = GC; cp += 5; } if (!strncmp(cp, "#=RF", 4)) { aln->markup = GC; } if (!strncmp(cp, "#=GR", 4)) { aln->markup = GR; cp += 5; } if (stripCoordTokensOn) stripCoordTokens(cp); /* Name */ strncpy(aln->name, cp, MAXNAMESIZE); aln->name[MAXNAMESIZE] = 0; if (cq = strchr(aln->name, ' ')) *cq = 0; /* Add Start and End coords */ if (stripCoordTokensOn && (IN_FORMAT != RAW) ) sscanf(strchr(cp, ' ') + 1, "%d%d", &aln->start, &aln->end); if (aln->markup == GR) { /* Add GR markup names to name #=GR O83071 192 246 SA 999887756453524252..55152525....36463774777.....948472782969685958 ->name = O83071_SA */ if (IN_FORMAT == MUL) { sscanf(cp + strlen(aln->name), "%d%d%s", &aln->start, &aln->end, GRfeat); } else { sscanf(cp + strlen(aln->name), "%s", GRfeat); } if (strlen(aln->name)+strlen(GRfeat)+1 > MAXNAMESIZE) messout("Too long name or/and feature name"); strcat(aln->name, " "); strncat(aln->name, GRfeat, MAXNAMESIZE-strlen(aln->name)); } } /* ReadMul * parses an alignment file in mul or selex format and puts it in the Align array * * Assumes header contains ' ' or '#' * * Alignment can have one of the following formats: * CSW_DROME VTHIKIQNNGDFFDLYGGEKFATLP * CSW_DROME 51 75 VTHIKIQNNGDFFDLYGGEKFATLP * CSW_DROME 51 75 VTHIKIQNNGDFFDLYGGEKFATLP P29349 * KFES_MOUSE/458-539 .........WYHGAIPW.....AEVAELLT........HTGDFLVRESQG * */ static void readMul(FILE *pipe) { char *cp, *cq, ch; int len, i, alnstart; BOOL new; stack = stackReCreate (stack, 100); AnnStack = stackReCreate (AnnStack, 100); /* Parse header to check for MSF or Fasta format */ while (!feof (pipe)) { if (!isspace(ch = fgetc(pipe))) { if (ch == '>') { ungetc(ch, pipe); readFastaAln(pipe); return; } else break; } else if (ch != '\n') if (!fgets (line, MAXLENGTH, pipe)) break; if (strstr(line, "MSF:") && strstr(line, "Type:") && strstr(line, "Check:")) { readMSF(pipe); return; } } if (!feof(pipe)) ungetc(ch, pipe); /* Read raw alignment into stack *******************************/ alnstart = MAXLENGTH; while (!feof (pipe)) { if (!fgets (line, MAXLENGTH, pipe) || (unsigned char)*line == (unsigned char)EOF ) break; /* EOF checking to make acedb calling work */ if (!strncmp(line, "PileUp", 6)) { readMSF(pipe); return; } if (cp = strchr(line, '\n')) *cp = 0; len = strlen (line); /* Sequences */ if (len && *line != '#' && strcmp(line, "//")) { if (!(cq = strchr(line, ' '))) fatal("No spacer between name and sequence in %s", line); /* Find leftmost start of alignment */ for (i = 0; line[i] && line[i] != ' '; i++); for (; line[i] && !isAlign(line[i]); i++); if (i < alnstart) alnstart = i; /* Remove optional accession number at end of alignment */ /* FOR PRODOM STYLE ALIGNMENTS W/ ACCESSION TO THE RIGHT - MAYBE MAKE OPTIONAL This way it's incompatible with alignments with ' ' gapcharacters */ for (; line[i] && isAlign(line[i]); i++); line[i] = 0; pushText(stack, line); } /* Markup line */ else if (!strncmp(line, "#=GF ", 5) || !strncmp(line, "#=GS ", 5)) { pushText(AnnStack, line); } else if (!strncmp(line, "#=GC ", 5) || !strncmp(line, "#=GR ", 5) || !strncmp(line, "#=RF ", 5)) { pushText(stack, line); } /* Match Footer */ else if (!strncmp(line, "# matchFooter", 13)) { matchFooter = 1; break; } } /* Read alignment from Stack into Align array ***************************************/ /* * First pass - Find number of unique sequence names and the total sequence lengths */ stackCursor(stack, 0); while (cp = stackNextText(stack)) { parseMulLine(cp, &aln); /* Sequence length */ len = strlen(cp+alnstart); if (arrayFind(Align, &aln, &ip, (void*)alphaorder)) arrp(Align, ip, ALN)->len += len; else { aln.len = len; aln.nr = ++nseq; arrayInsert(Align, &aln, (void*)alphaorder); /*printf("\nInserted %s %6d %6d %s %d\n", aln.name, aln.start, aln.end, cp+alnstart, len); fflush(stdout);*/ } } /* Find maximum length of alignment; allocate it */ for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (alnp->len > maxLen) maxLen = alnp->len; } for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); alnp->seq = messalloc(maxLen+1); } /* * Second pass - allocate and read in sequence data */ stackCursor(stack, 0); while (cp = stackNextText(stack)) { parseMulLine(cp, &aln); if (arrayFind(Align, &aln, &i, (void*)alphaorder)) { alnp = arrp(Align, i, ALN); strcat(alnp->seq, cp+alnstart); } else { messout("Cannot find back %s %d %d seq=%s.\n", aln.name, aln.start, aln.end, aln.seq); } } /* Parse sequence attributes from #=GS lines */ stackCursor(AnnStack, 0); while (cp = stackNextText(AnnStack)) { char *namep, /* Seqname */ name[MAXNAMESIZE+1],/* Seqname */ *labelp, /* Label (OS) */ *valuep; /* Organism (c. elegans) */ if (strncmp(cp, "#=GS", 4)) continue; strcpy(name, cp+4); namep = name; while (*namep == ' ') namep++; labelp = namep; while (*labelp != ' ') labelp++; *labelp = 0; /* Terminate namep */ labelp++; /* Walk across terminator */ while (*labelp == ' ') labelp++; valuep = labelp; while (*valuep != ' ') valuep++; while (*valuep == ' ') valuep++; if (!strncasecmp(labelp, "LO", 2)) { int i, colornr; for (colornr = -1, i = 0; i < NUM_TRUECOLORS; i++) if (!strcasecmp(colorNames[i], valuep)) colornr = i; if (colornr == -1) { printf("Unrecognized color: %s\n", valuep); colornr = 0; } str2aln(namep, &aln); /* Find the corresponding sequence */ if (!arrayFind(Align, &aln, &i, (void*)alphaorder)) { messout("Cannot find %s %d %d in alignment", aln.name, aln.start, aln.end); continue; } alnp = arrp(Align, i, ALN); alnp->color = colornr; } else if (!strncasecmp(labelp, OrganismLabel, 2)) { /* Add organism to table of organisms. This is necessary to make all sequences of the same organism point to the same place and to make a non-redundant list of present organisms */ /* Find string in permanent stack */ if (!(valuep = strstr(cp, valuep))) { messout("Failed to parse organism properly"); continue; } aln.organism = valuep; /* Find string in permanent stack */ arrayInsert(organismArr, &aln, (void*)organism_order); if (strchr(cp, '/') && strchr(cp, '-')) { str2aln(namep, &aln); /* Find the corresponding sequence */ if (!arrayFind(Align, &aln, &i, (void*)alphaorder)) { messout("Cannot find %s %d %d in alignment", aln.name, aln.start, aln.end); continue; } alnp = arrp(Align, i, ALN); /* Store pointer to unique organism in ALN struct */ aln.organism = valuep; arrayFind(organismArr, &aln, &ip, (void*)organism_order); alnp->organism = arrp(organismArr, ip, ALN)->organism; } else { /* Organism specified for entire sequences. Find all ALN instances of this sequences. */ for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (!strcmp(alnp->name, namep)) { /* Store pointer to unique organism in ALN struct */ aln.organism = valuep; arrayFind(organismArr, &aln, &ip, (void*)organism_order); alnp->organism = arrp(organismArr, ip, ALN)->organism; } } } } } /* For debugging * / for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); printf("\n%-10s %4d %4d %s %d %s\n", alnp->name, alnp->start, alnp->end, alnp->seq, alnp->len, alnp->organism); } for (i=0; i < arrayMax(organismArr); i++) printf("%s\n", arrp(organismArr, i, ALN)->organism); */ if (!nseq || !maxLen) fatal("Unable to read sequence data"); strcpy(saveFormat, MulStr); } static void readScores(char *filename) { char line[MAXLENGTH+1], linecp[MAXLENGTH+1], *cp; FILE *file; ALN aln; int scoreLen; if (!(file = fopen(filename, "r"))) fatal("Cannot open file %s", filename); while (!feof (file)) { if (!fgets (line, MAXLENGTH, file)) break; strcpy(linecp, line); resetALN(&aln); if (!(cp = strtok(line, " "))) fatal("Error parsing score file %s.\nLine: %s", filename, linecp); if (!(sscanf(cp, "%f", &aln.score))) fatal("Error parsing score file %s - bad score.\nLine: %s", filename, linecp); if (!(cp = strtok(0, "/"))) fatal("Error parsing score file %s.\nLine: %s", filename, linecp); strncpy(aln.name, cp, MAXNAMESIZE); aln.name[MAXNAMESIZE] = 0; if (!(cp = strtok(0, "-"))) fatal("Error parsing score file %s.\nLine: %s", filename, linecp); if (!(aln.start = atoi(cp))) fatal("Error parsing score file %s - no start coordinate.\nLine: %s", filename, linecp); if (!(cp = strtok(0, "\n"))) fatal("Error parsing score file %s.\nLine: %s", filename, linecp); if (!(aln.end = atoi(cp))) fatal("Error parsing score file %s - no end coordinate.\nLine: %s", filename, linecp); if (!alignFind(Align, &aln, &ip)) { /* printf("Warning: %s/%d-%d (score %.1f) not found in alignment\n", aln.name, aln.start, aln.end, aln.score);*/ } else { arrp(Align, ip, ALN)->score = aln.score; scoreLen = strlen(messprintf("%.1f", aln.score)); if (scoreLen > maxScoreLen) maxScoreLen = scoreLen; } } fclose(file); displayScores = 1; } static void checkAlignment() { int i, g, cres, nres, tmp; maxNameLen = maxStartLen = maxEndLen = 0; for (i = 0; i < nseq; i++) { alnp = arrp(Align, i, ALN); if (!alnp->markup) { /* Count residues */ for (cres = g = 0; g < maxLen; g++) { if (isAlign(alnp->seq[g]) && !isGap(alnp->seq[g])) cres++; } if (!alnp->start) { /* No coords provided - reconstruct them */ alnp->start = 1; alnp->end = cres; } else { /* Check if provided coordinates are correct */ nres = abs(alnp->end - alnp->start) + 1; if ( nres != cres) fprintf(stderr, "Found wrong number of residues in %s/%d-%d: %d instead of %d\n", alnp->name, alnp->start, alnp->end, cres, nres); } } /* Find max string length of name, start, and end */ if (maxNameLen < strlen(alnp->name)) maxNameLen = strlen(alnp->name); if (maxStartLen < (tmp = strlen(messprintf("%d", alnp->start)))) maxStartLen = tmp; if (maxEndLen < (tmp = strlen(messprintf("%d", alnp->end)))) maxEndLen = tmp; } } void argvAdd(int *argc, char ***argv, char *s) { char **v; int i; v = (char **)malloc(sizeof(char *)*(*argc+2)); /* Copy argv */ for (i = 0; i < (*argc); i++) v[i] = (*argv)[i]; /* Add s */ v[*argc] = (char *)malloc(strlen(s)+1); strcpy(v[*argc], s); (*argc)++; /* Terminator - ANSI C standard */ v[*argc] = 0; /* free(*argv); */ /* Too risky */ *argv = v; } int main(int argc, char **argv) { FILE *file, *pipe; char command[MAXLINE+1], *scoreFile = 0, *readMatchFile = 0, *colorCodesFile = 0, *markupColorCodesFile = 0, *cp, *output_format = 0, *optargc; int i, dx, dy, pw, ph, /* pixel width and height */ output_probs = 0, init_tree = 0, only_tree = 0, show_ann = 0; float makeNRinit = 0.0, init_rmEmptyColumns = 0.0, init_rmGappySeqs = 0.0, cw; /* character width of initial alignment */ /* extern int printOnePage;*/ int optc; extern int optind; extern char *optarg; char *optstring="ab:Ccgl:L:m:n:O:o:PpQ:q:RrS:s:T:t:uz:"; static char *cc_date = #if defined(__DATE__) __DATE__ #else "" #endif ; char *usage; static char usageText[] = "\ \n\ Belvu - View multiple alignments in good-looking colours.\n\ \n\ Usage: belvu [options] |- [X options]\n\ \n\ |- = alignment file or pipe.\n\ \n\ Options:\n\ -c Print Conservation table.\n\ -l Load color code file.\n\ -L Load Markup color code file.\n\ -m Read file with matching sequence segments.\n\ -r Read alignment in 'raw' format (Name sequence).\n\ -R Do not parse coordinates when reading alignment.\n\ -o Write alignment or tree to stdout in this format and exit.\n\ Valid formats: MSF, Mul(Stockholm), Selex, \n\ FastaAlign, Fasta, tree.\n\ -n Make non-redundant to %identity at startup.\n\ -Q Remove columns more gappy than .\n\ -q Remove sequences more gappy than .\n\ -P Remove partial sequences at startup.\n\ -C Don't write coordinates to saved file.\n\ -z Separator char between name and coordinates in saved file.\n\ -a Show alignment annotations on screen (Stockholm format only).\n\ -p Output random model probabilites for HMMER.\n\ (Based on all residues.)\n\ -S Sort sequences in this order.\n\ a -> alphabetically\n\ o -> by Swissprot organism, alphabetically\n\ s -> by score\n\ n -> by Neighbor-joining tree\n\ u -> by UPGMA tree\n\ S -> by similarity to first sequence\n\ i -> by identity to first sequence\n\ -s Read in file of scores.\n\ -T Start up with tree calculated by method:\n\ n -> Neighbor-joining\n\ u -> UPGMA\n\ N -> Neighbor-joining, only show tree\n\ U -> UPGMA, only show tree\n\ c -> Don't color tree by organism\n\ o -> Don't display sequence coordinates in tree\n\ s -> Use Storm & Sonnhammer distance correction\n\ r -> Use uncorrected distances\n\ -b <#> Print out # bootstrapped trees and exit\n\ (Negative value -> display bootstrap trees on screen)\n\ -O