CREATE OR REPLACE PROCEDURE calc_ic50(pRunsetCode IN NUMBER) AS /* NOTE: disable trigger PRE_CALC_Q_AFTER_INSERT. */ /* pRunsetCode - parameters are the same for whole runset. */ /* pLow = fixed or initial value of the LOW parameter. */ /* pSlope = fixed or initial value of the SLOPE parameter. */ /* pIC50 = fixed or initial value of the IC50 parameter. */ /* pHigh = fixed or initial value of the HIGH parameter. */ /* pAlist1 = 1 (variable LOW parameter; if = 0, fixed) */ /* pAlist2 = 1 (variable SLOPE parameter; if = 0, fixed) */ /* pAlist3 = 1 (variable IC50 parameter; if = 0, fixed) */ /* pAlist4 = 1 (variable HIGH parameter; if = 0, fixed) */ /* pTolerance = percentage of the range; default is 5 = 5%; cannot be 0. */ /* pConvergence = ratio of decrease that counts as convergence, */ /* can be from 0.01 to 0.000001, default is 0.001 */ pLow NUMBER := NULL; pSlope NUMBER := NULL; pIC50 NUMBER := NULL; pHigh NUMBER := NULL; /* here the pAlist parameters refer to the Runset_Parameters table */ /* where 0 means the parameter is not fixed, and so is fitted. */ pAlist1 NUMBER := 0; pAlist2 NUMBER := 0; pAlist3 NUMBER := 0; pAlist4 NUMBER := 0; pTolerance NUMBER := 5; pConvergence NUMBER := 0.001; pIterLimit NUMBER := 1000; vIter NUMBER; vConcentration NUMBER; vResponse1 NUMBER; vStdDev NUMBER; vNPoints NUMBER; vCount1 NUMBER; vCount2 NUMBER; vMaxSeries NUMBER; /* used to calculate number of variable, fitted parameters. */ vNFit NUMBER := 0; vConvergence NUMBER; vALambda NUMBER; vChiSq NUMBER; vCount NUMBER; vIndex NUMBER; vIndex2 NUMBER; vMinDelta NUMBER; vNumRows NUMBER; vOldChiSq NUMBER; vResponseAtMinC NUMBER; vResponseAtMaxC NUMBER; vTemp NUMBER; vTemp2 NUMBER; vTemp3 NUMBER; vA IC50_ENGINE.NumTabType; vConc IC50_ENGINE.NumTabType; vLista IC50_ENGINE.NumTabType; vResponse IC50_ENGINE.NumTabType; vResults IC50_ENGINE.NumTabType; vSig IC50_ENGINE.NumTabType; vAlfa IC50_ENGINE.MultiDNumTabType; vCovar IC50_ENGINE.MultiDNumTabType; /* for IC50, the total number of parameters is 4. */ sMa CONSTANT NUMBER := 4; /* select each row of AMPOINT table created for the runset, valid or not. */ /* must contain a sample (control series have reference samples). */ CURSOR cDataSet IS SELECT * FROM ampoint WHERE runset_code = pRunsetCode AND sample_code IS NOT NULL; /* the next line is in case you want to exclude reference samples */ /* AND reference_flag IS NULL or reference_flag <> 'T'); */ /* select each concentration in the series for an mpoint code. */ /* base the series on the first replicate series. */ /* (all replicate series will have the same number of concentrations */ /* and the same dilutions). */ CURSOR cDilutions (pMpointCode IN NUMBER) IS SELECT * FROM awell WHERE mpoint_code = pMpointCode AND replicant_number = 1; /* pick all fields for all rows about an mpoint_code from PRE_CALC_DATA */ /* and order them. */ CURSOR cPreCalcData(pMpointCode IN NUMBER) IS SELECT * FROM pre_calc_data WHERE mpoint_code = pMpointCode ORDER BY concentration; BEGIN /* check the Runset_Parameters table for any parameter values. */ SELECT COUNT(RUNSET_CODE) INTO vCount2 from RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; IF vCount2 = 1 THEN SELECT ACTUAL_IC50_HIGH INTO pHigh FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT ACTUAL_IC50_LOW INTO pLow FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT ACTUAL_IC50_SLOPE INTO pSlope FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT ACTUAL_IC50_IC50 INTO pIC50 FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT FIX_ACTUAL_IC50_HIGH INTO pAlist4 FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT FIX_ACTUAL_IC50_LOW INTO pAlist1 FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT FIX_ACTUAL_IC50_SLOPE INTO pAlist2 FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; SELECT FIX_ACTUAL_IC50_IC50 INTO pAlist3 FROM RUNSET_PARAMETERS WHERE RUNSET_CODE = pRunsetCode; END IF; /* Runset_Parameters table marks fixed parameters as 1, but the algorithm */ /* marks fitted parameters as 1, so the pAlist values must be reversed */ /* since from now on the pAlist parameters refer to the algorithm. */ vTemp := pAlist1; pAlist1 := 0; IF vTemp = 0 OR vTemp IS NULL THEN pAlist1 := 1; END IF; vTemp := pAlist2; pAlist2 := 0; IF vTemp = 0 OR vTemp IS NULL THEN pAlist2 := 1; END IF; vTemp := pAlist3; pAlist3 := 0; IF vTemp = 0 OR vTemp IS NULL THEN pAlist3 := 1; END IF; vTemp := pAlist4; pAlist4 := 0; IF vTemp = 0 OR vTemp IS NULL THEN pAlist4 := 1; END IF; /* find the number of variable parameters. */ IF pAlist1 = 1 THEN vNFit := vNFit + 1; END IF; IF pAlist2 = 1 THEN vNFit := vNFit + 1; END IF; IF pAlist3 = 1 THEN vNFit := vNFit + 1; END IF; IF pAlist4 = 1 THEN vNFit := vNFit + 1; END IF; /* select each concentration series with a sample code, in the runset. */ <> FOR vDataSet IN cDataSet LOOP /* set values in AMPOINT to NULL. */ UPDATE ampoint SET fit1_value = NULL, fit2_value = NULL, fit3_value = NULL, fit4_value = NULL, data_value = NULL, n_points = NULL, stat1_value = NULL, stat2_value = NULL WHERE mpoint_code = vDataSet.mpoint_code; /* for other SCREEN validity purposes, enter a calc_value into each awell */ UPDATE AWELL SET calc_value = raw_value WHERE mpoint_code = vDataSet.mpoint_code; /* find the valid LOW and HIGH data values for the mpoint code. */ SELECT MAX(raw_value), MIN(raw_value) /* you could work with calc_values instead, if you calculate them first */ /* SELECT MAX(calc_value), MIN(calc_value) */ INTO vMaxSeries, vCount1 FROM awell WHERE mpoint_code = vDataSet.mpoint_code AND bogus_flag IS NULL; /* set number of concentration positions to 0. */ vNPoints := 0; /* check each concentration in the series, for the sample. */ <> FOR vDilutions IN cDilutions(vDataSet.mpoint_code) LOOP /* find average conc and raw value over all replicates at a concentration,*/ /* and the raw value standard deviation there, */ /* for current series and position in series, */ /* for valid data points. */ /* do not use any negative or 0 concentrations. */ /* assumes averaging of all data values at a concentration. */ /* (to use all valid data points independently, add the AND statement: */ /* AND replicant_number = vDilutions.replicant_number */ /* before "AND bogus_flag IS NULL;". */ vConcentration := 0; vResponse1 := 0; vStddev := 0; SELECT AVG(concentration), AVG(raw_value), STDDEV(raw_value) /* you could work with calc_values instead, if you calculate them first */ /* SELECT AVG(concentration), AVG(calc_value), STDDEV(calc_value) */ INTO vConcentration, vResponse1, vStddev FROM AWELL WHERE mpoint_code = vDilutions.mpoint_code AND series_order = vDilutions.series_order AND series_number = vDilutions.series_number AND concentration > 0 AND bogus_flag IS NULL; /* check if any valid points were at this concentration. */ /* if no data points, skip to next concentration. */ IF (vConcentration > 0 AND vConcentration IS NOT NULL) THEN /* standard deviation of actual values is not good for measurement error */ /* but it is commented out here in case a user wants to use it. */ /* use entered error pTolerance (from actual use or machine instructions).*/ /* if no value is entered SCREEN will use a value = 5% of range of values.*/ /* Also different errors for different concentrations or */ /* responses could be coded here. */ /* IF (vStddev IS NULL OR vStddev = 0) THEN */ /* vStddev := 0.05 * (vMaxSeries - vCount1); */ /* END IF; */ vStddev := pTolerance * (vMaxSeries - vCount1) / 100; IF pTolerance > 100 THEN vStddev := 5 * (vMaxSeries - vCount1) / 100; END IF; IF (pTolerance IS NULL OR pTolerance <= 0) THEN vStddev := 5 * (vMaxSeries - vCount1) / 100; END IF; /* add current values to a data row in temp table. */ INSERT INTO PRE_CALC_DATA (MPOINT_CODE, CONCENTRATION, RESPONSE, STD_DEV) VALUES (vDilutions.mpoint_code, vConcentration, vResponse1, vStdDev); /* increase number of concentration positions by 1. */ vNPoints := vNPoints + 1; END IF; END LOOP dilutions_loop; /* make sure number of concentration positions is more than the */ /* number of variable parameters, so curve can be defined. */ IF vNPoints > vNFit THEN /* insert the current mpoint code row into a row of the temp table. */ INSERT INTO PRE_CALC_QUEUE (MPOINT_CODE, N_POINTS, LISTA_1, LISTA_2, LISTA_3, LISTA_4, N_FIT, LOW_VAL, SLOPE, IC50, HIGH_VAL) VALUES (vDataSet.mpoint_code, vNPoints, pAlist1, pAlist2, pAlist3, pAlist4, vNFit, pLow, pSlope, pIC50, pHigh ); /* the TRIGGER Pre_Calc_Q_After_Insert must be disabled in the database, */ /* using "ALTER TRIGGER Pre_Calc_Q_After_Insert DISABLE;" */ /* before using this calculation. */ /* Call the standard IC50 algorithm. */ /* set vLista(1), vLista(2), vLista(3), vLista(4). */ vLista(1) := pAlist1; vLista(2) := pAlist2; vLista(3) := pAlist3; vLista(4) := pAlist4; /* use the given HIGH activity value for the runset. */ vA(4) := pHigh; /* If no initial value was given, find an initial HIGH value for mpoint. */ /* if the entered low value is higher than the entered high value, sub. */ IF (pHigh IS NULL OR pHigh < pLow) THEN SELECT MAX(response) INTO vA(4) FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code; END IF; /* use the given LOW activity value for the runset. */ vA(1) := pLow; /* If no initial value was given, find an initial LOW value for mpoint. */ /* if the entered low value is higher than the entered high value, sub. */ IF (pLow IS NULL OR pHigh < pLow) THEN SELECT MIN(response) INTO vA(1) FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code; END IF; /* find maximum amd minimum concentration. */ /* initialize IC50 to the geometric mean of the two, if not given. */ SELECT MAX(concentration) INTO vIndex FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code; SELECT MIN(concentration) INTO vMinDelta FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code; vA(3) := pIC50; IF (pIC50 IS NULL OR pIC50 <= 0) THEN vA(3) := SQRT(vIndex * vMinDelta); END IF; /* initialize SLOPE, if not given. */ vA(2) := pSlope; /* first find midpoint of responses. */ vTemp := (vA(4) + vA(1)) / 2; /* find responses at maximum and minimum concentrations. */ SELECT response INTO vResponseAtMaxC FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code AND concentration = vIndex; SELECT response INTO vResponseAtMinC FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code AND concentration = vMinDelta; /* find trend of slope, ascending or descending. */ /* SLOPE is not the actual slope of curve but a parameter of the slope. */ /* it varies between -10 and +10. */ /* if origin is at 0 response and 0 concentration: */ /* flat line slope= 0,descending slope SLOPE> 0,ascending slope SLOPE < 0.*/ /* a steeper curve has a greater slope in absolute value. */ IF (pSlope IS NULL OR pSlope = 0 OR pSlope > 10 OR pSlope < -10 OR (pSlope > -0.01 AND pSlope < 0.01)) THEN vA(2) := 1; IF vResponseAtMaxC > vResponseAtMinC THEN vA(2) := - vA(2); END IF; END IF; /* initialize index. */ vIndex := 1; /* record activity value, concentration and expected error for points */ /* with values in PRE_CALC_DATA, in order of concentration, for mpoint. */ FOR cPreCalcDataRec IN cPreCalcData(vDataSet.mpoint_code) LOOP vConc(vIndex) := cPreCalcDataRec.concentration; vResponse(vIndex) := cPreCalcDataRec.response; vSig(vIndex) := cPreCalcDataRec.std_dev; vIndex := vIndex + 1; END LOOP; /* record number of rows in cPreCalcDataRec. */ vNumRows := vIndex - 1; /* set lambda to -1 initially. */ vALambda := -1.0; /* use MRQMIN to find the lambda, chi square, */ /* LOW, SLOPE, IC50, HIGH, and pAlist1, pAlist2, pAlist3, and pAlist4, */ /* as well as the covariance matrix COVAR and alpha matrix ALFA. */ /* MRQMIN calls MRQCOF to create the chi square, and to create the */ /* ALFA matrix (inverse matrix) and BETA matrix (solution vector). */ /* MRQCOF calls the IC50 function to use its mathematical function,*/ /* and the derivatives with respect to the parameters, */ /* and the actual data points to find the chi square. */ /* The IC50 function is the standard one to give a sigmoid curve. */ /* MRQMIN calls GAUSSJ to create an inverse matrix from a original matrix */ /* and return both matrices. */ /* GAUSSJ uses the MAP2D and SWAP functions, which are not part of */ /* the algorithm, just technical ways to make SQL code work */ /* MRQMIN calls COVSRT to create the COVAR covariant matrix. */ /* COVAR uses the MAP2D function. */ vTemp := mrqmin(vConc, vResponse, vSig, vNumRows, vA, vLista, sMa, vCovar, vAlfa, vChiSq, vALambda); /* check if IC50 < =0 or HIGH < LOW or SLOPE is way off. */ IF (vA(3) <= 0 OR vA(4) < vA(1) OR vA(2) > 10 OR vA(2) < -10 OR (vA(2) > -0.01 AND vA(2) < 0.01)) THEN vTemp := 1; END IF; /* if MRQMIN = 0, valid. */ IF vTemp = 0 THEN /* begin regular test. */ vIndex := 0; vTemp2 := 0; vIter := 1; /* record the current chi square value. */ vOldChiSq := vChiSq; WHILE (vIndex = 0 AND vTemp2 = 0 AND vIter < pIterLimit) LOOP /* test. */ vTemp2 := mrqmin(vConc, vResponse, vSig, vNumRows, vA, vLista, sMa, vCovar, vAlfa, vChiSq, vALambda); vIter := vIter + 1; /* check if IC50 < =0 or HIGH < LOW or SLOPE is way off. */ IF (vA(3) <= 0 OR vA(4) < vA(1) OR vA(2) > 10 OR vA(2) < -10 OR (vA(2) > -0.01 AND vA(2) < 0.01)) THEN vTemp2 := 1; END IF; /* MRQMIN sets new chi square to the old value if new is more than old. */ /* loop more if new chi square is more than old. */ /* go on if new is less than old only by less than the convergence ratio. */ /* pConvergence can range from 0.01 to 0.000001. */ /* vOldChiSq cannot be 0. */ vConvergence := pConvergence; IF pConvergence < 0.000001 THEN vConvergence := 0.000001; END IF; IF pConvergence > 0.01 THEN vConvergence := 0.01; END IF; IF vChiSq < vOldChiSq THEN IF vChiSq / vOldChiSq > 1 - vConvergence THEN vIndex := 1; END IF; vOldChiSq := vChiSq; END IF; END LOOP; /* continue if calculation is still valid. */ IF (vTemp2 = 0 AND vIter < 1000) THEN /* when new chi square is a little less than old, set lambda = 0. */ vALambda := 0.0; /* test. */ vTemp3 := mrqmin(vConc, vResponse, vSig, vNumRows, vA, vLista, sMa, vCovar, vAlfa, vChiSq, vALambda); /* If it was valid before, it will still be valid. */ UPDATE ampoint SET fit1_value = vA(1), fit2_value = vA(4), fit3_value = vA(2), /* fit4_value = , */ data_value = vA(3), n_points = vNumRows, stat1_value = vChiSq, stat2_value = vALambda WHERE mpoint_code = vDataSet.mpoint_code; END IF; END IF; END IF; /* remove all rows from PRE_CALC_DATA. */ DELETE FROM pre_calc_data WHERE mpoint_code = vDataSet.mpoint_code; /* remove all rows from PRE_CALC_QUEUE. */ DELETE FROM pre_calc_queue WHERE mpoint_code = vDataSet.mpoint_code; END LOOP dataset_loop; COMMIT; END; / SHOW ERRORS PROCEDURE CALC_IC50