1 /* @source embmol *************************************************************
2 **
3 ** Routines for molecular weight matching.
4 **
5 ** @author Copyright (c) 1999 Alan Bleasby
6 ** @version $Revision: 1.23 $
7 ** @modified $Date: 2012/07/14 14:52:40 $ by $Author: rice $
8 ** @@
9 **
10 ** This library is free software; you can redistribute it and/or
11 ** modify it under the terms of the GNU Lesser General Public
12 ** License as published by the Free Software Foundation; either
13 ** version 2.1 of the License, or (at your option) any later version.
14 **
15 ** This library is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 ** Lesser General Public License for more details.
19 **
20 ** You should have received a copy of the GNU Lesser General Public
21 ** License along with this library; if not, write to the Free Software
22 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
23 ** MA 02110-1301, USA.
24 **
25 ******************************************************************************/
26
27
28 #include "ajlib.h"
29
30 #include "embmol.h"
31 #include "embprop.h"
32
33 #include "ajarr.h"
34
35 #include <string.h>
36 #include <stdlib.h>
37 #include <limits.h>
38 #include <stdio.h>
39 #include <math.h>
40
41
42
43 static ajint embMolFragSort(const void* a, const void* b);
44
45
46
47
48 /* @func embMolGetFrags *******************************************************
49 **
50 ** Create a sorted list of molwt fragments
51 **
52 ** @param [r] thys [const AjPStr] sequence
53 ** @param [r] rno [ajint] 1=Trypsin 2=LysC 3=ArgC 4=AspN 5=V8b 6=V8p
54 ** 7=Chy 8=CNBr
55 ** @param [r] mwdata [EmbPPropMolwt const *] molecular weight data
56 ** @param [r] mono [AjBool] true for monoisotopic data
57 ** @param [w] l [AjPList*] list for results
58 ** @return [ajint] number of fragments
59 **
60 ** @release 1.13.0
61 ******************************************************************************/
62
embMolGetFrags(const AjPStr thys,ajint rno,EmbPPropMolwt const * mwdata,AjBool mono,AjPList * l)63 ajint embMolGetFrags(const AjPStr thys, ajint rno, EmbPPropMolwt const *mwdata,
64 AjBool mono, AjPList *l)
65 {
66 static struct enz
67 {
68 const char *ename;
69 const char *residues;
70 const char *type;
71 const char *partial;
72 } zyme[]=
73 {
74 {"Trypsin","KR","CC","KRIFL"},
75 {"Lys-C","K","C",""},
76 {"Arg-C","R","C",""},
77 {"Asp-N","D","N",""},
78 {"V8b","E","C","KR"},
79 {"V8p","DE","CC",""},
80 {"Chymotrypsin","FYWLM","CCCCC",""},
81 {"CNBr","M","C",""}
82 };
83
84 EmbPMolFrag frag = NULL;
85 EmbPMolFrag *ptr = NULL;
86
87 ajint len;
88 ajint pos;
89 const char *p;
90
91 static AjPInt defcut =NULL;
92 ajint defcnt;
93
94 ajint beg;
95 ajint end;
96 ajint i;
97 double mw;
98
99
100 if(!defcut)
101 defcut=ajIntNew();
102
103 --rno;
104
105 len = ajStrGetLen(thys);
106 p = ajStrGetPtr(thys);
107
108 defcnt=0;
109
110 /* Positions of complete digest cuts */
111 for(pos=0;pos<len;++pos)
112 {
113 if(!strchr(zyme[rno].residues,(ajint)p[pos]))
114 continue;
115
116 if(len==pos+1)
117 continue;
118
119 if(p[pos+1]=='P' && rno!=3 && rno!=7)
120 continue;
121
122 if(rno==4 && p[pos+1]=='E')
123 continue;
124
125 ajIntPut(&defcut,defcnt++,pos);
126 }
127
128
129 /* Molwts of definite cuts */
130 beg = 0;
131 for(i=0;i<defcnt;++i)
132 {
133 end = ajIntGet(defcut,i);
134
135 if(strchr(zyme[rno].type,(ajint)'N'))
136 --end;
137
138 mw = embPropCalcMolwt(p,beg,end, mwdata, mono);
139
140 if(rno==7)
141 mw -= (double)(17.0079 + 31.095);
142
143 AJNEW0(frag);
144 frag->begin = beg+1;
145 frag->end = end+1;
146 frag->mwt = mw;
147 ajListPush(*l,(void *)frag);
148 beg = end+1;
149 }
150
151 if(defcnt)
152 {
153 mw = embPropCalcMolwt(p,beg,len-1,mwdata,mono);
154
155 if(rno==7)
156 mw -= (double)(17.0079 + 31.095);
157
158 AJNEW0(frag);
159 frag->begin = beg+1;
160 frag->end = len;
161 frag->mwt = mw;
162 ajListPush(*l,(void *)frag);
163 }
164
165 /* Overlaps */
166 if(defcnt)
167 {
168 ajListReverse(*l);
169 ajListToarray(*l,(void ***)&ptr);
170
171 for(i=0;i<defcnt-1;++i)
172 {
173 beg = ptr[i]->begin;
174 end = ptr[i+1]->end;
175 AJNEW0(frag);
176 frag->begin = beg;
177 frag->end = end;
178 mw = embPropCalcMolwt(p,beg-1,end-1,mwdata,mono);
179 frag->mwt = mw + EMBMOLPARDISP;
180 ajListPush(*l,(void *)frag);
181 }
182
183 AJFREE(ptr);
184 }
185
186
187 ajListSort(*l, &embMolFragSort);
188 ajIntDel(&defcut);
189
190 return (ajuint) ajListGetLength(*l);
191 }
192
193
194
195
196 /* @funcstatic embMolFragSort *************************************************
197 **
198 ** Sort routine for molwt fragments
199 **
200 ** @param [r] a [const void*] EmbPMolFrag pointer
201 ** @param [r] b [const void*] EmbPMolFrag pointer
202 **
203 ** @return [ajint] molwt difference
204 **
205 ** @release 1.5.0
206 ******************************************************************************/
207
embMolFragSort(const void * a,const void * b)208 static ajint embMolFragSort(const void* a, const void* b)
209 {
210 return (ajint)((*(EmbPMolFrag const *)a)->mwt -
211 (*(EmbPMolFrag const *)b)->mwt);
212 }
213