/* File: fmapcontrol.c * Author: Richard Durbin (rd@sanger.ac.uk) * Copyright (C) J Thierry-Mieg and R Durbin, 1992 * ------------------------------------------------------------------- * Acedb is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * or see the on-line version at http://www.gnu.org/copyleft/gpl.txt * ------------------------------------------------------------------- * This file is part of the ACEDB genome database package, written by * Richard Durbin (Sanger Centre, UK) rd@sanger.ac.uk, and * Jean Thierry-Mieg (CRBM du CNRS, France) mieg@kaa.cnrs-mop.fr * * Description: display choice control for fmap * Exported functions: * HISTORY: * Last edited: Oct 2 10:02 2002 (edgrif) * * Mar 13 22:55 2002 (rd): fixed fMapFindSpan() to use sMapInverseMap() not SEGs and removed fMapFindZoneFather(), replacing calls with fMapFindSpan() * * Apr 3 14:33 2000 (edgrif): On behalf of Lincoln I have: * "Reintegrated some of the GFF dumping features introduced in 4.7l." * * Mar 30 15:42 2000 (edgrif): Added code for "no redrawing" of columns while * selecting them. * * Oct 5 16:30 1999 (fw): new methodcache system * * May 25 15:15 1999 (edgrif): Put column reuse into main menu. * * May 21 22:29 1999 (edgrif): Add support for retaining column selections on map reuse. * * Oct 21 11:09 1998 (edgrif): Corrected func. proto for fMapFollowVirtualMultiTrace * a whole load of fmap functions are nulled for ACEMBLY way down the file. * * Jul 27 10:36 1998 (edgrif): Add calls to fMapIntialise() for externally * called routines. * * Jul 16 10:06 1998 (edgrif): Introduce private header fmap_.h * * Jun 24 14:17 1998 (edgrif): Fix initialisation bug to make sure isDone * is set in fmapInitialise. Also remove methodInitialise from fmapInitialise, * this is not needed, method does its own intialisation. * * Dec 21 09:00 1995 (mieg) * * Jun 21 00:18 1995 (rd): fused with Jean - remaining discrepancies in * readMethod(s) sections - probably he found a bug here with 1, 2 etc. * and with HOMOL unshading in fMapSelect - I am quite confident here. * Created: Sat Jul 25 20:19:11 1992 (rd) * CVS info: $Id: fmapcontrol.c,v 1.253 2002/11/06 18:58:15 srk Exp $ *------------------------------------------------------------------- */ #undef FMAP_DEBUG #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /*************** globals within fMap Package *******************/ /******* addresses of these two are used as unique identifiers */ /* find the FeatureMap struct on the active graph */ magic_t GRAPH2FeatureMap_ASSOC = "FeatureMap"; /* verify a generic pointer to be a FeatureMap */ magic_t FeatureMap_MAGIC = "FeatureMap"; /* Can be set to TRUE to turn on debugging. */ BOOL fMapDebug = FALSE ; /* TEMP for testing stuff for James... */ int fMapCutCoords_G = FALSE ; /********************** local functions***********************/ #ifdef USE_SMAP static BOOL reportErrors; static consMapDNAErrorReturn callback(KEY key, int position); #endif static void fMapDrag(float y); static void fMapSelectBox (FeatureMap look, int box, double x, double y); static void fMapFollow (FeatureMap look, double x, double y) ; static void fMapCentre (FeatureMap look, KEY key, KEY from) ; static void fMapPick (int box, double x, double y) ; static void fMapCompHelp (void) ; static void fMapReverse (void) ; static BOOL fMapDoRecalculate (FeatureMap look, int start, int stop) ; static void fMapZoneKeySet (void) ; static void fMapKbd (int k) ; static FeatureMap fMapCreateLook (KEY seq, int start, int stop, BOOL reuse, FmapDisplayData *display_data); static BOOL setStartStop(KEY seq, int *start_inout, int *stop_inout, KEY *key_out, FeatureMap look) ; static void fMapDefaultColumns (MAP map) ; static BOOL checkSequenceLength(KEY seq) ; static void addToSet (ACEIN command_in, DICT *dict) ; static void dumpstats(FeatureMap look) ; /* One day these should be in menu code. */ static int getPositionInMenu(MenuBaseFunction func, MENUOPT *menu) ; static BOOL removeFromMenu(MenuBaseFunction func, MENUOPT *menu) ; static void toggleTextInMenu(MENUOPT *m, MENUFUNCTION menu_func, char *menu_text) ; /* prototypes for menu funcs */ static void setDisplayPreserve(void) ; static void toggleColumnButtons (void) ; static void toggleReportMismatch(MENUITEM unused) ; static void togglePreserveColumns(MENUITEM unused) ; static void toggleRedrawColumns(MENUITEM unused) ; static void toggleCutSelection(MENUITEM unused) ; static void zoom1 (void) ; static void zoom10 (void) ; static void zoom100(void) ; static void zoom1000 (void) ; static void fMapWhole (void) ; static void fMapMenuToggleTranslation (void) ; static void fMapMenuReverseComplement (void) ; static void fMapMenuRecalculate (void) ; static void fMapMenuComplementSurPlace (void) ; static void fMapMenuComplement (void) ; static void fMapMenuReverseComplement (void) ; static void fMapMenuToggleHeader (void); /* "Hide Header" */ static void fMapMenuDumpSeq (void); /* "Export Sequence" */ static void fMapMenuDumpSegs (void) ; /* "Export Features" */ static void geneStats (void) ; /* "Statistics" */ extern void emblDump (void) ; /* "EMBL dump" */ static void fMapMenuReadFastaFile (void); /* "Import Sequence" */ /* functions that manipulate the given FeatureMap look */ static void fMapToggleDna (FeatureMap look); static BOOL fMapToggleTranslation (FeatureMap look); static void fMapToggleHeader (FeatureMap look); #ifdef ACEMBLY static void fMapTraceJoin (void) ; #endif /* Strings that must be toggle in menu items, a better solution would be */ /* toggle buttons for many items which I will do sometime. */ #define PRESERVE_COLS "Preserve Current Columns" #define UN_PRESERVE_COLS "Do not Preserve Current Columns" #define NOREDRAW_COLS "No Automatic Redraw of Columns" #define REDRAW_COLS "Automatic Redraw of Columns" #define REPORT_MISMATCHES "Report DNA mismatches" #define NO_REPORT_MISMATCHES "Don't report DNA mismatches" #define CUT_DNA "Put DNA in cut buffer" #define CUT_COORDS "Put Object Coords in cut buffer" /********************** locals *******************************/ static KEY _Feature ; static KEY _Tm, _Temporary ; static KEY _EMBL_feature ; KEY _Arg1_URL_prefix ; /* really ugly to make global */ KEY _Arg1_URL_suffix ; KEY _Arg2_URL_prefix ; KEY _Arg2_URL_suffix ; static KEY _Spliced_cDNA, _Transcribed_gene ; static KEY _VTranscribed_gene = 0, _VIST = 0 ; /* mhmp 25.09.98 */ static int weight, seqDragBox, nbTrans ; static FeatureMap selectedfMap = NULL ; /************************************************************/ /* Default map column settings for fmap. */ static FeatureMapColSettingsStruct defaultMapColumns[] = { /* all column position 100, -90 ... should be different */ #define STL_STATUS_0 {-100.0, TRUE, "Locator", fMapShowMiniSeq}, {-90.0, TRUE, "Sequences & ends", fMapShowCanonical}, #ifdef STL_STATUS /* mike holman status column */ {-89.0, TRUE, "Cosmids by group", fMapShowOrigin}, #endif {-2.1, FALSE, "Up Gene Translation", fMapShowUpGeneTranslation}, {-1.9, TRUE, "-Confirmed introns", fMapShowSoloConfirmed}, {-0.1, TRUE, "Restriction map", fMapShowCptSites}, #ifdef ACEMBLY /* Contig bar supersedes locator */ {0.0, FALSE, "Summary bar", fMapShowSummary}, #else {0.0, TRUE, "Summary bar", fMapShowSummary}, #endif {0.1, TRUE, "Scale", fMapShowScale}, {1.9, TRUE, "Confirmed introns", fMapShowSoloConfirmed}, {3.0, TRUE, "EMBL features", fMapShowEmblFeatures}, #ifdef STL_STATUS {3.1, TRUE, "Cosmid Status", fMapShowStatus}, #endif {3.2, FALSE, "CDS Lines", fMapShowCDSLines}, {3.25, FALSE, "CDS Boxes", fMapShowCDSBoxes}, {3.3, TRUE, "Alleles", fMapShowAlleles}, {3.4, FALSE, "cDNAs", fMapShowcDNAs}, {3.5, FALSE, "Gene Names", fMapShowGeneNames}, {3.7, TRUE, "Assembly Tags", fMapShowAssemblyTags}, {3.8, TRUE, "Oligos", fMapOspShowOligo}, {3.82, TRUE, "Oligo_pairs", fMapOspShowOligoPairs}, /* isFrame starts if either of next 2 are On */ {4.0, FALSE, "3 Frame Translation", fMapShow3FrameTranslation}, {4.05, FALSE, "ORF's", fMapShowORF}, {4.1, TRUE, "Coding Frame", fMapShowCoding}, /* only shows if isFrame */ {4.2, FALSE, "ATG", fMapShowATG}, /* frame dependent stuff ends */ /* {4.98, FALSE, "Gene Translation", fMapShowGeneTranslation}, */ {4.99, FALSE, "Down Gene Translation", fMapShowDownGeneTranslation}, #ifdef ACEMBLY {5.5, TRUE, "Alignements", fMapShowAlignments}, {5.6, TRUE, "Previous Contigs", fMapShowPreviousContigs}, {5.7, TRUE, "Contigs", fMapShowContig}, {5.8, TRUE, "Trace Assembly", fMapShowTraceAssembly}, {5.9, TRUE, "Multiplets", fMapShowVirtualMultiplets}, #endif {6.0, FALSE, "Coords", fMapShowCoords}, {6.1, FALSE, "DNA Sequence", fMapShowDNA}, #ifdef ACEMBLY {6.2, FALSE, "Assembly DNA", fMapShowVirtualDna}, #endif {6.5, FALSE, "Brief Identifications", fMapShowBriefID}, {7.0, TRUE, "Text Features", fMapShowText}, {0.0, FALSE, NULL, NULL} } ; /************************************************************/ /* Main fmap control buttons. */ /* */ /* The buttons... */ static MENUOPT buttonOpts[] = { #ifndef ACEMBLY {toggleColumnButtons, "Columns"}, #endif /* ACEMBLY */ { mapZoomIn, "Zoom In.."}, { mapZoomOut, "Zoom Out.."}, #ifdef ACEMBLY /* dna... called in fmaptrace */ { fMapMenuReverseComplement, "Complement"}, { dnaAnalyse, "Analysis.."}, { fMapTraceJoin, "Assembly.."}, #else /* not ACEMBLY */ { fMapClear, "Clear"}, { fMapMenuReverseComplement, "Rev-Comp.."}, { fMapMenuToggleDna, "DNA.."}, { dnaAnalyse, "Analysis.."}, #endif /* not ACEMBLY */ { fMapMenuAddGfSegs, "GeneFind.."}, {0, 0}} ; #define NUM_BUTTON_BOXES UtArraySize(buttonOpts) /* The button menus... */ static MENUOPT fMapZoomOpts[] = { { zoom1, "1 bp/line"}, { zoom10, "10 bp/line"}, { zoom100, "100 bp/line"}, { zoom1000, "1000 bp/line"}, { fMapWhole, "Whole"}, { 0, 0 } } ; static MENUOPT fMapDnaOpts[] = { { fMapMenuToggleDna, "DNA On-Off" }, { fMapMenuToggleTranslation, "Translate Gene On-Off" }, { fMapSetDnaWidth, "Set DNA Width" }, { fMapToggleAutoWidth, "Toggle Auto-DNA-width" }, { 0, 0 } } ; static MENUOPT fMapAnalysisOpts[] = { { menuSpacer, "Built in tools"}, { dnaAnalyse, " Restriction Analysis, codon usage etc." }, { fMapGelDisplay, " Artificial gels" }, { fMapMenuAddGfSegs, " Genefinder (Green)" }, { menuSpacer, "External tools"}, { fMapBlastInit, " BLAST: search external database (Altschul & al)" }, { fMapMenesInit, " BioMotif: search a complex motif (Menessier)" }, { fMapOspInit, " OSP: Oligos Selection Program (Hillier)" }, #ifdef ACEMBLY { dnaCptFindRepeats, " Repeats: Printout repeat structure (Parsons)" }, { abiFixFinish, " Finish: select new reads, StLouis specific (Marth)" }, #endif { 0, 0} } ; static MENUOPT fMapComplementOpts[] = { { fMapMenuReverseComplement, "Reverse-Complement" }, { fMapReverse, "Reverse" }, { fMapMenuComplement, "Complement" }, { fMapMenuComplementSurPlace, "Complement in place" }, { fMapCompHelp, "Help-Explanation" }, { 0, 0} } ; /********************************************************************/ /*********** main menu ************/ /* */ static MENUOPT fMapMenu[] = { { graphDestroy, "Quit" }, { help, "Help" }, { graphPrint, "Print Screen" }, { mapPrint, "Print Whole Map" }, { setDisplayPreserve, "Preserve" }, { MENU_FUNC_CAST togglePreserveColumns, UN_PRESERVE_COLS}, { MENU_FUNC_CAST toggleRedrawColumns, NOREDRAW_COLS}, { MENU_FUNC_CAST toggleReportMismatch, NO_REPORT_MISMATCHES}, { MENU_FUNC_CAST toggleCutSelection, CUT_COORDS}, { menuSpacer, "" }, { fMapMenuRecalculate, "Recalculate" }, { fMapMenuToggleHeader, "Hide Header" }, { menuSpacer, "" }, { fMapMenuDumpSeq, "Export Sequence" }, { fMapMenuDumpSegs, "Export Features" }, { fMapExportTranslations, "Export translations" }, { fMapZoneKeySet, "Make a key set from the active zone" }, { emblDump, "EMBL dump" }, { fMapMenuReadFastaFile, "Import Sequence" }, { menuSpacer, "" }, #ifdef ACEMBLY { mapColControl, "Columns"} , #else /* !ACEMBLY */ { geneStats, "Statistics" }, { fMapAddSegments, "Read Segments" }, { fMapClearSegments, "Clear Segments" }, #endif /* !ACEMBLY */ { 0, 0 } } ; #define FMAPMENU_SIZE sizeof(fMapMenu) char* fMapSegTypeName[] = { "MASTER", "ASSEMBLY_TAG", "SEQUENCE", "SEQUENCE_UP", "CDS", "CDS_UP", "INTRON", "INTRON_UP", "EXON", "EXON_UP", "EMBL_FEATURE", "EMBL_FEATURE_UP", "FEATURE", "FEATURE_UP", "ATG", "ATG_UP", "SPLICE3", "SPLICE3_UP", "SPLICE5", "SPLICE5_UP", "CODING", "CODING_UP", "TRANSCRIPT", "TRANSCRIPT_UP", "SPLICED_cDNA", "SPLICED_cDNA_UP", "VIRTUAL_SUB_SEQUENCE_TAG", "VIRTUAL_SUB_SEQUENCE_TAG_UP", "VIRTUAL_MULTIPLET_TAG", "VIRTUAL_MULTIPLET_TAG_UP", "VIRTUAL_ALIGNED_TAG", "VIRTUAL_ALIGNED_TAG_UP", "VIRTUAL_PREVIOUS_CONTIG_TAG", "VIRTUAL_PREVIOUS_CONTIG_TAG_UP", "HOMOL_GAP", "HOMOL_GAP_UP", "HOMOL", "HOMOL_UP", "PRIMER", "PRIMER_UP", "OLIGO", "OLIGO_UP", "OLIGO_PAIR", "OLIGO_PAIR_UP", "TRANS_SEQ", "TRANS_SEQ_UP", "ALLELE", "ALLELE_UP", "VISIBLE", "DNA_SEQ", "PEP_SEQ", "ORF", "VIRTUAL_PARENT_SEQUENCE_TAG", "VIRTUAL_CONTIG_TAG", "CLONE_END" } ; /**************************************************************************/ /* Verify that a pointer actually points to a FeatureMap struct. */ void uVerifyFeatureMap(FeatureMap fmap, char *caller, char *filename, int line_num) { char *err_site = NULL ; err_site = hprintf(0, "(called from %s() in %s at line %d)", caller, filename, line_num) ; if (!fmap) messcrash("verifyFeatureMap() received NULL FeatureMap pointer %s.", err_site); if (fmap && fmap->magic != &FeatureMap_MAGIC) messcrash("verifyFeatureMap() received non-magic FeatureMap pointer %s.", err_site); messfree(err_site) ; return ; } /* get the FeatureMap struct for the current graph and verify the pointer */ FeatureMap uCurrentFeatureMap(char *caller, char *filename, int line_num) { FeatureMap fmap; char *err_site = NULL ; err_site = hprintf(0, "(called from %s() in %s at line %d)", caller, filename, line_num) ; if (!graphAssFind (&GRAPH2FeatureMap_ASSOC,&fmap)) messcrash("currentFeatureMap() could not find FeatureMap on graph %s.", err_site); uVerifyFeatureMap(fmap, caller, filename, line_num) ; messfree(err_site) ; return fmap; } /* currentFeatureMap */ #ifdef FMAP_DEBUG void fMapPrintLook(FeatureMap look, char *prompt) { if (prompt) printf("------- %s ---------\n", prompt); printf("pleaseRecompute = %d\n", look->pleaseRecompute); printf("start = %d, stop = %d\n", look->start, look->stop); printf("min = %d, max = %d\n", look->min, look->max); printf("centre = %f\n", look->map->centre); } #endif /************** default column specification ************/ static void fMapDefaultColumns (MAP map) { int i = 0 ; while (defaultMapColumns[i].name != NULL && defaultMapColumns[i].func != NULL) { mapInsertCol(map, defaultMapColumns[i].priority, defaultMapColumns[i].isOn, defaultMapColumns[i].name, defaultMapColumns[i].func) ; i++ ; } return ; } /*************** tags and their initialisation ***************/ /* This routine is called whenever an fmap interface routine is called */ /* to make sure that the fmap package is initialised, initialisation is */ /* only done on the first call. */ void fMapInitialise (void) { static int isDone = FALSE ; KEY key ; if (!isDone) { isDone = TRUE ; _Feature = str2tag("Feature") ; _EMBL_feature = str2tag("EMBL_feature") ; _Arg1_URL_prefix = str2tag("Arg1_URL_prefix"); _Arg1_URL_suffix = str2tag("Arg1_URL_suffix") ; _Arg2_URL_prefix = str2tag("Arg2_URL_prefix") ; _Arg2_URL_suffix = str2tag("Arg2_URL_suffix") ; _Spliced_cDNA = str2tag ("Spliced_cDNA") ; _Transcribed_gene = str2tag ("Transcribed_gene") ; _Tm = str2tag ("Tm") ; _Temporary = str2tag ("Temporary") ; if (lexword2key ("Transcribed_gene", &key, _VMainClasses)) _VTranscribed_gene = KEYKEY(key) ; if (lexword2key ("IST", &key, _VMainClasses)) _VIST = KEYKEY(key) ; /* THIS IS GRIM, SURELY THIS SHOULD NOT BE SET HERE ???? */ /* Actually I now know why this is here...the zooming all gets screwed up initially * if its not....this is because all the topmargin gubbins is tied up with * the calculations of size of fmap etc. etc.....horribleness.... * the topMargin is needed BEFORE all the top margin stuff is * drawn, but its then incorrect and this means the magnification is incorrect * which means loads of stuff is in the wrong place. The solution is really to * put all the button stuff in a separate window......*/ #ifdef ACEMBLY topMargin = 6.5 ; #else topMargin = 4.5 ; /* fixes minor bug from Sam */ #endif } return ; } /***************************************************/ /******** Communication toolKit ********************/ static void fMapUnselect (void) { if (selectedfMap) { if (selectedfMap->magic != &FeatureMap_MAGIC) messcrash ("fMapUnselect() received non-magic FeatureMap pointer"); if (selectedfMap->selectBox) { Graph hold = graphActive () ; graphActivate (selectedfMap->graph) ; graphBoxDraw (selectedfMap->selectBox, WHITE, WHITE) ; graphActivate (hold) ; } selectedfMap = NULL ; } return ; } /***************/ static void fMapSelect (FeatureMap look) { if (selectedfMap && selectedfMap != look) fMapUnselect() ; selectedfMap = look ; if (look->selectBox) graphBoxDraw (look->selectBox,BLACK,LIGHTRED) ; return ; } /***************/ /* can be called with NULL adresses for undesired parameters */ BOOL fMapActive(Array *dnap, Array *dnaColp, KEY *seqKeyp, FeatureMap* lookp) { BOOL mapIsActive ; fMapInitialise() ; if (selectedfMap && graphExists(selectedfMap->graph)) { if (selectedfMap->magic != &FeatureMap_MAGIC) messcrash("fMapActive has corrupt selectedfMap->magic"); if (lookp) *lookp = selectedfMap ; if (seqKeyp) *seqKeyp = selectedfMap->seqKey ; if (dnap) *dnap = selectedfMap->dna ; if (dnaColp) *dnaColp = selectedfMap->colors ; mapIsActive = TRUE ; } else { selectedfMap = 0 ; if (lookp) *lookp = 0 ; if (seqKeyp) *seqKeyp = 0 ; if (dnap) *dnap = 0 ; if (dnaColp) *dnaColp = 0 ; mapIsActive = FALSE ; } return mapIsActive ; } Graph fMapActiveGraph (void) { FeatureMap fmap; fMapInitialise() ; return (fMapActive (0,0,0,&fmap) ? fmap->graph : 0) ; } BOOL fMapActivateGraph (void) { FeatureMap fmap; return (fMapActive (0,0,0,&fmap) && graphActivate (fmap->graph)) ; } MAP fMapGetMap (FeatureMap look) /* to get private member */ { return look->map ; } /*****************************************************************/ static void zoom1 (void) { FeatureMap look = currentFeatureMap("zoom1") ; look->map->mag = 1 ; fMapDraw(look, 0) ; } static void zoom10 (void) { FeatureMap look = currentFeatureMap("zoom10") ; look->map->mag = 0.1 ; fMapDraw(look, 0) ; } static void zoom100 (void) { FeatureMap look = currentFeatureMap("zoom100") ; look->map->mag = 0.01 ; fMapDraw(look, 0) ; } static void zoom1000 (void) { FeatureMap look = currentFeatureMap("zoom1000") ; look->map->mag = 0.001 ; fMapDraw(look, 0) ; } static void fMapWhole (void) /* simply show whole object */ { FeatureMap look = currentFeatureMap("fMapWhole") ; if (!look->length) return ; look->map->centre = look->length / 2 ; look->map->mag = (mapGraphHeight - topMargin - 5) / (1.05 * look->length) ; fMapDraw (look, 0) ; return ; } #ifdef ACEMBLY extern MENUOPT fMapAcemblyOpts[] ; static void fMapTraceJoin (void) { FeatureMap look = currentFeatureMap("fMapTraceJoin") ; messout("This is a menu. Please use the right mouse button.") ; } #else /* not ACEMBLY */ static void toggleColumnButtons (void) { FeatureMap look = currentFeatureMap("toggleColumnButtons") ; look->flag ^= FLAG_COLUMN_BUTTONS ; fMapDraw (look, 0) ; } #endif /* not ACEMBLY */ /* Set fmap so that it is not reused again, in this case the menu items for */ /* "Preserve" and "Preserve Cols" are now useless as the fmap will stay */ /* preserved for ever more, so they are removed from the main menu. */ static void setDisplayPreserve(void) { FeatureMap look = currentFeatureMap("setDisplayPreserve") ; BOOL found ; found = removeFromMenu( MENU_FUNC_CAST togglePreserveColumns, look->mainMenu) ; found = removeFromMenu(setDisplayPreserve, look->mainMenu) ; displayPreserve() ; fMapDraw (look, 0) ; return ; } /* Turns on/off whether DNA mismatches get reported or not, needed because */ /* sometimes user knows that there are lots of mismatches and doesn't want */ /* to see them because they are fixing up a link object or some such. */ static void toggleReportMismatch(MENUITEM unused) { FeatureMap look = currentFeatureMap("toggleReportMismatch") ; MENUOPT *m = look->mainMenu; look->reportDNAMismatches = !(look->reportDNAMismatches) ; toggleTextInMenu(m, toggleReportMismatch, look->reportDNAMismatches ? NO_REPORT_MISMATCHES : REPORT_MISMATCHES) ; graphMenu (look->mainMenu) ; return; } /* Turns on/off whether the whole fmap gets redrawn each time user clicks on */ /* a column to turn it on/off or not. */ /* */ static void toggleRedrawColumns(MENUITEM unused) { FeatureMap look = currentFeatureMap("toggleRedrawColumns") ; MENUOPT *m = look->mainMenu; look->redrawColumns = !(look->redrawColumns) ; toggleTextInMenu(m, toggleRedrawColumns, look->redrawColumns ? NOREDRAW_COLS : REDRAW_COLS) ; graphMenu (look->mainMenu) ; return; } /* Changes what fmap puts in the PRIMARY selection buffer, by default its */ /* DNA but can be toggled to Coordinate information. */ /* */ static void toggleCutSelection(MENUITEM unused) { FeatureMap look = currentFeatureMap("toggleCutSelection") ; MENUOPT *m = look->mainMenu; char *menu_txt = NULL ; if (look->cut_data == FMAPCUT_DNA) { look->cut_data = FMAPCUT_COORDS ; menu_txt = CUT_DNA ; } else { look->cut_data = FMAPCUT_DNA ; menu_txt = CUT_COORDS ; } toggleTextInMenu(m, toggleCutSelection, menu_txt) ; graphMenu (look->mainMenu) ; return; } /* Note that toggling the preserve columns options does NOT affect the */ /* current column settings, it only affects whether the _next_ fmap shown */ /* will keep the same columns. */ static void togglePreserveColumns(MENUITEM unused) { FeatureMap look = currentFeatureMap("togglePreserveColumns") ; MENUOPT *m = look->mainMenu; look->preserveColumns = !(look->preserveColumns) ; toggleTextInMenu(m, togglePreserveColumns, look->preserveColumns ? UN_PRESERVE_COLS : PRESERVE_COLS) ; graphMenu (look->mainMenu) ; return; } /* toggle menu text - pretty ugly this way, but easiest way using the MENUOPT system (new menu system not final yet) */ static void toggleTextInMenu(MENUOPT *m, MENUFUNCTION menu_func, char *menu_text) { while (m->text && m->f) { if (m->f == (MenuBaseFunction)menu_func) { m->text = menu_text ; break ; } ++m; } return ; } /****************************/ static void fMapShiftOrigin(void) { KEY key ; int i,x0 ; FeatureMap look = currentFeatureMap("fMapShiftOrigin") ; fMapSelect(look) ; if (lexword2key (look->originBuf, &key, _VSequence)) { for (i = 1 ; i < arrayMax(look->segs) ; ++i) if (arrp(look->segs,i,SEG)->key == key) { look->origin = arrp(look->segs,i,SEG)->x1 ; goto done ; } messout ("Sorry, object %s is not being displayed in this graph", name(key)) ; return ; } else if (sscanf(look->originBuf, "%d", &x0)) { look->origin += x0 - 1 ; /* -1 Plato dixit */ strcpy (look->originBuf, "1") ; } else { messout ("Give either a sequence name or a position in the " "current units.") ; return ; } done: fMapDraw(look,0) ; /* since i want to rewrite the coordinate system */ } static BOOL fMapSetOriginFromKey (FeatureMap look, KEY key) { int i ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) if (arrp(look->segs,i,SEG)->key == key) { if (look->flag & FLAG_REVERSE) /* top of object is base 1 */ look->origin = look->length - arrp(look->segs,i,SEG)->x2 - 1 ; else look->origin = arrp(look->segs,i,SEG)->x1 ; return TRUE ; } return FALSE ; } static void fMapDoPickOrigin (KEY key) { FeatureMap look = currentFeatureMap("fMapDoPickOrigin") ; fMapSelect(look) ; if (fMapSetOriginFromKey (look, key)) fMapDraw(look,0) ; /* since i want to rewrite the coordinate system */ else messout ("Sorry, %s is not being displayed in this window", name(key)) ; } static void fMapOriginButton (void) { graphRegister (MESSAGE_DESTROY, displayUnBlock) ; /*mhmp 15.05.98*/ displayBlock (fMapDoPickOrigin, "Pick an object on this window\n" "to set the first base of its sequence as base 1.\n" "Remove this message to cancel.") ; } /**********************************/ static BOOL setZone (FeatureMap look, int newMin, int newMax) { BOOL res = FALSE ; look->zoneMin = newMin ; if (look->zoneMin < 0) { look->zoneMin = 0 ; res = TRUE ; } look->zoneMax = newMax ; if (look->zoneMax > look->length) { look->zoneMax = look->length ; res = TRUE ; } return res ; } static void setZoneUserCoords (FeatureMap look, int x0, int x1) { int tmp ; if (look->flag & FLAG_REVERSE) { tmp = x0 ; x0 = look->length - look->origin - x1 ; x1 = look->length - look->origin - tmp + 1 ; } else { x0 = x0 + look->origin - 1 ; x1 = x1 + look->origin ; } if (x0 > x1) { tmp = x0 ; x0 = x1 ; x1 = tmp ;} if (setZone (look, x0, x1)) messout ("The requested zone extends beyond the sequence " "and is being clipped") ; if (look->flag & FLAG_REVERSE) strncpy (look->zoneBuf, messprintf ("%d %d", COORD(look, look->zoneMax)+1, COORD(look, look->zoneMin)), 15) ; else strncpy (look->zoneBuf, messprintf ("%d %d", COORD(look, look->zoneMin), COORD(look, look->zoneMax)-1), 15) ; if (look->zoneBox) graphBoxDraw (look->zoneBox, -1, -1) ; fMapReDrawDNA (look) ; } void fMapSetZoneEntry (char *zoneBuf) { int x0, x1 ; FeatureMap look = currentFeatureMap("fMapSetZoneEntry") ; fMapSelect(look) ; if (sscanf(zoneBuf,"%d %d", &x0, &x1) != 2) { messout ("The zone must be a pair of integers") ; return ; } setZoneUserCoords (look, x0, x1) ; /* mhmp 02.06.98 & 08.09.98 */ graphEntryDisable () ; graphRegister (KEYBOARD, fMapKbd) ; } static void fMapDoSetZone (KEY key) { int i ; FeatureMap look = currentFeatureMap("fMapDoSetZone") ; fMapSelect(look) ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) if (arrp(look->segs,i,SEG)->key == key) { setZone (look, arrp(look->segs,i,SEG)->x1, arrp(look->segs,i,SEG)->x2 + 1) ; fMapDraw(look,0) ; return ; } messout ("Sorry, sequence %s is not being displayed in this window", name(key)) ; } static void fMapSetZoneButton (void) { graphRegister (MESSAGE_DESTROY, displayUnBlock) ; /*mhmp 15.05.98*/ displayBlock (fMapDoSetZone, "Pick an object in this window\n" "to set the limits of the active zone.\n" "Remove this message to cancel.") ; } /**********************************/ BOOL fMapFindZone (FeatureMap look, int *min, int *max, int *origin) { verifyFeatureMap(look, fMapFindZone) ; if (min) *min = look->zoneMin ; if (max) *max = look->zoneMax ; if (origin) *origin = look->origin ; return TRUE ; } /***********************************/ /* Note that there are two functions here : fMapFindSpan - find enclosing sequence fMapFindWriteableSpan - like fMapFindSpan, but sequence found must be either Genomic_canonical or have the Link tag set - this ensures that it is suitable for adding extra data to. */ /* Allowable results must be class sequence, not have exons and must be Genomic_canonical or have the Link tag set. */ static BOOL pred(KEY key) { if (class(key) != _VSequence) return FALSE; if (bIndexTag(key, str2tag("Source_exons"))) return FALSE; return TRUE; } /* fmap coords run 0 -> (length - 1), smap coords run 1 -> length, so we * must correct input and output coords here. */ BOOL fMapFindSpan(FeatureMap look, KEY *key_out, int *start_out, int *end_out) { BOOL result ; int start, end ; messAssert(key_out != NULL && start_out != NULL && end_out != NULL) ; start = *start_out + 1 ; end = *end_out + 1 ; if ((result = sMapFindSpan(look->smap, key_out, &start, &end, pred))) { *start_out = start - 1 ; *end_out = end - 1 ; } return result ; } static BOOL predWrite(KEY key) { if (class(key) != _VSequence) return FALSE; if (bIndexTag(key, str2tag("Source_exons"))) return FALSE; if (bIndexTag(key, str2tag("Genomic_canonical"))) return TRUE; if (bIndexTag(key, str2tag("Link"))) return TRUE; return FALSE; } /* fmap coords run 0 -> (length - 1), smap coords run 1 -> length, so we * must correct input and output coords here. */ BOOL fMapFindSpanWriteable(FeatureMap look, KEY *key_out, int *start_out, int *end_out) { BOOL result ; int start, end ; verifyFeatureMap(look, fMapFindSpanWriteable) ; messAssert(key_out != NULL && start_out != NULL && end_out != NULL) ; start = *start_out + 1 ; end = *end_out + 1 ; if ((result = sMapFindSpan(look->smap, key_out, &start, &end, predWrite))) { if (start > end) { *start_out = start - 1 ; *end_out = end - 1 ; } else { *start_out = start + 1 ; *end_out = end + 1 ; } } return result ; } /***********************************/ /* Returns the bottom position of the header. */ static float drawHeader(FeatureMap look) { int box, mainButs ; float bottom = 0 ; look->selectBox = graphBoxStart () ; graphText ("Selected DNA",1.0,.5) ; graphBoxEnd () ; if (look == selectedfMap) graphBoxDraw (look->selectBox,BLACK,LIGHTRED) ; else graphBoxDraw (look->selectBox,WHITE,WHITE) ; graphButton ("Origin:", fMapOriginButton, 16.5, .4) ; if (!strlen(look->originBuf)) strcpy (look->originBuf,"1") ; look->originBox = graphTextEntry (look->originBuf, 14, 25, .5, TEXT_ENTRY_FUNC_CAST fMapShiftOrigin) ; graphButton ("Active zone:", fMapSetZoneButton, 41, .4) ; if (look->flag & FLAG_REVERSE) strncpy (look->zoneBuf, messprintf ("%d %d", COORD(look, look->zoneMax)+1, COORD(look, look->zoneMin)), 15) ; else strncpy (look->zoneBuf, messprintf ("%d %d", COORD(look, look->zoneMin), COORD(look, look->zoneMax)-1), 15) ; look->zoneBox = graphTextScrollEntry (look->zoneBuf, 23, 14, 54.5, .5, fMapSetZoneEntry) ; /* mhmp 02.06.98 + 08.09.98 */ graphEntryDisable() ; graphRegister (KEYBOARD, fMapKbd) ; if (look->flag & FLAG_COMPLEMENT) graphText ("COMPLEMENT", 70, 0) ; if (look->flag & FLAG_COMPLEMENT_SUR_PLACE) { graphText ("COMP_DNA", 81, 0) ; graphText ("IN_PLACE", 81, 1) ; } if (((look->flag & FLAG_COMPLEMENT) && !(look->flag & FLAG_REVERSE)) || (!(look->flag & FLAG_COMPLEMENT) && (look->flag & FLAG_REVERSE))) graphText ("REVERSED", 70, 1) ; look->segBox = graphBoxStart() ; graphTextPtr (look->segTextBuf, 2, 2, 126) ; graphBoxEnd() ; graphBoxDraw (look->segBox, BLACK, PALEBLUE) ; /* default is to put dna in cut buffer, alternative is coord data. */ if (look->cut_data == FMAPCUT_COORDS) graphPostBuffer(look->segTextBuf) ; /* Draw fmap main control buttons and find out how far down they come. */ mainButs = graphBoxStart() ; box = graphButtons(buttonOpts, 1, 3.2, mapGraphWidth) ; graphBoxEnd() ; graphBoxDim(mainButs, 0, 0, 0, &bottom) ; /* menus for main buttons. */ #ifdef ACEMBLY graphBoxMenu (box + getPositionInMenu(mapZoomIn, buttonOpts), fMapZoomOpts) ; graphBoxMenu (box + getPositionInMenu(mapZoomOut, buttonOpts), fMapZoomOpts) ; graphBoxMenu (box + getPositionInMenu(dnaAnalyse, buttonOpts), fMapAnalysisOpts) ; graphBoxMenu (box + getPositionInMenu(fMapTraceJoin, buttonOpts), fMapAcemblyOpts) ; if (look->isGeneFinder) graphBoxMenu (box + getPositionInMenu(fMapMenuAddGfSegs, buttonOpts), fMapGeneOpts) ; else graphBoxClear (box + getPositionInMenu(fMapMenuAddGfSegs, buttonOpts)) ; #else mapColMenu (box) ; graphBoxMenu (box + getPositionInMenu(mapZoomIn, buttonOpts), fMapZoomOpts) ; graphBoxMenu (box + getPositionInMenu(mapZoomOut, buttonOpts), fMapZoomOpts) ; graphBoxMenu (box + getPositionInMenu(fMapMenuReverseComplement, buttonOpts), fMapComplementOpts) ; graphBoxSetMenu (box + getPositionInMenu(fMapClear, buttonOpts), FALSE); /* no menu */ graphBoxMenu (box + getPositionInMenu(fMapMenuToggleDna, buttonOpts), fMapDnaOpts) ; graphBoxMenu (box + getPositionInMenu(dnaAnalyse, buttonOpts), fMapAnalysisOpts) ; graphBoxMenu (box + getPositionInMenu(fMapMenuAddGfSegs, buttonOpts), fMapGeneOpts) ; #endif look->minLiveBox = 1 + graphBoxStart() ; graphBoxEnd() ; return(bottom) ; } /********************************************************************/ /* The menu stuff is all being changed, really I need a menu function to do */ /* these things for me...... */ /* This routine just returns the position (starting at 0) of the menu item */ /* with a given menu function. It saves us from having to hard code a load of*/ /* menu positions. */ /* If it can't find the item it returns -1 */ static int getPositionInMenu(MenuBaseFunction func, MENUOPT *menu) { int position = -1 ; int i = 0 ; while(position == -1 && menu[i].f != NULL) { if (menu[i].f == func) position = i ; else i++ ; } return position ; } /* This function removes an item from a menu by finding it and shuffling up */ /* all the following items by one place. */ /* If it can't find the item it returns FALSE. */ static BOOL removeFromMenu(MenuBaseFunction func, MENUOPT *menu) { BOOL found = FALSE ; int i = 0 ; /* Find the position of the item. */ while(found == FALSE && menu[i].f != NULL) { if (func == menu[i].f) found = TRUE ; else i++ ; } /* Shuffle up the entries. */ if (found == TRUE) { while(menu[i].text != NULL) { menu[i].f = menu[i + 1].f ; menu[i].text = menu[i + 1].text ; i++ ; } } return found ; } /********************************************************************/ /* menu functions that manipulate the fMap look of the active graph */ /********************************************************************/ void fMapMenuToggleDna (void) /* public - also called from fmaptrace.c */ { FeatureMap look = currentFeatureMap("fMapMenuToggleDNA") ; fMapToggleDna (look) ; fMapDraw (look,0) ; } static void fMapMenuToggleTranslation (void) { FeatureMap look = currentFeatureMap("fMapMenuToggleTranslation") ; if (fMapToggleTranslation(look)) fMapDraw(look,0) ; } static void fMapMenuToggleHeader (void) { FeatureMap look = currentFeatureMap("fMapMenuToggleHeader") ; fMapToggleHeader(look); fMapDraw (look, 0) ; } void fMapMenuAddGfSegs (void) /* public - also called from fmapgene.c */ { FeatureMap look = currentFeatureMap("fMapMenuAddGfSegs") ; fMapAddGfSegs(look); /* includes fMapDraw() */ } static void fMapMenuDumpSeq (void) { dnaAnalyse() ; dnacptDontUseKeySet () ; dnacptFastaDump() ; } static void fMapMenuReadFastaFile (void) { KEY key ; FILE *fil ; int c ; ACEIN fi = 0 ; BOOL ok = FALSE ; static char directory[DIR_BUFFER_SIZE] = ""; static char filename[FIL_BUFFER_SIZE] = ""; char *unused ; if (!(fil = filqueryopen (directory, filename, "", "r", "DNA sequence file in plain or FASTA format"))) return ; c = getc(fil) ; ungetc(c,fil) ; /* peek at first character */ filclose (fil); fi = aceInCreateFromFile (messprintf("%s/%s", directory, filename), "r", "", 0); if (c == '>') aceInCard (fi) ; /* throw out title line */ lexaddkey ("temp_seq", &key, _VDNA) ; if (dnaParse(fi, key, &unused) == PARSEFUNC_OK) ok = TRUE ; aceInDestroy (fi) ; if (ok) { lexaddkey ("temp_seq", &key, _VSequence) ; /* now display it */ display (key, 0, "FMAP") ; } return; } /* fMapMenuReadFastaFile */ /**********************************************************/ static void fMapToggleDna (FeatureMap look) { mapColSetByName (look->map, "DNA Sequence", 2) ; /* toggle */ mapColSetByName (look->map, "Coords", /* set to same state */ mapColSetByName (look->map, "DNA Sequence", -1)) ; return; } /* fMapToggleDna */ static BOOL fMapToggleTranslation (FeatureMap look) { int min, max ; Array cds, index ; if (fMapGetCDS (look, look->translateGene, &cds, &index)) { min = arr(index,0,int) - look->min ; max = arr(index,arrayMax(index)-1,int) - look->min ; if (min > max) mapColSetByName (look->map, "Up Gene Translation", 2) ; else mapColSetByName (look->map, "Down Gene Translation", 2) ; return TRUE; } messout ("Select a gene to translate with the " "pulldown menu on the genes") ; return FALSE; } /* fMapToggleTranslation */ static void fMapToggleHeader (FeatureMap look) { MENUOPT *m = look->mainMenu; look->flag ^= FLAG_HIDE_HEADER ; /* toggle menu text - pretty ugly this way, but easiest way * using the MENUOPT system (new menu system not final yet) */ while (m->text && m->f) { if (m->f == fMapMenuToggleHeader) { if (look->flag & FLAG_HIDE_HEADER) m->text = "Show Header"; else m->text = "Hide Header"; } ++m; } return; } /* fMapToggleHeader */ /****************************************************************/ void fMapPleaseRecompute (FeatureMap look) /* set look->pleaseRecompute */ { if (!look) return ; fMapInitialise() ; if (look->magic == &FeatureMap_MAGIC && graphExists (look->graph)) look->pleaseRecompute = TRUE ; /* trace.c may called using a stale look, i do not want to crash on that, the really correct solution would be inside fmapdestroy to deregister this look from the trace,c code, i do not think it is worth the effort */ } /**********************************************************************/ static int fromTrace = 0, fromTracePos = 0 ; /* specialised entry point */ void fMapDrawFromTrace (FeatureMap look, KEY from, int xx, int type) { fMapInitialise() ; fromTrace = type ; fromTracePos = xx ; fMapDraw (look, from) ; } static void fMapDrag(float y) { FeatureMap look = currentFeatureMap("drag"); if (look && look->segBox) { sprintf(look->segTextBuf, "%d", (int)GRAPH2MAP(fMapGetMap(look),y)); graphBoxDraw (look->segBox, BLACK, PALEBLUE) ; } } void fMapDraw (FeatureMap look, KEY from) { SEG *seg = 0 ; char *cp ; int i, fromIndex = 0 ; float absMag, x1, x2 ; fMapInitialise() ; if (look && look->map && graphActivate(look->graph)) graphPop() ; else return ; /* If user toggled a column on/off and they don't want the whole fmap */ /* redrawn, then just redraw the toggled column button. */ if (mapColButPressed(look->map) && !(look->redrawColumns)) { /* Only redraw the column button if column buttons are displayed. */ if (look->flag & FLAG_COLUMN_BUTTONS) mapRedrawColBut(look->map) ; return ; } absMag = (look->map->mag > 0) ? look->map->mag : -look->map->mag ; chrono("fMapDraw") ; restart: if (!class(from) && look->activeBox && look->activeBox < arrayMax(look->boxIndex) && arr(look->boxIndex, look->activeBox, int) < arrayMax(look->segs) && BOXSEG(look->activeBox)->type != DNA_SEQ) { fromIndex = arr(look->boxIndex, look->activeBox, int) ; from = BOXSEG(look->activeBox)->parent ; } /* clear all DNA highlighting */ if (arrayExists (look->colors)) { cp = arrp(look->colors, 0, char); for( i=0; icolors); i++) *cp++ &= ~(TINT_HIGHLIGHT1 | TINT_HIGHLIGHT2); } graphWhere(&x1, 0, &x2, 0); graphClear () ; graphGoto((x2+x1)/2.0, 0.0); look->summaryBox = 0 ; graphFitBounds (&mapGraphWidth, &mapGraphHeight) ; /* Need to set up our very own look menu here.... */ graphMenu (look->mainMenu) ; /* if(!look->isOspSelecting) */ graphRegister (PICK,(GraphFunc) fMapPick) ; arrayMax(look->segs) = look->lastTrueSeg ; look->boxIndex = arrayReCreate(look->boxIndex,200,int) ; look->visibleIndices = arrayReCreate(look->visibleIndices,64,int) ; if (look->flag & FLAG_HIDE_HEADER) { look->selectBox = look->originBox = look->zoneBox = look->segBox = 0 ; look->minLiveBox = 1 ; topMargin = 0 ; } else { topMargin = drawHeader (look) ; memset (look->segTextBuf, 0, 63) ; } halfGraphHeight = 0.5 * (mapGraphHeight - topMargin) ; if (class(from) == _VCalcul) { i = KEYKEY (from) ; look->map->centre = i ; from = 0 ; } else if (fromTrace && arrayExists(look->segs)) { for (i = 0, seg = arrp(look->segs,0,SEG) ; i < arrayMax(look->segs) ; ++i, ++seg) if (from == seg->key) { look->map->centre = seg->x1 + fromTracePos ; break ; } if (i == arrayMax(look->segs)) look->pleaseRecompute = 1 ; else if (fromTrace == 1) { from = 0 ; fromIndex = 0 ; } /* no highlighting */ } fromTrace = 0 ; /* I guess this might happen as a result of one of the database reads from */ /* above but should never happen as a result of user interaction. */ /* mhmp 05.02.98 */ if (look->map->centre < look->map->min) { #ifdef FMAP_DEBUG printf("Clamping centre.\n"); #endif look->map->centre = look->map->min ; } if (look->map->centre > look->map->max) { look->map->centre = look->map->max ; #ifdef FMAP_DEBUG printf("Clamping centre.\n"); #endif } look->min = look->map->centre - (halfGraphHeight-2.5)/absMag ; look->max = look->map->centre + (halfGraphHeight-2.5)/absMag ; if (look->pleaseRecompute != -1) { #ifdef FMAP_DEBUG fMapPrintLook(look, "In fMapDraw"); #endif if (look->flag & FLAG_COMPLEMENT) { if (look->pleaseRecompute || (look->min < 0 && look->stop < look->fullLength-1) || (look->max > look->length && look->start)) { look->pleaseRecompute = -1 ; if (fMapDoRecalculate (look, look->min, look->max)) goto restart ; else { messerror("Could not recalculate, Fmap will be removed.") ; graphDestroy() ; return ; } } } else { if (look->pleaseRecompute || (look->min < 0 && look->start) || (look->max > look->length && look->stop < look->fullLength-1)) { look->pleaseRecompute = -1 ; if (fMapDoRecalculate (look, look->min, look->max)) goto restart ; else { messerror("Could not recalculate, Fmap will be removed.") ; graphDestroy() ; return ; } } } } look->pleaseRecompute = FALSE ; if (look->min < 0) look->min = 0 ; while (look->min % 3) /* establish standard frame */ ++look->min ; if (look->max >= look->length) look->max = look->length ; look->map->fixFrame = &look->min ; /* fiddle for column drawing */ seqDragBox = 0 ;/* mhmp 07.10.98 */ /* The guts of column drawing, this map routine call draws all the columns in the fmap. */ graphTextBounds(mapDrawColumns (look->map), 0) ; /* Display column control buttons...this should be in a separate window. */ if (look->flag & FLAG_COLUMN_BUTTONS) { int first_box; first_box = mapDrawColumnButtons (look->map) ; /* switch off background menu over these column buttons */ for (i = 0; i < arrayMax(look->map->cols); ++i) graphBoxMenu(first_box + i, NULL); } look->activeBox = 0 ; if (fromIndex) for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; ++i) { if (arr(look->boxIndex, i, int) == fromIndex) { fMapSelectBox (look, i, 0, 0) ; break ; } } else if (class(from)) for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; ++i) { if (BOXSEG(i)->key == from) { fMapSelectBox (look, i, 0, 0) ; break ; } } chrono("f-graphRedraw") ; graphRedraw() ; fMapSelect(look) ; chronoReturn() ; return ; } /* fmapDraw */ /************************************************************/ /****************** Registered routines *********************/ void fMapDestroy (FeatureMap look) { verifyFeatureMap(look, fMapDestroy) ; handleDestroy (look->handle) ; arrayDestroy (look->dna) ; arrayDestroy (look->dnaR) ; arrayDestroy (look->sites) ; /* sites come from dnacpt.c */ stackDestroy(look->siteNames) ; arrayDestroy (look->multiplets) ; arrayDestroy (look->oligopairs) ; mapDestroy (look->map) ; /* If this fmap is the "current" (i.e. non-preserved fmap), then remove it from the display * systems list of "current" displays. */ if (look->graph == displayGetGraph("FMAP")) displaySetGraph("FMAP", GRAPH_NULL) ; if (graphActivate(look->graph)) { graphAssRemove (&GRAPH2FeatureMap_ASSOC) ; graphAssRemove (&MAP2LOOK_ASSOC) ; graphRegister (DESTROY, 0) ; /* deregister */ } fMapTraceDestroy (look) ; fMapOspDestroy (look) ; /* Sometimes fmap is called to display a temporarily created object, e.g. */ /* cDNA display. So now we are finished we need to delete the object. */ if (look->destroy_obj) { OBJ obj ; if ((obj = bsUpdate(look->seqKey)) != NULL) bsKill(obj) ; } if (graphActivate (look->selectGraph)) /* from fmapgene.c */ graphDestroy () ; if (look == selectedfMap) selectedfMap = 0 ; look->magic = 0 ; messfree (look) ; return; } /* fMapDestroy */ static void fMapDestroyCallback (void) /* callback */ { FeatureMap look = currentFeatureMap("fMapDestroy") ; fMapDestroy (look) ; } /*************************************************************/ /* Sets the start and stop coords for the selected sequence and returns the root * parent sequence. * Can be called in two ways: * if look is NULL, then just works off the start/end coords and the key, * this happens when we first create an fmap of a sequence. * else tries to keep the active zone etc. in scope for recalculations of. * an existing fmap. * Note that the start/stop can be (n,0) or (0,n) meaning "the sequence from * "n" to the end. */ static BOOL setStartStop(KEY seq, int *start_inout, int *stop_inout, KEY *key_out, FeatureMap look) { BOOL result = FALSE ; int start, stop, length ; BOOL swop ; KEY *key_ptr ; start = *start_inout ; stop = *stop_inout ; if (key_out != NULL) key_ptr = key_out ; else key_ptr = NULL ; /* Only check coords if they are _not_ in form (n,0), in this form they can be passed * straight through to sMapTreeRoot(). */ if (!(start == 0) && !(stop == 0)) { swop = (start > stop) ; /* ascending order for easier arithmetic. */ if (swop) { int tmp ; tmp = start, start = stop, stop = tmp ; } if (!look) { /* If we are constructing a new look, make it a decent length otherwise user may be * given a ridiculously tiny map. */ length = stop - start + 1 ; if (length < 10000) length = 10000 ; start -= length ; stop += length ; } else { /* Logic is: keep length the same or bigger, if start/stop are within an * already calculated bit of sequence then just keep all coords the same, * otherwise slide coords up to be centred around new start/stop. */ length = stop - start + 1 ; if (length < look->length) length = look->length ; if ((start >= (look->start + 1) && stop <= (look->stop + 1))) { start = look->start + 1 ; stop = look->stop + 1 ; } else { int extend = (length - (stop - start + 1)) / 2 ; start -= extend ; stop += extend ; } } /* smap coords start at 1 so make sure they do, no need to clamp the stop value, this is * done by sMapTreeRoot() */ if (start < 1) start = 1 ; if (swop) { int tmp ; tmp = start, start = stop, stop = tmp ; } } /* Ok now find the actual start/stop for the sequence. */ if (sMapTreeRoot(seq, start, stop, key_ptr, &start, &stop)) { result = TRUE ; *start_inout = start ; *stop_inout = stop ; } return result ; } static BOOL fMapDoRecalculate(FeatureMap look, int start, int stop) { BOOL result ; int diff, oldMax ; BOOL isComplement = (look->flag & FLAG_COMPLEMENT) != 0 ; /* NB. At this point start and stop are in the co-ordinate system if the calculated DNA. Below they get moved into the co-ordinate system of the root object. The same variables are used to store both of these. Magic. */ #ifdef FMAP_DEBUG fMapPrintLook(look, "Before recalc."); #endif arrayDestroy (look->dna) ; arrayDestroy (look->dnaR) ; fMapTraceDestroy (look) ; fMapOspDestroy (look) ; #ifdef USE_SMAP sMapDestroy(look->smap); #endif /* NB start and stop arguments are in the co-ordinate system of the current fMap: This code converts them to the co-ordinate system of the root object. */ if (isComplement) { start = look->stop + 2 - start ; stop = look->stop + 2 - stop ; } else { start += look->start ; stop += look->start ; } #ifdef USE_SMAP if (!setStartStop(look->seqKey, &start, &stop, NULL, look)) { messerror("Could not find start/stop of " "key: %s, start: %d stop: %d", name(look->seqKey), start, stop) ; return FALSE ; } look->smap = sMapCreate(look->handle, look->seqKey, start, stop, NULL); reportErrors = look->reportDNAMismatches; if (!(look->dna = sMapDNA(look->smap, 0, callback))) #else if (!(look->dna = fMapFindDNA(look->seqKey, &start, &stop, look->reportDNAMismatches))) #endif { /* If we wander completely away from the DNA, sMapDNA returns NULL. Since we need something to be able to darw a screen fake it up here. This is not, pretty, but it works. */ look->dna = arrayHandleCreate(abs(stop - start) + 1, unsigned char, 0); arrayMax(look->dna) = abs(stop - start); } /* sMapDNA might not get full range. */ look->length = arrayMax(look->dna) ; stop = start + look->length + 1; look->colors = arrayReCreate (look->colors, arrayMax(look->dna), char) ; arrayMax(look->colors) = arrayMax(look->dna) ; oldMax = look->map->max; look->map->max = arrayMax(look->dna) ; if (isComplement) diff = (start-1) - look->stop ; else diff = look->start - (start-1) ; look->map->centre += diff ; look->zoneMin +=diff ; look->zoneMax +=diff ; { /* Change offsets in choosen and antichosen arrays */ Array vA = arrayCreate(40, char *); Array v1A = arrayCreate(40, char *); char *v, *v1; int i; #ifdef FMAP_DEBUG printf("Munging assocs.\n"); #endif v = 0; while (assNext (look->chosen, &v, &v1)) { array(vA, arrayMax(vA), char *) = v; array(v1A, arrayMax(v1A), char *) = v1; #ifdef FMAP_DEBUG printf("Chosen: %s: %d,%d\n", fMapSegTypeName[UNHASH_TYPE(v)], assInt(v1), UNHASH_X2(v1)); #endif } look->chosen = assReCreate(look->chosen); for (i = 0 ; i < arrayMax(vA) ; ++i) { v = array(vA, i, char *); v1 = array(v1A, i, char *); assInsert(look->chosen, HASH(UNHASH_TYPE(v), UNHASH_X2(v) + diff), assVoid(assInt(v1)+diff)); } #ifdef FMAP_DEBUG v = 0; while (assNext (look->chosen, &v, &v1)) printf("Chosen-out: %s: %d,%d\n", fMapSegTypeName[UNHASH_TYPE(v)], assInt(v1), UNHASH_X2(v1)); #endif arrayMax(vA) = 0; arrayMax(v1A) = 0; v = 0; while (assNext (look->antiChosen, &v, &v1)) { array(vA, arrayMax(vA), char *) = v; array(v1A, arrayMax(v1A), char *) = v1; } look->antiChosen = assReCreate(look->antiChosen); for (i = 0 ; i < arrayMax(vA) ; ++i) { v = array(vA, i, char *); v1 = array(v1A, i, char *); assInsert(look->antiChosen, HASH(UNHASH_TYPE(v), UNHASH_X2(v) + diff), assVoid(assInt(v1)+diff)); } arrayDestroy(vA); arrayDestroy(v1A); } if (look->flag & FLAG_REVERSE) look->origin = look->length + 1 - look->origin ; look->length = arrayMax(look->dna) ; look->origin += diff ; if (look->flag & FLAG_REVERSE) look->origin = look->length + 1 - look->origin ; if (isComplement) { look->start = stop-1 ; look->stop = start-1 ; } else { look->start = start-1 ; look->stop = stop-1 ; } if (!fMapConvert(look)) result = FALSE ; else { look->activeBox = 0 ; seqDragBox = 0 ; /* since segs array reorganised */ result = TRUE ; } #ifdef FMAP_DEBUG fMapPrintLook(look, "After recalc."); #endif return result ; } static void fMapMenuRecalculate (void) { int start, stop ; float absMag ; KEY curr ; FeatureMap look = currentFeatureMap("fMapMenuRecalculate") ; if (look->activeBox) { curr = BOXSEG(look->activeBox)->parent ; if (!class(curr)) curr = 0 ; /* avoids key being a tag */ } else curr = 0 ; absMag = (look->map->mag > 0) ? look->map->mag : -look->map->mag ; start = look->map->centre - (halfGraphHeight-2)/absMag ; stop = look->map->centre + (halfGraphHeight-2)/absMag ; if (fMapDoRecalculate (look, start, stop)) fMapDraw (look, curr) ; else { messerror("Could not recalculate, Fmap will be removed.") ; graphDestroy() ; } return; } /* fMapMenuRecalculate */ /***********************************/ static double oldx, oldy, oldx2, oldy2, oldx3, oldx4 ; static float newx2, newy2 , xdeb, xfin, xlimit ; static float newx, newy ; /*mhmp 22.04.98 */ static int oldbox ; static FeatureMap oldlook ; static void fMapLeftDrag (double x, double y) { graphXorLine (0, oldy, mapGraphWidth, oldy) ; oldy = y ; graphXorLine (0, y, mapGraphWidth, y) ; } static void frameDNA (void) { FeatureMap look = currentFeatureMap("frameDNA") ; if (look->flag & FLAG_REVERSE) { graphXorLine (newx2, newy2, xfin, newy2) ; graphXorLine (xfin, newy2, xfin, oldy2+1) ; graphXorLine (xfin, oldy2+1, oldx2, oldy2+1) ; graphXorLine (oldx2, oldy2+1, oldx2, oldy2) ; graphXorLine (oldx2, oldy2, xdeb, oldy2) ; graphXorLine (xdeb, oldy2, xdeb, newy2-1) ; graphXorLine (xdeb, newy2-1, newx2, newy2-1) ; graphXorLine (newx2, newy2-1, newx2, newy2) ; } else { graphXorLine (oldx2, oldy2, xfin, oldy2) ; graphXorLine (xfin, oldy2, xfin, newy2-1) ; graphXorLine (xfin, newy2-1, newx2, newy2-1) ; graphXorLine (newx2, newy2-1, newx2, newy2) ; graphXorLine (newx2, newy2, xdeb, newy2) ; graphXorLine (xdeb, newy2, xdeb, oldy2+1) ; graphXorLine (xdeb, oldy2+1, oldx2, oldy2+1) ; graphXorLine (oldx2, oldy2+1, oldx2, oldy2) ; } } static void fMapLeftDrag2 (double x, double y) { float newx1, newy1, xx1, xx2, yy1, yy2 ; int box ; frameDNA() ; oldx3 = oldx3 + x - oldx4 ; box =graphBoxAt (x, y, &newx1, &newy1) ; if (x < xlimit && (box > oldbox || (box == oldbox && newx1 >= (int)oldx))) { seqDragBox = box ; newx = newx1 ; newy = newy1 ; graphBoxDim (seqDragBox, &xx1, &yy1, &xx2, &yy2) ; newy2 = (int)newy + yy1 + 1; newx2 = (int)newx + xx1 + 1; } oldx4 = x ; frameDNA() ; } static void fMapLeftDrag3 (double x, double y) { float newx1, newy1, xx1, xx2, yy1, yy2 ; int box ; frameDNA() ; box =graphBoxAt (x, y, &newx1, &newy1) ; if (x < xlimit && (box > oldbox || (box == oldbox && newx1 >= (int)oldx))) { seqDragBox = box ; newx = newx1 ; newy = newy1 ; graphBoxDim (seqDragBox,&xx1, &yy1, &xx2, &yy2) ; /* xdeb = xx1 ; xfin = xx2 ;*/ newy2 = (int)newy + yy1 + 1; newx2 = (int)newx + xx1 + 1; } frameDNA() ; } static void fMapLeftUp (double x, double y) { float pos, pos1, min = 1000000 ;/* mhmp 30.04*/ int i ; SEG *seg ; KEY keymin = 0 ; FeatureMap look = currentFeatureMap("fMapLeftUp") ; graphXorLine (0, oldy, mapGraphWidth, oldy) ; graphRegister (LEFT_DRAG, 0) ; graphRegister (LEFT_UP, 0) ; if (x < (look->map->thumb.x - 1.5)) /* show closest clone or locus! */ { pos1 = WHOLE2MAP(look->map, y) ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) { seg = arrp(look->segs,i,SEG) ; if (class(seg->key) == _VClone || class(seg->key) == _VLocus) { pos = (seg->x1 + seg->x2)*0.5 ; if (abs (pos - pos1) < min) { min = abs (pos - pos1) ; keymin = seg->key ;} } } if (keymin) { if (class(keymin) == _VClone) display (keymin, 0, "PMAP") ; else display (keymin, 0, "GMAP") ; } } return; } /* fMapLeftUp */ static void fMapLeftUp2 (double x, double y)/*mhmp 22.04.98 */ { frameDNA() ; graphRegister (LEFT_UP, 0) ; graphRegister (LEFT_DRAG, 0) ; fMapSelectBox (oldlook, oldbox, oldx, oldy) ; } /***************************************************************/ static void fMapPick (int box, double x, double y) { float x1,x2,y1,y2 ; int xxx, ii ; FeatureMap look = currentFeatureMap("fMapPick") ; if (look != selectedfMap) fMapSelect(look) ; if (!box) /* cursor line */ { graphBoxDim (box, &x1, &y1, &x2, &y2) ; y += y1 ; oldy = y ; graphColor (BLACK) ; graphXorLine (0, y, mapGraphWidth, y) ; graphRegister (LEFT_DRAG, (GraphFunc)fMapLeftDrag) ; graphRegister (LEFT_UP, (GraphFunc)fMapLeftUp) ; } else if (box == look->map->thumb.box) graphBoxDrag (look->map->thumb.box, mapThumbDrag) ; else if (box == look->originBox) graphTextEntry (look->originBuf,0,0,0,0) ; else if (box == look->zoneBox) graphTextEntry (look->zoneBuf,0,0,0,0) ; else if (box >= look->minLiveBox && box != look->summaryBox) { /* Make sure the box is not one of the command buttons, or the yellow */ /* summary box. */ if (box == look->activeBox && BOXSEG(box)->type != DNA_SEQ && BOXSEG(box)->type != TRANS_SEQ && BOXSEG(box)->type != TRANS_SEQ_UP && BOXSEG(box)->type != PEP_SEQ) { fMapFollow (look, x, y) ; } else if (BOXSEG(box)->type == DNA_SEQ || BOXSEG(box)->type == PEP_SEQ) { graphBoxDim (box,&x1, &y1, &x2, &y2) ; graphColor (BLACK) ; oldx = x ; oldy = y; oldlook = look ; oldbox = box; oldx2 = (int)x + x1 + 0.0001; /* mhmp 006.10.98 pb d'arrondi */ oldy2 = (int)y + y1; newx2 = oldx2 +1; newy2 = oldy2; xfin = x2 ; xlimit = x2 ; xdeb = x1 ; frameDNA() ; oldx3 = oldx ; oldx4 = oldx2 ; fMapLeftDrag3 (oldx2, oldy2) ; graphRegister (LEFT_DRAG, (GraphFunc)fMapLeftDrag3) ; graphRegister (LEFT_UP, (GraphFunc)fMapLeftUp2) ; } else if ( BOXSEG(box)->type == TRANS_SEQ || BOXSEG(box)->type == TRANS_SEQ_UP) { if (look->flag & FLAG_REVERSE) { xxx = COORD (look, BOXSEG(box)->x1) - 3*(int)x ; if (xxx < array(look->maxDna,arrayMax(look->maxDna) - 1 , int)) return ; for (ii=1; iimaxDna); ii++) if (xxx > array(look->maxDna, ii, int)) { if (xxx > array(look->minDna, ii-1, int)) return ; break ; } } else { xxx = COORD (look, BOXSEG(box)->x1) + 3*(int)x ; if (xxx > array(look->maxDna,arrayMax(look->maxDna) - 1 , int)) return ; for (ii=1; iimaxDna); ii++) if (xxx < array(look->maxDna, ii, int)) { if (xxx < array(look->minDna, ii-1, int)) return ; break ; } } graphBoxDim (box,&x1, &y1, &x2, &y2) ; graphColor (BLACK) ; oldx = x ; oldy = y; oldlook = look ; oldbox = box; oldx2 = (int)x + x1 + 0.0001; /* mhmp 006.10.98 pb d'arrondi */ oldy2 = (int)y + y1; newx2 = oldx2 +1; newy2 = oldy2; xfin = x2 ; xlimit = x2 ; xdeb = x1 ; frameDNA() ; oldx3 = oldx ; oldx4 = oldx2 ; fMapLeftDrag2 (oldx2, oldy2) ; graphRegister (LEFT_DRAG, (GraphFunc)fMapLeftDrag2) ; graphRegister (LEFT_UP, (GraphFunc)fMapLeftUp2) ; } else /* not active and not sequence */ fMapSelectBox (look, box, x, y) ; } return; } /* fMapPick */ /***********************************/ static void fMapKbd (int k) { int box ; FeatureMap look = currentFeatureMap("fMapKbd") ; if(look != selectedfMap) fMapSelect(look) ; if (!look->activeBox) return ; box = look->activeBox ; switch (k) { case UP_KEY : if (box > look->minLiveBox) --box ; break ; case DOWN_KEY : if (box < arrayMax(look->boxIndex) - 1) ++box ; break ; } fMapSelectBox (look, box, 0, 0) ; if (look->activeBox != box) fMapFollow (look,0.,0.) ; return; } /* fMapKbd */ /************** little function to act as a callback ***********/ static void fMapRefresh (void) { FeatureMap look = currentFeatureMap("fMapRefresh") ; fMapDraw (look, 0) ; } /*********************************************************/ /*********************************************************/ #ifndef USE__SMAP /* Check the length of the sequence to be viewed, main use of this is to */ /* give the user a chance to stop if they have clicked on a huge sequence. */ /* The technique is to look at any of the tags that give a length that are */ /* present in the object they clicked on. I make the assumption that if none */ /* of these are present then the object simply cannot be shown because it */ /* has no position in the sequence. This may turn out to be wrong, perhaps */ /* I should be looking in the objects parents etc. */ /* */ /* Threshold for asking user if they want to continue is currently set at */ /* a MBase. */ /* */ static BOOL checkSequenceLength(KEY seq) { BOOL result = FALSE ; OBJ obj = NULL ; Array positions = NULL ; BOOL found_position = FALSE ; int key_index, start_index, end_index, array_incr ; BSunit *u ; int i ; int min, max ; enum {MAX_SEQ_SPAN = 1048576} ; if (!seq || !(obj = bsCreate(seq))) return result ; positions = arrayCreate(256, BSunit) ; /* First look for Source_exons as these will set the initial span (code */ /* may make this bigger later but this doesn't matter), if can't find that */ /* then look for S_Child tags (i.e. the Smap system), if can't find that */ /* look for Subsequence tags (i.e. the Sequence system). */ if (bsFindTag(obj, _Source_Exons) && bsFlatten(obj, 2, positions)) { found_position = TRUE ; key_index = -1 ; start_index = 0 ; end_index = 1 ; array_incr = 2 ; } else if (bsFindTag (obj, str2tag ("S_Child")) && bsGetArray (obj, str2tag ("S_Child"), positions, 4)) { found_position = TRUE ; key_index = 1 ; start_index = 2 ; end_index = 3 ; array_incr = 4 ; } #ifdef ACEDB4 else if (bsFindTag(obj, _Subsequence) && bsFlatten(obj, 3, positions)) { found_position = TRUE ; key_index = 0 ; start_index = 1 ; end_index = 2 ; array_incr = 3 ; } #endif /* ACEDB4 */ /* If we find some positional information from which to get the sequence */ /* length then get the max/min from it. */ if (found_position) { min = INT_MAX ; /* Assume all sequence coords in */ max = 0 ; /* int range. */ for (i = 0 ; i < arrayMax(positions) ; i += array_incr) { u = arrp(positions, i, BSunit) ; if (key_index > 0 && !u[key_index].k) { messerror ("Tag with missing subsequence in %s", name(seq)) ; continue ; } if (!u[start_index].i || !u[end_index].i) { messerror ("Coords of subsequence %s missing in %s", name(u[key_index].k), name(seq)) ; continue ; } if (u[start_index].i < u[end_index].i) { if (u[start_index].i < min) min = u[start_index].i ; if (u[end_index].i > max) max = u[end_index].i ; } else { if (u[start_index].i > max) max = u[start_index].i ; if (u[end_index].i < min) min = u[end_index].i ; } } if ((max - min) > MAX_SEQ_SPAN) result = messQuery("You have asked to view a sequence that is %d bases long," " do you really want to continue ?", (max - min)) ; else result = TRUE ; } else if (bsFindTag(obj, _DNA)) result = TRUE; arrayDestroy(positions) ; bsDestroy(obj); return result ; } /* This routine attempts to find the position of the given sequence within */ /* its parent, its parents parent and so on. */ /* */ /* I use an associator to remember which sequences I've been to before in an */ /* attempt to prevent this routine from looping endlessly if somewhere a */ /* parent points back to a sequence below it. */ /* */ static BOOL findStartStop (KEY *seq, int *start, int *stop) { OBJ Seq ; int x1, x2 ; KEY parent ; static Array s_children = NULL ; Associator sequence_assoc = NULL ; /* Used to remember visited sequences. */ void *unused = NULL ; /* We don't care about the actual value stored in the assoc. */ if (!*seq || !(Seq = bsCreate(*seq))) return FALSE ; sequence_assoc = assCreate() ; while (Seq) /* recurse up */ { if (bsFindTag (Seq, str2tag("S_Parent")) && /* tag2 system */ bsGetKeyTags (Seq, _bsRight, 0) && /* a tag */ bsGetKey (Seq, _bsRight, &parent)) { int i ; bsDestroy (Seq) ; /* find the child in the parent */ s_children = arrayReCreate (s_children, 128, BSunit) ; if (!(Seq = bsCreate(parent)) || !bsGetArray (Seq, str2tag("S_Child"), s_children, 4)) break ; /* stop the recursion */ for (i = 0 ; i < arrayMax(s_children) ; i += 4) if (arrp(s_children, i+1, BSunit)->k == *seq) break ; if (i >= arrayMax(s_children)) break ; /* stop the recursion */ x1 = arrp(s_children, i+2, BSunit)->i ; x2 = arrp(s_children, i+3, BSunit)->i ; if (x1 < x2) { if (!*stop) { *start = x1 ; *stop = x2 ; } else { *start += x1-1 ; *stop += x1-1 ; } } else { if (!*stop) { *start = x2 ; *stop = x1 ; } else { int tmp = *start ; *start = x1 - *stop + 1 ; *stop = x1 - tmp + 1 ; } } *seq = parent ; } #ifdef ACEDB4 else if (bsGetKey (Seq, _Source, &parent)) { bsDestroy (Seq) ; /* Try to insert key for new sequence, if this fails then we've */ /* seen this sequence before. */ if (!assInsert(sequence_assoc, assVoid(parent), unused)) { messerror("loop in parent/child sequence objects, current object is %s", nameWithClass(parent)) ; break ; } /* find the child in the parent */ if (!(Seq = bsCreate(parent)) || !bsFindKey (Seq, _Subsequence, *seq) || !bsGetData (Seq, _bsRight, _Int, &x1) || !bsGetData (Seq, _bsRight, _Int, &x2)) break ; if (x1 < x2) { if (!*stop) { *start = x1 ; *stop = x2 ; } else { *start += x1-1 ; *stop += x1-1 ; } } else { if (!*stop) { *start = x2 ; *stop = x1 ; } else { int tmp = *start ; *start = x1 - *stop + 1 ; *stop = x1 - tmp + 1 ; } } *seq = parent ; } #endif /* ACEDB4 */ else break ; } if (Seq) bsDestroy (Seq) ; assDestroy(sequence_assoc) ; return (*start != 1 || *stop != 0) ; } #endif /* !USE_SMAP */ /*********************************************************/ #ifdef USE_SMAP static consMapDNAErrorReturn callback(KEY key, int position) { BOOL OK; if (!reportErrors) return sMapErrorReturnContinue; OK = messQuery("Mismatch error at %d in %s.\n" "Do you wish to see further errors?", position, name(key)); return OK ? sMapErrorReturnContinue : sMapErrorReturnSilent; } #endif static FeatureMap fMapCreateLook(KEY seq, int start, int stop, BOOL reuse, FmapDisplayData *display_data) { FeatureMap look ; #ifndef USE_SMAP FeatureMap curr_look = NULL ; #endif int startOrig = start ; KEY seqOrig = seq ; /* do again - no problem if done already in fMapDisplay */ #ifdef USE_SMAP if (!setStartStop(seq, &start, &stop, &seq, NULL)) { messerror("Could not find start/stop of " "key: %s, start: %d stop: %d", name(seq), start, stop) ; return NULL ; } #else findStartStop (&seq, &start, &stop) ; #endif look = (FeatureMap) messalloc (sizeof(struct FeatureMapStruct)) ; look->magic = &FeatureMap_MAGIC ; look->handle = handleCreate () ; look->segsHandle = handleHandleCreate (look->handle) ; look->seqKey = seq ; look->seqOrig = seqOrig ; #ifdef USE_SMAP look->smap = sMapCreate(look->handle, look->seqKey, start, stop, NULL); #endif /* Need to look at some old settings. */ #ifndef USE_SMAP if (reuse) curr_look = currentFeatureMap("fMapCreateLook") ; if (reuse) /* Must come before fMapFindDNA */ look->reportDNAMismatches = curr_look->reportDNAMismatches ; else look->reportDNAMismatches = TRUE ; #else /* This can go once all old dna stuff is removed. */ look->reportDNAMismatches = FALSE ; #endif /* test code.... */ if (display_data && display_data->features) { look->feature_set = sMapFeaturesSet(look->handle, display_data->features, display_data->include_features) ; } #ifndef USE_SMAP look->zoneMin = start ; look->zoneMax = stop + 1 ; #endif #ifdef USE_SMAP reportErrors = look->reportDNAMismatches; if (!(look->dna = sMapDNA(look->smap, 0, callback))) #else if (!(look->dna = fMapFindDNA(seq, &start, &stop, look->reportDNAMismatches))) #endif { messfree (look->handle) ; messfree (look) ; return NULL ; } #ifdef USE_SMAP look->zoneMin = 0 ; look->zoneMax = stop - start + 1; look->origin = 0 ; stop = start + arrayMax(look->dna) + 1; #else /* This coordinate fiddling was done because fMapFindDNA() (see above) alters start/stop once more. */ look->zoneMin -= start ; look->zoneMax -= start ; look->origin = look->zoneMin - startOrig + 1 ; #endif look->colors = arrayHandleCreate (arrayMax(look->dna), char, look->handle) ; arrayMax(look->colors) = arrayMax(look->dna) ; fMapClearDNA (look) ; /* clear colors array to _WHITE */ look->chosen = assHandleCreate (look->handle) ; look->antiChosen = assHandleCreate (look->handle) ; look->length = arrayMax(look->dna) ; #ifdef USE_SMAP look->fullLength = sMapLength (seq) ; #else look->fullLength = sequenceLength (seq) ; #endif look->start = start-1; look->stop = stop-1; look->featDict = dictHandleCreate (1000, look->handle) ; dictAdd (look->featDict, "null", 0) ; /* so real entries are non-zero */ look->homolInfo = arrayHandleCreate (256, HOMOLINFO, look->handle) ; look->seqInfo = arrayHandleCreate (32, SEQINFO, look->handle) ; look->homolBury = bitSetCreate (256, look->handle) ; look->homolFromTable = bitSetCreate (256, look->handle) ; look->visibleIndices = arrayHandleCreate (64, int, look->handle) ; look->preserveColumns = TRUE ; look->redrawColumns = TRUE ; if (display_data) look->destroy_obj = display_data->destroy_obj ; else look->destroy_obj = FALSE ; look->mcache = methodCacheCreate (look->handle); /* must create map before convert, because it adds columns */ look->map = mapCreate (fMapRefresh) ; /* map redraw func */ look->map->drag = fMapDrag; /* no graph perhaps, so do not attach to graph! */ fMapDefaultColumns (look->map) ; if (!fMapConvert(look)) { messfree (look->handle) ; messfree (look) ; return NULL ; } look->boxIndex = arrayHandleCreate (64, int, look->handle) ; look->activeBox = 0 ; seqDragBox = 0 ; /* Some items of the main menu are specific to an individual look */ /* so we must copy the main menu so we can alter it. */ look->mainMenu = handleAlloc(NULL, look->handle, FMAPMENU_SIZE) ; memcpy(look->mainMenu, fMapMenu, FMAPMENU_SIZE) ; /* Some menu items need to be set according to current program */ /* state, e.g. cut data state can be set from command line. */ if (fMapCutCoords_G) /* call will reverse what we set here. */ look->cut_data = FMAPCUT_COORDS ; else look->cut_data = FMAPCUT_DNA ; toggleTextInMenu(look->mainMenu, toggleCutSelection, (look->cut_data == FMAPCUT_DNA) ? CUT_COORDS : CUT_DNA) ; look->flag = 0 ; look->flag |= FLAG_AUTO_WIDTH ; /* I honestly cannot understand why you would take that away */ return look ; } /* fMapCreateLook */ /************************************************************/ static void fMapAttachToGraph (FeatureMap look) { graphAssociate (&GRAPH2FeatureMap_ASSOC, look) ; /* attach look for fmap package */ graphAssociate (&MAP2LOOK_ASSOC, look) ; /* attach look for map package */ } /************************************************************/ /* this is being used a DisplayFunc called by display() */ BOOL fMapDisplay (KEY key, KEY from, BOOL isOldGraph, void *app_data) { FmapDisplayData *create_data = (FmapDisplayData *)app_data ; int start, stop ; FeatureMap look, curr_look = NULL ; KEY seq = 0 ; OBJ obj ; BOOL shiftOrigin = FALSE ; BOOL reuse ; fMapInitialise() ; if (isOldGraph) curr_look = currentFeatureMap("fMapDisplay") ; if (class(key) == _VSequence || bIndexTag (key, str2tag("SMAP"))) seq = key ; else if (class(key) == _VDNA) lexReClass (key, &seq, _VSequence) ; else if ((obj = bsCreate (key))) { if ((!bsGetKey (obj, _Sequence, &seq)) && bsGetKey (obj, str2tag("Genomic_sequence"), &seq)) shiftOrigin = TRUE ; bsDestroy (obj) ; } if (key != seq) { from = key ; key = seq ; } if (!seq || iskey(seq) != 2) /* iskey tests a full object */ return FALSE ; /* Crude check of size of object to be displayed, gives user a chance to */ /* stop if they accidentally try to display a huge sequence. */ if (!isGifDisplay && !checkSequenceLength(seq)) { return FALSE ; } /* recurse up to find topmost link - must do here rather than wait for fMapCreateLook() so we can check if we have it already */ #ifdef USE_SMAP sMapTreeRoot (seq, 1, 0, &seq, &start, &stop) ; /* smap keeps coords in the original order but fmap needs them to be in */ /* ascending order. */ if (start > stop) { int tmp ; tmp = start ; start = stop ; stop = tmp ; } #else start = 1 ; stop = 0 ; /* go to end of sequence */ findStartStop (&seq, &start, &stop) ; #endif /* Is the requested part of the fMap already on display ? */ if (isOldGraph && graphAssFind (&GRAPH2FeatureMap_ASSOC, &look) && look->seqKey == seq && arrayExists (look->dna) && start-1 >= look->start && stop-1 <= look->stop) { look->seqOrig = key ; fMapCentre (look, key, from) ; if (shiftOrigin) fMapSetOriginFromKey (look, from) ; graphRetitle (name(from ? from : key)) ; fMapDraw (look, key) ; return TRUE ; } /* Reuse the current fmap graph ? */ if (isOldGraph && graphAssFind (&GRAPH2FeatureMap_ASSOC, 0)) reuse = TRUE ; else reuse = FALSE ; /* make a new look now */ look = fMapCreateLook (seq, start, stop, reuse, create_data) ; if (!look) { /* The drawing would be meaningless or perhaps take too long.*/ return FALSE ; } /* If the fmap is being reused the default is now to keep the */ /* currently selected active columns the same and to display */ /* the column buttons if they were being displayed. */ /* Otherwise column settings will revert to the default and the */ /* column buttons will not be displayed. */ if (reuse) { if (curr_look->preserveColumns == TRUE) { look->preserveColumns = TRUE ; mapColCopySettings(curr_look->map, look->map) ; if (curr_look->flag & FLAG_COLUMN_BUTTONS) look->flag |= FLAG_COLUMN_BUTTONS ; } else { look->preserveColumns = FALSE ; look->flag &= !FLAG_COLUMN_BUTTONS ; } /* Set menu texts for settings preserved between fmap recalcs. */ toggleTextInMenu(look->mainMenu, togglePreserveColumns, look->preserveColumns ? UN_PRESERVE_COLS : PRESERVE_COLS) ; /* This needs to be cleared up and integrated with new smapdna stuff... */ toggleTextInMenu(look->mainMenu, toggleReportMismatch, (look->reportDNAMismatches ? NO_REPORT_MISMATCHES : REPORT_MISMATCHES)) ; } if (reuse) { KEY kk = from ? from : key ; if (class(kk) == _VSequence || class(kk) == _VTranscribed_gene) graphRetitle(name (kk)) ; else if (from || look->seqKey) graphRetitle(name(from ? from : look->seqKey)) ; fMapDestroyCallback () ; /* kill look currently attached */ graphRegister (DESTROY, fMapDestroyCallback) ; /* reregister */ graphClear () ; } else /* open a new window */ { displayCreate ("FMAP"); graphRetitle (name(from ? from : key)) ; graphRegister (RESIZE, fMapRefresh) ; graphRegister (DESTROY, fMapDestroyCallback) ; graphRegister (KEYBOARD, fMapKbd) ; graphRegister (MESSAGE_DESTROY, displayUnBlock) ; } look->seqOrig = key ; look->map->min = 0 ; look->map->max = look->length ; look->graph = graphActive() ; fMapAttachToGraph (look) ; mapAttachToGraph (look->map) ; fMapcDNAInit (look) ; fMapCentre (look, key, from) ; /* centre on key if possible and choose magnification */ if (shiftOrigin) fMapSetOriginFromKey (look, from) ; if (!getenv("ACEDB_PROJECT") && (look->zoneMax - look->zoneMin < 1000)) { mapColSetByName(look->map, "DNA Sequence", TRUE) ; mapColSetByName(look->map, "Coords", TRUE) ; } #ifdef FMAP_DEBUG fMapPrintLook(look, "fMapDisplay"); #endif fMapDraw (look, key) ; fMapSelect (look) ; return TRUE ; } /* fMapDisplay */ /************************************************/ /************************************************/ /* also called from fmapgene.c */ int fMapOrder (void *a, void *b) { SEG *seg1 = (SEG*)a, *seg2 = (SEG*)b ; int diff ; diff = seg1->type - seg2->type ; if (diff > 0) return 1 ; else if (diff < 0) return -1 ; if ((seg1->type | 0x1) == SPLICED_cDNA_UP) { if (!(seg1->source * seg2->source)) return seg1->source ? -1 : 1 ; diff = seg1->source - seg2->source ; if (diff) return diff ; diff = seg1->parent - seg2->parent ; if (diff) return diff ; { int x1, x2, y1, y2 ; x1 = seg1->data.i >> 14 & 0x3FFF ; x2 = seg1->data.i & 0x3FFF ; y1 = seg2->data.i >> 14 & 0x3FFF ; y2 = seg2->data.i & 0x3FFF ; if ((x1 - x2) * (y1 - y2) < 0) return x1 < x2 ? -1 : 1 ; } } diff = seg1->x1 - seg2->x1 ; if (diff > 0) return 1 ; else if (diff < 0) return -1 ; diff = seg1->x2 - seg2->x2 ; if (diff > 0) return 1 ; else if (diff < 0) return -1 ; return 0 ; } /*****************************************************************/ /* need key and from for Calcul positioning */ static void fMapCentre (FeatureMap look, KEY key, KEY from) { int i, orig = -999999 ; float nx, ny ; SEG *seg, *seg1 ; BOOL fromSeq = (class(from) == _VSequence || (_VTranscribed_gene && class(from) == _VTranscribed_gene) || (_VIST && class(from) == _VIST) || (_VProtein && class(from) == _VProtein)) ; /* so that protein homol get centered */ graphFitBounds (&mapGraphWidth,&mapGraphHeight) ; halfGraphHeight = 0.5 * (mapGraphHeight - topMargin) ; nx = mapGraphWidth ; ny = mapGraphHeight - topMargin - 5 ; /* WHY IS IT "- 5" ???????? */ seg = arrp(look->segs, 0 , SEG) ; if (key) for (i = 1 ; i < arrayMax (look->segs) ; ++i) { seg1 = arrp (look->segs, i, SEG) ; if (fromSeq && seg1->key == key) orig = seg->x1 ; if ((!fromSeq && seg1->key == key) || (from && seg1->key == from)) { seg = seg1 ; break ; } } i = KEYKEY(from)*10 ; if (class(from) == _VCalcul && i < seg->x2 - seg->x1 + 1) { if (seg->type == SEQUENCE_UP) look->map->centre = seg->x2 - i ; else look->map->centre = seg->x1 + i ; /* Go to rather high magnification look->map->mag = (ny-5)/ 5000.0 ; */ } else look->map->centre = 0.5 * (seg->x1 + seg->x2) ; /* shouldn't the divisor be * 1.05 here like other mags ??? e.g. fMapWhole ????? */ look->map->mag = ny/(seg->x2 - seg->x1 + 1.0) ; #ifdef ED_G_NEVER_INCLUDE_THIS_CODE look->map->mag = ny / ((seg->x2 - seg->x1) * 1.05) ; #endif /* ED_G_NEVER_INCLUDE_THIS_CODE */ look->origin = orig != -999999 ? orig : seg->x1 ; setZone (look, seg->x1, seg->x2 + 1) ; if (fromSeq) look->map->mag /= 1.4 ; if (look->map->mag <= 0) /* safeguard */ look->map->mag = 0.05 ; return ; } /*****************************************************************/ static void processMethod (FeatureMap look, KEY key, float offset, MapColDrawFunc func) { METHOD *meth = NULL; float p ; if (!(meth = methodCacheGet(look->mcache, key, ""))) /* Note the string (now "") will eventually come from the * aceDiff-Text in the view object for this FMAP */ return; /* Some methods need to not be displayed so we don't even put them in a */ /* column....don't know if this is correct really. */ if (meth->no_display) return ; if (meth->priority) p = meth->priority ; else if (meth->flags & METHOD_FRAME_SENSITIVE) p = 4 + offset ; else p = 5 + offset ; mapInsertCol (look->map, p, TRUE, name(key), func) ; if (meth->flags & METHOD_SHOW_UP_STRAND) mapInsertCol (look->map, -p, TRUE, messprintf("-%s", name(key)), func) ; return; } /* processMethod */ void fMapProcessMethods (FeatureMap look) { int i; KEY key; if (!look->map) return ; for (i = 0 ; i < arrayMax(look->homolInfo) ; ++i) { key = arrp(look->homolInfo, i, HOMOLINFO)->method ; if (key) processMethod (look, key, 0.5, fMapShowHomol) ; } for (i = 0 ; i < arrayMax(look->seqInfo) ; ++i) { key = arrp(look->seqInfo, i, SEQINFO)->method ; if (key) processMethod (look, key, -3.0, fMapShowSequence) ; } for (i = 0 ; i < arrayMax(look->segs) ; ++i) { key = arrp(look->segs, i, SEG)->key ; if (class(key) == _VMethod) /* set up methods */ switch (arrp(look->segs, i, SEG)->type) { case FEATURE: case FEATURE_UP: processMethod (look, key, 0.7, fMapShowFeature) ; break ; case SPLICE5: case SPLICE5_UP: case SPLICE3: case SPLICE3_UP: processMethod (look, key, 0.9, fMapShowSplices) ; break ; case TRANSCRIPT: case TRANSCRIPT_UP: processMethod (look, key, 2.1, fMapcDNAShowTranscript) ; break ; case SPLICED_cDNA: case SPLICED_cDNA_UP: processMethod (look, key, 2.2, fMapcDNAShowSplicedcDNA) ; break ; case ATG: case ATG_UP: break ; /* do nothing yet */ default: messerror ("Method key %s for unknown seg type %s", name(key), fMapSegTypeName[arrp(look->segs, i, SEG)->type]) ; } } return; } /* fMapProcessMethods */ /*****************************/ static BOOL pcDown ; static int pcStart, pcEnd ; static Array pcArray = 0 ; static Array pcExons = NULL ; static void posConvertInit (SEG *seg, Array exons, BOOL isDown) { int i ; pcDown = isDown ; pcStart = seg->x1 - 1 ; pcEnd = seg->x2 + 1 ; /* I think this is a useless array... */ pcArray = arrayReCreate(pcArray,20,int) ; for (i = 0 ; i < arrayMax(exons)-2 ; i+=2) { array(pcArray,i,int) = arr(exons,i+1,BSunit).i ; array(pcArray,i+1,int) = arr(exons,i+2,BSunit).i - arr(exons,i+1,BSunit).i - 1 ; } /* Here is my array of just the exons... */ pcExons = arrayReCreate(pcExons, 20, int) ; for (i = 0 ; i < arrayMax(exons) ; i+=2) { array(pcExons, i, int) = arr(exons, i, BSunit).i ; array(pcExons, i+1, int) = arr(exons, i+1, BSunit).i ; } return ; } /* Given a position in source exon coords, returns a position in the spliced */ /* exon coords. */ /* If given the value zero, will return the end position. */ /* If given a position outside of the exons, returns -1 */ static int posConvertGetSplicePos(int pos) { int position = -1 ; int spliced_length ; int exon_length ; BOOL pos_found ; int i ; int arr1, arr2 ; for (i = 0, spliced_length = 0, pos_found = FALSE ; i < arrayMax(pcExons) && pos_found == FALSE ; i += 2) { arr1 = arr(pcExons, i, int) ; arr2 = arr(pcExons, i+1, int) ; exon_length = arr(pcExons, i+1, int) - arr(pcExons, i, int) + 1 ; spliced_length += exon_length ; if (pos == 0) position = spliced_length ; else if (pos <= spliced_length) { position = arr(pcExons, i, int) + (exon_length - (spliced_length - pos) - 1) ; pos_found = TRUE ; } } if (pos != 0 && pos_found == FALSE) position = -1 ; return position ; } static void posConvert (SEG *seg, int pos1, int pos2) { int i, arr1, arr2 ; /* Are there any exons ? */ if (arrayMax(pcArray)) { /* It would seem that these two loops are only called for CDS mapping, */ /* this is unfortunate because they are incorrect for this purpose... */ /* So here I add my own code to do the job.... */ /* If it turns out that I am right and this is all wrong then I should fix */ /* up the array in posConvertInit to hold more useful numbers... */ /* */ if (seg->type == CDS || seg->type == CDS_UP) { pos1 = posConvertGetSplicePos(pos1) ; pos2 = posConvertGetSplicePos(pos2) ; } else { for (i = 0 ; i < arrayMax(pcArray) ; i += 2) { arr1 = arr(pcArray, i, int) ; arr2 = arr(pcArray, i+1, int) ; if (pos1 > arr(pcArray, i, int)) pos1 += arr(pcArray, i+1, int) ; else break ; } for (i = 0 ; i < arrayMax(pcArray) ; i += 2) { arr1 = arr(pcArray, i, int) ; arr2 = arr(pcArray, i+1, int) ; if (pos2 > arr(pcArray, i, int)) pos2 += arr(pcArray, i+1, int) ; else break ; } } } /* Most segs go straight here because they do not have some array of */ /* source exons.... */ if (pcDown) { seg->x1 = pos1 + pcStart ; seg->x2 = pos2 + pcStart ; } else { seg->x1 = pcEnd - pos2 ; seg->x2 = pcEnd - pos1 ; } return ; } /*****************************/ /* selective for CALCULATED segs */ static void addOldSegs (Array segs, Array oldSegs, MethodCache mcache) { SEG *seg, *oldMaster, *newMaster ; int i, j, length, offset ; METHOD *meth; oldMaster = arrp(oldSegs,0,SEG) ; newMaster = arrp(segs,0,SEG) ; if (newMaster->key != oldMaster->key) return ; length = newMaster->x2 - newMaster->x1 + 1 ; offset = newMaster->x1 - oldMaster->x1 ; j = arrayMax(segs) ; for (i = 1 ; i < arrayMax(oldSegs) ; ++i) { seg = arrp(oldSegs, i, SEG) ; switch (seg->type) { case FEATURE: case FEATURE_UP: case ATG: case ATG_UP: case SPLICE5: case SPLICE5_UP: case SPLICE3: case SPLICE3_UP: if (class(seg->key) != _VMethod) messcrash ("Non-method key in addOldSegs") ; meth = methodCacheGet(mcache, seg->key, ""); if (meth && meth->flags & METHOD_CALCULATED) { seg->x1 -= offset ; seg->x2 -= offset ; if (seg->x1 >= 0 && seg->x2 < length) array(segs, j++, SEG) = *seg ; } break ; default: /* needed for picky compiler */ break ; } } return; } /* addOldSegs */ KEY defaultSubsequenceKey (char* name, int colour, float priority, BOOL isStrand) { KEY key ; OBJ obj ; lexaddkey (name, &key, _VMethod) ; /* mieg, dec 27 */ if (iskey(key)!=2) /* i.e. no object */ if ((obj = bsUpdate (key))) /* make it */ { float width = 2 ; KEY colourKey = _WHITE + colour ; if (bsAddTag (obj, str2tag("Colour"))) { bsPushObj (obj) ; /* get into sub-model #Colour */ bsAddTag (obj, colourKey) ; bsGoto (obj, 0) ; } bsAddData (obj, str2tag("Width"), _Float, &width) ; bsAddData (obj, str2tag("Right_priority"), _Float, &priority) ; if (isStrand) bsAddTag (obj, str2tag("Show_up_strand")) ; bsSave (obj) ; } return key ; } BOOL fMapFindCoding (FeatureMap look) { BOOL result = TRUE ; int i, j, j1 ; Array index = NULL ; SEG *seg ; KEY parent ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) { /* User can interrupt by pressing F4 */ if (messIsInterruptCalled()) { if (messQuery("Finding coding - do you really want to interrupt fmap ?")) { if (index != NULL) arrayDestroy(index) ; result = FALSE ; break ; } else messResetInterrupt() ; } seg = arrp(look->segs,i,SEG) ; if (seg->type == CDS && (parent = seg->parent) && fMapGetCDS (look, parent, 0, &index)) { j1 = 0 ; for (j = 0 ; j < arrayMax(index) ; ++j) if (j == arrayMax(index) - 1 || arr(index,j+1,int) > arr(index,j,int)+1) /* boundary */ { seg = arrayp (look->segs, arrayMax(look->segs), SEG) ; seg->key = _CDS ; seg->type = CODING ; seg->x1 = arr(index, j1, int) ; seg->x2 = arr(index, j, int) ; seg->parent = parent ; seg->data.i = j1 ; j1 = j+1 ; } arrayDestroy (index) ; index = NULL ; } else if (seg->type == CDS_UP && (parent = seg->parent) && fMapGetCDS (look, parent, 0, &index)) { j1 = 0 ; for (j = 0 ; j < arrayMax(index) ; ++j) if (j == arrayMax(index) - 1 || arr(index,j+1,int) < arr(index,j,int)-1) /* boundary */ { seg = arrayp (look->segs, arrayMax(look->segs), SEG) ; seg->key = _CDS ; seg->type = CODING_UP ; seg->x1 = arr(index,j,int) ; seg->x2 = arr(index,j1,int) ; seg->parent = parent ; seg->data.i = j1 ; j1 = j+1 ; } arrayDestroy (index) ; index = NULL ; } /* remove duplicate INTRON segs */ if (seg->type == INTRON || seg->type == INTRON_UP) { SEG *seg2 = seg ; int flags = 0 ; BOOL existsParent = FALSE ; SEQINFO *sinf ; while (i < arrayMax(look->segs)-1 && seg->type == seg2->type && seg->x1 == seg2->x1 && seg->x2 == seg2->x2) { sinf = arrp(look->seqInfo, seg->data.i, SEQINFO) ; flags |= sinf->flags ; if (seg->parent) existsParent = TRUE ; ++i ; ++seg ; } for (--i, --seg ; seg2 <= seg ; ++seg2) if (seg2->parent || !existsParent) { sinf = arrp(look->seqInfo, seg2->data.i, SEQINFO) ; if (sinf->flags != flags) { sinf = arrayp(look->seqInfo, arrayMax(look->seqInfo), SEQINFO) ; *sinf = arr(look->seqInfo, seg2->data.i, SEQINFO) ; sinf->flags = flags ; seg2->data.i = arrayMax(look->seqInfo)-1 ; } existsParent = TRUE ; } else { *seg2 = array(look->segs,arrayMax(look->segs)-1,SEG) ; --arrayMax(look->segs) ; } } } return result ; } /************************************************/ /* NOTE, this small routine makes use of the new smap convert code to */ /* make the array of segs. Eventually this routine should merge/replace */ /* fMapConvert(). If you look at fMapConvert() you will see that that it */ /* calls this routine and then returns immediately making the rest of */ /* fMapConvert redundant. */ /* */ /* */ BOOL fMapSMapConvert(FeatureMap look, Array oldSegs) { BOOL result = FALSE ; /* This routine now encapsulates the convert. */ result = sMapFeaturesMake(look->smap, look->seqKey, look->start, look->stop, look->length, oldSegs, (look->flag & FLAG_COMPLEMENT), look->mcache, look->dna, look->feature_set, look->featDict, look->segs, look->seqInfo, look->homolInfo, look->homolBury, look->homolFromTable, look->segsHandle) ; if (result) { /* The comment about these 2 routines not working is deeply distressing */ /* I'm not sure if its true or what.... */ fMapProcessMethods (look) ; /* always do this - to pick up edits */ arraySort(look->segs, fMapOrder) ; /* these 2 routines no longer work after reverse complement because RC calls sort * they should be written differently */ fMapOspFindOligoPairs (look) ; /* why AFTER sorting ? not sure, but probably ok */ look->lastTrueSeg = arrayMax(look->segs) ; if (fMapDebug) dumpstats(look) ; } return result; } /* This routine needs to be merged with the above routine to give a new */ /* fmapconvert routine. */ /* */ BOOL fMapConvert(FeatureMap look) { KEY key, parent ; Array segs = NULL ; Array oldSegs = NULL ; STORE_HANDLE oldSegsHandle ; BSunit *u ; SEG *seg, *seg1, *seg2 ; OBJ obj ; int i, nsegs, iseg, pos1, pos2, tmp ; BOOL isDown, isExons ; SEQINFO *sinf ; KEY M_Transposon, M_SPLICED_cDNA, M_RNA ; KEY M_TRANSCRIPT, M_Pseudogene, M_Coding ; static Array units = 0 ; BITSET_MAKE_BITFIELD /* Add bitField array */ static Array s_children = 0 ; int cds_min, cds_max ; BOOL result = FALSE ; /* do a little subroutine that will replace this and be called from */ /* posconvert...much cleaner....sigh... */ /* Need POS_TO_SEG1 since posConvert only valid after posConvertInit */ #define POS_TO_SEG1 if (isDown) \ { seg1->x1 = seg->x1 + pos1 - 1 ; \ seg1->x2 = seg->x1 + pos2 - 1 ; \ } \ else \ { seg1->x1 = seg->x2 - pos2 + 1 ; \ seg1->x2 = seg->x2 - pos1 + 1 ; \ } chrono("fMapConvert") ; fMapInitialise () ; M_Coding = defaultSubsequenceKey ("Coding", BLUE, 2.0, TRUE) ; M_RNA = defaultSubsequenceKey ("RNA", GREEN, 2.0, TRUE) ; M_Pseudogene = defaultSubsequenceKey ("Pseudogene", DARKGRAY, 2.0, TRUE) ; M_Transposon = defaultSubsequenceKey ("Transposon", DARKGRAY, 2.05, FALSE) ; M_TRANSCRIPT = defaultSubsequenceKey ("TRANSCRIPT", DARKGRAY, 2.1, FALSE) ; M_SPLICED_cDNA = defaultSubsequenceKey ("SPLICED_cDNA", DARKGRAY, 2.2, FALSE) ; if (!look || !look->seqKey) messcrash ("No sequence for fMapConvert") ; if (look->pleaseRecompute > 0) /* -1 should be untouched */ look->pleaseRecompute = FALSE ; if (!(obj = bsCreate (look->seqKey))) { messout ("No data for sequence %s", name(look->seqKey)) ; chronoReturn() ; return FALSE ; } bsDestroy (obj) ; messStatus ("Recalculating DNA window") ; /* Initialise handles, arrays and bitsets. */ oldSegs = look->segs ; oldSegsHandle = look->segsHandle ; look->segsHandle = handleHandleCreate (look->handle) ; segs = look->segs = arrayHandleCreate(oldSegs ? arrayMax(oldSegs) : 128, SEG, look->segsHandle) ; if (oldSegs) arrayMax(oldSegs) = look->lastTrueSeg ; units = arrayReCreate (units, 256, BSunit) ; for (i=0; ihomolInfo); i++) arrayDestroy(arr(look->homolInfo, i, HOMOLINFO).gaps); look->homolInfo = arrayReCreate (look->homolInfo, 256, HOMOLINFO) ; look->seqInfo = arrayReCreate (look->seqInfo, 32, SEQINFO) ; look->homolBury = bitSetReCreate (look->homolBury, 256) ; look->homolFromTable = bitSetReCreate (look->homolFromTable, 256) ; nsegs = 0 ; /* Main seg index. */ #ifdef USE_SMAP /* If smap is being used then the rest of fMapConvert() becomes redundant. */ result = fMapSMapConvert(look, oldSegs); if (oldSegs) { handleDestroy(oldSegsHandle) ; /* kills old segs and related stuff */ oldSegs = NULL ; } return result ; #endif /* first seg is master */ seg1 = arrayp(segs,nsegs++,SEG) ; seg1->parent = 0 ; /* important */ seg1->source = 0 ; seg1->key = look->seqKey ; seg1->x1 = look->start ; seg1->x2 = look->stop ; seg1->type = MASTER ; if (fMapDebug) printf ("Converting %s: %d %d\n", name(seg1->key), seg1->x1, seg1->x2) ; /* second is self */ seg2 = arrayp(segs,nsegs++,SEG) ; seg2->key = seg2->parent = look->seqKey ; seg2->source = 0 ; seg2->x1 = -look->start ; seg2->x2 = look->fullLength - 1 - look->start ; seg2->type = SEQUENCE ; seg2->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg2->data.i, SEQINFO) ; /* recurses through all sequences since they get added in turn */ /* N.B. iseg = 1 means that we start with the top sequence object referred */ /* to by look->seqKey. */ for (iseg = 1 ; iseg < nsegs ; ++iseg) { /* User can interrupt by pressing F4 */ if (messIsInterruptCalled()) { if (messQuery("Constructing virtual sequence - do you really want to interrupt fmap ?")) { result = FALSE ; goto abort ; } else messResetInterrupt() ; } /* only interested in sequences */ seg = arrp(segs,iseg,SEG) ; if (seg->type == SEQUENCE || seg->type == TRANSCRIPT) isDown = TRUE ; else if (seg->type == SEQUENCE_UP || seg->type == TRANSCRIPT_UP) isDown = FALSE ; else continue ; if (!(obj = bsCreate (seg->key))) continue ; /* Find all the subsequences using the new S_Child model. */ /* */ /* first create subsequences (new tag2 system)*/ if (bsFindTag (obj, str2tag ("S_Child"))) s_children = arrayReCreate (s_children, 128, BSunit) ; if (bsFindTag (obj, str2tag ("S_Child")) && bsGetArray (obj, str2tag ("S_Child"), s_children, 4)) { for (i = 0 ; i < arrayMax(s_children) ; i += 4) { u = arrp(s_children,i,BSunit) ; if (!u[1].k) continue ; if (look->flag & FLAG_HIDE_SMALL_CONTIGS) { pos1 = u[2].i ; pos2 = u[3].i ; /* jump if smaller than 1 kb */ if ((pos1 - pos2) * ( pos1 - pos2) < 1000000) continue ; } seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->key = seg1->parent = u[1].k ; seg1->source = seg->key ; seg1->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg1->data.i, SEQINFO) ; if (!u[2].i || !u[3].i) { messerror ("Coords of subsequence %s missing in %s", name(seg1->key), name(seg->key)) ; --nsegs ; } else { pos1 = u[2].i ; pos2 = u[3].i ; if (isDown) seg1->type = (pos2 > pos1) ? SEQUENCE : SEQUENCE_UP ; else seg1->type = (pos1 > pos2) ? SEQUENCE : SEQUENCE_UP ; if (pos1 > pos2) { tmp = pos2 ; pos2 = pos1 ; pos1 = tmp ;} POS_TO_SEG1 ; if (seg1->x2 < 0 || seg1->x1 >= look->length) --nsegs ; /* ignore this subsequence */ else { OBJ obj1 = bsCreate(seg1->key) ; if (obj1 && bsFindTag(obj1, _Assembled_from)) { seg1->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg1->data.i, SEQINFO) ; sinf->flags |= SEQ_VIRTUAL_ERRORS ; } bsDestroy (obj1) ; } } } } #ifdef ACEDB4 /* Shouldn't the above and this code be mutually exclusive ?? */ /* Find all the subsequences using the old Subsequence tags. */ /* */ if (bsFindTag (obj, _Subsequence) && bsFlatten (obj, 3, units)) { for (i = 0 ; i < arrayMax(units) ; i += 3) { u = arrp(units,i,BSunit) ; if (!u[0].k) continue ; if (look->flag & FLAG_HIDE_SMALL_CONTIGS) { pos1 = u[1].i ; pos2 = u[2].i ; /* jump if smaller than 1 kb */ if ((pos1 - pos2) * ( pos1 - pos2) < 1000000) continue ; } if (arrp(segs,iseg,SEG)->key == u[0].k) { /* about to add another seg for itself, this would go round in circles * until we run out of memory (fw-990511, see SANgc04317) */ messerror("Error : %s is subsequence of itself", name(u[0].k)); continue; } seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->key = seg1->parent = u[0].k ; seg1->source = seg->key ; seg1->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg1->data.i, SEQINFO) ; if (!u[1].i || !u[2].i) { messerror ("Coords of subsequence %s missing in %s", name(seg1->key), name(seg->key)) ; --nsegs ; } else { pos1 = u[1].i ; pos2 = u[2].i ; if (isDown) seg1->type = (pos2 > pos1) ? SEQUENCE : SEQUENCE_UP ; else seg1->type = (pos1 > pos2) ? SEQUENCE : SEQUENCE_UP ; if (pos1 > pos2) { tmp = pos2 ; pos2 = pos1 ; pos1 = tmp ;} POS_TO_SEG1 ; if (seg1->x2 < 0 || seg1->x1 >= look->length) --nsegs ; /* ignore this subsequence */ else { OBJ obj1 = bsCreate(seg1->key) ; if (obj1 && bsFindTag(obj1, _Assembled_from)) { seg1->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg1->data.i, SEQINFO) ; sinf->flags |= SEQ_VIRTUAL_ERRORS ; } bsDestroy (obj1) ; } } } } #endif /* Used for exons/introns and other bits where there are multiple segs */ /* per sequence object. */ parent = seg->key ; /* Set up the "extras" info, this used when we have a set of segs */ /* e.g. exons) where they all share some info., e.g. a method. */ sinf = arrp(look->seqInfo, seg->data.i, SEQINFO) ; if (bsFindTag (obj, str2tag("Genomic_canonical")) || bsFindTag (obj, str2tag("Genomic"))) sinf->flags |= SEQ_CANONICAL ; /* Set up the drawing style for this feature... */ /* set a default drawing style according to this features type. */ if (bsFindTag (obj, _Pseudogene)) sinf->method = M_Pseudogene ; else if (bsFindTag (obj, _Transposon)) sinf->method = M_Transposon ; else if (bsFindTag (obj, _tRNA) || bsFindTag (obj, _rRNA) || bsFindTag (obj, _snRNA)) sinf->method = M_RNA ; else if (bsFindTag (obj, _CDS)) sinf->method = M_Coding ; else sinf->method = 0 ; /* check if the default has been overridden by an explicit style in the sequence. */ if (bsGetKey (obj, str2tag("Method"), &sinf->method)) if (bsGetData (obj, _bsRight, _Float, &sinf->score)) sinf->flags |= SEQ_SCORE ; /* source exons */ pos1 = 0 ; if (bsFindTag (obj, _Source_Exons) && bsFlatten (obj, 2, units)) { dnaExonsSort (units) ; for (i = 0 ; i < arrayMax(units) ; i += 2) { pos2 = arr(units, i, BSunit).i ; if (pos1) { seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; POS_TO_SEG1 ; seg1->key = 0 ; seg1->parent = seg1->source = parent ; seg1->data.i = seg->data.i ; /* index in seqInfo */ seg1->x1++ ; /* Because the bsTree gives the coordinates of the exon */ seg1->x2-- ; seg1->type = (isDown) ? INTRON : INTRON_UP ; if (fMapDebug) printf("INTRON: %d --> %d\n", seg1->x1, seg1->x2) ; } pos1 = pos2 ; pos2 = arr(units, i+1, BSunit).i ; if (pos2) { seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; POS_TO_SEG1 ; seg1->key = 0 ; seg1->parent = seg1->source = parent ; seg1->data.i = seg->data.i ; /* index in seqInfo */ seg1->type = (isDown) ? EXON : EXON_UP ; if (fMapDebug) printf("EXON: %d --> %d\n", seg1->x1, seg1->x2) ; } else messerror ("Missing pos2 for exon in %s", name(seg->key)) ; pos1 = pos2 ; } sinf->flags |= SEQ_EXONS ; isExons = TRUE ; } else { arrayMax (units) = 0 ; array (units, 0, BSunit).i = 1 ; array (units, 1, BSunit).i = seg->x2 - seg->x1 + 1 ; isExons = FALSE ; } /* Now we have the exons, we can set up our position finding code for */ /* features that are mapped on to the exons. */ posConvertInit (seg, units, isDown) ; /* CDS */ /* OK, pay attention, when the CDS gets drawn (fMapShowSequence), the */ /* CDS seg must come _before_ the exon segs, this is guaranteed by */ /* the numerical value of its seg->type (fmap_.h) and the arraySort at */ /* the end of this routine. */ if (bsFindTag (obj, _CDS)) { seg1 = arrayp(segs, nsegs++, SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->key = _CDS ; seg1->type = isDown ? CDS : CDS_UP ; seg1->parent = seg1->source = parent ; /* CDS defaults to the whole of the spliced dna. */ cds_min = pos1 = 1 ; cds_max = pos2 = posConvertGetSplicePos(0) ; /* zero means get "end" position */ /* If the CDS start or stop positions are set then we will need to */ /* convert positions for finding the spliced DNA, e.g. for later */ /* protein translation. */ if (bsGetData(obj, _bsRight, _Int, &cds_min)) { bsGetData(obj, _bsRight, _Int, &cds_max) ; if (cds_min < pos1 || cds_max > pos2) { messerror("The start/end of the CDS in object %s" " is outside of the spliced DNA coordinates for that object" " (spliced DNA start = %d & end = %d," " CDS start = %d & end = %d)." " Make sure the CDS positions have been specified" " in spliced DNA, not Source_exon, coordinates." " CDS start/end positions have been reset to start/end of" " spliced DNA.", name(parent), pos1, pos2, cds_min, cds_max) ; } else { pos1 = cds_min ; pos2 = cds_max ; } posConvert(seg1, pos1, pos2) ; } else { seg1->x1 = seg->x1 - 1 + pos1 ; seg1->x2 = seg->x2 ; } /* We need to record position and other CDS information for later */ /* drawing of the CDS. */ sinf->flags |= SEQ_CDS ; sinf->cds_info.cds_seg = *seg1 ; sinf->cds_info.cds_only = TRUE ; if (fMapDebug) printf("CDS: %d --> %d\n\n", seg1->x1, seg1->x2) ; if (!isExons) /* must have EXON to make CODING */ { seg2 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg2 = *seg ; seg2->source = parent ; seg2->key = 0 ; seg2->type = isDown ? EXON : EXON_UP ; sinf->flags |= SEQ_EXONS ; } } /* Tag signals that this object is incomplete, there is more upstream */ /* somewhere. Tag can be followed by a frame shift to correct any */ /* resultant reading frame error (0 means ignore setting, default). */ /* If there is a CDS, then it just start at position 1 for this tag to */ /* be valid. */ sinf->start_not_found = 0 ; if (bsFindTag (obj, _Start_not_found)) { if (sinf->cds_info.cds_only && (cds_min != 1)) { messerror("Data inconsistency: the Start_not_found tag is set for object %s," " but the CDS begins at position %d instead of at 1," " tag will be ignored", name(parent), cds_min) ; } else { /* Record any value found in seqinfo for later use. */ if (bsGetData(obj, _bsRight, _Int, &sinf->start_not_found)) { if (sinf->start_not_found < 1 || sinf->start_not_found > 3) { messerror("Data inconsistency: the Start_not_found tag for object %s" " should be given the value 1, 2 or 3 but has been" " set to %d," " tag will be ignored", name(parent), sinf->start_not_found) ; sinf->start_not_found = 0 ; } } } } /* remarks, clones etc */ if (bsFindTag (obj, _Visible) && bsFlatten (obj, 2, units)) { for (i = 0 ; i < arrayMax(units) ; i += 2) { /* RD 980414 checks to prevent NULL KEY in display moved here from Jean's checks in fMapShowVisible() */ if (!iskey (arr(units,i,BSunit).k) || class(arr(units,i,BSunit).k)) continue ; /* should be a tag */ if (!iskey (arr(units,i+1,BSunit).k)) continue ; /* there is a problem if the model says: Visible foo Text a hack because iskey() could be true for a text pointer should store these as different type, e.g. VISIBLE_TEXT */ seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; seg1->data.k = arr(units, i, BSunit).k ; seg1->key = arr(units, i+1, BSunit).k ; seg1->type = VISIBLE ; } } /* assembly tags */ if (bsFindTag (obj, _Assembly_tags) && bsFlatten (obj, 4, units)) for (i = 0 ; i < arrayMax(units)/4 ; ++i) { int length; char *string; char *type_text, *text_text; seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; /* get start and end of tag */ pos1 = arr(units, 4*i+1, BSunit).i; pos2 = arr(units, 4*i+2, BSunit).i; /* Convert co-ords of pos1, pos2 and stuff into seg as x1 and x2 */ posConvert(seg1, pos1, pos2); /* Key is the sort of assemly tag */ seg1->type = ASSEMBLY_TAG ; seg1->key = _Assembly_tags; /* Get the text */ length = 0; type_text = arr(units, 4*i, BSunit).s; if (type_text) length += strlen(type_text); text_text = arr(units, 4*i+3, BSunit).s; if (text_text) length += strlen(text_text); seg1->data.s = string = halloc(length+3, look->segsHandle); string[0] = '\0'; if (type_text) strcpy(string, type_text) ; strcat(string, ": "); if (text_text) strcat(string, text_text); seg1->parent = 0; /* assembly tags don't need parents */ } /* primers (directed sequencing */ /* mieg */ if (bsFindTag (obj, _Primer) && bsFlatten (obj, 1, units)) for (i = 0 ; i < arrayMax(units) ; ++i) { seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; seg1->key = arr(units,i, BSunit).k ; seg1->type = PRIMER ; } /* virtual sequence SCF */ /* mieg */ if (bsFindTag (obj, _SCF_File)) { seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; seg1->type = VIRTUAL_SUB_SEQUENCE_TAG ; if (!isDown) { if (seg->x1 > seg->x2) /* never happens i suppose */ { seg1->x1 = seg->x2 ; seg1->x2 = seg->x1 ; } seg1->type = VIRTUAL_SUB_SEQUENCE_TAG_UP ; } } /* Transcript, mieg */ pos1 = 0 ; if (bsFindTag (obj, _Transcribed_gene) && bsFlatten (obj, 3, units)) { for (i = 0 ; i < arrayMax(units) ; i += 3) { u = arrp(units,i,BSunit) ; if (!u[0].k) continue ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->key = seg1->parent = seg1->source = u[0].k ; seg1->data.i = 0 ; if (!u[1].i || !u[2].i) { /* messerror ("Coords of transcript %s missing in %s", name(seg1->key), name(seg->key)) ; */ --nsegs ; } else { pos1 = u[1].i ; pos2 = u[2].i ; if (isDown) seg1->type = (pos2 > pos1) ? TRANSCRIPT : TRANSCRIPT_UP ; else seg1->type = (pos1 > pos2) ? TRANSCRIPT_UP : TRANSCRIPT ; if (pos1 > pos2) { tmp = pos2 ; pos2 = pos1 ; pos1 = tmp ;} POS_TO_SEG1 ; if (seg1->x2 < 0 || seg1->x1 >= look->length) --nsegs ; /* ignore this transcript */ else { KEYSET tmp = queryKey (seg1->key, ">Annotations Fully_edited") ; if (keySetMax (tmp)) seg1->data.i = 0x1 ; keySetDestroy (tmp) ; tmp = queryKey(seg1->key, "Transpliced_to = SL1") ; if (keySetMax (tmp)) seg1->data.i |= 2 ; keySetDestroy (tmp) ; tmp = queryKey(seg1->key, "Transpliced_to = SL2") ; if (keySetMax (tmp)) seg1->data.i |= 4 ; keySetDestroy (tmp) ; if (seg1->data.i >= 6) seg1->data.i ^= 0x0F ; /* leave bit 1, zero bit 2 and 4 , set bit 8 */ } } } } /* transcript pieces, mieg */ pos1 = 0 ; if ((seg->type | 0x1) == TRANSCRIPT_UP && bsFindTag (obj, str2tag("Splicing")) && bsFlatten (obj, 4, units)) { seg1 = arrayp (segs,nsegs++,SEG) ; seg1->key = M_TRANSCRIPT ; seg1->type = (isDown) ? TRANSCRIPT : TRANSCRIPT_UP ; for (i = 0 ; i < arrayMax(units) ; i += 4) { char *cp ; if (arrayMax(units) < i + 3) continue ; pos1 = arr(units, i , BSunit).i ; pos2 = arr(units, i + 1, BSunit).i ; seg1 = arrayp (segs,nsegs++,SEG) ; POS_TO_SEG1 ; seg1->key = arr(units, i + 2, BSunit).k ; seg1->parent = seg->key ; seg1->source = parent ; seg1->type = (isDown) ? TRANSCRIPT : TRANSCRIPT_UP ; cp = arr(units, i + 3, BSunit).s ; if (cp && *cp && (seg1->key == str2tag("Intron") || seg1->key == str2tag("Alternative_Intron"))) if (!lexword2key(cp, &(seg1->data.k), 0)) lexaddkey ("Other", &seg1->data.k,0) ; } } /* spliced cDNA , mieg */ pos1 = 0 ; if ((seg->type | 0x1) == TRANSCRIPT_UP && bsFindTag (obj, _Assembled_from) && bsFlatten (obj, 5, units)) { seg1 = arrayp (segs,nsegs++,SEG) ; seg1->key = M_SPLICED_cDNA ; seg1->type = (isDown) ? SPLICED_cDNA : SPLICED_cDNA_UP ; for (i = 0 ; i < arrayMax(units) ; i += 5) { if (arrayMax(units) < i + 5) continue ; pos1 = arr(units, i + 0, BSunit).i ; pos2 = arr(units, i + 1, BSunit).i ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; POS_TO_SEG1 ; seg1->key = arr(units, i+2, BSunit).k ; seg1->parent = seg1->key ; seg1->source = parent ; pos1 = arr(units, i + 3, BSunit).i ; pos2 = arr(units, i + 4, BSunit).i ; /* reserve bit 31 and 30 for polyAinfo */ /* reserve bit 29 for ambiguous, 28 free and 30 for polyAinfo */ seg1->data.i = ((pos1 & 0x3FFF) << 14) + (pos2 & 0x3FFF) ; cDNAAlignInit () ; fMapcDNADoFillData(seg1) ; /* implies objCreate, unfortunate, but needed for ordering */ seg1->type = (isDown) ? SPLICED_cDNA : SPLICED_cDNA_UP ; seg1->source = seg->x1 ; } } /* virtual sequence assembly tags, */ /* mieg */ if ((seg->type | 0x1) == SEQUENCE_UP && bsFindTag (obj, _Assembled_from) && bsFlatten (obj, 3, units)) for (i = 0 ; i < arrayMax(units)/3 ; ++i) { int seg1index = nsegs ; seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; /* get start and end of tag */ pos1 = arr(units, 3*i+1, BSunit).i; pos2 = arr(units, 3*i+2, BSunit).i; /* Convert co-ords of pos1, pos2 and stuff into seg as x1 and x2 */ posConvert(seg1, pos1, pos2); if (seg1->type == SEQUENCE) { if (seg1->x1 > seg1->x2) seg1->type = VIRTUAL_SUB_SEQUENCE_TAG_UP ; else seg1->type = VIRTUAL_SUB_SEQUENCE_TAG ; } else /* sequence_up */ { if (seg1->x1 < seg1->x2) seg1->type = VIRTUAL_SUB_SEQUENCE_TAG_UP ; else seg1->type = VIRTUAL_SUB_SEQUENCE_TAG ; } if (seg1->x1 > seg1->x2) { int tmp = seg1->x1 ; seg1->x1 = seg1->x2 ; seg1->x2 = tmp ; } /* Key is the subsequence */ seg1->key = arr(units, 3*i, BSunit).k ; seg1->parent = 0 ; seg1->data.i = 0 ; { OBJ obj1 ; KEY _Proposed, _Stolen, _New_Read ; lexaddkey ("Proposed", &_Proposed, _VSystem) ; lexaddkey ("Stolen", &_Stolen, _VSystem) ; lexaddkey ("New_Read", &_New_Read, _VSystem) ; if ((obj1 = bsCreate(seg1->key))) { KEY key1, myCol ; if (bsFindTag(obj1, _Stolen)) seg1->data.i |= 1 << 24 ; if (bsFindTag(obj1, _Proposed)) seg1->data.i |= 1 << 25 ; if (bsFindTag(obj1, _Vector)) seg1->data.i |= 1 << 26 ; if (bsFindTag(obj1, _New_Read)) seg1->data.i |= 1 << 27 ; if (bsFindTag(obj1, _Significant_bases)) seg1->data.i |= 1 << 31 ; if (bsGetKeyTags(obj1, _Colour, &myCol)) seg1->data.i |= ((myCol - _WHITE + WHITE) & 0x3f) << 16 ; #ifdef ACEMBLY fMapColorSeg (seg1) ; #endif /* primers (directed sequencing */ /* mieg */ if (bsGetKey (obj1, _Primer, &key1)) do { seg2 = arrayp(segs,nsegs++,SEG) ; *seg2 = arr(segs,seg1index,SEG) ; seg2->key = key1 ; seg2->parent = key1 ; /* so both half co color */ seg2->data.i = 0 ; seg2->type = PRIMER | (seg1->type & 1) ; } while (bsGetKey (obj, _bsDown, &key1)) ; bsDestroy (obj1) ; } } } /* virtual sequence assembly tags, dummy seq */ /* ulrich */ if (bsFindTag (obj, _Previous_contig) && bsFlatten (obj, 3, units)) for (i = 0 ; i < arrayMax(units)/3 ; ++i) { seg1 = arrayp(segs, nsegs++, SEG) ; seg = arrp(segs, iseg, SEG) ; *seg1 = *seg ; seg1->source = parent ; /* get start and end of tag */ pos1 = arr(units, 3*i+1, BSunit).i; pos2 = arr(units, 3*i+2, BSunit).i; /* Convert co-ords of pos1, pos2 and stuff into seg as x1 and x2 */ posConvert(seg1, pos1, pos2); if (seg1->type == SEQUENCE) { if (seg1->x1 > seg1->x2) seg1->type = VIRTUAL_PREVIOUS_CONTIG_TAG_UP ; else seg1->type = VIRTUAL_PREVIOUS_CONTIG_TAG ; } else /* sequence_up */ { if (seg1->x1 < seg1->x2) seg1->type = VIRTUAL_PREVIOUS_CONTIG_TAG_UP ; else seg1->type = VIRTUAL_PREVIOUS_CONTIG_TAG ; } if (seg1->x1 > seg1->x2) { int tmp = seg1->x1 ; seg1->x1 = seg1->x2 ; seg1->x2 = tmp ; } /* Key is the subsequence */ seg1->key = arr(units, 3*i, BSunit).k ; seg1->parent = 0 ; seg1->data.i = 0 ; } /* virtual alignment */ /* mieg */ if (bsFindTag (obj, _Aligned) && bsFlatten (obj, 3, units)) for (i = 0 ; i < arrayMax(units)/3 ; ++i) { seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; /* get start and end of tag */ pos1 = arr(units, 3*i+1, BSunit).i; pos2 = arr(units, 3*i+2, BSunit).i; /* Convert co-ords of pos1, pos2 and stuff into seg as x1 and x2 */ posConvert(seg1, pos1, pos2); if (seg1->type == SEQUENCE) { if (seg1->x1 > seg1->x2) seg1->type = VIRTUAL_ALIGNED_TAG_UP ; else seg1->type = VIRTUAL_ALIGNED_TAG ; } else /* sequence_up */ { if (seg1->x1 < seg1->x2) seg1->type = VIRTUAL_ALIGNED_TAG_UP ; else seg1->type = VIRTUAL_ALIGNED_TAG ; } if (seg1->x1 > seg1->x2) { int tmp = seg1->x1 ; seg1->x1 = seg1->x2 ; seg1->x2 = tmp ; } /* Key is the subsequence */ seg1->key = arr(units, 3*i, BSunit).k ; seg1->parent = 0 ; seg1->data.i = 0 ; } /* virtual sequence assembly tags, */ /* mieg */ if (bsGetKey (obj, _Assembled_into, &key)) do { seg1 = arrayp(segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; *seg1 = *seg ; seg1->source = parent ; seg1->key = key ; seg1->type = VIRTUAL_PARENT_SEQUENCE_TAG ; seg1->parent = seg->key ; seg1->data.i = 0 ; { OBJ obj1 ; KEY Proposed ; lexaddkey ("Proposed", &Proposed, _VSystem) ; if ((obj1 = bsCreate(seg1->key))) { if (bsFindTag(obj1, Proposed)) seg1->data.i = -1 ; bsDestroy (obj1) ; } } seg1->x1 = seg->x1 ; seg1->x2 = seg->x2 ; } while (bsGetKey (obj, _bsDown, &key)) ; #define POS_PROCESS(Z,Z_UP) \ if (isDown) seg1->type = (pos2 > pos1) ? Z : Z_UP ; \ else seg1->type = (pos1 > pos2) ? Z : Z_UP ; \ if (pos1 > pos2) { tmp = pos2 ; pos2 = pos1 ; pos1 = tmp ;} \ posConvert (seg1, pos1, pos2) ; /* Oligos, mieg */ if (bsFindTag (obj, _Oligo) && bsFlatten (obj, 3, units)) for (i = 0 ; i < arrayMax (units) ; i += 3) { OBJ obj1 ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = arr(units,i, BSunit).k ; seg1->parent = seg1->key ; seg1->data.i = 0 ; pos1 = arr(units,i+1, BSunit).i ; pos2 = arr(units,i+2, BSunit).i ; if (!pos1 || !pos2) fMapOspPositionOligo (look, seg, seg1->key, &pos1, &pos2) ; if (!pos1 || !pos2) continue ; POS_PROCESS(OLIGO,OLIGO_UP) ; if ((obj1 = bsCreate(seg1->key))) { float xf ; int xi ; if (bsGetData(obj1, _Tm, _Float, &xf)) { xi = 10 * xf ; seg1->data.i = (xi & 0x3FF) << 16 ; } if (bsGetData(obj1, _Score, _Float, &xf)) { xi = xf ; seg1->data.i |= (xi & 0xfff) ; } if (!bsFindTag(obj1, _Temporary)) seg1->data.i |= 0x4000 ; /* old */ bsDestroy (obj1) ; } } /* Alleles */ if (bsFindTag (obj, _Allele) && bsFlatten (obj, 4, units)) for (i = 0 ; i < arrayMax (units) ; i += 4) { char *string ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = arr(units,i, BSunit).k ; seg1->parent = 0 ; seg1->data.i = 0 ; pos1 = arr(units,i+1, BSunit).i ; pos2 = arr(units,i+2, BSunit).i ; if (pos1 || pos2) seg1->parent = seg1->key ; /* flag to draw it */ else seg1->parent = seg->key ; /* couple with parent */ POS_PROCESS(ALLELE,ALLELE_UP) ; if ((string = arr(units,i+3,BSunit).s)) seg1->data.s = strnew (string, look->segsHandle) ; else seg1->data.s = 0 ; } /* EMBL features */ if (bsFindTag(obj, _EMBL_feature) && bsFlatten(obj, 4, units)) for (i = 0 ; i < arrayMax (units) ; i += 4) { char *string ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->parent = nsegs ; seg1->key = arr(units,i, BSunit).k ; pos1 = arr(units,i+1, BSunit).i ; pos2 = arr(units,i+2, BSunit).i ; POS_PROCESS(EMBL_FEATURE, EMBL_FEATURE_UP) ; /* 4'th column: can't distinguish KEY from Text securely assume if iskey() that it is a key, e.g. ?Text */ if (iskey (key = arr(units,i+3,BSunit).k)) seg1->data.s = strnew (name(key), look->segsHandle) ; else if ((string = arr(units,i+3,BSunit).s)) seg1->data.s = strnew (string, look->segsHandle) ; else seg1->data.s = 0 ; } /* Homologies */ /* Model looks like: Homol Float Int UNIQUE Int Int UNIQUE Int #Homo_info #Homol_info Align Int UNIQUE Int UNIQUE Int // not needed if full alignment is ungapped // one row per ungapped alignment section // first int is position in query, second in target // third int is length - only necessary when gaps in // both sequences align. mark1 cursors down mark2 cursors down mark3 cursors down Float mark4 cursors down Int mark5 cursors down Int */ { HOMOLINFO *hinf ; BSMARK mark1 = 0, mark2 = 0, mark3 = 0, mark4 = 0, mark5 = 0, mark6 = 0; KEY method, homol, tag2; float score; int target1, target2; Array ga = arrayCreate(16, BSunit); if (bsGetKeyTags(obj, _Homol, &tag2)) do { mark1 = bsMark(obj, mark1); if (bsGetKey(obj, _bsRight, &homol)) do { mark2 = bsMark(obj, mark2); if (bsGetKey(obj, _bsRight, &method)) do { mark3 = bsMark(obj, mark3); if (bsGetData(obj, _bsRight, _Float, &score)) do { mark4 = bsMark(obj, mark4); if (bsGetData(obj, _bsRight, _Int, &pos1)) do { mark5 = bsMark(obj, mark5); if (bsGetData(obj, _bsRight, _Int, &pos2) && bsGetData(obj, _bsRight, _Int, &target1)) do { mark6 = bsMark(obj, mark6); if (bsGetData(obj, _bsRight, _Int, &target2)) { char *feature = name(homol) ; int pos1save = pos1, pos2save = pos2; seg1 = arrayp (segs,nsegs++,SEG); seg1->source = parent ; seg1->key = seg1->parent = homol; /* target */ seg1->data.i = arrayMax(look->homolInfo) ; hinf = arrayp(look->homolInfo, seg1->data.i, HOMOLINFO) ; hinf->method = method; hinf->score = score; bsPushObj(obj); /* into Homo_info */ hinf->gaps = sMapLocalMap(obj, pos1, pos2, look->handle); if (target1 > target2) { int tmp = pos1; pos1 = pos2; pos2 = tmp; hinf->x1 = target2; hinf->x2 = target1; } else { hinf->x1 = target1; hinf->x2 = target2; } POS_PROCESS(HOMOL,HOMOL_UP) ; pos1 = pos1save; pos2 = pos2save; /* in case they got flipped */ } bsGoto(obj, mark6); } while (bsGetData(obj, _bsDown, _Int, &target1)); bsGoto(obj, mark5); } while (bsGetData(obj, _bsDown, _Int, &pos1)); bsGoto(obj, mark4); } while (bsGetData(obj, _bsDown, _Float, &score)); bsGoto(obj, mark3); } while (bsGetKey(obj, _bsDown, &method)); bsGoto(obj, mark2); } while (bsGetKey(obj, _bsDown, &homol)); bsGoto(obj, mark1); } while (bsGetKeyTags(obj, _bsDown, &tag2)); arrayDestroy(ga); bsMarkFree(mark1); bsMarkFree(mark2); bsMarkFree(mark3); bsMarkFree(mark4); bsMarkFree(mark5); bsMarkFree(mark6); } /* MatchTable - more homologies */ { Array ma ; MATCH *m ; HOMOLINFO *hinf ; if (bsGetKey (obj, str2tag("Match_table"), &key) && (ma = arrayGet (key, MATCH, matchFormat))) { FLAGSET fBury = flag ("Match", "Bury") ; for (i = 0 ; i < arrayMax(ma) ; ++i) { m = arrp (ma, i, MATCH) ; if (class(m->meth) != _VMethod || !m->q1 || !m->q2) continue ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = key ; /* the match table key */ seg1->key = seg1->parent = m->key ; seg1->data.i = arrayMax(look->homolInfo) ; hinf = arrayp(look->homolInfo, seg1->data.i, HOMOLINFO) ; hinf->method = m->meth ; hinf->score = m->score ; hinf->gaps = 0; bitSet (look->homolFromTable, seg1->data.i) ; if (m->t1 > m->t2) { pos1 = m->q2 ; pos2 = m->q1 ; hinf->x1 = m->t2 ; hinf->x2 = m->t1 ; } else { pos1 = m->q1 ; pos2 = m->q2 ; hinf->x1 = m->t1 ; hinf->x2 = m->t2 ; } POS_PROCESS(HOMOL,HOMOL_UP) ; if (m->flags & fBury) bitSet (look->homolBury, seg1->data.i) ; } arrayDestroy (ma) ; } } /* necessary to ensure array long enough */ bitExtend (look->homolBury, arrayMax(look->homolInfo)) ; bitExtend (look->homolFromTable, arrayMax(look->homolInfo)) ; /* Features */ if (bsFindTag(obj, _Feature) && bsFlatten(obj, 5, units)) for (i = 0 ; i < arrayMax(units) ; i += 5) { u = arrayp(units, i, BSunit) ; if (class(u[0].k) != _VMethod || !u[1].i || !u[2].i) continue ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = u[0].k ; pos1 = u[1].i ; pos2 = u[2].i ; POS_PROCESS(FEATURE,FEATURE_UP) ; seg1->data.f = u[3].f ; if (u[4].s) { dictAdd (look->featDict, u[4].s, &tmp) ; seg1->parent = -tmp ; } else seg1->parent = 0 ; } /* Splice sites */ { KEY _Predicted_5 = str2tag("Predicted_5") ; KEY _Predicted_3 = str2tag("Predicted_3") ; if (bsFindTag(obj, str2tag("Splices")) && bsFlatten(obj, 5, units)) for (i = 0 ; i < arrayMax(units) ; i += 5) { u = arrayp(units, i, BSunit) ; if (u[0].k != _Predicted_5 && u[0].k != _Predicted_3) continue ; if (class(u[1].k) != _VMethod || !u[2].i || !u[3].i) continue ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = u[1].k ; pos1 = u[2].i ; pos2 = u[3].i ; if (u[0].k == _Predicted_5) { POS_PROCESS(SPLICE5,SPLICE5_UP) ; } else { POS_PROCESS(SPLICE3,SPLICE3_UP) ; } seg1->data.f = u[4].f ; } } /* clone end information */ { int z ; KEY tag ; for (z = 0 ; z < 2 ; ++z) { tag = z ? _Clone_left_end : _Clone_right_end ; if (bsFindTag(obj, tag) && bsFlatten(obj,2,units)) for (i = 0 ; i < arrayMax(units) ; i += 2) { u = arrayp(units, i, BSunit) ; seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = u[0].k ; /* Clone */ seg1->type = CLONE_END ; pos1 = u[1].i ; pos2 = u[1].i ; posConvert(seg1, pos1, pos2) ; seg1->data.k = tag ; lexaddkey (name(seg1->key), &key, _VSequence) ; seg1->parent = key ; } } } /* confirmed introns: add all now, and remove duplicates with those from genes later in fMapFindCoding() */ { KEY _Confirmed_intron = str2tag("Confirmed_intron") ; KEY _EST = str2tag("EST") ; KEY _cDNA = str2tag("cDNA") ; KEY _Homology = str2tag("Homology") ; KEY _UTR = str2tag("UTR") ; KEY _False = str2tag("False"); if (bsFindTag (obj, _Confirmed_intron) && bsFlatten (obj, 3, units)) for (i = 0 ; i < arrayMax(units) ; i += 3) { u = arrayp(units, i, BSunit) ; if (!u[0].i || !u[1].i) continue ; if (!i || u[-3].i != u[0].i || u[-2].i != u[1].i) { seg1 = arrayp (segs,nsegs++,SEG) ; seg = arrp(segs,iseg,SEG) ; seg1->source = parent ; seg1->key = 0 ; pos1 = u[0].i ; pos2 = u[1].i ; seg1->data.i = 0 ; POS_PROCESS(INTRON,INTRON_UP) ; } seg1->data.i = arrayMax(look->seqInfo) ; sinf = arrayp (look->seqInfo, seg1->data.i, SEQINFO) ; if (u[2].k == _EST) sinf->flags |= SEQ_CONFIRM_EST ; else if (u[2].k == _cDNA) sinf->flags |= SEQ_CONFIRM_CDNA ; else if (u[2].k == _Homology) sinf->flags |= SEQ_CONFIRM_HOMOL ; else if (u[2].k == _UTR) sinf->flags |= SEQ_CONFIRM_UTR ; else if (u[2].k == _False) sinf->flags |= SEQ_CONFIRM_FALSE ; else sinf->flags |= SEQ_CONFIRM_UNKNOWN ; } } bsDestroy (obj) ; } /* END OF SEGS LOOP....PHEW... */ if (look->flag & FLAG_COMPLEMENT) { int top = look->length - 1 ; seg = arrp (look->segs, 1, SEG) ; for (i = 1 ; i < arrayMax(segs) ; ++i, ++seg) { tmp = seg->x1 ; seg->x1 = top - seg->x2 ; seg->x2 = top - tmp ; if (seg->type >= SEQUENCE && seg->type <= ALLELE_UP) seg->type ^= 1 ; } } if (oldSegs) { arrayMax(oldSegs) = look->lastTrueSeg ; /* still the old value */ addOldSegs (segs, oldSegs, look->mcache) ; handleDestroy (oldSegsHandle) ; /* kills old segs and related stuff */ oldSegs = NULL ; } arraySort (segs, fMapOrder) ; /* must sort for FindCoding and */ /* RemoveSelfHomol */ if (!(result = fMapFindCoding(look))) /* make CODING segs */ goto abort ; fMapRemoveSelfHomol (look) ; fMapTraceFindMultiplets (look) ; /* make virtual multiplets */ fMapProcessMethods (look) ; /* always do this - to pick up edits */ arraySort (segs, fMapOrder) ; /* these 2 routines no longer work after reverse complement because RC calls sort * they should be written differently */ fMapOspFindOligoPairs (look) ; /* why AFTER sorting ? not sure, but probably ok */ look->lastTrueSeg = arrayMax(segs) ; result = TRUE ; /* belt and braces. */ abort: /* Currently the caller is responsible for cleaning up fmap resources */ /* acquired by this routine, don't know if this should continue. */ /* But we do need to clean up oldsegs... */ if (oldSegs) { handleDestroy (oldSegsHandle) ; /* kills old segs and related stuff */ oldSegs = NULL ; } chronoReturn() ; { int i ; SEG *seg ; for (i = 0 ; i < arrayMax(look->segs) ; i++) { seg = arrp(look->segs, i, SEG) ; } } return result ; } /*****************************************************************/ /*****************************************************************/ static char *title (KEY key) { OBJ obj ; char *result ; KEY textKey ; result = name(key) ; if ((obj = bsCreate(key))) { if (bsGetKey (obj, _Title, &textKey)) result = name(textKey) ; bsDestroy (obj) ; } return result ; } /**********************************************************/ void fMapReportLine (FeatureMap look, SEG *seg, BOOL isGIF, float x) { char *cp ; HOMOLINFO *hinf ; SEQINFO *sinf ; METHOD *meth ; int i, i2, c1, ii ; SEG *seg2 ; memset (look->segTextBuf, 0, 63) ; /* object name */ if (isGIF) *look->segTextBuf = 0 ; else if (iskey(seg->key)) { strncpy (look->segTextBuf, name(seg->key), 40) ; look->segTextBuf[40] = '\0'; } else if (seg->parent && iskey(seg->parent)) { strncpy (look->segTextBuf, name(seg->parent), 40) ; look->segTextBuf[40] = '\0'; } else *look->segTextBuf = 0 ; /* coordinates */ if (look->flag & FLAG_REVERSE) strcat (look->segTextBuf, messprintf (" %d %d (%d) ", COORD(look, seg->x2), COORD(look, seg->x1), 1 - COORD(look, seg->x2) + COORD(look, seg->x1))) ; else strcat (look->segTextBuf, messprintf (" %d %d (%d) ", COORD(look, seg->x1), COORD(look, seg->x2), 1 - COORD(look, seg->x1) + COORD(look, seg->x2))) ; /* description */ if (seqDragBox) seg2 = BOXSEG (seqDragBox) ; else seg2 = seg ; switch (seg->type) { case SEQUENCE: case SEQUENCE_UP: case EXON: case EXON_UP: case INTRON: case INTRON_UP: sinf = arrp (look->seqInfo, seg->data.i, SEQINFO) ; if (sinf->method) strncat (look->segTextBuf, name(sinf->method), 40) ; if (sinf->flags & SEQ_SCORE) strncat (look->segTextBuf, messprintf ("%.1f", sinf->score), 10) ; if (sinf->flags & SEQ_CONFIRMED) { strcat (look->segTextBuf, "Confirmed") ; if (sinf->flags & SEQ_CONFIRM_CDNA) strcat (look->segTextBuf, " by_cDNA") ; if (sinf->flags & SEQ_CONFIRM_EST) strcat (look->segTextBuf, " by_EST") ; if (sinf->flags & SEQ_CONFIRM_HOMOL) strcat (look->segTextBuf, " by_homology") ; if (sinf->flags & SEQ_CONFIRM_UTR) strcat (look->segTextBuf, " in_UTR") ; if (sinf->flags & SEQ_CONFIRM_FALSE) strcat (look->segTextBuf, " as_false") ; } break ; case DNA_SEQ: i = seg->x1 + (int)x ; i2 = seg2->x1 + (int)newx ; if (look->flag & FLAG_REVERSE) strcat (look->segTextBuf, messprintf ("Selection %d ---> %d (%d)", COORD(look, i2), COORD(look, i), 1 - COORD(look, i2)+ COORD(look, i))) ; else strcat (look->segTextBuf, messprintf ("Selection %d ---> %d (%d)", COORD(look, i), COORD(look, i2), 1 + COORD(look, i2)- COORD(look, i))) ; if (look->gf.cum && i >= look->gf.min && i < look->gf.max) strncat (look->segTextBuf, messprintf (" cumulative coding score %.2f", look->gf.cum[i - look->gf.min]), 64) ; break ; case PEP_SEQ: i = seg->x1 + 3*(int)x ; i2 = seg2->x1 + 3*(int)newx ; if ((cp = pepName[(int) codon (arrp(look->dna, i, char))])) { if (look->flag & FLAG_REVERSE) strncat (look->segTextBuf, messprintf ("%s at %d to %d Selection weight: %d",cp, COORD(look, i+2), COORD(look,i), weight), 64) ; else strncat (look->segTextBuf, messprintf ("%s at %d to %d Selection weight: %d",cp, COORD(look, i), COORD(look,i+2), weight), 64) ; } break ; case TRANS_SEQ: i = seg->x1 + 3*(int)x ; c1 = COORD(look, i) ; if ((cp = pepName[(int) codon (arrp(look->dna, i, char))])) { if (look->flag & FLAG_REVERSE) { if(!seqDragBox) nbTrans = 1 ; strncat (look->segTextBuf, messprintf ("%s at %d ", cp, c1 - nbTrans + 1), 64) ; if (seqDragBox && c1 <= array(look->minDna, 0, int)) { for (ii=1; iimaxDna); ii++) if (c1 > array(look->maxDna, ii, int) && c1 <= array(look->minDna, ii- 1, int)) { c1 = (array(look->decalCoord, ii-1, int) - c1)/3 ; break ; } strcat (look->segTextBuf, messprintf (" Selection %d ---> %d (%d) weight: %d", c1, c1 + nbTrans - 1, nbTrans, weight)) ; } } else { strncat (look->segTextBuf, messprintf ("%s at %d ", cp, c1), 64) ; if (seqDragBox && c1 >= array(look->minDna, 0, int)) { for (ii=1; iimaxDna); ii++) if (c1 < array(look->maxDna, ii, int) && c1 >= array(look->minDna, ii- 1, int)) { c1 = (c1 - array(look->decalCoord, ii-1, int))/3 ; break ; } strcat (look->segTextBuf, messprintf ( " Selection %d ---> %d (%d) weight: %d", c1, c1 + nbTrans - 1, nbTrans, weight)) ; } } } break ; case TRANS_SEQ_UP: i = seg->x1 + 3*(int)x ; c1 = COORD(look, i) ; if ((cp = pepName[(int) reverseCodon (arrp(look->dna, i-2, char))])) { if (look->flag & FLAG_REVERSE) { if(!seqDragBox) nbTrans = 1 ; strncat (look->segTextBuf, messprintf ("up-strand %s at %d ", cp, c1 + nbTrans - 1), 64) ; if (seqDragBox && c1 <= array(look->minDna, 0, int)) { for (ii=1; iimaxDna); ii++) if (c1 > array(look->maxDna, ii, int) && c1 <= array(look->minDna, ii- 1, int) ) { c1 = (c1 - array(look->decalCoord, ii-1, int))/3 ; break ; } strcat (look->segTextBuf, messprintf (" Selection %d ---> %d (%d) weight: %d", c1, c1 - nbTrans + 1, nbTrans, weight)) ; } } else { strncat (look->segTextBuf, messprintf ("up-strand %s at %d ", cp, c1), 64) ; if (seqDragBox && c1 >= array(look->minDna, 0, int)) { for (ii=1; iimaxDna); ii++) if (c1 < array(look->maxDna, ii, int) && c1 >= array(look->minDna, ii- 1, int)) { c1 = (array(look->decalCoord, ii-1, int) - c1)/3 ; break ; } strcat (look->segTextBuf, messprintf ( " Selection %d ---> %d (%d) weight: %d", c1, c1 - nbTrans + 1, nbTrans, weight)) ; } } } break ; case ORF: strncat (look->segTextBuf, messprintf ("Frame %d", seg->data.i), 64) ; break ; case SPLICE3: case SPLICE5: case ATG: strncat (look->segTextBuf, fMapSegTypeName[seg->type], 32) ; strncat (look->segTextBuf, messprintf (" Score %f", seg->data.f), 32) ; break ; case ASSEMBLY_TAG: strncat (look->segTextBuf, messprintf ("%s", seg->data.s), 64); break; case ALLELE: case ALLELE_UP: case EMBL_FEATURE: case EMBL_FEATURE_UP: if (seg->data.s) strncat (look->segTextBuf, messprintf ("%s", seg->data.s), 64) ; break ; case HOMOL: case HOMOL_UP: hinf = arrp(look->homolInfo, seg->data.i, HOMOLINFO) ; meth = methodCacheGet(look->mcache, hinf->method, "") ; if (meth->flags & METHOD_BLASTN) strncat (look->segTextBuf, messprintf("%s %.1f %d%% (%d - %d) ", name(hinf->method), hinf->score, (int) (100 * (4 + hinf->score/(seg->x2-seg->x1+1)) / 9.0), hinf->x1, hinf->x2), 64) ; else strncat (look->segTextBuf, messprintf("%s %.1f (%d - %d) ", name(hinf->method), hinf->score, hinf->x1, hinf->x2), 64) ; if (!isGIF) /* title() requires DB access */ strncat (look->segTextBuf, title(seg->key), 127 - strlen(look->segTextBuf)) ; break ; case FEATURE: case FEATURE_UP: meth = methodCacheGet(look->mcache, seg->key, "") ; if (meth->flags & METHOD_PERCENT) strcat (look->segTextBuf, messprintf ("%.0f%% ", seg->data.f)) ; else if (meth->flags & METHOD_SCORE) strcat (look->segTextBuf, messprintf ("%.3g ", seg->data.f)) ; if (seg->parent) strncat (look->segTextBuf, dictName (look->featDict, -seg->parent), 58) ; break ; case CODING: if (look->gf.cum && seg->x1 >= look->gf.min && seg->x2 < look->gf.max) strncat (look->segTextBuf, messprintf ("Coding score %.2f", look->gf.cum[seg->x2 + 1 - look->gf.min] - look->gf.cum[seg->x1 - look->gf.min]), 64) ; break ; case SPLICED_cDNA: case SPLICED_cDNA_UP: fMapcDNAReportLine (look->segTextBuf,seg,64) ; break ; case TRANSCRIPT: case TRANSCRIPT_UP: fMapcDNAReportTranscript (look->segTextBuf,seg,64) ; break ; case VIRTUAL_MULTIPLET_TAG: fMapSelectVirtualMultiplet (look, seg) ; break ; case CLONE_END: strncat (look->segTextBuf, name(seg->data.k), 64) ; /* left/right */ break ; case OLIGO: case OLIGO_UP: i = (seg->data.i >> 16) & 0x3ff ; if (i > 0) strcat (look->segTextBuf, messprintf (" Tm: %3.1f ", (float)(i/10.0))) ; i = seg->data.i & 0xfff ; if (i > 0) strcat (look->segTextBuf, messprintf (" Score: %d ", i)) ; break ; case OLIGO_PAIR: case OLIGO_PAIR_UP: fMapOspSelectOligoPair (look, seg) ; break ; case VISIBLE: if (iskey (seg->data.k)) strcat (look->segTextBuf, messprintf (" (%s) ", name(seg->data.k))) ; break ; default: /* needed for picky compiler */ break ; } } /**********************************************************/ static void fMapSelectBox (FeatureMap look, int box, double x, double y) { int i, i0, max , imax, w=0 ; SEG *seg, *seg2, *seg3; KEY parent ; unsigned char *cp ; char *buf, *cq, *cq1 ; HOMOLINFO *hinf ; METHOD *meth ; static int finDna = 0 ; int ii, c1 ; BOOL isProt ; int nExon=0, tframe, oldTframe=0, diffTframe ; BOOL isTiret, calcWeight ; BOOL doComplementSurPlace = (look->flag & FLAG_COMPLEMENT_SUR_PLACE) ; if (box >= arrayMax(look->boxIndex) || box < look->minLiveBox) return ; if (look->activeBox) { seg = BOXSEG(look->activeBox) ; parent = seg->parent ; for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; ++i) { seg2 = BOXSEG(i) ; if (seg2 == seg || (parent && seg2->parent == parent)) switch (seg2->type) { case DNA_SEQ: case PEP_SEQ: case TRANS_SEQ: case TRANS_SEQ_UP: break ; case HOMOL: case HOMOL_UP: hinf = arrp(look->homolInfo, seg2->data.i, HOMOLINFO) ; meth = methodCacheGet(look->mcache, hinf->method, "") ; graphBoxDraw (i, -1, meth->colour) ; break ; case FEATURE: case FEATURE_UP: meth = methodCacheGet(look->mcache, seg2->key, "") ; graphBoxDraw (i, -1, meth->colour) ; break ; case ASSEMBLY_TAG: { int color ; unsigned char *c = (unsigned char*) seg2->data.s; switch(c[0] ^ c[1] ^ c[2] ^ c[3]) { case 'a'^'n'^'n'^'o': color = BLACK; break; case 'c'^'o'^'m'^'m': color = BLUE; break; case 'o'^'l'^'i'^'g': color = YELLOW; break; case 'c'^'o'^'m'^'p': color = RED; break; case 's'^'t'^'o'^'p': color = RED; break; case 'r'^'e'^'p'^'e': color = GREEN; break; case 'c'^'o'^'s'^'m': color = LIGHTGRAY; break; case 'A'^'l'^'u'^' ': color = GREEN; break; case 'i'^'g'^'n'^'o': color = LIGHTGRAY; break; default: color = LIGHTBLUE; break; } graphBoxDraw (i, BLACK, color) ; } break; case INTRON: case INTRON_UP: if (arrp(look->seqInfo,seg2->data.i,SEQINFO)->flags & SEQ_CONFIRMED) graphBoxDraw (i, -1, (arrp(look->seqInfo,seg2->data.i,SEQINFO)->flags & SEQ_CONFIRM_UTR ? YELLOW : GREEN)) ; else graphBoxDraw (i, -1, WHITE) ; break ; case OLIGO: case OLIGO_UP: if (seg2->data.i & 0x4000) /* old */ graphBoxDraw (i, -1, ORANGE) ; else if (seg2->data.i & 0x8000) /* selected */ graphBoxDraw (i, -1, GREEN) ; else graphBoxDraw (i, -1, WHITE) ; break ; default: /* Includes exon drawing. */ if (assFind (look->chosen, SEG_HASH(seg2), 0)) graphBoxDraw (i, -1, GREEN) ; else if (assFind (look->antiChosen, SEG_HASH(seg2), 0)) graphBoxDraw (i, -1, LIGHTGREEN) ; else graphBoxDraw (i, -1, WHITE) ; break ; } } i = seg->x1 ; /* mhmp 11.06.98 + 08.09.98 */ if (look->isRetro) i -=2 ; if (i <= 0) i = 0 ; if (!finDna) finDna = (seg->x2 >= look->length) ? look->length - 1 : seg->x2 ; if (look->colors) { cp = arrp(look->colors, i, unsigned char) ; for ( ; i <= finDna && i < arrayMax(look->colors) ; ++i) *cp++ &= ~(TINT_HIGHLIGHT1 | TINT_HIGHLIGHT2); } for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; ++i) { seg2 = BOXSEG(i) ; if ((seg2->type == DNA_SEQ || seg2->type == PEP_SEQ || seg2->type == TRANS_SEQ || seg2->type == TRANS_SEQ_UP) && seg2->x1 <= finDna && seg2->x2 >= seg->x1) graphBoxDraw (i, -1, -1) ; } } look->activeBox = box ; seg = BOXSEG(look->activeBox) ; if (seqDragBox) seg3 = BOXSEG(seqDragBox) ; else seg3 = seg ; /* record exact dna postition for gf add splice site functions */ if (seg->type == DNA_SEQ) look->DNAcoord = seg->x1 + (int)x; i = (seg->x1 > 0) ? seg->x1 : 0 ; max = (seg->x2 >= look->length) ? look->length - 1 : seg->x2 ; finDna = max ; if (seg->type == DNA_SEQ && seg3->type == DNA_SEQ) { i = i + (int)x ; i0 = i ; imax = seg3->x1 + (int)newx ; if (i0 <= imax) /* mhmp 27.04.98 A VOIR: i0 <--> imax*/ { finDna = (seg3->x2 >= look->length) ? look->length - 1 : seg3->x2 ; cp = arrp(look->colors, i, unsigned char) ; for ( ; i <= imax ; ++i) *cp++ |= TINT_HIGHLIGHT1; /* mhmp 27.04.98 dna en lignes de 50 */ buf = messalloc((imax - i0 + 2) *51/50) ; cq = buf ; cq1 = arrp(look->dna, i0, char) ; i = imax ; while(i-- >= i0 ) { *cq++ = doComplementSurPlace ? dnaDecodeChar[(int)complementBase [(int)*cq1++]] : dnaDecodeChar[(int)*cq1++] ; if ( !((imax - i) % 50) ) *cq++ = '\n' ; } *cq++ = 0 ; graphPostBuffer(buf) ; messfree(buf) ; } } else if (seg->type == PEP_SEQ || seg->type == TRANS_SEQ || seg->type == TRANS_SEQ_UP) { i = i + 3*(int)x ; i0 = i ; imax = seg3->x1 + 3*(int)newx ; finDna = seg3->x2 + 2 ; /* mhmp 07.10.98 */ finDna = (finDna >= look->length) ? look->length - 1 : finDna ; cp = arrp(look->colors, i, unsigned char) ; if (i0 <= imax) /* mhmp 27.04.98 A VOIR: i0 <--> imax*/ { if (i0 < imax) for ( ; i <= imax + 2 ; ++i) *cp++ |= TINT_HIGHLIGHT1; else { cp = arrp(look->colors, seg->x1 + 3*(int)x, unsigned char) ; if (seg->type == TRANS_SEQ_UP) /* mhmp 10.06.98 */ { *cp-- |= TINT_HIGHLIGHT1; *cp-- |= TINT_HIGHLIGHT2; *cp-- |= TINT_HIGHLIGHT2; look->isRetro = TRUE ; } else { *cp++ |= TINT_HIGHLIGHT1; *cp++ |= TINT_HIGHLIGHT2; *cp++ |= TINT_HIGHLIGHT2; look->isRetro = FALSE ; } } /* mhmp 27.04.98 proteines en lignes de 50 */ buf = messalloc((imax - i0 + 2) *51/50) ; cq = buf ; cq1 = arrp(look->dna, i0, char) ; i = i0 ; c1 = COORD (look, i); if (seg->type != PEP_SEQ) { for (ii=1; iimaxDna); ii++) if (look->flag & FLAG_REVERSE) { if (c1 > array(look->maxDna, ii, int) && c1 <= array(look->minDna, ii- 1, int)) { nExon = ii - 1; oldTframe = array(look->tframe, ii- 1, int) ; break ; } } else { if (c1 < array(look->maxDna, ii, int) && c1 >= array(look->minDna, ii- 1, int)) { nExon = ii - 1; oldTframe = array(look->tframe, ii- 1, int) ; break ; } } isTiret = FALSE ; nbTrans = 0 ; calcWeight = TRUE ; weight = 18 ; while(i <= imax ) { isProt = FALSE ; c1 = COORD (look, i); for (ii = nExon; iimaxDna); ii++) if (look->flag & FLAG_REVERSE) { if (c1 > array(look->maxDna, ii, int) && c1 <= array(look->minDna, ii- 1, int)) { isProt = TRUE ; isTiret = FALSE ; break ; } } else if (c1 < array(look->maxDna, ii, int) && c1 >= array(look->minDna, ii- 1, int)) { isProt = TRUE ; isTiret = FALSE ; break ; } if (isProt) { if (seg->type == TRANS_SEQ) *cq++ = doComplementSurPlace ? antiCodon (cq1) : codon(cq1) ; else /* TRANS_SEQ_UP */ *cq++ = doComplementSurPlace ? codon (cq1) : antiCodon (cq1) ; nbTrans++ ; if (calcWeight) { if (seg->type == TRANS_SEQ) w = doComplementSurPlace ? molecularWeight[(int) antiCodon (cq1)] : molecularWeight[(int) codon(cq1)] ; else /* TRANS_SEQ_UP */ w = doComplementSurPlace ? molecularWeight[(int) codon(cq1)] : molecularWeight[(int) antiCodon (cq1)] ; } if (w < 0) calcWeight = FALSE ; else if (w == 0) { weight = 0 ; calcWeight = FALSE ; } else weight = weight + w - 18 ; } else if (!isTiret) { nExon++ ; tframe = array(look->tframe, nExon, int) ; diffTframe = (tframe - oldTframe)%3 ; if (diffTframe < 0) diffTframe +=3 ; i += diffTframe ; cq1 += diffTframe ; oldTframe = tframe ; isTiret = TRUE ; } cq1 += 3 ; i += 3 ; if ( !(nbTrans % 50) ) *cq++ = '\n' ; } *cq++ = 0 ; graphPostBuffer(buf) ; messfree(buf) ; } else /* PEPTIDE */ { i = imax ; while(i >= i0 ) { *cq++ = doComplementSurPlace ? antiCodon (cq1) : codon(cq1) ; cq1 += 3 ; i -= 3 ; if ( !((imax - i) % 50) ) *cq++ = '\n' ; } *cq++ = 0 ; graphPostBuffer(buf) ; messfree(buf) ; weight = 18 ; cq1 = arrp(look->dna, i0, char) ; i = imax ; while(i >= i0 ) { w = doComplementSurPlace ? molecularWeight[(int) antiCodon (cq1)] : molecularWeight[(int) codon(cq1)] ; if (w < 0) break ; else if (w == 0) { weight = 0 ; break ; } else weight = weight + w - 18 ; cq1 += 3 ; i -= 3 ; } } }/* endif i0 */ } /* endif seg->type */ else { cp = arrp(look->colors, i, unsigned char) ; for ( ; i <= max ; ++i) *cp++ |= TINT_HIGHLIGHT1; } /* write blue report line */ if (!(look->flag & FLAG_HIDE_HEADER)) { fMapReportLine (look, seg, FALSE, x) ; graphBoxDraw (look->segBox, BLACK, PALEBLUE) ; /* default is to put dna in cut buffer, alternative is coord data. */ if (look->cut_data == FMAPCUT_COORDS) graphPostBuffer(look->segTextBuf) ; } for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; ++i) if (i != look->map->thumb.box) { seg2 = BOXSEG(i) ; if ((seg2->type == DNA_SEQ || seg2->type == PEP_SEQ || seg2->type == TRANS_SEQ || seg2->type == TRANS_SEQ_UP) && seg2->x1 <= finDna && seg2->x2 >= seg->x1) graphBoxDraw (i, -1, -1) ; } /* colour siblings */ if ((parent = seg->parent)) for (i = look->minLiveBox ; i < arrayMax(look->boxIndex) ; i++) { if (i == look->map->thumb.box) continue; seg2 = BOXSEG(i) ; if (seg2->parent == parent) switch (seg2->type) { case DNA_SEQ: case PEP_SEQ: case TRANS_SEQ: case TRANS_SEQ_UP: break ; case HOMOL: case HOMOL_UP: graphBoxDraw (i, -1, PALERED) ; break ; default: if (assFind (look->chosen, SEG_HASH(seg2), 0)) graphBoxDraw (i, -1, CYAN) ; else if (assFind (look->antiChosen, SEG_HASH(seg2), 0)) graphBoxDraw (i, -1, DARKRED) ; else graphBoxDraw (i, -1, PALEBLUE) ; break ; } } switch (seg->type) /* colour box itself over sibling colour */ { case DNA_SEQ: case PEP_SEQ: case TRANS_SEQ: case TRANS_SEQ_UP: break ; case VIRTUAL_SUB_SEQUENCE_TAG: case VIRTUAL_SUB_SEQUENCE_TAG_UP: graphBoxDraw (box, -1, PALERED) ; default: if (assFind (look->chosen, SEG_HASH(seg), 0) || assFind (look->antiChosen, SEG_HASH(seg), 0)) graphBoxDraw (box, -1, MAGENTA) ; else graphBoxDraw (box, -1, PALERED) ; break ; } } /**************************************/ static void fMapFollow (FeatureMap look, double x, double y) { SEG* seg = BOXSEG(look->activeBox) ; KEY key = seg->key ; int table = class(key) ; if (!table || table == _VText) { key = seg->parent ; table = class(key) ; } if ((seg->type | 0x1) == SPLICED_cDNA_UP || (seg->type | 0x1) == TRANSCRIPT_UP ) { if (fMapcDNAFollow (look->activeBox)) return ; } else if ((seg->type | 0x1) == TRANSCRIPT_UP) { display (seg->parent, look->seqKey, "TREE") ; return ; } else if (seg->type == VIRTUAL_SUB_SEQUENCE_TAG || seg->type == VIRTUAL_SUB_SEQUENCE_TAG_UP || seg->type == VIRTUAL_PARENT_SEQUENCE_TAG || seg->type == VIRTUAL_CONTIG_TAG) { if (fMapFollowVirtualMultiTrace (look->activeBox)) return ; } if (table == _VSequence) display (key, look->seqKey, "TREE") ; else if (table) display (key, look->seqKey, 0) ; } /*************************************************************/ /************* a couple of general utilities *****************/ BOOL fMapFindSegBounds (FeatureMap look, SegType type, int *min, int *max) { for (*min = 1 ; *min < arrayMax(look->segs) ; ++*min) if (arrp(look->segs, *min, SEG)->type == type) break ; for (*max = *min ; *max < arrayMax(look->segs) ; ++*max) if (arrp(look->segs, *max, SEG)->type != type) break ; return (*min < arrayMax(look->segs)) ; } float fMapSquash (float value, float midpoint) { if (value < 0) value = 0 ; value /= midpoint ; return (1 - 1/(value + 1)) ; } /*************************************************************/ /**************** Reverse complement code ********************/ /* "...where angels fear to tread..." */ /* */ /* */ static void fMapCompHelp (void) { graphMessage ( "Reverse-Complement: does what you expect (I hope!)\n" "Complement: keep coordinates, view opposite strand,\n" " which runs bottom to top 5' to 3'.\n" "Reverse: reverse coordinates, view the same strand,\n" " which now runs bottom to top 5' to 3'.\n" " In both the last two cases the DNA display is still\n" " read 5' to 3' left to right. Doing both Complement\n" " and Reverse is equivalent to Reverse-Complement.\n" "Complement DNA in place: Don't change coords or\n" " direction. DNA display swaps G-C, A-T. Now\n" " the sequence reads 3' to 5' left to right.") ; } static void fMapMenuComplementSurPlace (void) { FeatureMap look = currentFeatureMap("fMapMenuComplementSurPlace") ; look->flag ^= FLAG_COMPLEMENT_SUR_PLACE ; fMapDraw (look, 0) ; } static void fMapReverse (void) { FeatureMap look = currentFeatureMap("fMapReverse") ; look->map->mag = -look->map->mag ; look->flag ^= FLAG_REVERSE ; look->origin = look->length + 1 - look->origin ; fMapDraw (look, 0) ; } /* This is the guts of the reverse complement, you should note that the code */ /* only knows about the seg coords and only operates on those. If you have */ /* other features that have coordinates stored else where you need to add */ /* code to translate them here (see CDS stuff below). */ /* */ void fMapRC (FeatureMap look) { int top = look->length - 1 ; int i, tmp ; SEG *seg ; char *ci, *cj, ctmp ; float *ftmp ; SEQINFO *sinf ; sMapFeaturesComplement(look->length, look->segs, look->seqInfo) ; tmp = arrayMax(look->segs) ; arrayMax(look->segs) = look->lastTrueSeg ; /* avoid mixing temporary and permanent segs */ arraySort (look->segs, fMapOrder) ; arrayMax(look->segs) = tmp ; if (look->dnaR) { Array tmp = look->dna ; look->dna = look->dnaR ; look->dnaR = tmp ; } else { ci = arrp(look->dna, 0, char) ; cj = arrp(look->dna, top, char) ; while (ci <= cj) { ctmp = *ci ; *ci++ = complementBase[(int)*cj] ; *cj-- = complementBase[(int)ctmp] ; } } if (look->gf.cum) { tmp = look->gf.max ; look->gf.max = top + 1 - look->gf.min ; look->gf.min = top + 1 - tmp ; ftmp = look->gf.cum ; look->gf.cum = look->gf.revCum ; /* already reversed */ look->gf.revCum = ftmp ; } look->map->centre = top - look->map->centre ; tmp = look->zoneMin ; look->zoneMin = look->length - look->zoneMax ; look->zoneMax = look->length - tmp ; look->origin = look->length + 1 - look->origin ; look->flag ^= FLAG_COMPLEMENT ; fMapClearDNA (look) ; /* The sMap must be the same way round and the displayed sequence, so we re-create it here. It would be much easier to re-calculate from scratch to reverse-complement. */ sMapDestroy(look->smap); if (look->flag & FLAG_COMPLEMENT) look->smap = sMapCreate(look->handle, look->seqKey, look->stop-1, look->start-1, NULL); else look->smap = sMapCreate(look->handle, look->seqKey, look->start-1, look->stop-1, NULL); } static void fMapCompToggleTranslation (FeatureMap look) { if ((mapColSetByName (look->map, "Up Gene Translation", -1)) || (mapColSetByName (look->map, "Down Gene Translation", -1))) { mapColSetByName (look->map, "Up Gene Translation", 2) ; mapColSetByName (look->map, "Down Gene Translation", 2) ; } } static void fMapMenuComplement (void) { FeatureMap look = currentFeatureMap("fMapMenuComplement"); fMapCompToggleTranslation (look) ; fMapRC (look) ; fMapReverse () ; /* does a lookDraw() */ } static void fMapMenuReverseComplement (void) { FeatureMap look = currentFeatureMap("fMapMenuReverseComplement") ; fMapCompToggleTranslation (look) ; fMapRC (look) ; fMapDraw (look, 0) ; } void fMapDoReverseComplement (FeatureMap look, BOOL complement) { Graph old = graphActive() ; BOOL iscpl = look->flag & FLAG_COMPLEMENT ; if (complement != iscpl) { fMapRC (look) ; fMapDraw (look, 0) ; } graphActivate(old) ; } /**************************************************************/ /********* major routine to dump out in GFF format ************/ static char *cleanNewlines (char *s, STORE_HANDLE handle) { int n = 0 ; char *cp, *cq, *copy ; for (cp = s ; *cp ; ++cp) if (*cp == '\n') ++n ; if (!n) return s ; copy = halloc (cp-s+n+1, handle) ; for (cp = s, cq = copy ; *cp ; ++cp) if (*cp == '\n') { *cq++ = '\\' ; *cq++ = 'n' ; } else *cq++ = *cp ; *cq = 0 ; return copy ; } /***********************************/ static char strandOfKey(FeatureMap look, KEY key, int rootReversed) { /* Determine direction of key in Dump. If the dump is of a sequence which is reversed wrt to the root of the fMap, rootReversed == 1) */ SMapKeyInfo *info = sMapKeyInfo(look->smap, key); BOOL isReverse; if (!info) return '.'; isReverse = sMapIsReverse(info, 1, 0); if (isReverse) return rootReversed == 0 ? '-' : '+'; return rootReversed == 1 ? '-' : '+'; } static BOOL fMapDumpGFF (FeatureMap look, int version, KEY *refSeq, int *offPtr, BOOL isList, DICT *sourceSet, DICT *featSet, ACEOUT gff_out) { int i ; SEG *seg ; HOMOLINFO *hinf ; SEQINFO *sinf ; METHOD *meth ; KEY seqKey, sourceKey ; int x, y, type, offset = 0 ; int chopped ; /* only used for version 1 */ int tmp = 0, reversed = 0 ; /* used by AcePerl */ char *featName, *sourceName, *tagText = 0 ; float score ; char strand, frame ; BOOL flipped ; BOOL isScore ; STORE_HANDLE handle = handleCreate() ; Associator key2sinf = assHandleCreate (handle) ; Associator key2source = assHandleCreate (handle) ; Associator key2feat = assHandleCreate (handle) ; Array stats = arrayHandleCreate (64, int, handle) ; DICT *listSet = 0 ; /* first establish seqKey and offset */ seqKey = 0 ; if (refSeq && !*refSeq && version > 1) *refSeq = look->seqOrig ; if (refSeq && *refSeq) { for (i = 1 ; i < arrayMax(look->segs) ; ++i) { seg = arrp(look->segs, i, SEG) ; if ((seg->type == SEQUENCE || seg->type == SEQUENCE_UP) && seg->key == *refSeq) break ; } if (i < arrayMax(look->segs)) #if 0 { if (seg->type == SEQUENCE) { seqKey = *refSeq ; offset = - seg->x1 ; } else if (seg->type == SEQUENCE_UP) { messout ("Can't GFF dump from a reversed refseq for now. Sorry.") ; return FALSE ; } } #endif { if (seg->type == SEQUENCE || seg->type == SEQUENCE_UP) { seqKey = *refSeq ; offset = - seg->x1 ; reversed = seg->type == SEQUENCE_UP ? 1 : 0; } } else { messout ("Can't find reference sequence %s", name (*refSeq)) ; return FALSE ; } } if (!seqKey) /* find minimal spanning sequence */ { x = look->zoneMin ; y = look->zoneMax ; seqKey = 0 ; if (!fMapFindSpan(look, &seqKey, &x, &y)) seqKey = look->seqKey ; if (x > y) /* what if reversed? */ { messout ("Can't GFF dump from reversed sequences for now. Sorry.") ; return FALSE ; } offset = x - look->zoneMin ; } if (refSeq) *refSeq = seqKey ; if (offPtr) *offPtr = offset ; if (!isList) { aceOutPrint (gff_out, "##gff-version %d\n", version) ; aceOutPrint (gff_out, "##date %s\n", timeShow (timeParse ("today"))) ; aceOutPrint (gff_out, "##sequence-region %s %d %d %s\n", name(seqKey), look->zoneMin+1+offset, look->zoneMax+offset, reversed ? "(reversed)" : "") ; } if (isList) listSet = dictHandleCreate (64, handle) ; for (i = 0 ; i < arrayMax(look->segs) ; ++i) { seg = arrp(look->segs, i, SEG) ; if (seg->x1 > look->zoneMax || seg->x2 < look->zoneMin) continue ; if (seg->type == MASTER || seg->type == VISIBLE || seg->type == CDS || seg->type == CDS_UP || seg->type == DNA_SEQ || seg->type == PEP_SEQ || seg->type == ORF || seg->type == TRANS_SEQ || seg->type == TRANS_SEQ_UP) continue ; /* Set the strand of the feature. In the seg this is usually indicated by * the "_UP", but some segs have been flipped from one strand to the other * for display (usually for EST 3' reads), we flip them back to their original * orientation here. Otherwise people don't know where they are for the GFF. */ if (seg->type & 0x01 && (seg->type >= SEQUENCE && seg->type <= ALLELE_UP)) { type = seg->type - 1 ; x = seg->x2+1 ; y = seg->x1+1 ; } else { type = seg->type ; x = seg->x1+1 ; y = seg->x2+1 ; } if (seg->flags & SEGFLAG_FLIPPED) { int tmp = x ; x = y ; y = tmp ; } chopped = 0 ; if (x < y) { flipped = FALSE ; strand = '+' ; /* Note that this might be overridden. */ if (version == 1 && x <= look->zoneMin) chopped = look->zoneMin + 1 - x ; } else if (x > y) { int tmp = x ; x = y ; y = tmp ; flipped = TRUE ; strand = '-' ; if (version == 1 && y > look->zoneMax) chopped = y - look->zoneMax ; } else /* x = y */ { strand = strandOfKey(look, seg->source, reversed); if (strand == '.') strand = strandOfKey(look, seg->parent, reversed); } if (version == 1) /* clip */ { if (x <= look->zoneMin) x = look->zoneMin + 1 ; if (y > look->zoneMax) y = look->zoneMax ; } x += offset ; y += offset ; frame = '.' ; /* other defaults */ isScore = (version == 1) ? TRUE : FALSE ; score = 0.0 ; sourceName = (version == 1 || version == 2) /* mieg, why not in version2 ? */ ? fMapSegTypeName[type] : 0 ; featName = 0 ; sourceKey = 0 ; sourceName = ""; /* LS otherwise everything ends up being a SEQUENCE */ switch (type) { case EXON: case INTRON: case SEQUENCE: sinf = arrp(look->seqInfo, seg->data.i, SEQINFO) ; if ( seg->key && (version == 2) ) { featName = className(seg->key); } else featName = (type == SEQUENCE) ? "sequence" : (type == EXON) ? "exon" : "intron" ; if (sinf->method) { sourceKey = sinf->method ; } if (type == SEQUENCE) { assInsert (key2sinf, assVoid(seg->key), sinf) ; if (sinf->flags & SEQ_SCORE) { score = sinf->score ; isScore = TRUE ; } } break ; case CODING: if (version == 1) featName = "coding_exon" ; else featName = "CDS" ; assFind (key2sinf, assVoid(seg->parent), &sinf) ; if (sinf->method) sourceKey = sinf->method ; switch (seg->data.i%3) { case 0: frame = '0' ; break ; case 1: frame = '2' ; break ; case 2: frame = '1' ; break ; } break ; case HOMOL: hinf = arrp(look->homolInfo, seg->data.i, HOMOLINFO) ; score = hinf->score ; sourceKey = hinf->method ; featName = "similarity" ; meth = methodCacheGet(look->mcache, hinf->method, "") ; if (!(meth->flags & METHOD_STRAND_SENSITIVE)) strand = '.' ; if (meth->flags & METHOD_SCORE) isScore = TRUE ; if (meth->flags & METHOD_FRAME_SENSITIVE) frame = '0' ; break ; case FEATURE: case ATG: case SPLICE5: case SPLICE3: sourceKey = seg->key ; if (type == SPLICE3) featName = "splice3" ; else if (type == SPLICE5) featName = "splice5" ; else if (type == ATG) featName = "atg" ; else if (version == 1) featName = "feature" ; else featName = "misc_feature" ; meth = methodCacheGet(look->mcache, seg->key, "") ; if (meth->flags & METHOD_SCORE) { score = seg->data.f ; isScore = TRUE ; } if (!(meth->flags & METHOD_STRAND_SENSITIVE)) strand = '.' ; if (meth->flags & METHOD_FRAME_SENSITIVE) frame = '0' ; break ; case ASSEMBLY_TAG: tagText = strnew (seg->data.s, handle) ; sourceName = "assembly_tag" ; featName = tagText ; { char *cp ; for (cp = tagText ; *cp && *cp != ':' && *cp != ' ' && *cp != '\t' && *cp != '\n' ; ++cp) ; if (*cp) { *cp++ = 0 ; while (*cp == ' ' || *cp == '\t' || *cp == '\n') ++cp ; } tagText = cp ; } break ; case ALLELE: if (version == 1) sourceName = "variation" ; if (seg->data.s) { if (*seg->data.s == '-') featName = "deletion" ; else { char *cp = seg->data.s ; featName = "variation" ; while (*cp) if (!strchr ("AGCTagct", *cp++)) featName = (version == 1) ? "insertion_site" : "insertion" ; } } break ; case EMBL_FEATURE: featName = name(seg->key) ; break ; case CLONE_END: featName = name(seg->data.k) ; strand = '.' ; break ; default: ; } /* finalise source and feat fields */ if (sourceKey) { if (assFind (key2source, assVoid(sourceKey), &sourceName)) assFind (key2feat, assVoid(sourceKey), &featName) ; else { OBJ obj = bsCreate (sourceKey) ; sourceName = name(sourceKey) ; if (obj) { bsGetData (obj, str2tag("GFF_source"), _Text, &sourceName) ; if (bsGetData (obj, str2tag("GFF_feature"), _Text, &featName)) { featName = strnew (featName, handle) ; assInsert (key2feat, assVoid(sourceKey), featName) ; } } sourceName = strnew(sourceName, handle); assInsert (key2source, assVoid(sourceKey), sourceName) ; if (obj) /* only after strnew(sourceName) !!! */ bsDestroy(obj); } } if (!sourceName) /* only possible if version > 1 */ sourceName = "unknown" ; if (!featName) /* from method or from seg->type above */ { if (version == 1) featName = sourceName ; else featName = fMapSegTypeName[type] ; } if (frame >= '0' && frame <= '2' && chopped) { frame = frame + 3 - (chopped % 3) ; if (frame == '3') frame = '0' ; if (frame == '4') frame = '1' ; if (frame == '5') frame = '2' ; } if (sourceSet && featSet && (dictMax (sourceSet) || dictMax (featSet)) && !dictFind (sourceSet, sourceName, 0) && !dictFind (featSet, featName, 0)) continue ; if (isList) /* just accumulate stats */ { int k ; dictAdd (listSet, messprintf("%s\t%s",sourceName,featName), &k) ; ++array(stats, k, int) ; continue ; /* move on to next line */ } /* !isList from here on */ /* write the main part of the line */ /* LS/AcePerl: fixup reversed reference sequence */ if (reversed) { tmp = look->zoneMax + offset + 1 - y; y = look->zoneMax + offset + 1 - x; x = tmp; if ( strand == '+' ) strand = '-'; else if ( strand == '-' ) strand = '+'; } if (isScore) aceOutPrint (gff_out, "%s\t%s\t%s\t%d\t%d\t%g\t%c\t%c", name(seqKey), sourceName, featName, x, y, score, strand, frame) ; else aceOutPrint (gff_out, "%s\t%s\t%s\t%d\t%d\t.\t%c\t%c", name(seqKey), sourceName, featName, x, y, strand, frame) ; /* extras */ switch (type) { case SEQUENCE: case EXON: case CODING: if (seg->parent) { if (version == 1) aceOutPrint (gff_out, "\t%s", name(seg->parent)) ; else aceOutPrint (gff_out, "\t%s \"%s\"", className(seg->parent), name(seg->parent)) ; } break ; case INTRON: sinf = arrp(look->seqInfo, seg->data.i, SEQINFO) ; if (version == 1) { if (seg->parent) { /* RD 010129 added class for version 2 like EXON */ if (version == 1) aceOutPrint (gff_out, "\t%s", name(seg->parent)) ; else aceOutPrint (gff_out, "\t%s \"%s\"", className(seg->parent), name(seg->parent)) ; if (sinf->flags & SEQ_CONFIRMED) aceOutPrint (gff_out, " Confirmed") ; } else if (sinf->flags & SEQ_CONFIRMED) aceOutPrint (gff_out, "\tConfirmed") ; if (sinf->flags & SEQ_CONFIRM_EST) aceOutPrint (gff_out, " by_EST") ; if (sinf->flags & SEQ_CONFIRM_CDNA) aceOutPrint (gff_out, " by_cDNA") ; if (sinf->flags & SEQ_CONFIRM_HOMOL) aceOutPrint (gff_out, " by_homology") ; if (sinf->flags & SEQ_CONFIRM_UTR) aceOutPrint (gff_out, " in_UTR") ; if (sinf->flags & SEQ_CONFIRM_FALSE) aceOutPrint (gff_out, " as_false") ; } else { BOOL isData = FALSE ; if (seg->parent) { aceOutPrint (gff_out, "\tSequence \"%s\"", name(seg->parent)) ; isData = TRUE ; } if (sinf->flags & SEQ_CONFIRM_EST) { aceOutPrint (gff_out, "%s", isData ? " ; " : "\t") ; isData = TRUE ; aceOutPrint (gff_out, "Confirmed_by_EST") ; } if (sinf->flags & SEQ_CONFIRM_CDNA) { aceOutPrint (gff_out, "%s", isData ? " ; " : "\t") ; isData = TRUE ; aceOutPrint (gff_out, "Confirmed_by_cDNA") ; } if (sinf->flags & SEQ_CONFIRM_HOMOL) { aceOutPrint (gff_out, "%s", isData ? " ; " : "\t") ; isData = TRUE ; aceOutPrint (gff_out, "Confirmed_by_homology") ; } if (sinf->flags & SEQ_CONFIRM_UTR) { aceOutPrint (gff_out, "%s", isData ? " ; " : "\t") ; isData = TRUE ; aceOutPrint (gff_out, "Confirmed_in_UTR") ; } if (sinf->flags & SEQ_CONFIRM_FALSE) { aceOutPrint (gff_out, "%s", isData ? " ; " : "\t") ; isData = TRUE ; aceOutPrint (gff_out, "Confirmed_as_FALSE") ; } } break ; case HOMOL: if (version == 1) aceOutPrint (gff_out, "\t%s:%s", className(seg->key), name(seg->key)) ; else aceOutPrint (gff_out, "\tTarget \"%s:%s\"", className(seg->key), name(seg->key)) ; hinf = arrp(look->homolInfo, seg->data.i, HOMOLINFO) ; if (flipped) aceOutPrint (gff_out, " %d %d", hinf->x2, hinf->x1) ; else aceOutPrint (gff_out, " %d %d", hinf->x1, hinf->x2) ; break ; case FEATURE: case SPLICE5: case SPLICE3: case ATG: { BOOL isTab = FALSE ; if (seg->parent) { if (version == 1) aceOutPrint (gff_out, "\t%s", cleanNewlines(dictName (look->featDict, -seg->parent), handle)) ; else aceOutPrint (gff_out, "\tNote \"%s\"", cleanNewlines(dictName (look->featDict, -seg->parent), handle)) ; isTab = TRUE ; } if (assFind (look->chosen, SEG_HASH(seg), 0)) aceOutPrint (gff_out, "%sSelected", isTab ? " ; " : "\t") ; else if (assFind (look->antiChosen, SEG_HASH(seg), 0)) aceOutPrint (gff_out, "%sAntiselected", isTab ? " ; " : "\t") ; } break ; case ASSEMBLY_TAG: if (tagText && *tagText) { if (version == 1) aceOutPrint (gff_out, "\t%s", cleanNewlines(tagText, handle)) ; else aceOutPrint (gff_out, "\tNote \"%s\"", cleanNewlines(tagText, handle)) ; } break ; case ALLELE: if (version == 1) aceOutPrint (gff_out, "\t%s", name(seg->key)) ; else aceOutPrint (gff_out, "\tAllele \"%s\"", name(seg->key)) ; if (seg->data.s && *seg->data.s != '-') { if (version == 1) aceOutPrint (gff_out, "\t%s", cleanNewlines(seg->data.s, handle)) ; else if (strcmp (featName, "variation") == 0) aceOutPrint (gff_out, " ; Variant \"%s\"", cleanNewlines(seg->data.s, handle)) ; else aceOutPrint (gff_out, " ; Insert \"%s\"", cleanNewlines(seg->data.s, handle)) ; } break ; case EMBL_FEATURE: if (seg->data.s) { if (version == 1) aceOutPrint (gff_out, "\t%s", cleanNewlines(seg->data.s, handle)) ; else aceOutPrint (gff_out, "\tNote \"%s\"", cleanNewlines(seg->data.s, handle)) ; } break ; case CLONE_END: if (version == 1) aceOutPrint (gff_out, "\t%s", name(seg->key)) ; else aceOutPrint (gff_out, "\tClone \"%s\"", name(seg->key)) ; break ; default: ; } aceOutPrint (gff_out, "\n") ; } if (isList) for (i = 0 ; i < arrayMax(stats) ; ++i) aceOutPrint (gff_out, "%s\t%d\n", dictName (listSet, i), arr(stats,i,int)) ; handleDestroy (handle) ; return TRUE ; } /* fMapDumpGFF */ /********* entry point to fMapDumpGFF for fmap menu *********/ static void fMapMenuDumpSegs (void) { static char directory[DIR_BUFFER_SIZE] = ""; static char filename[FIL_BUFFER_SIZE] = ""; ACEOUT gff_out; FeatureMap look = currentFeatureMap("fMapMenuDumpSegs") ; gff_out = aceOutCreateToChooser ("File to write features into", directory, filename, "gff", "w", 0); if (gff_out) { fMapDumpGFF (look, 2, 0, 0, 0, 0, 0, gff_out); aceOutDestroy (gff_out); } return; } /* fMapMenuDumpSegs */ /********* entry point to fMapDumpGFF from dnacpt.c *********/ void fMapDumpSegsKeySet (KEYSET kSet, ACEOUT gff_out) { int i ; #ifdef USE_SMAP KEY key; #endif FeatureMap look ; fMapInitialise() ; look = (FeatureMap) messalloc (sizeof(struct FeatureMapStruct)) ; look->origin = 0 ; look->zoneMin = 0 ; for (i = 0 ; i < keySetMax(kSet) ; ++i) { #ifdef USE_SMAP key = keySet(kSet, i) ; if (!sMapTreeRoot (key, 1, 0, &look->seqKey, &look->start, &look->stop)) { if (look->stop = sMapLength (key)) { look->seqKey = key ; look->start = 1 ; } else continue ; /* can't get limits */ } #else look->seqKey = keySet(kSet, i) ; look->start = 1 ; look->stop = 0 ; /* find whole sequence */ if (!findStartStop (&look->seqKey, &look->start, &look->stop)) { Array dna ; if ((dna = dnaGet (look->seqKey))) /* try DNA */ { look->start = 1 ; look->stop = arrayMax(dna) ; arrayDestroy (dna) ; } else continue ; /* can't get limits */ } #endif --look->start ; --look->stop ; look->fullLength = look->length = look->stop - look->start + 1 ; look->zoneMax = look->length - 1 ; if (fMapConvert(look)) fMapDumpGFF(look, 2, 0, 0, 0, 0, 0, gff_out) ; arrayDestroy (look->segs) ; /* THIS LOOKS LIKE A MEMORY HOLE TO ME, WE NEED BETTER CLEANUP... we should at least do a messfree (look->handle) ; and probably much else besides, there needs to be a proper routine, separate from fmapdestroy() to do this clear up. */ } messfree (look) ; return; } /* fMapDumpSegsKeySet */ /****************/ static void fMapDumpAlign (FeatureMap look, BOOL isPep, ACEOUT dump_out) { int i, ii, bufMax ; char *cp, *cq ; Array dna = 0 ; SEG *seg ; HOMOLINFO *hinf = 0 ; KEY seqKey = look->seqOrig ; int a1, a2, x1, x2, x, y ; char *buf = 0 ; BOOL first = TRUE ; STORE_HANDLE handle = handleCreate() ; x = look->zoneMin ; y = look->zoneMax ; bufMax = y - x + 1 ; if (isPep) bufMax /= 3 ; messfree(buf) ; buf = messalloc (bufMax + 1) ; buf[bufMax] = 0 ; /* extras */ for (ii = 0 ; ii < arrayMax(look->segs) ; ++ii) { seg = arrp(look->segs, ii, SEG) ; if (seg->x1 > look->zoneMax || seg->x2 < look->zoneMin) continue ; switch (seg->type) { case HOMOL: if (first) { first = FALSE ; dna = arrayCopy (look->dna) ; x = look->zoneMin ; y = look->zoneMax ; memset (buf, (int)'.', bufMax) ; aceOutPrint (dump_out, "%s/%d-%d ", name(seqKey), isPep ? (x + 1 - look->origin)/3 : x + 1 - look->origin , isPep ? (y - look->origin)/3 : y - look->origin) ; if (x > 0 && x < y && y < arrayMax(dna)) { for (i = x , cp = arrp(dna, i, char), cq = buf ; i < y ; i++, cp++, cq++) if (isPep) { *cq = codon(cp) ; i +=2 ; cp += 2 ; } else *cq = dnaDecodeChar[(int)*cp] ; } aceOutPrint (dump_out, "%s %s\n", buf, name(seqKey)) ; arrayDestroy (dna) ; } dna = dnaGet (seg->key) ; if (dna && arrayMax(dna)) { hinf = arrp(look->homolInfo, seg->data.i, HOMOLINFO) ; if (hinf->x1 > hinf->x2) /* flipped) */ { } else { x1 = hinf->x1 ; x2 = hinf->x2 ; a1 = seg->x1 ; a2 = seg->x2 ; i = look->zoneMin - a1 ; if (i > 0) { a1 += i ; x1 += i ; } i = seg->x2 - look->zoneMax ; if (i > 0) { a2 -= i ; x2 -= i ; } if (x1 >= 0 && x1 < x2 && x2 < arrayMax(dna)) { aceOutPrint (dump_out, "%s/", name(seg->key)) ; aceOutPrint (dump_out, "%d-%d ", isPep ? x1/3 : x1, isPep ? x2/3 : x2 ) ; memset (buf, (int)'.', bufMax) ; for (i = a1 , cp = arrp(dna, x1 - 1, char) , cq = (isPep ? buf + (a1 - look->zoneMin)/3 : buf + a1 - look->zoneMin) ; i <= a2 ; i++, cp++, cq++) if (isPep) { *cq = codon(cp) ; i += 2 ; cp += 2 ; } else *cq = dnaDecodeChar[(int)*cp] ; aceOutPrint (dump_out, "%s %s\n",buf, name(seg->key)) ; } } arrayDestroy (dna) ; } break ; default: ; } } messfree (buf) ; handleDestroy (handle) ; return; } /* fMapDumpAlign */ /*************************************/ static void geneStats (void) { #if !defined(MACINTOSH) int i ; SEG *seg ; int len = 0 ; int nIntron = 0 ; int nExon = 0 ; int nCDS = 0 ; int intronLength = 0 ; int exonLength = 0 ; int cdsLength = 0 ; int nLoci = 0 ; int ncDNA = 0 ; int nMatch = 0 ; char *n ; FeatureMap look = currentFeatureMap("geneStats") ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) { seg = arrp(look->segs, i, SEG) ; n = name (seg->parent) ; if (n[strlen(n)-1] == 'a') continue ; switch (seg->type) { case INTRON: case INTRON_UP: ++nIntron ; intronLength += (seg->x2 - seg->x1 + 1) ; break ; case EXON: case EXON_UP: ++nExon ; exonLength += (seg->x2 - seg->x1 + 1) ; break ; case CDS: case CDS_UP: ++nCDS ; cdsLength += (seg->x2 - seg->x1 + 1) ; printf (" %s\n", name(seg->parent)) ; break ; case VISIBLE: printf (" %s: %s %s\n", name(seg->parent), name(seg->data.k), name(seg->key)) ; if (seg->data.k == _Locus) ++nLoci ; if (seg->data.k == _Matching_cDNA) ++ncDNA ; if (seg->data.k == _Brief_identification) ++nMatch ; break ; default: /* needed for picky compiler */ break ; } } for (i = 0, n = arrp(look->dna,0,char) ; i < arrayMax(look->dna) ; ++i) if (*n++) ++len ; printf ("%d bp sequenced DNA from %d in LINK\n", len, look->length) ; printf ("%d introns, total %d bp\n", nIntron, intronLength) ; printf ("%d exons, total %d bp\n", nExon, exonLength) ; printf ("%d CDS (predicted genes), total %d bp\n", nCDS, cdsLength) ; printf ("%d database matches\n", nMatch) ; printf ("%d cDNAs\n", ncDNA) ; printf ("%d Genetic loci\n", nLoci) ; #else messout ("Sorry, this does nothing on a macintosh " "(you aren't missing much - it was a quick hack)") ; #endif } #ifndef ACEMBLY /* trace, assembly stuff */ void fMapTraceDestroy (FeatureMap look) { ;} void fMapShowTraceAssembly (LOOK look, float *offset) { ; } BOOL fMapFollowVirtualMultiTrace (int box) { return FALSE ; } void abiFixFinish (void) { ; } void fMapGelDisplay (void) { gelDisplay () ; } void fMapTraceFindMultiplets (FeatureMap look) { ; } void fMapSelectVirtualMultiplet (FeatureMap look, SEG *seg) { ; } #endif /****************** giface hooks *****************/ void fMapGifsMap(ACEIN command_in, ACEOUT result_out) { char *cp, *cp1; BOOL gotCoords = FALSE, dump = FALSE; int x1, x2, y1, y2, nx1, nx2; STORE_HANDLE handle = handleCreate(); KEY to=0, from=0; SMap* smap = NULL; SMapKeyInfo *info; SMapStatus status; if (!command_in || !result_out) messcrash("Bad args to fMapGifsMap") ; while ((cp = strnew(aceInWord(command_in), handle))) { if (strcmp(cp, "-dump") == 0) dump = TRUE; else if (strcmp(cp, "-coords") == 0) { if (!aceInInt (command_in, &x1) || !aceInInt (command_in, &x2) || (x2 == x1)) goto usage ; gotCoords = TRUE; } else if (strcmp(cp, "-to") == 0) { if (!(cp1 = strnew(aceInWord(command_in), handle))) goto usage ; if (!lexClassKey(cp1, &to)) { aceOutPrint (result_out, "// gif smap error: invalid to object\n") ; goto abort ; } } else if (strcmp(cp, "-from") == 0) { if (!(cp1 = strnew(aceInWord(command_in), handle))) goto usage ; if (!lexClassKey(cp1, &from)) { aceOutPrint (result_out, "// gif smap error: invalid from object\n") ; goto abort ; } } else goto usage ; } if (from == 0) { aceOutPrint (result_out, "// gif smap error: must have -from\n") ; goto abort ; } if (dump) { if (!(smap = sMapCreate(handle, from, 1, sMapLength(from), NULL))) { aceOutPrint (result_out, "// gif smap error: Cannot create smap\n") ; goto abort ; } sMapDump(smap, result_out); goto abort; } if (!gotCoords) { x1 = 1; x2 = sMapLength(from); } if (to == 0) sMapTreeRoot(from, x1, x2, &to, NULL, NULL); if (!(smap = sMapCreate(handle, to, 1, sMapLength(to), NULL))) { aceOutPrint (result_out, "// gif smap error: Cannot create smap\n") ; goto abort ; } if (!(info = sMapKeyInfo(smap, from))) { aceOutPrint (result_out, "// gif smap error: -from object not in sMap\n") ; goto abort ; } status = sMapMap(info, x1, x2, &y1, &y2, &nx1, &nx2); aceOutPrint(result_out, "SMAP %s", className(to)); aceOutPrint(result_out, ":%s", name(to)); aceOutPrint(result_out, " %d %d %d %d ", y1, y2, nx1, nx2); if (status == SMAP_STATUS_PERFECT_MAP) aceOutPrint(result_out, "%s", sMapErrorString(status)); else { int bit = 1; char *sep=""; while (bit) { char *string; if ((status & bit) && (string = sMapErrorString(status&bit))) { aceOutPrint(result_out, "%s%s", sep, string); sep = ","; } bit = bit<<1; } } aceOutPrint(result_out, "\n"); handleDestroy(handle); return; usage: aceOutPrint (result_out, "// gif smap error: usage: " "smap -coords x1 x2 -from class:object -to class:object\n") ; abort: handleDestroy(handle); } /* Construct a sequence but don't display it. */ FeatureMap fMapGifGet (ACEIN command_in, ACEOUT result_out, FeatureMap look_in) { KEY key ; int x1, x2 ; char *cp, *word = NULL ; FeatureMap look = NULL ; if (!command_in || !result_out) messcrash("Bad args to fMapGifGet") ; fMapInitialise() ; if (look_in) { fMapDestroy (look_in) ; /* could check for possible reuse */ } x1 = 1 ; x2 = 0 ; /* default for whole sequence */ while ((cp = aceInWord(command_in))) { if (strcmp(cp, "-coords") == 0) { if (!aceInInt (command_in, &x1) || !aceInInt (command_in, &x2) || (x2 == x1)) goto usage ; } else word = strnew(cp, 0) ; } if (!word) /* try to get from graph */ { if (fMapActive(0, 0, 0, &look)) return look ; else { aceOutPrint (result_out, "// gif seqget error: no sequence attached to active window\n") ; return NULL ; } } if (!lexword2key (word, &key, _VSequence)) { aceOutPrint (result_out, "// gif seqget error: Sequence %s not known\n", word) ; goto usage ; } if (x1 < x2 || (x1 == 1 && x2 == 0)) look = fMapCreateLook (key, x1, x2, FALSE, NULL) ; else look = fMapCreateLook (key, x2, x1, FALSE, NULL) ; /* set the zone */ if (look) { if (x1 == 1 && x2 == 0) { int i ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) { if (arrp(look->segs,i,SEG)->key == key) { setZone(look, arrp(look->segs,i,SEG)->x1, arrp(look->segs,i,SEG)->x2 + 1) ; break ; } } } else if (x1 > x2) fMapRC (look) ; } if (!look) aceOutPrint (result_out, "// gif seqget error: no coordinate data for %s\n", name(key)) ; if (word) messfree(word) ; return look ; usage: if (word) messfree(word) ; aceOutPrint (result_out, "// gif seqget error: usage: seqget [sequence [-coords x1 x2]]\n") ; return NULL ; } /* fMapGifGet */ /* this is my version of fMapGifGet() which allows for selection of which */ /* sets of features will be included. */ /* */ FeatureMap fMapEdGifGet(ACEIN command_in, ACEOUT result_out, FeatureMap look) { KEY key ; int x1, x2 ; char *word ; STORE_HANDLE handle = NULL ; DICT *features = NULL ; FmapDisplayData *fmap_data = NULL ; BOOL already_features, inclusive ; fMapInitialise() ; if (look) fMapDestroy (look) ; /* could check for possible reuse */ if (!(word = aceInWord(command_in))) /* try to get from graph */ { if (fMapActive (0, 0, 0, &look)) return look ; else { aceOutPrint (result_out, "// gif seqget error: no sequence attached to active window\n") ; return 0 ; } } handle = handleCreate() ; if (!lexword2key (word, &key, _VSequence)) { aceOutPrint (result_out, "// gif seqget error: Sequence %s not known\n", word) ; goto usage ; } already_features = FALSE ; x1 = 1 ; x2 = 0 ; /* default for whole sequence */ while (aceInCheck (command_in, "w")) { word = aceInWord(command_in) ; if (strcmp (word, "-coords") == 0) { if (!aceInInt (command_in, &x1) || !aceInInt (command_in, &x2) || (x2 == x1)) goto usage ; } else if (strcmp (word, "-feature") == 0 && aceInCheck (command_in, "w")) { if (already_features) goto usage ; features = dictHandleCreate(16, handle) ; addToSet(command_in, features) ; inclusive = TRUE ; already_features = TRUE ; } else if (strcmp (word, "+feature") == 0 && aceInCheck (command_in, "w")) { if (already_features) goto usage ; features = dictHandleCreate(16, handle) ; addToSet(command_in, features) ; inclusive = FALSE ; already_features = TRUE ; } else goto usage ; } if (features) { fmap_data = (FmapDisplayData *)halloc(sizeof(FmapDisplayData), handle) ; fmap_data->destroy_obj = FALSE ; fmap_data->include_features = inclusive ; fmap_data->features = features ; } if (x1 < x2 || (x1 == 1 && x2 == 0)) look = fMapCreateLook(key, x1, x2, FALSE, fmap_data) ; else look = fMapCreateLook(key, x2, x1, FALSE, fmap_data) ; if (!look) { aceOutPrint (result_out, "// gif seqget error: no coordinate data for %s\n", name(key)) ; if (handle) handleDestroy (handle) ; return NULL ; } /* set the zone */ if (x1 == 1 && x2 == 0) { int i ; for (i = 1 ; i < arrayMax(look->segs) ; ++i) { if (arrp(look->segs,i,SEG)->key == key) { setZone(look, arrp(look->segs,i,SEG)->x1, arrp(look->segs,i,SEG)->x2 + 1) ; break ; } } } else if (x1 > x2) fMapRC (look) ; if (handle) handleDestroy (handle) ; return look ; usage: aceOutPrint (result_out, "// gif fexseqget error: usage: " "fexseqget [sequence [-coords x1 x2][-feature feat1|feat2]]\n") ; if (handle) handleDestroy (handle) ; return NULL ; } /**************************************************************/ void fMapGifDisplay (ACEIN command_in, ACEOUT result_out, FeatureMap look) { int x1 = 0, x2 = 0 ; char *word ; Graph oldGraph = graphActive() ; char *title = NULL ; BOOL new = FALSE ; if (!look || !command_in || !result_out) messcrash("Bad args to fMapGifDisplay") ; while ((word = aceInWord(command_in))) { if (strcmp(word, "-visible_coords") == 0) { if (!aceInInt (command_in, &x1) || !aceInInt (command_in, &x2) || (x2 == x1)) goto usage ; if (x1 > x2) { int t = x1 ; x1 = x2 ; x2 = t ; } } else if (strcmp(word, "-new") == 0) new = TRUE ; else if (strcmp(word, "-title") == 0) { if ((word = aceInWord(command_in))) title = strnew(word, 0) ; else goto usage ; } else goto usage ; } if (!look->graph || !graphActivate(look->graph)) { Graph curr_fmap = GRAPH_NULL ; curr_fmap = displayGetGraph("FMAP") ; if (curr_fmap && !graphActivate(curr_fmap)) messcrash("Could not activate current fmap, graph id: %d.", curr_fmap) ; if (!new && curr_fmap) { /* Replace existing current fmap with new one. */ fMapDestroyCallback() ; /* kill look currently attached */ graphRegister(DESTROY, fMapDestroyCallback) ; /* reregister */ graphClear() ; look->graph = curr_fmap ; displaySetGraph("FMAP", look->graph) ; /* reset fmap in display system, was * removed by fMapDestroyCallback() */ } else { /* Keep existing fmap if a new one is required. */ if (new && curr_fmap) displayPreserve() ; /* Create a new fmap graph */ look->graph = displayCreate("FMAP"); graphActivate(look->graph) ; graphRegister(RESIZE, fMapRefresh) ; graphRegister(DESTROY, fMapDestroyCallback) ; graphRegister(KEYBOARD, fMapKbd) ; graphRegister(MESSAGE_DESTROY, displayUnBlock) ; /* If there was no current graph, we need to set it up in "display" system. */ if (!curr_fmap) displaySetGraph("FMAP", look->graph) ; } if (!title) title = name(look->seqOrig) ; fMapAttachToGraph (look) ; mapAttachToGraph (look->map) ; look->map->min = 0 ; look->map->max = look->length ; if (!x1 && !x2) /* default to show the active zone */ { if (look->flag & FLAG_REVERSE) { x1 = COORD(look, look->zoneMax) + 1 ; x2 = COORD(look, look->zoneMin) ; } else { x1 = COORD(look, look->zoneMin) ; x2 = COORD(look, look->zoneMax) - 1 ; } } if (!getenv("ACEDB_PROJECT") && (look->zoneMax - look->zoneMin < 1000)) { mapColSetByName(look->map, "DNA Sequence", TRUE) ; mapColSetByName(look->map, "Coords", TRUE) ; } } if (title) graphRetitle(title) ; graphFitBounds (&mapGraphWidth,&mapGraphHeight) ; if (x1 || x2) /* show this region */ { halfGraphHeight = 0.5*(mapGraphHeight - topMargin) ; look->map->centre = look->origin + (x1+x2) / 2.0 - 1 ; look->map->mag = (mapGraphHeight - topMargin - 5) / (1.05 * (x2 - x1)) ; if (x1 > x2) { look->map->mag = -look->map->mag ; fMapRC (look) ; } } fMapDraw (look, 0) ; fMapSelect (look) ; { float d = (mapGraphHeight - topMargin - 5.0) / (look->map->mag * 1.05) ; float xx2 = look->map->centre + 0.5*d ; x1 = (int)(xx2 - d + 0.5) ; x2 = (int)(xx2 + 0.5) ; aceOutPrint(result_out, "// FMAP_COORDS %s %d %d", freeprotect (name(look->seqKey)), x1, x2) ; aceOutPrint(result_out, " -visible_coords %d %d\n", COORD(look,x1), COORD(look, x2)) ; } /* clumsy...uuggghhh */ goto end ; usage: aceOutPrint (result_out, "// gif seqdisp error: usage: seqdisp [-visible_coords x1 x2]\n") ; end: graphActivate (oldGraph) ; return ; } /* fMapGifDisplay */ /**************************************************************/ void fMapGifActions (ACEIN command_in, ACEOUT result_out, FeatureMap look) { char *word; if (!look) return; while (aceInCheck (command_in, "w")) { word = aceInWord (command_in) ; if (strcmp (word, "-dna") == 0) { fMapToggleDna (look) ; fMapDraw (look,0) ; } else if (strcmp (word, "-gf_features") == 0) { fMapAddGfSegs(look) ; /* includes fMapDraw() */ } else if (strcmp (word, "-hide_header") == 0) { fMapToggleHeader(look); fMapDraw (look, 0) ; } else if (strcmp (word, "-rev_comp") == 0) { fMapCompToggleTranslation (look) ; fMapRC (look) ; fMapDraw (look, 0) ; } else goto usage ; } return ; usage: aceOutPrint (result_out, "// gif seqactions error: usage: SEQACTIONS [-dna] [-gf_features] [-hide_header] [-rev_comp]\n") ; return; } /* fMapGifActions */ /**************************************************************/ void fMapGifColumns (ACEIN command_in, ACEOUT result_out, FeatureMap look) { char *word; if (!look || !look->map) return ; while (aceInCheck (command_in, "w")) { word = aceInWord (command_in) ; if (strcmp (word, "-on") == 0) { word = aceInWord (command_in) ; if (!word) goto usage ; if (!mapColSetByName (look->map, word, TRUE)) aceOutPrint (result_out, "// gif seqcolumns error: unknown column %s", word); } else if (strcmp (word, "-off") == 0) { word = aceInWord (command_in) ; if (!word) goto usage ; if (!mapColSetByName (look->map, word, FALSE)) aceOutPrint (result_out, "// gif seqcolumns error: unknown column %s", word); } else if (strcmp (word, "-list") == 0) { STORE_HANDLE hh = handleCreate(); Array colList = mapColGetList(look->map, hh); int i; BOOL status; char *name; for (i = 0; i < arrayMax(colList); i++) { name = array(colList, i, char*); status = mapColSetByName(look->map, name, -1); aceOutPrint (result_out, "%s \"%s\"\n", status ? "+" : "-", name); } handleDestroy(hh); } else goto usage ; } return ; usage: aceOutPrint (result_out, "// gif seqcolumns error: usage: SEQCOLUMNS {-on name} {-off name} {-list}\n") ; return; } /* fMapGifColumns */ /**************************************************************/ void fMapGifAlign (ACEIN command_in, ACEOUT result_out, FeatureMap look) { int x1, x2, pepFrame = 0 ; char *word ; BOOL isPep = FALSE ; ACEOUT fo = 0, dump_out = 0; if (!look) return; /* by default output goes to same place as normal command results */ dump_out = aceOutCopy (result_out, 0); while (aceInCheck (command_in, "w")) { word = aceInWord (command_in) ; if (strcmp (word, "-coords") == 0 && aceInInt (command_in, &x1) && aceInInt (command_in, &x2) && (x2 != x1)) { setZoneUserCoords (look, x1, x2) ; } else if (strcmp (word, "-file") == 0 && aceInCheck (command_in, "w") && (fo = aceOutCreateToFile (aceInPath (command_in), "w", 0))) { /* replace default output with file-output */ aceOutDestroy (dump_out); dump_out = fo; } else if (strcmp (word, "-peptide") == 0) isPep = TRUE ; else goto usage ; } pepFrame = (pepFrame + 999) % 3 ; fMapDumpAlign (look, isPep, dump_out) ; aceOutDestroy (dump_out); return ; usage: aceOutDestroy (dump_out); aceOutPrint (result_out, "// gif seqalign error: usage: SEQALIGN [-coords x1 x2] [-file fname] \n") ; return; } /* fMapGifAlign */ /**************************************************************/ static void addToSet (ACEIN command_in, DICT *dict) /* little utility for parsing below */ { char cutter, *cutset, *word ; if (aceInStep (command_in, '"')) cutset = "|\"" ; else cutset = "| \t" ; while ((word = aceInWordCut (command_in, cutset, &cutter))) { dictAdd (dict, word, 0) ; if (cutter != '|') break ; } return; } /* addToSet */ /******************/ void fMapGifFeatures (ACEIN command_in, ACEOUT result_out, FeatureMap look) { KEY refseq = 0 ; int x1, x2, offset ; int version = 2 ; char *word ; STORE_HANDLE handle = handleCreate() ; DICT *sourceSet = dictHandleCreate (16, handle) ; DICT *featSet = dictHandleCreate (16, handle) ; BOOL isList = FALSE ; BOOL result ; ACEOUT fo = 0, dump_out = 0; if (!look) return; /* by default output goes to same place as normal command results */ dump_out = aceOutCopy (result_out, 0); while (aceInCheck (command_in, "w")) { word = aceInWord (command_in) ; if (strcmp (word, "-coords") == 0 && aceInInt (command_in, &x1) && aceInInt (command_in, &x2) && (x2 != x1)) { setZoneUserCoords (look, x1, x2) ; } else if (strcmp (word, "-file") == 0 && aceInCheck (command_in, "w") && (fo = aceOutCreateToFile (aceInPath (command_in), "w", 0))) { /* replace default output with file output */ aceOutDestroy (dump_out); dump_out = fo; } else if (strcmp (word, "-refseq") == 0 && aceInCheck (command_in, "w") && (lexword2key (aceInWord (command_in), &refseq, _VSequence))) { ; } else if (strcmp (word, "-version") == 0 && aceInInt (command_in, &version)) { ; } else if (strcmp (word, "-list") == 0) { isList = TRUE ; } else if (strcmp (word, "-source") == 0 && aceInCheck (command_in, "w")) { addToSet (command_in, sourceSet) ; } else if (strcmp (word, "-feature") == 0 && aceInCheck (command_in, "w")) { addToSet (command_in, featSet) ; } else goto usage ; } result = fMapDumpGFF(look, version, &refseq, &offset, isList, sourceSet, featSet, dump_out) ; aceOutDestroy (dump_out); if (result) aceOutPrint (result_out, "// FMAP_FEATURES %s %d %d\n", freeprotect (name(refseq)), look->zoneMin+1+offset, look->zoneMax+offset) ; handleDestroy (handle) ; return ; usage: aceOutDestroy (dump_out); aceOutPrint (result_out, "// gif seqfeatures error: usage: SEQFEATURES [-coords x1 x2] ") ; aceOutPrint (result_out, "[-file fname] [-refseq sequence] [-version 1|2] ") ; aceOutPrint (result_out, "[-list] [-source source(s)] [-feature feature(s)]\n") ; handleDestroy (handle) ; return ; } /* fMapGifFeatures */ /**************************************************************/ void fMapGifDNA (ACEIN command_in, ACEOUT result_out, FeatureMap look) { KEY key ; int x1, x2 ; ACEOUT fo = 0, dump_out = 0; char *word ; if (!look) return; /* by default output goes to same place as normal command results */ dump_out = aceOutCopy (result_out, 0); while (aceInCheck (command_in, "w")) { word = aceInWord (command_in) ; if (strcmp (word, "-coords") == 0 && aceInInt (command_in, &x1) && aceInInt (command_in, &x2) && (x2 != x1)) { setZoneUserCoords (look, x1, x2) ; } else if (strcmp (word, "-file") == 0 && aceInCheck (command_in, "w") && (fo = aceOutCreateToFile (aceInPath (command_in), "w", 0))) { /* write to file instead of default output */ aceOutDestroy (dump_out); dump_out = fo; } else goto usage ; } key = 0 ; x1 = look->zoneMin ; x2 = (look->zoneMax - 1) ; if (!fMapFindSpan(look, &key, &x1, &x2)) key = look->seqKey ; /* fMapFindSpan returns 0-based coords and we need to report in 1-based coords. */ x1++, x2++ ; /* If we extend this to do spliced and dna stuff then * WE NEED TO CHECK OF look->seq is same as look->seqorig and if not, we should output * look->seqorig for the spliced dna _only_.....*/ dnaDumpFastA (look->dna, look->zoneMin, look->zoneMax-1, messprintf("%s/%d-%d", name(key), x1, x2), dump_out) ; aceOutPrint (result_out, "// FMAP_DNA %s %d %d\n", freeprotect (name(key)), x1, x2) ; aceOutDestroy (dump_out); return; usage: aceOutDestroy (dump_out); aceOutPrint (result_out, "// gif seqdna error: usage: SEQDNA [-coords x1 x2] [-file fname]\n") ; return; } /* fMapGifDNA */ /**************************************************************/ void fMapGifRecalculate(ACEIN command_in, ACEOUT result_out, FeatureMap look) { Graph oldGraph = graphActive() ; Graph display_graph = GRAPH_NULL ; /* Actually look is almost irrelevant...I think ??? */ if (!command_in || !result_out) messcrash("Bad args to fMapGifRecalculate") ; if (!look) { if ((display_graph = displayGetGraph("FMAP"))) { graphActivate(display_graph) ; graphAssFind(&GRAPH2FeatureMap_ASSOC, &look) ; } } if (look) { fMapPleaseRecompute(look) ; fMapDraw (look, 0) ; } else aceOutPrint (result_out, "// gif seqrecalc error: could not find fmap to recalculate\n") ; graphActivate (oldGraph) ; return ; } /**************************************************************/ /**************************************************************/ /**************************************************************/ static void fMapZoneKeySet (void) { SEG *seg ; int i, n = 0 ; KEYSET ks = keySetCreate () ; FeatureMap look = currentFeatureMap("fMapZoneKeySet") ; for (i = 0 ; i < arrayMax(look->segs) ; ++i) { seg = arrp(look->segs, i, SEG) ; if (seg->x1 > look->zoneMax || seg->x2 < look->zoneMin) continue ; if (class(seg->key)) keySet(ks, n++) = seg->key ; } keySetSort(ks) ; keySetCompress (ks) ; if (!keySetMax(ks)) { messout ("The active zone does not contain any key") ; keySetDestroy (ks) ; } else keySetNewDisplay (ks, messprintf("%s %d %d", name(look->seqKey),look->zoneMin - look->origin + 1, look->zoneMax - look->origin)) ; } /* Routine to dump out stats for a look. * The intention is augment this routine to produce more comprehensive and * selectable stats to aid with debugging/performance testing etc. */ static void dumpstats(FeatureMap look) { int seg_totals[CLONE_END + 1] = {0} ; int i ; SEG *seg ; int seg_size = sizeof(SEG) ; float megabyte = 1048576.0 ; for (i = 0 ; i < arrayMax(look->segs) ; i++) { seg = arrp(look->segs,i,SEG) ; seg_totals[seg->type]++ ; } printf("\n\t\tStats for FMap of %s\n" , name(look->seqKey)) ; for (i = 0 ; i <= CLONE_END ; i++) { if (seg_totals[i]) { printf("%15s: %15d\t\t(%2d%%)\n", fMapSegTypeName[i], seg_totals[i], (seg_totals[i] * 100)/arrayMax(look->segs)) ; } } printf("%15s: %15d\t\t(%.1fMB)\n", "Total segs", arrayMax(look->segs), (((float)seg_size * (float)arrayMax(look->segs)) / megabyte)) ; return ; } /****************** end of file *****************/