PageRenderTime 61ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 1ms

/code_comm/SeqFuncs.c

http://research-code-base-animesh.googlecode.com/
C | 1221 lines | 715 code | 146 blank | 360 comment | 357 complexity | 4f1d097f6c9a33be26bb695adc201b9a MD5 | raw file

Large files files are truncated, but you can click here to view the full file

    

  1. /* tacg - a command line tool for the analysis of nucleic acids and protein  */

    2. 

  2. /* Copyright Š 1994-2001 Harry J Mangalam, tacg Informatics

    3. 

  3. (mangalam@tacgi.com | mangalam@home.com | 949 856 2847) */

    4. 

  4. 

    5. 

  5. /* $Id: SeqFuncs.c,v 1.13 2001/11/14 01:18:09 hjm Exp $  */

    6. 

  6. 

    7. 

  7. /* The use of this software (except that by Harald T. Alvestrand, which is described in 'udping.c')

    8. 

  8.    and that by James Knight, which is described in 'seqio.c' is bound by the notice that appears 

    9. 

  9.    in the file 'tacg.h' which should accompany this file.  In the event that 'tacg.h' is not bundled 

    10. 

  10.    with this file, please contact the author.

    11. 

  11. */

    12. 

  12. 

    13. 

  13. #include <stdio.h>

    14. 

  14. #include <ctype.h> 

    15. 

  15. #include <string.h>

    16. 

  16. #include <stdlib.h>

    17. 

  17. #include <sys/stat.h>

    18. 

  18. #include <unistd.h>

    19. 

  19. #include <regex.h>

    20. 

  20.   

    21. 

  21. #include "seqio.h"

    22. 

  22. #include "tacg.h" 

    23. 

  23. 

    24. 

  24. /* takes care of:

    25. 

  25. 	- reading rules from the file (path validated via SearchPaths()), 

    26. 

  26.    - parsing the file into the Rules struct

    27. 

  27.    - validating that the rules are formed correctly (parens match up, logicals are correct)

    28. 

  28.    - making sure that each of the subpattern names has been loaded in RE (have to hash the

    29. 

  29.    	names as in satisfying the -x flag in ReadEnz().

    30. 

  30.    - eventually, create --rules 'thisrule,thatrule,mars,promoter, etc'

    31. 

  31.      (can be specified w/o --rulefile in which case it will read the deafult 'rules.data'.)

    32. 

  32.      as opposed to current --rule for specifying a rule at a time from the commandline.

    33. 

  33.      (which should write out that rule to the tacg.patterns file for incoporation into the 

    34. 

  34.      main rules file.  require that it be named as well. 

    35. 

  35.      ie: --rule 'Name,((()()())(())),width'

    36. 

  36.      

    37. 

  37.      - open the file, getting path name from F.RuleFile

    38. 

  38.      - read in the strings, concatting as we go, then validate the logic 

    39. 

  39.          and stuff all the bits into the Rule_struct

    40. 

  40.      - at each iter, make sure there's enough mem

    41. 

  41.      - when finished with the file read, and all the bits are in Rule_struct, go thru

    42. 

  42.      		Rules[].Name and check against the names in RE[].E_nam.

    43. 

  43.      - make 

    44. 

  44.  */    

    45. 

  45. int tacg_Eval_Rules_from_File(struct SE_struct SelEnz[MAX_SEL_ENZ],  char datestamp[80]) {

    46. 

  46.                               /*   long ProtoNameHash[], int NumREs, */

    47. 

  47. 	

    48. 

  48.    FILE *fpRF; 

    49. 

  49.    char junk[10], *ctemp1, ct, *rulestring, *packed_RS, *rstr;

    50. 

  50.    short NameCheck[RE_NAME_TAB_SIZ]; /* indices from realhash(PatName) to check against multiple identicals  */

    51. 

  51.    int SEi=0;

    52. 

  52.    long j, EOpRS=0,

    53. 

  53.         window_size = 0,

    54. 

  54.         NumRules=0, 

    55. 

  55.    	  rstr_len, linenum=0,

    56. 

  56.         RS_Cnt = 0, /* RuleStruct Counter */

    57. 

  57.         EOAM_RS = INIT_NUM_RULES; /* End Of Allocated Memory for the Rule_Struct; allocated in tacg */

    58. 

  58.         

    59. 

  59.    short cont = 0; /* continue indicator */

    60. 

  60. 

    61. 

  61.    /* alloc the Rule_struct Rules */

    62. 

  62.    if ((Rules = calloc(EOAM_RS, sizeof(*Rules))) == NULL) {

    63. 

  63.       BadMem("Can't init mem for Rules", 1);

    64. 

  64.    }

    65. 

  65. 

    66. 

  66. 	if ((ctemp1 = calloc(512,sizeof(char))) == NULL) { BadMem("Can't init mem for 'ctemp1'", 1); }

    67. 

  67. 	if ((rulestring = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rulestring'", 1); }

    68. 

  68. 	if ((packed_RS = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'packed_RS'", 1); }

    69. 

  69. 	if ((rstr = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rstr'", 1); }

    70. 

  70.    memset(NameCheck,0, sizeof(short)*RE_NAME_TAB_SIZ); /* make sure everything is set to 0 */

    71. 

  71.    

    72. 

  72.    if ((fpRF = fopen(F.RuleFile,"r")) == NULL) {  /* or the standard/optional one from SetFlags() */

    73. 

  73.    	fprintf(stderr,"Cannot open the Rules file \"%s\" for reading!!\n", F.RuleFile); /* DON'T HIDE behind -V */

    74. 

  74.       exit(1); /* print an error and die gracefully */

    75. 

  75.    } else if (F.Verbose > 0) { fprintf(stderr,"Rule File %s opened OK for reading\n", F.RuleFile); }

    76. 

  76.    

    77. 

  77. 	/*  Scan thru comments in rebase file to separator ("..")   */

    78. 

  78.    junk[0] = 'k';   /* make sure the strcmp below fails the first time */

    79. 

  79.    while ((feof(fpRF) == 0) && (strncmp(junk, "..",2) != 0))  { 

    80. 

  80.       fscanf (fpRF, "%2s", junk); ct = 'm'; /* and read to end of line - this needs to be able to deal */

    81. 

  81.       while (ct != '\n' && ct != '\r' && (feof(fpRF) == 0)) {  /* with end of line conditions for Mac, Win, and unix - what have I missed? */

    82. 

  82.          ct = fgetc(fpRF); linenum++;

    83. 

  83.       }

    84. 

  84.    }

    85. 

  85. 

    86. 

  86. /* we're at "..", so go into file and start slurping directly data into struct

    87. 

  87.    except for RE_rawsite which has to be filtered a couple ways before being acceptable */

    88. 

  88.    

    89. 

  89.    ct = 'm';

    90. 

  90.    while (fgets(rstr, 512, fpRF) && (NumRules <= MAX_NUM_RULES)) { 

    91. 

  91.       /* check the size of the rule struct to see if it's going to run into the end of the alloc'ed mem */

    92. 

  92.       if (RS_Cnt == EOAM_RS) { /* then we need to alloc more mem for the Rules struct */

    93. 

  93.           EOAM_RS +=  RULES_INCR;

    94. 

  94.           if ((Rules = realloc(Rules, sizeof(*Rules)*EOAM_RS)) == NULL) {

    95. 

  95.             BadMem("Failed to get more mem for Rules.\n", 1);

    96. 

  96.           }  

    97. 

  97.       } 

    98. 

  98.       

    99. 

  99.       /* now the whole line's in rstr; use strsep to parse it, make decisions re: parsing, skipping */

    100. 

  100.       j = stripwhitespace(rstr,0,0);

    101. 

  101.       if (rstr[0] != ';'  && (strlen(rstr) != 0)) {  /* if the line hasn't been commented out, get all the bits */

    102. 

  102.          /* REMEMBER!  strsep chews the string to pieces and changes the start of the string */

    103. 

  103.       if (j>0 && F.Verbose > 0) {fprintf(stderr,"DEBUG: stripped %ld whitespace characters from RuleName\n", j);} 

    104. 

  104.          

    105. 

  105.          strcpy(packed_RS, rstr);  /* copy over to pass to tacg_SlidingWindow - may not be needed */

    106. 

  106.          

    107. 

  107. 

    108. 

  108.          /* now have to break on ',' to get all the rulestring in 1 shot 

    109. 

  109.          	- also have to detect if the last char is a continuation char '\' */

    110. 

  110.             

    111. 

  111.          /* if last char is '\', then rulestring is incomplete and there won't be a window term  */

    112. 

  112.          cont = 0;

    113. 

  113.          EOpRS = strlen(packed_RS) - 1; /* */

    114. 

  114.          if (packed_RS[EOpRS] == '\\') { /* then there's a continuation char */

    115. 

  115.          	cont = 1; 

    116. 

  116.             /* keep on going until there's no more continuation lines */

    117. 

  117.             while ((cont == 1) && fgets(rstr, 512, fpRF)) {  /*     */

    118. 

  118.                j = stripwhitespace(rstr,0,0); /* no whitespaces, eols */

    119. 

  119.             	rstr_len = strlen(rstr); /* length excluding terminating \n */

    120. 

  120.                EOpRS = strlen(packed_RS);

    121. 

  121.                if (EOpRS + rstr_len >= MAX_LEN_RULESTRING) {

    122. 

  122.                   fprintf(stderr, "Total length of the rulestring exceeded; change MAX_LEN_RULESTRING \n"

    123. 

  123.                   "or check your continuation characters for the rule %s\n", ctemp1);

    124. 

  124.                   exit (1);

    125. 

  125.                } 

    126. 

  126.                /* rstr should be the right size to copy mod the possible final \ and \n */

    127. 

  127.             	strcpy(packed_RS+EOpRS-1,rstr); /* copy all but the last char over */

    128. 

  128. 

    129. 

  129. /*                fgets(rstr, 512, fpRF); */

    130. 

  130. 

    131. 

  131.                

    132. 

  132.                fprintf(stderr, "\nend char = %s\n", rstr+(strlen(rstr)-3) );

    133. 

  133.                if (rstr[(strlen(rstr)-1)] == '\\' ){ cont = 1;}   

    134. 

  134.                else                                { cont = 0;}   

    135. 

  135.             }

    136. 

  136.             fprintf(stderr, "\n\n\npacked_RS = %s\n\n", packed_RS);

    137. 

  137. 	         if ((rstr = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rstr'", 1); }

    138. 

  138.             strcpy(rstr,packed_RS);   

    139. 

  139.             

    140. 

  140.          }  

    141. 

  141.          /* coming out of the continuation loop, packed_RS contains the entire Name,RuleString,Window concat

    142. 

  142.             that the whole thing should look like. rstr has only the last line of a continued rule, so it needs to 

    143. 

  143.             be synced, so copy packed_RS to rstr again */

    144. 

  144.                                  

    145. 

  145.          ctemp1 = strsep(&rstr, ","); /* using ctemp1 (=Name) for congruence with previous code  */

    146. 

  146.          if (strlen(ctemp1) > 10) {

    147. 

  147.             fprintf(stderr, "Name %s is longer than allowed (%d) - check it!\n", ctemp1, MAX_PAT_NAME_LEN);

    148. 

  148.             exit(1);

    149. 

  149.          } 

    150. 

  150.             

    151. 

  151.             

    152. 

  152.          /* so now packed_RS has the whole packed rulestring in it at this point; rstr has everything 

    153. 

  153.             BUT the RuleName in it. So now we finish off the line, so break on ',' to get the rulestring

    154. 

  154.             and ONLY THEN break on the ',' to separate the rule term from the sliding window term 

    155. 

  155.             (since already have processed the Name) */  

    156. 

  156.          if (strlen(rstr) >= MAX_LEN_RULESTRING) {

    157. 

  157.             fprintf(stderr, "Total length of the rulestring exceeded; change MAX_LEN_RULESTRING \n"

    158. 

  158.                             "or check your continuation characters for the rule %s\n", ctemp1);

    159. 

  159.             exit (1);

    160. 

  160.          }

    161. 

  161.          

    162. 

  162.          /* break on the comma so we can ignore the window term - this probably needs some more error-checking */

    163. 

  163.          memset(rulestring, 0, MAX_LEN_RULESTRING);

    164. 

  164.          rulestring = strsep(&rstr,",");  /* rulestring now contains ONLY the packed rulestring */

    165. 

  165.          window_size = (long)atoi(strsep(&rstr,","));

    166. 

  166. 

    167. 

  167.          if (window_size < 1) {  /* now test whether windowsize & rulestring are valid  */

    168. 

  168.             fprintf (stderr, "Window size for %s is less than 1, which is silly.\n", ctemp1);

    169. 

  169.             exit (1);

    170. 

  170.          }

    171. 

  171.          /* Are the component pattern Names valid? ie is this name represented in RE by being loaded 

    172. 

  172.             from the RE database?  need to hash the Names of the individual patterns and see if they 

    173. 

  173.             hit any of the Names that were loaded from whatever rebase was being used.

    174. 

  174.             Soooooo, go thru the rulestring and hash all words breaking on ()|^& to cut out the 'Name:min:Max' stanzas */

    175. 

  175.          tacg_Fill_SelEnz_fr_RuleString(SelEnz, &SEi, rulestring, NameCheck, ctemp1);

    176. 

  176.          

    177. 

  177. 

    178. 

  178. 

    179. 

  179.          /* stuff ALL the bits (Name, rulestring, window_size) into the struct (the validity of the rule itself 

    180. 

  180.             tested in the evaluation code) */

    181. 

  181.          strcpy(Rules[RS_Cnt].Name, ctemp1);

    182. 

  182.          if ((Rules[RS_Cnt].RuleString = calloc((strlen(rulestring)+1),sizeof(char))) == NULL) { 

    183. 

  183.             BadMem("Can't init mem for 'Rules[RS_Cnt].RuleString'", 1); 

    184. 

  184.          }

    185. 

  185.          strcpy(Rules[RS_Cnt].RuleString, rulestring);

    186. 

  186.          Rules[RS_Cnt].Window = window_size;         

    187. 

  187. 

    188. 

  188.          if ((rstr = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rstr'", 1); }

    189. 

  189.          RS_Cnt++;

    190. 

  190. 		} else {

    191. 

  191. /*          fprintf(stderr, "\nDEBUG: line ignored: %s\n", rstr); */

    192. 

  192.       }

    193. 

  193.    }

    194. 

  194.    F.Xplicit = -1; /* notify ReadEnz */

    195. 

  195.    F.Rules = RS_Cnt-1; /* does this have to be decr by one 1st?? */

    196. 

  196.    return SEi;

    197. 

  197. }

    198. 

  198. 

    199. 

  199. 

    200. 

  200. /* ***************************************************************************************************** */

    201. 

  201. 

    202. 

  202. 

    203. 

  203. /* tacg_Fill_SelEnz_fr_RuleString() takes the SelEnz sruct and the rulestring (and a few other housekeeping vars 

    204. 

  204.    and parses the rulestring passed to it to populate SelEnz and do a namecheck to make sure that the names being 

    205. 

  205.    parsed out of the rulestring haven't been seen before and aren't otherwise objectionable      

    206. 

  206.    takes as args:

    207. 

  207.        SelEnz, 

    208. 

  208.        the starting SEi (as an int * as it's changed in the fn())

    209. 

  209.        the rulestring (copied in the function so that the original rulestring is untouched)

    210. 

  210.        NameCheck (has to survive fn() to check for multiple identicals over multiple rulestrings,

    211. 

  211.           as well as to check for multiple identicals in the same rulestring (A&B|A&C)

    212. 

  212.        ctemp1 - Name of the Rule (for error messages)

    213. 

  213.     and back comes 

    214. 

  214.        the (more) populated SelEnz

    215. 

  215.        and the ending SEi.  

    216. 

  216. */

    217. 

  217. 

    218. 

  218. void tacg_Fill_SelEnz_fr_RuleString(struct SE_struct SelEnz[MAX_SEL_ENZ], int *SEi, 

    219. 

  219.             char *rulestring, short NameCheck[RE_NAME_TAB_SIZ],  char *ctemp1) {

    220. 

  220.    char *Name_mM, *rstr, *PatName = "            ",

    221. 

  221.    *amin = "           ", 

    222. 

  222.    *aMax = "           ";

    223. 

  223.    

    224. 

  224.    int ind;

    225. 

  225.  	if ((Name_mM = calloc(33,sizeof(char))) == NULL) { BadMem("Can't init mem for 'Name_mM'", 1); } 

    226. 

  226. 	if ((rstr = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rstr'", 1); }

    227. 

  227.    strcpy(rstr,rulestring); /* re-using rstr; mem should be fine after re-calloc above */

    228. 

  228.    Name_mM[0] = '\0'; /*  to pass on 1st pass */

    229. 

  229.    while (Name_mM[0] == '\0') { 

    230. 

  230.    Name_mM = strsep(&rstr, "()|^&"); } /* this should break out a 'Name:#:#' set */

    231. 

  231. 

    232. 

  232.    while (Name_mM != 0L) {

    233. 

  233.       if (Name_mM[0] != '\0') {

    234. 

  234.          PatName = strsep(&Name_mM, ":");

    235. 

  235.          ind = (int) realhash(PatName, RE_NAME_TAB_SIZ);

    236. 

  236. /*          fprintf(stderr, "DEBUG: realhash '%s' = %d\n", PatName, ind); */

    237. 

  237.          if (strlen(PatName) < 11) {

    238. 

  238.             if (NameCheck[ind] == 0) { 

    239. 

  239.                strcpy(SelEnz[*SEi].PName,PatName);

    240. 

  240.                NameCheck[ind] = 1; /* now mark it to make sure that we don't take it again */

    241. 

  241.                if ((amin = strsep(&Name_mM, ":")) == '\0') { 

    242. 

  242.                  SelEnz[*SEi].min =   0; /* RE[ii].min =  */

    243. 

  243.                } else { 

    244. 

  244.                   SelEnz[*SEi].min =  atoi(amin); /* RE[ii].min =  */

    245. 

  245.                }

    246. 

  246.                if ((aMax = strsep(&Name_mM, ":")) == '\0') { 

    247. 

  247.                   SelEnz[*SEi].Max = 32000; /* RE[ii].Max = */

    248. 

  248.                } else { 

    249. 

  249.                   SelEnz[*SEi].Max =  atoi(aMax); /* RE[ii].Max =   */

    250. 

  250.                }

    251. 

  251.                (*SEi)++; 

    252. 

  252. 

    253. 

  253.             } else { /* else it's already been seen in another rule or sub-rule and we don't care */

    254. 

  254. /*                fprintf(stderr, "DEBUG: The pattern named %s has already been seen previously; continuing .. \n", PatName); */

    255. 

  255.             }

    256. 

  256.          } else { 

    257. 

  257.             fprintf (stderr, "\nIn Rule named %s, Pattern Name %s is too long (must be <=10 chars)\n", 

    258. 

  258.             ctemp1, PatName); 

    259. 

  259.             exit(1);

    260. 

  260.          }

    261. 

  261.       }

    262. 

  262.       Name_mM = strsep(&rstr, "()|^&");  /* get the next one */

    263. 

  263.    }/* @ end of loop, all constituent patterns are loaded into SelEnz or the app has exited on an err */

    264. 

  264. }

    265. 

  265. 

    266. 

  266. 

    267. 

  267. /* ***************************************************************************************************** */

    268. 

  268. 

    269. 

  269. 

    270. 

  270. 

    271. 

  271. 

    272. 

  272. 

    273. 

  273. /* RE_Filter returns 1 if the RE indexed by Proto[REi] matches the tests for

    274. 

  274. 	Magnitude

    275. 

  275.    Number of Cuts

    276. 

  276.    Cost

    277. 

  277.    Overlap

    278. 

  278.    etc that are can be used to fine-tune the selection of REs */

    279. 

  279. int RE_Filters_OK(int REi) {  /* REi should be the ProtoCutters index, as that is the only RE el that we have to check */

    280. 

  280. 	int OK=0;

    281. 

  281. /*    fprintf(stdout, "\n REi = %d, RE name = %s \n", REi, RE[REi].E_nam); */

    282. 

  282.    if (((RE[REi].E_mag   >= (int)F.Mag)       &&    /* E_mag must be >= to -n flag */

    283. 

  283.         (REi >= F.RE_ST)                              &&    /* Don't EVER want dam/dcm sites printed */

    284. 

  284.         (RE[REi].E_Ncuts >= (int)F.Min)       &&    /* Ncuts >= -m, */

    285. 

  285.         (RE[REi].E_Ncuts <= (int)F.Max)       &&    /* Ncuts <= -M */

    286. 

  286.         (RE[REi].Cost    >= (int)F.Cost)      && 

    287. 

  287. /* the dam/dcm part of this COULD be rolled in if the RWC (aka base_real_pos) was included as an optional 

    288. 

  288. 	argument (or if -1, ignore it), but right now it doesn't seem to be used in the same way, so I'll let

    289. 

  289.    the code stay but commented out for now */

    290. 

  290. /*        (DamDcmOK(base_real_pos, REi) == 1)   &&     print only RE's !masked by dam/dcm methylases */

    291. 

  291.         

    292. 

  292.     (F.Overhang == 1))      || /* if want all, only 1st has to be true to short circuit */

    293. 

  293.        ((RE[REi].E_olap > 0) && (F.Overhang == 5)) || /* if olap < 0 && overhang = 5' */

    294. 

  294.        ((RE[REi].E_olap < 0) && (F.Overhang == 3)) || /* if olap > 0 && overhang = 3' */

    295. 

  295.        ((RE[REi].E_olap ==0) && (F.Overhang == 0)))  /* if olap = 0 && overhang = blunt */ 

    296. 

  296. 	   { OK = 1; } 

    297. 

  297.    else { OK = 0; }

    298. 

  298.    return OK;

    299. 

  299. }

    300. 

  300. 

    301. 

  301. /* if either the --dam or --dcm flags have been used, DamDcmOK() checks to see 

    302. 

  302. 	if the REi under consideration should be printed - has it been masked or is 

    303. 

  303.    it OK to print it.  What it does: passed the RE index of the enzyme, it

    304. 

  304.    checks to see if that RE is possibly affected and if so, whether the 'raw'

    305. 

  305.    hit recorded in DD.Dat is present in RE[].Sites.  The index to the

    306. 

  306.    last-checked position in RE[].Sites is kept in: 

    307. 

  307.    RE[].Sites[RE[].Sites[0]]       as 'RE[].Sites[0]' points to the end of the array.

    308. 

  308.    In this way, only have to incr by 1 to get to the next check point.

    309. 

  309.    RWC = Real World Coordinates

    310. 

  310.    REi = RE index

    311. 

  311.    Hmmm - this could be the prototype for an overall RE filter to be used in a number of cases..

    312. 

  312. */

    313. 

  313. 

    314. 

  314. /* Clone() goes thru RE[] and figures out which REs match the conditions of the option flag, where

    315. 

  315.    --clone '#_#,#x#,#x#,#_#' (up to MAX_NUM_CLONE_RANGES).  '_' indicates that the range cannot be cut;

    316. 

  316.    'x' indicates that it can be cut.  

    317. 

  317.    There will probably be some interactions between --clone and other options, but

    318. 

  318.    for now, finish the basic functionality and sort out the interactions later.  In fact,

    319. 

  319.    RE[].. should not be a problem as it will be populated only with sites that are requested, 

    320. 

  320.    so that if only a subset or special file of REs are used, only they will be used for 

    321. 

  321.    for the clone() analysis.

    322. 

  322.    This fn() is called AFTER the digest is done in the normal way and should only 

    323. 

  323.    need the number of REs (NumREs - not modified by this so it doesn;t ave to be 

    324. 

  324.    passed as a pointer) as a max to go thru to determine the ranges that agree with the 

    325. 

  325.    ranges specified.

    326. 

  326. */

    327. 

  327. /* should really be running thru Protos[] from 0->ProtoCutters */

    328. 

  328. void Clone(int ProtoCutters, int *Protos) {

    329. 

  329. 	long i, c, j, m, p, r, 

    330. 

  330.    Best[F.Clone+2][ProtoCutters], 	/* Best = array that stores and sorts the best REs for a set of cloning criteria

    331. 

  331.    												last el is for storing summary truths to check if ALL rules have been met */

    332. 

  332.    BadProtos[ProtoCutters];  	/* BadProtos keep track of Protos that map inside a forbidden zone and therefore have to 

    333. 

  333.    										be removed from contention after the fact.  If a Protos index shows up in BadProtos,

    334. 

  334.                                  it cannot be a possible cloning RE (unless cloner wants to try a partial digest.. */

    335. 

  335.    short N_Crits;

    336. 

  336.    int LadderProtos[ProtoCutters];

    337. 

  337.    

    338. 

  338.    N_Crits =  F.Clone; /*  simplicity; F.Clone stores the # of ranges in the --clone option string */ 

    339. 

  339.    for (i=0; i<= F.Clone; i++) {

    340. 

  340.    	memset(Best[i],0,sizeof(int)*ProtoCutters);

    341. 

  341.    }

    342. 

  342.    memset(BadProtos,0,sizeof(long)*ProtoCutters);

    343. 

  343.    memset(LadderProtos,0,sizeof(int)*ProtoCutters);

    344. 

  344.    

    345. 

  345.    /* Need additional rules to elim REs.  

    346. 

  346.    	If an RE hits in the forbidden zone, that RE is then added to a FORBIDDEN array to be checked post-processing.

    347. 

  347.    	now #x# = MUST HIT range  Should there be an option to set prefs for this? Or just show how many

    348. 

  348.       	there are in that zone?

    349. 

  349.       How should REs that hit in zones NOT explicitly marked as allowed (x) or forbidden (_) be treated?  In typical 

    350. 

  350.       	cloning, don't want lots of extra cuts as it increases ends to contend with.   

    351. 

  351.    */

    352. 

  352.    

    353. 

  353.    for (r=0; r<=N_Crits; r++) { /* for each range parsed out in SetFlags() */

    354. 

  354.    	for (i=0; i<ProtoCutters; i++) {

    355. 

  355. 	      p = Protos[i]; /* p is the alias for the RE prototype index, so RE[p].etc will hit only prototypes */

    356. 

  356.          if (Clone_S[r].to_cut == 1) { /* if this is a range where we REQUIRE hits... */  

    357. 

  357.             m = Clone_S[r].matching_REs[1]; /* [1] is the 1st free position of the array */

    358. 

  358. 				for (j=1; j<RE[p].Sites[0]; j++) { /* RE[p].Sites[0] holds the index of the 1st FREE position in the Sites array */

    359. 

  359.                if (RE[p].Sites[j] >= Clone_S[r].begin &&     /* site must match BOTH conditions */

    360. 

  360.                    RE[p].Sites[j] <= Clone_S[r].end ) {

    361. 

  361.                  	/* do we need a per-range log for positives & negatives? */

    362. 

  362.                	Clone_S[r].matching_REs[m++] = i;  		/* incr m, store the Proto index */

    363. 

  363.                   Clone_S[r].matching_REs[1] = m;      	/* keep track of the 1st free el as well */

    364. 

  364.                   /* and check for realloc conditions */

    365. 

  365.                   if ((m+1) == Clone_S[r].matching_REs[0]) { realloc_Clone_matching_REs(m, r); }

    366. 

  366.                } /* else we silently ignore the bad hits */

    367. 

  367.             }

    368. 

  368.             Clone_S[r].matching_REs[1] = m; /* and then store the position back into the array */

    369. 

  369.          } else { /* this is a range we where we will NOT TOLERATE hits */

    370. 

  370.             m = Clone_S[r].matching_REs[1]; /* [1] is the 1st free position of the array  */

    371. 

  371.             for (j=1; j<RE[p].Sites[0]; j++) { /* RE[p].Sites[0] holds the index of the 1st FREE position in the Sites array */

    372. 

  372.             	if (RE[p].Sites[j] <= Clone_S[r].begin ||      /* site can match either condition */

    373. 

  373.                    RE[p].Sites[j] >= Clone_S[r].end ) {

    374. 

  374.                	Clone_S[r].matching_REs[m++] = i;  /* incr m, store the Proto index */

    375. 

  375.                   Clone_S[r].matching_REs[1] = m;            /* keep track of the 1st free el as well */

    376. 

  376.                   /* and check for realloc conditions  */

    377. 

  377.                   if ((m+1) == Clone_S[r].matching_REs[0]) {realloc_Clone_matching_REs(m, r); }

    378. 

  378.                } else { /* else we have to track the REs to mask out the results by adding their index to BadProtos[] */

    379. 

  379.                	BadProtos[i]++; /* if an el is > 0, it's EVIL */

    380. 

  380.                }

    381. 

  381.             }

    382. 

  382.             Clone_S[r].matching_REs[1] = m; /* and then store the position back into the array */                     	

    383. 

  383.          }

    384. 

  384.       }

    385. 

  385.    }

    386. 

  386.    

    387. 

  387.       /* now we've run thru all the criteria for judging whether the REs map into the ranges 

    388. 

  388.    	now have to output the data in some usable form.  1st, are there any REs that match ALL requirements?

    389. 

  389.       2nd output all ranges with all matching REs, sorted alphabetically.  Should --clone automatically set

    390. 

  390.       --sites, so can output all the sites at the same time??  Hmmm...

    391. 

  391.    */

    392. 

  392.    /* 1st are there any that match ALL crits? for all crits, are there any indices that appear in all of them??  

    393. 

  393.    	set up array of length N_Crits and increment each RE index as go thru all the ranges.  Also has the effect of

    394. 

  394.       being able to sort those REs that match the MOST crits as well. */

    395. 

  395.    for (i=0; i<=N_Crits; i++) {   /*  cycle thru all the crits */

    396. 

  396.    	for (j=2; j<Clone_S[i].matching_REs[1]; j++) { /* cycle thru all the prototype REs that match for each crit */

    397. 

  397.       	Best[i][Clone_S[i].matching_REs[j]]++; 	/* this will count the # of times a proto came up positive in the criteria  */

    398. 

  398.       }

    399. 

  399.    } /*  so now have all the REs counted. */

    400. 

  400.    

    401. 

  401.    /* and print the results in some kind of visually appealing way. */

    402. 

  402.    /* 1st print out any REs that met ALL crits.. */

    403. 

  403.   	fprintf(stdout,"\n\n== Here's the Clone output for the %d Criteria:\n", (F.Clone+1));

    404. 

  404. 	for (i=0; i<=F.Clone; i++) { fprintf(stdout, "    Rule %ld:  %s\n", i+1, Clone_S[i].range_rule); }

    405. 

  405. 

    406. 

  406.    for (j=0; j <= F.Clone; j++) {  /* sum the rules results to generate a summary */

    407. 

  407.    	for (i=0; i<ProtoCutters; i++) {

    408. 

  408.    	   if (Best[j][i] > 0 && BadProtos[i] == 0) {

    409. 

  409.          	Best[F.Clone+1][i]++; /* for each rule passed, incr the truth element once */

    410. 

  410.          }

    411. 

  411.       }

    412. 

  412.    }

    413. 

  413.    /* and finally print out the Protos that match ALL the rules BUT this should probably be mod to allow selection and

    414. 

  414.    	filtering by the standard rules -o -n -m -M --cost, etc but later... this will require a test fn() that tests all the

    415. 

  415.       Protos against all the various restrictions, but it will be applicable to at least 4 cases */

    416. 

  416.    fprintf(stdout,"\nREs that met ALL criteria:\n");

    417. 

  417. /* 

    418. 

  418. 	for (i=0; i<ProtoCutters; i++) {   

    419. 

  419.       if (Best[F.Clone+1][i] == F.Clone+1) {

    420. 

  420.    	   fprintf(stdout,"\n  RE[Protos[i]].E_nam);

    421. 

  421.       }

    422. 

  422.    } 

    423. 

  423. */

    424. 

  424.    i = c = m = 0;

    425. 

  425.    while (i<ProtoCutters) {

    426. 

  426.       if (m >= F.Width) { m = 0; fprintf(stdout, "\n"); }

    427. 

  427.    	if (Best[F.Clone+1][i] == F.Clone+1 ) { 

    428. 

  428.    		fprintf(stdout,"%10s", RE[Protos[i]].E_nam);

    429. 

  429.          LadderProtos[i] = Protos[i];

    430. 

  430.          c++; m += 10;

    431. 

  431.       }

    432. 

  432.       i++;

    433. 

  433.    } 

    434. 

  434.    if (c == 0) { fprintf(stdout,"There were no REs that met ALL criteria\n"); }

    435. 

  435. 	else {

    436. 

  436.    	PrintGelLadderMap(i, F.Seq_Len, LadderProtos, 2);

    437. 

  437.    }

    438. 

  438.    memset(LadderProtos,0,sizeof(int)*ProtoCutters);

    439. 

  439.    

    440. 

  440.    fprintf(stdout,"\n\nResults (listed by Rules) that met SOME criteria & did not violate a hits-forbidden range.");

    441. 

  441.    for (j=0; j<=F.Clone; j++) {

    442. 

  442.    	fprintf(stdout, "\n\n--- For Rule %ld (%s), the following REs matched:\n\n", j+1, Clone_S[j].range_rule);

    443. 

  443.    	i = c = m = 0;

    444. 

  444.       while (i<ProtoCutters) {

    445. 

  445.       	if (m >= F.Width) { m = 0; fprintf(stdout, "\n"); }

    446. 

  446.    	   if (Best[j][i] > 0 && BadProtos[i] == 0) { 

    447. 

  447.    		   fprintf(stdout,"%10s", RE[Protos[i]].E_nam);

    448. 

  448.             LadderProtos[i] = Protos[i];

    449. 

  449.             c++; m += 10; 

    450. 

  450.          }

    451. 

  451.       	i++;

    452. 

  452.       }

    453. 

  453.       if (c == 0) { fprintf(stdout,"There were no REs that met EVEN SOME of the criteria\n"); }

    454. 

  454. 	   else {

    455. 

  455.    	   PrintGelLadderMap(i, F.Seq_Len, LadderProtos, 2);

    456. 

  456.       }

    457. 

  457.       memset(LadderProtos,0,sizeof(int)*ProtoCutters);

    458. 

  458.    }

    459. 

  459.    

    460. 

  460.    fprintf(stdout,"\n\n--- Results (listed by RE) that met SOME criteria & did not violate a hits-forbidden range\n");

    461. 

  461.    fprintf(stdout,"\n   RE Name");

    462. 

  462.    for (j=0; j<=F.Clone; j++) { fprintf(stdout,"  Rule %ld", j+1); }

    463. 

  463.    fprintf(stdout,"\n");

    464. 

  464.    for (i=0; i<ProtoCutters; i++) {

    465. 

  465.    	if (BadProtos[i] == 0) { 

    466. 

  466.    	   fprintf(stdout,"\n%10s", RE[Protos[i]].E_nam);

    467. 

  467.    	   for (j=0; j<=F.Clone; j++) {

    468. 

  468.    	      if (Best[j][i] > 0) { 

    469. 

  469.    		      fprintf(stdout,"    X   ");

    470. 

  470.             } else {

    471. 

  471.             	fprintf(stdout,"        ");

    472. 

  472.             }

    473. 

  473.          }

    474. 

  474.       }

    475. 

  475.    }

    476. 

  476.    fprintf(stdout,"\n\n");

    477. 

  477. };

    478. 

  478. 

    479. 

  479. 

    480. 

  480. void realloc_Clone_matching_REs(long current, long index) {

    481. 

  481. 	long newsize;

    482. 

  482.    /* and check for realloc conditions */

    483. 

  483.    if ((current+1) == Clone_S[index].matching_REs[0]) { /* if we're nearly up to the alloc'ed mem, need to realloc */

    484. 

  484.       newsize = Clone_S[index].matching_REs[0] + CLONE_INCR_MRE; /* this is how big the new size should be  */

    485. 

  485.       if ((Clone_S[index].matching_REs = (long *) realloc(Clone_S[index].matching_REs, sizeof(long) * newsize)) == NULL) {

    486. 

  486.          BadMem("Can't realloc mem for Clone_S[].matching_REs", 1); 

    487. 

  487.       } else { Clone_S[index].matching_REs[0] = newsize; }   /* new alloc range */

    488. 

  488.    }

    489. 

  489. };

    490. 

  490.    

    491. 

  491. 

    492. 

  492. /* the example function included below is a mock-up to show how to include your own function.

    493. 

  493. 	Obviously, there are 10s of functions included with this source code and some of them may be

    494. 

  494.    immediately applicable to your work.  Feel free to bash them into a useful state for yourself.

    495. 

  495.    If you think that you've created a widely useful function that you would like to have included 

    496. 

  496.    with tacg, please give me a shout, and I'll include a pointer to it or provide other information

    497. 

  497.    about it on my web site.	  

    498. 

  498.    

    499. 

  499.    void Example(void){

    500. 

  500. 

    501. 

  501. 	}   

    502. 

  502. */

    503. 

  503. 

    504. 

  504. 

    505. 

  505. /* LinearMap() is simply a functionization of the Linear Mapping chunk that

    506. 

  506. 	should have been functionized long ago */

    507. 

  507.    

    508. 

  508. void LinearMap(long seq_len, int basesPerLine, int NumREs, char *sequence, int max_okline, 

    509. 

  509. 				char Codons[][][], int codon_Table, int *Protos) {

    510. 

  510. 

    511. 

  511. 	/* Declarations */

    512. 

  512.    int i, d, k, mm=0, rc, n_letters=0, Xn=0, okline, HTML, spacer, spacermod, stanzaspace,

    513. 

  513.    	   min_okline, tic_line, Cur_RE, block_Cnt, block_repeat, block_cut_pos, lastline = O_SEQ_LINE+2, hits,

    514. 

  514.    		ok2wr[MAX_OUTPUT_LINES], p_width, abcp, in_OL=-1, seq_offset, tripclick=0;

    515. 

  515. 

    516. 

  516.    long  *BackCheck, /* BackCheck[] checks proto entries in the Linear map output to remove duplicates */

    517. 

  517.          DCnt, Cnt_reset=0, base_real_pos, numstart; 

    518. 

  518.          					

    519. 

  519.    char  ctemp1[30], *s,  Prot_out[MAX_BASES_PER_LINE],

    520. 

  520.    		O_txt[MAX_OUTPUT_LINES][O_LMARG+O_RMARG+MAX_BASES_PER_LINE];

    521. 

  521.    

    522. 

  522. /*    char Codons[MAX_ORGS][N_CODONS][7]; */

    523. 

  523.    char os = '~';  /* started out as '<' to mark other strand, but causes problems with HTML generation */

    524. 

  524. 

    525. 

  525.    BackCheck = (long *)calloc((size_t)NumREs+2, sizeof(long));  /* Needs to be '+2' as combination cut will add one and 1 extra */

    526. 

  526.    if (BackCheck == NULL) BadMem("Can't init mem for BackCheck", 1);

    527. 

  527. 

    528. 

  528.    s = (char *) calloc(256, sizeof(char)); 

    529. 

  529.    if (s == NULL) BadMem("failed to calloc mem for 's'", 0);

    530. 

  530. 	memset(O_txt,' ',(size_t)((O_LMARG+O_RMARG+MAX_BASES_PER_LINE)*MAX_OUTPUT_LINES));  /* set all of O_txt to blanks - works!*/

    531. 

  531. 	p_width = O_LMARG + O_RMARG + basesPerLine;	  

    532. 

  532.    

    533. 

  533. 	/* Routine to print out the text of the digestion, with names going where they should go

    534. 

  534. 	a la Strider - have to keep track of how long the names are and where an allowable place 

    535. 

  535. 	to write the name - */

    536. 

  536. 

    537. 

  537. 	/* In brief, the following function, does this: 

    538. 

  538. 	It reads thru DD->Dat (the array that keeps track of the linear sequence of cuts in the sequence) 

    539. 

  539. 	and matches the position of the enzymes to the current block being printed, taking into account 

    540. 

  540. 	the cut offset of the RE from the 1st base of the recognition sequence.  The tricky part is that 

    541. 

  541. 	because of this latter functionality, it has to check the sequence that will be the *following* 

    542. 

  542. 	block to see if there are any RE recognition seqs within the 1st few bases that will cause the actual 

    543. 

  543. 	cut site to be in the *previous* block.  Doing this seems to take a great amount of the cpu time - can 

    544. 

  544. 	probably be optimized a good deal more, but it's still >10 faster than GCG ;). */

    545. 

  545. 	

    546. 

  546.    HTML = F.HTML;         /* simplicity */

    547. 

  547.    n_letters = F.Lab1or3; /* simplicity */

    548. 

  548.    Xn = F.XLinFr;         /* simplicity */

    549. 

  549.    numstart = F.NumStart;  /* simplicity */

    550. 

  550.    

    551. 

  551. 	fprintf(stdout, "\n\n");

    552. 

  552. 	if (HTML > -1) fprintf(stdout, "<H3><A NAME=\"linmap\">");

    553. 

  553. 	fprintf(stdout, "== Linear Map of Sequence:\n"); /* Write the header */ 

    554. 

  554. 	if (HTML > -1) fprintf(stdout, "</A></H3>\n");

    555. 

  555. 

    556. 

  556. 	/* inits, etc */ 

    557. 

  557. 	if (seq_len % basesPerLine == 0)  block_repeat = (int) seq_len/basesPerLine; /*  # times thru the compose/print cycle */

    558. 

  558. 	else  block_repeat = (int) seq_len/basesPerLine + 1;  

    559. 

  559. 	DCnt=2;  /* start the DD->Dat counter at the beginning of REAL DATA [0]=Cnt; [1]=Siz */

    560. 

  560. 	base_real_pos = DD->Dat[DCnt];

    561. 

  561.   

    562. 

  562.   /* spacer dictates the line spacing between the 1st 3 and 2nd 3 translation stanzas; if 6 there

    563. 

  563.       will be no space between the 2 stanzas, if 7 there will be 1 space; could make this another

    564. 

  564.       option but sheesh... */

    565. 

  565.   spacer = 7;  spacermod = 0;

    566. 

  566.   if (F.Strands == 1) { spacermod++; }   

    567. 

  567.   if (F.Tics == 0) { spacermod++;}

    568. 

  568.   spacer = spacer - spacermod;

    569. 

  569.   if (spacer == 7) stanzaspace = spacermod + 1; 

    570. 

  570. 

    571. 

  571. 	for (block_Cnt=1; block_Cnt <= block_repeat; block_Cnt++) {

    572. 

  572.       /* hits is set to 1 if there are hits on the line; if no hits, decr the top buffer line to output 

    573. 

  573.       	so that there is a constant 1 blank line spacer between stanzas */

    574. 

  574.    	hits = -1;  

    575. 

  575. 		/* if there was an OVERLAP failure previously, fix it */ 

    576. 

  576. 		if (in_OL == 1) {

    577. 

  577.          memset(BackCheck,0,sizeof(long)*(NumREs+2)); /* clear BackCheck[] as it's filled with stale data */

    578. 

  578.          DCnt = Cnt_reset;  /* and reset DCnt to that pointer */ 

    579. 

  579. 		}

    580. 

  580.       in_OL = -1;  /* and reset 'in OVERLAP' marker for new block */

    581. 

  581. 		base_real_pos = DD->Dat[DCnt];  /* Make sure that ' base_real_pos' is set */

    582. 

  582.  		okline = min_okline = O_SEQ_LINE-2; /* set the min (highest on page) line of buffer to print out */

    583. 

  583. 

    584. 

  584. 		/* Set up the seq and its rev compl in the output buffer */

    585. 

  585. 		memset(ok2wr,0,sizeof(int)*MAX_OUTPUT_LINES);  /* set/reset all of ok2wr to 0's for the new block*/

    586. 

  586. 		/* next line was deleted in a brief fork of tacg.c - note this if get unexlained segfs */

    587. 

  587.       memset(BackCheck, 0, sizeof(long)*NumREs); /* reset the BackCheck[] */

    588. 

  588. 

    589. 

  589. 		/*write the correct numbering in Left margin  */

    590. 

  590.     if (F.NumStart < 1 && (((block_Cnt-1)*basesPerLine) + 1 + F.NumStart) == 0) {numstart++; tripclick++;}

    591. 

  591.     sprintf((char *)s, "%7ld",(((block_Cnt-1)*basesPerLine) + 1 + numstart)); /* numstart is just a + offset */

    592. 

  592.     memcpy(&O_txt[O_SEQ_LINE][0], s, 7);

    593. 

  593. 		/* and in the right margin - combine with the above bit? */

    594. 

  594.     if (F.NumStart < 0 && (tripclick || (((block_Cnt)*basesPerLine + numstart)) > 0)) {

    595. 

  595.       if (!tripclick) {tripclick++; numstart++;}

    596. 

  596.     }

    597. 

  597. 		sprintf((char *)s, "%7ld",(((block_Cnt)*basesPerLine + numstart)));

    598. 

  598. 		memcpy(&O_txt[O_SEQ_LINE][O_LMARG+basesPerLine], s, 7);

    599. 

  599. 		/* and drop in the sequence directly after it. */

    600. 

  600. 

    601. 

  601. 		if (block_Cnt != block_repeat) k = basesPerLine; /* test for end of sequence to end gracefully */ 

    602. 

  602. 		else k = (int)seq_len - (block_Cnt-1)*basesPerLine;

    603. 

  603. 		if ((F.SeqBeg != 1 || F.NumStart != 0)  &&

    604. 

  604.         (F.Strands != 1 || F.Tics != 0)) { /* if subsequence or dif #ing, print the #s of the original seq too */

    605. 

  605. 			/* write the correct numbering in Left margin  */

    606. 

  606. 			sprintf((char *)s, "%7d",((int)(((block_Cnt-1)*basesPerLine)+F.SeqBeg))); /* these 2 lines could be combined...? */

    607. 

  607. 			memcpy(&O_txt[O_SEQ_LINE+1][0], s, 7);

    608. 

  608. 			/* and in the right margin - combine with the above bit? */

    609. 

  609. 			sprintf((char *)s, "%7d",(int)(((block_Cnt)*basesPerLine - 1 + F.SeqBeg)));

    610. 

  610. 			memcpy(&O_txt[O_SEQ_LINE+1][O_LMARG+basesPerLine], s, 7);

    611. 

  611. 		}

    612. 

  612. 

    613. 

  613. 		memcpy(&O_txt[O_SEQ_LINE][O_LMARG], &sequence[(block_Cnt-1)*basesPerLine+BASE_OVERLAP+1], (size_t)k); 

    614. 

  614. 

    615. 

  615.       if (F.Strands == 2) { /* if we want the bottom strand as well */

    616. 

  616. 			Rev_Compl(sequence+BASE_OVERLAP+1+((block_Cnt-1) * basesPerLine), s, (long)k);

    617. 

  617. 			/* and then plunk the converted sequence into O_txt with a similar statement */

    618. 

  618. 			memcpy(&O_txt[O_SEQ_LINE+1][10], s, (size_t)k);

    619. 

  619.       }

    620. 

  620. 		if (F.Tics == 1) {  /* if we want tics  */

    621. 

  621.          /* write out the minor and major tics below the reverse complement */

    622. 

  622. 			tic_line = O_SEQ_LINE+2;

    623. 

  623.          /* if there's only 1 strand AND we want tics, have to move them up 1 line */

    624. 

  624.          if (F.Strands == 1) tic_line -= 1; 

    625. 

  625. 			for (i=0;i<(int)(basesPerLine/10);i++) memcpy (&O_txt[tic_line][i*10+O_LMARG],"    ^    *",10);

    626. 

  626.       }

    627. 

  627. 

    628. 

  628. 		/* Following call to Translate has to calculate the # of bases at each call; otherwise end up with 

    629. 

  629. 			overrun condition at end; "k" in code above 

    630. 

  630.          This stanza doesn't have to change at all, even after the -T/-t change of 9.2.98 */

    631. 

  631. 

    632. 

  632.       rc = seq_len % 3;  /* rc is the correction factor for the trans, if not a perfect div by 3;

    633. 

  633.          if it is a perfect div by 3, then rc = 0, and there won't be any change */

    634. 

  634. 		if (n_letters != 0) {

    635. 

  635. 			for (mm=0;(mm<(int)Xn && mm<3); mm++) {

    636. 

  636. 			   if (k % 3 != 0) k = ((int) k /3)*3;   /* also added this code to Translate() for better modularity */

    637. 

  637. 			   /* this chunk doen't care about rc */

    638. 

  638.             seq_offset = ((block_Cnt-1)*basesPerLine+BASE_OVERLAP+1);

    639. 

  639. 			   Translate(sequence+seq_offset+mm, Prot_out, k, n_letters, Codons, codon_Table); 

    640. 

  640. 			   /* need to place the translation frame labels here before memcpy() call */

    641. 

  641.             d = 0;

    642. 

  642.             if (F.Strands == 1) { d += 1; }  

    643. 

  643.             if (F.Tics == 0) { d += 1; }

    644. 

  644.  			   sprintf(ctemp1, "%d", mm+1);   O_txt[O_SEQ_LINE+3+mm-d][0] = ctemp1[0]; 

    645. 

  645. 			   memcpy(&O_txt[O_SEQ_LINE+3+mm-d][O_LMARG+mm], &Prot_out, (size_t)k); /* and copy it into place */

    646. 

  646. 			   lastline = O_SEQ_LINE + 3 + mm - d; /* keep track of last line of Xl */

    647. 

  647.             

    648. 

  648.             if (Xn == 6) { /* but if it's 6 frames of Xlation, need to do oppo at the same time */

    649. 

  649.                /* but not EXACTLY the same - seq to be has to be shifted slightly to account for 

    650. 

  650.                   those seqs not exactly / 3 - should have to only add a correction to the expr below

    651. 

  651.                   AND add the same correction to the offset for printing it out, so this chunk DOES care

    652. 

  652.                   about rc */

    653. 

  653. 

    654. 

  654.                Anti_Par(sequence+seq_offset-mm+rc, s, k); /* get the *antipar* sequence */

    655. 

  655.                /* Translate() below shouldn't have to change */

    656. 

  656. 			      Translate (s, Prot_out, k, n_letters, Codons, codon_Table); /* Xlate it into N-letter code */

    657. 

  657. 			      if (n_letters == 1) Reverse (Prot_out); /* to match the 'backwards' DNA strand */

    658. 

  658. 			      else Triplet_Reverse (Prot_out);

    659. 

  659. 			      /* need to place the translation frame labels here before memcpy() call */

    660. 

  660. 			      sprintf(ctemp1, "%d", mm+4);   O_txt[O_SEQ_LINE+spacer+mm][0] = ctemp1[0];

    661. 

  661. 			      memcpy(&O_txt[O_SEQ_LINE+spacer+mm][O_LMARG-mm+rc], &Prot_out, (size_t)k); /* and copy it into place */

    662. 

  662.       			lastline = O_SEQ_LINE + spacer + mm; /* this is the last line of the Xl, so last to print */

    663. 

  663.             }

    664. 

  664. 			}

    665. 

  665. 		} else {

    666. 

  666.       d = 0;

    667. 

  667.       if (F.Strands == 1) { d += 1; }  

    668. 

  668.       if (F.Tics == 0) { d += 1; } 

    669. 

  669.       lastline = O_SEQ_LINE + 2 - d;

    670. 

  670.     }

    671. 

  671. 

    672. 

  672. 		/* This 'while' loop tests each text block for REs and also prints out those blocks that have no 

    673. 

  673. 		cuts in them at all but which have to be printed anyway....*/

    674. 

  674. 		while ((DD->Dat[DCnt] != -2) && (base_real_pos < ((block_Cnt * basesPerLine) + BASE_OVERLAP))) {

    675. 

  675. 			DCnt++; 

    676. 

  676. 			/* if we're not at the end and the real pos'n <  current block + the overlap.. */

    677. 

  677. 			/* !! Now DD->Dat[DCnt] should point to the corresponding  RE #  */

    678. 

  678. 

    679. 

  679. 			 /* if following test is true (they're !=) then it's safe to write it; */

    680. 

  680. 			 /* if it's false, then the same proto has already written to the same space */

    681. 

  681. 			if (BackCheck[labs(DD->Dat[DCnt])] != labs(DD->Dat[DCnt-1])) { /* if !=, do this, but 1st... */

    682. 

  682. 			   BackCheck[labs(DD->Dat[DCnt])] = labs(DD->Dat[DCnt-1]);  /*  ... set it = so that it won't do this one again */

    683. 

  683.              okline = O_SEQ_LINE-2; 

    684. 

  684. 

    685. 

  685.        /* This is where the RE filter could be inserted so that the DD data is ignored unless the RE that's 

    686. 

  686.          referenced matches the parameters set in the filters (can take this straight from the ladder-map 

    687. 

  687.          functions.  */

    688. 

  688.          /* direct from GelLad.. 

    689. 

  689.             these are the conditions that have to be passed in order to have the REs placed in the

    690. 

  690.             linear map. */

    691. 

  691.             /*  DD->Dat[DCnt] is still referencing a possibly neg #, so have to take the abs() */

    692. 

  692.             Cur_RE = labs((int)DD->Dat[DCnt]); 

    693. 

  693.             /*  this whole 'if' test thingie can and should be reduced to one extensible function call 

    694. 

  694.                (using the embedded DamDcmOK() as a template), but ... not yet... */

    695. 

  695.             if ((DamDcmOK(base_real_pos, Cur_RE) == 1)   &&  RE_Filters_OK(Cur_RE)) { /* if olap = 0 && overhang = blunt */ 

    696. 

  696.                /* this is the overall test for whether the RE will be printed out in the Linear Map or not  

    697. 

  697. 						if the RE DOESN'T meet the crit, it just gets ignored here, but will be printed in other outputs */

    698. 

  698. 

    699. 

  699. 			      block_cut_pos =  (int)base_real_pos - 1 - ((block_Cnt-1) * basesPerLine); /* see if it's this easy?! */

    700. 

  700. 			      /* Now locate it in the block, checking for under- and over-runs */

    701. 

  701. 			      if (block_cut_pos >= 0) { /* if the position is greater than the beginning of the block */

    702. 

  702. 			         if (block_cut_pos <= basesPerLine) { /* and if it's less than the end of the block... */

    703. 

  703. 			            abcp = O_LMARG + block_cut_pos - 1; /* abcp = adjusted block_cut_pos - used x times below */

    704. 

  704. 			            O_txt[O_SEQ_LINE-1][abcp] = '\\';   /* write the 'cut' indicator */

    705. 

  705.                      hits = 1;

    706. 

  706. 			            /* locate the position for writing the name */

    707. 

  707. 			            while (block_cut_pos < ok2wr[okline]) okline--; /* then go up to the lowest 'free' line */

    708. 

  708. 			            if (okline != O_SEQ_LINE-2)  { /* but if can't place the name on the lowest line.. */

    709. 

  709. 			               i = O_SEQ_LINE-1; k = 0;   /* check if there's any space lower between previous names */

    710. 

  710. 			               /* following 'while' must be ordered like it is - i will be decr only if k == 0 */

    711. 

  711. 			               while (k == 0 && --i != okline) { /* leaving 1 space before and after the name */

    712. 

  712. 			                  if (O_txt[i][abcp-1] == ' ' &&                /* this can be replaced with something */

    713. 

  713. 			                      O_txt[i][abcp+RE[Cur_RE].E_nam_l] == ' ' &&  /* more efficient later */

    714. 

  714. 			                      O_txt[i][abcp+3] == ' ' &&       /* must check for a space here... */               

    715. 

  715. 			                      O_txt[i][abcp+4] == ' ' )  k = 1;  /* as well as here */

    716. 

  716. 			                    /* i = vertical counter ~= okline, k = 1 when we find enuf space */

    717. 

  717. 			               }

    718. 

  718. 			               okline = i;

    719. 

  719. 			            }

    720. 

  720. 			            if (okline < 2) {

    721. 

  721. 			               fprintf(stderr, "!! Output Buffer Blown around base %ld - degeneracy caused too many Enz's\n"

    722. 

  722. 			                  "to be matched!! Increase the stringency of Enz selection w/ -n -o -M flags\n", base_real_pos);

    723. 

  723. 			               exit(0);

    724. 

  724. 			            }

    725. 

  725. 			            /* and memcpy the RE name from the struct to the output array */

    726. 

  726. 			            memcpy (&O_txt[okline][abcp],&RE[Cur_RE].E_nam, (size_t)RE[Cur_RE].E_nam_l);

    727. 

  727.                      if (DD->Dat[DCnt] < 0) { /* if the match is in the (-) orientation */

    728. 

  728.                         memcpy (&O_txt[okline][abcp+RE[Cur_RE].E_nam_l], &os, 1);

    729. 

  729.                      }

    730. 

  730. 

    731. 

  731. 			            /* and incr the ok2wr posn for that line by the length of the RE name + 1 */

    732. 

  732. 			            /* 'ok2wr[okline]' below is modified only if we didn't find any lower spaces */

    733. 

  733. 			            if (block_cut_pos+RE[Cur_RE].E_nam_l+1 > ok2wr[okline]) {

    734. 

  734. 			               ok2wr[okline] = block_cut_pos+RE[Cur_RE].E_nam_l + 1; 

    735. 

  735. 			            }

    736. 

  736. 			         } /* if (block_cut_pos < basesPerLine).. */

    737. 

  737. 			         else {  /* otherwise it's in the overlap area from the next block */

    738. 

  738. 			            if (in_OL == -1) {   /* if this is the 1st time in OVERLAP, mark it */

    739. 

  739. 			               in_OL = 1;        /* by setting in_OL to 1  */

    740. 

  740. 			               Cnt_reset = DCnt - 1;  /* Have to back up to 1st RE in the OVERLAP that did not 

    741. 

  741. 			                                                resolve into the current block */

    742. 

  742. 			            }   /* The above marker *should* back the DCnt to the position of the Cur_RE */ 

    743. 

  743. 			         }

    744. 

  744. 			      }     /* if (block_cut_pos >= 0) ... */

    745. 

  745. 			   }

    746. 

  746. 			   base_real_pos = DD->Dat[++DCnt]; /* make sure that this is set before the next loop */

    747. 

  747. 			   if (okline < min_okline) min_okline = okline;

    748. 

  748. 			} else {

    749. 

  749. 			   base_real_pos = DD->Dat[++DCnt]; /* this gonna fix it? */

    750. 

  750. 			}

    751. 

  751. 		}  /*  while (DD->Dat..... */

    752. 

  752. 

    753. 

  753. 		/* Now print out the block we've composed, from min_okline to MAX_OUTPUT_LINES */

    754. 

  754. 		/* adj to min_okline is nec to provide a 1 blank line between stanzas  */

    755. 

  755. 		for (i=min_okline-hits; i<= lastline; i++) fprintf(stdout, "%.*s\n",p_width,O_txt[i]);

    756. 

  756. 

    757. 

  757. 		/* And clean out the lines written into O_txt as well */

    758. 

  758. 		memset(O_txt,' ',(size_t)((O_LMARG + O_RMARG + MAX_BASES_PER_LINE) * MAX_OUTPUT_LINES));  /* set all of O_txt to blanks */

    759. 

  759. 	}   /* for (block_Cnt=1; block_Cnt<block_repeat; block_Cnt++) ... */

    760. 

  760. 	free(BackCheck); /* this line was also deleted from the bike fork of tacg.c at one point

    761. 

  761.                        along with the corresponding memset at ~ line 36 */

    762. 

  762. }

    763. 

  763. 

    764. 

  764. 

    765. 

  765. 

    766. 

  766.    

    767. 

  767. 

    768. 

  768. /* GenerateTOC() simply checks the status of all the major output headers and

    769. 

  769. 	generates the appro HTML to make up a Table of Contents for the results that 

    770. 

  770.    follow.  It does not actually check that they have completed correctly, but

    771. 

  771.    it does keep everything in one place.  Just a series of if()'s  */

    772. 

  772.    

    773. 

  773. void GenerateTOC(char *datestamp) {

    774. 

  774. 

    775. 

  775. 	/* print a header */

    776. 

  776.    fprintf(stdout, "<h2>Table of Results</h2>\n");

    777. 

  777.    /* maybe print some date, version, file, machine run on, domain run on? etc info */

    778. 

  778.    fprintf(stdout, "<b>tacg version %s; Date: %s </b><br>\n"

    779. 

  779.    "<i>Contact Harry Mangalam (mangalam@home.com) with comments, critiques</i><br>\n",TACG_VERSION, datestamp);

    780. 

  780.    

    781. 

  781.    /* Start the UNordered list */

    782. 

  782.    fprintf(stdout, "<ul>\n");

    783. 

  783.    

    784. 

  784.    if (F.Summary != 0) { /* this generates no-cuts and all-cuts */

    785. 

  785.    	fprintf(stdout, "<li><a href=\"#nocuts\">Patterns that DO NOT MAP to this sequence</a>\n");

    786. 

  786.    	fprintf(stdout, "<li><a href=\"#allcuts\">Total Number of Hits per Pattern</a>\n");

    787. 

  787.    }

    788. 

  788.    if (F.Ladder != 0) { /* generates ladder map and infrequent matches */

    789. 

  789.    	fprintf(stdout, "<li><a href=\"#ladmap\">Ladder Map of Pattern Matches</a>\n");

    790. 

  790.    	fprintf(stdout, "<li><a href=\"#summarycuts\">Summary Map of Infrequent Matches</a>\n");

    791. 

  791.    }

    792. 

  792.    if (F.GelLO != 0 || F.GelHI != 0) {

    793. 

  793.    	fprintf(stdout, "<li><a href=\"#gelmap\">Pseudo Gel Map</a>\n");

    794. 

  794.    }

    795. 

  795.    if (F.Sites != 0) {

    796. 

  796.    	fprintf(stdout, "<li><a href=\"#sitetab\">Sites by Pattern</a>\n");

    797. 

  797.    }

    798. 

  798.    if (F.Frags == 1) {

    799. 

  799.    	fprintf(stdout, "<li><a href=\"#fragunsort\">Fragment Sizes (UNSorted)</a>\n");

    800. 

  800.    }

    801. 

  801.    if (F.Frags == 2) {

    802. 

  802.    	fprintf(stdout, "<li><a href=\"#fragsort\">Fragment Sizes (Sorted)</a>\n");

    803. 

  803.    }

    804. 

  804.    if (F.ORFs != -1) {

    805. 

  805.    	fprintf(stdout, "<li><a href=\"#orfs\">Open Reading Frame Analysis</a>\n");

    806. 

  806.    }

    807. 

  807.    if (F.LinMap != 0) {

    808. 

  808.    	fprintf(stdout, "<li><a href=\"#linmap\">Linear Map</a>\n");

    809. 

  809.    }

    810. 

  810.    if (F.GrfBBin != 0) {

    811. 

  811.    	fprintf(stdout, "<li><a href=\"#graphics\">Graphics Plotting Data</a>\n");

    812. 

  812.    }

    813. 

  813.    if (F.Prox != 0) {

    814. 

  814.    	fprintf(stdout, "<li><a href=\"#pair_proximity\">Pairwise Proximity</a>\n");

    815. 

  815.    }

    816. 

  816.    if (F.Rules != 0) {

    817. 

  817.    	fprintf(stdout, "<li><a href=\"#rules_proximity\">Complex Rules-based Proximity</a>\n");

    818. 

  818.    }

    819. 

  819.    if (F.Hookey != -1) {

    820. 

  820.    	fprintf(stdout, "<li><a href=\"#hookeyunsort\">AFLP Analysis - UNSORTED</a>\n");

    821. 

  821.    	fprintf(stdout, "<li><a href=\"#hookeysort\">AFLP Analysis - SORTED</a>\n");

    822. 

  822.    }

    823. 

  823. 

    824. 

  824.    if (F.GrafMap != -1) {

    825. 

  825.    	if (F.PDF == 1) {

    826. 

  826.    		fprintf(stdout, "<li><a href=\"%s/tacg_Map.pdf \">PDF Plasmid Map</a>\n", F.tmppath);

    827. 

  827.       } else {

    828. 

  828.       	fprintf(stdout, "<li><a href=\"%s/tacg_Map.ps \">Postscript Plasmid Map</a>\n", F.tmppath);

    829. 

  829.       }

    830. 

  830.    }

    831. 

  831. /* spares below */

    832. 

  832. /*

    833. 

  833.    if (F. != ) {

    834. 

  834.    	fprintf(stdout, "<li><a href=\"#LABEL\">DESCRIPTION</a>\n");

    835. 

  835.    }

    836. 

  836.  */   

    837. 

  837.    /* below marks the end of the UNordered list and the beginning of the PRE-formatted output */

    838. 

  838.    fprintf(stdout, "</ul><HR noshade size=5><pre>\n"); 

    839. 

  839. /*    fprintf(stdout, "<>  </>\n"); */

    840. 

  840. 

    841. 

  841. }

    842. 

  842.    

    843. 

  843. 

    844. 

  844. 

    845. 

  845. /* if either the --dam or --dcm flags have been used, DamDcmOK() checks to see 

    846. 

  846. 	if the REi under consideration should be printed - has it been masked or is 

    847. 

  847.    it OK to print it.  What it does: passed the RE index of the enzyme, it

    848. 

  848.    checks to see if that RE is possibly affected and if so, whether the 'raw'

    849. 

  849.    hit recorded in DD.Dat is present in RE[].Sites.  The index to the

    850. 

  850.    last-checked position in RE[].Sites is kept in: 

    851. 

  851.    RE[].Sites[RE[].Sites[0]]       as 'RE[].Sites[0]' points to the end of the array.

    852. 

  852.    In this way, only have to incr by 1 to get to the next check point.

    853. 

  853.    RWC = Real World Coordinates

    854. 

  854.    REi = RE index

    855. 

  855.    Hmmm - this could be the prototype for an overall RE filter to be used in a number of cases..

    856. 

  856. */

    857. 

  857. int DamDcmOK(long RWC, int REi) {

    858. 

  858. 	int i, OK=0;

    859. 

  859.    /* 1st, check to see if this RE is sensitive at all and if not, skip it */

    860. 

  860.    /* also skip if it doesn't have ANY cuts at all */

    861. 

  861.    if ((RE[REi].E_Ncuts > 0) && ((F.dam == 1 && RE[REi].dam != 100) || (F.dcm == 1 && RE[REi].dcm != 100))) {

    862. 

  862.       if (RE[REi].Sites[RE[REi].Sites[0]] == 0) {

    863. 

  863.    	   RE[REi].Sites[RE[REi].Sites[0]] = 1; /* 1st time thru, set this to 1 (point it at the 1st real value */

    864. 

  864.       }

    865. 

  865.       OK = 0; 

    866. 

  866.       i = RE[REi].Sites[RE[REi].Sites[0]]; /* shorten: i = the last el to have been checked */

    867. 

  867.       while (RWC >= RE[REi].Sites[i] && OK == 0) {

    868. 

  868.       	if (RWC == RE[REi].Sites[i]) {

    869. 

  869.          	OK = 1; 

    870. 

  870.          }

    871. 

  871.          i++; 

    872. 

  872.       }

    873. 

  873.       RE[REi].Sites[RE[REi].Sites[0]] = i;

    874. 

  874.    } else OK = 1;

    875. 

  875.    if (RE[REi].E_Ncuts == 0) OK = 0; /* just in case */

    876. 

  876.    return OK;

    877. 

  877. }

    878. 

  878. 

    879. 

  879. /* Add_dam and Add_dcm just fill the RE[1] and RE[2] entries for these special 

    880. 

  880.    methylation enzymes. */

    881. 

  881. 

    882. 

  882. void Add_dam(int NumREs) {

    883. 

  883.    strcpy(RE[NumREs].E_nam,"dam\0");

    884. 

  884.    if ((RE[NumREs].E_raw_sit = calloc(5, sizeof(char))) == NULL) BadMem("Add_dam-1", 1);

    885. 

  885.    RE[NumREs].E_raw_sit = "GATC\0"; 

    886. 

  886.    if ((RE[NumREs].E_wsit[1] = calloc(5, sizeof(char))) == NULL) BadMem("Add_dam-2", 1);

    887. 

  887.    RE[NumREs].E_wsit[1] = "gatc\0"; 

    888. 

  888.    strcpy(RE[NumREs].E_hex[1],"gatcnn\0");

    889. 

  889.    RE[NumREs].E_tcut[1] = 0;

    890. 

  890.    RE[NumREs].E_olap = 0;

    891. 

  891.    RE[NumREs].E_len = 4;

    892. 

  892.    RE[NumREs].E_pal = 1;

    893. 

  893.    RE[NumREs].E_dgen = 16;

    894. 

  894.    RE[NumREs].E_mag = 4;

    895. 

  895.    RE[NumREs].proto = NumREs;

    896. 

  896.    RE[NumREs].Err = 0;

    897. 

  897.    RE[NumREs].E_nam_l = 3;

    898. 

  898. }

    899. 

  899. 

    900. 

  900. void Add_dcm(int NumREs) {

    901. 

  901.    strcpy(RE[NumREs].E_nam,"dcm\0");

    902. 

  902.    if ((RE[NumREs].E_raw_sit = calloc(6, sizeof(char))) == NULL) BadMem("Add_dcm-1", 1);

    903. 

  903.    RE[NumREs].E_raw_sit = "CCWGG\0"; 

    904. 

  904.    if ((RE[NumREs].E_wsit[1] = calloc(6, sizeof(char))) == NULL) BadMem("Add_dcm-2", 1);

    905. 

  905.    RE[NumREs].E_wsit[1] = "ccwgg\0"; 

    906. 

  906.    strcpy(RE[NumREs].E_hex[1],"ccwggn\0");

    907. 

  907.    RE[NumREs].E_tcut[1] = 0;

    908. 

  908.    RE[NumREs].E_olap = 0;

    909. 

  909.    RE[NumREs].E_len = 5;

    910. 

  910.    RE[NumREs].E_pal = 1;

    911. 

  911.    RE[NumREs].E_dgen = 8;

    912. 

  912.    RE[NumREs].E_mag = 4;

    913. 

  913.    RE[NumREs].proto = NumREs;

    914. 

  914.    RE[NumREs].Err = 0;

    915. 

  915.    RE[NumREs].E_nam_l = 3;

    916. 

  916. }

    917. 

  917. 

    918. 

  918. /***************************************************************************************************/

    919. 

  919. /* The following code performs the "pad the beginning sequence with the end and the end 

    920. 

  920.    with the beginning" regardless of whether it's the whole sequence or a subsequence, 

    921. 

  921.    using the vars 'begin'  and 'end'.  *Bare_Seq_Len is the the exactly "Real Sequence 

    922. 

  922.    Length" as determined from the SEQIO routines, no padding added or subtracted - 

    923. 

  923.    it's incremented to TotPaddedSeqLen in here and passed back to calling fn() via

    924. 

  924.    *Pad_Seq_Cnt  */

    925. 

  925. 

    926. 

  926. /* TotPaddedSeqLen    length of sequence + bracketing overlaps

    927. 

  927.    Bare_Seq_Len       actual number of bases read

    928. 

  928. */

    929. 

  929. /* will try to contruct this so that you give it the exact sequence you want padded, and it

    930. 

  930.    does all the calculations, work to alloc a new sequence space and do the padding, returning 

    931. 

  931.    exactly what you want */

    932. 

  932. 

    933. 

  933. char *PadEndForEnd(char *SeqIn, long *Pad_Seq_Len, long *Bare_Seq_Len) {

    934. 

  934.    long m, j, l, BareLen, TotPaddedSeqLen, FrPadSeqLen;

    935. 

  935. /*   int  f1=0, iseq; */

    936. 

  936.    char *sequence, b…
    937.

Large files files are truncated, but you can click here to view the full file