-
Notifications
You must be signed in to change notification settings - Fork 29
Experimental Artifact Suppression Parameters
I've added some new options to try to reign in artifacts. I'll give the fast how-to first and then the long explanation later.
In the Additional Deconvolution Parameters controls in regular UniDec, there are now 6 additional parameters for Z trimming, 4 types of suppression and 2 control parameters. The basic idea is to constrain the allowed charge state assignments for each data point. UniDec by default allows each data point to be assigned to any number of allowed charges, which is why the core data structure is a matrix with m/z on one dimension and z on the other. These parameters all constrain the Z dimension. All 4 types of suppression can be used in tandem with each other.
-
Suppression Top N: The default is 0, which turns this parameter off. By setting the value to >0, it limits the number of non-zero charge states that are allowed per data point. For example, a setting of 3 would allow 3 possible charge states to have positive intensity at each data point. A setting of 1 is the most extreme and will limit each data point to only one allowed charge state assignment. Higher settings than 3 probably won't do much, but you can try them out anyway. A value larger than the total number of possible charge states will be waste computer time.
-
Suppression Top X: The default is 0, which turns this parameter off. Setting the value to >0 will remove charge states that are below a threshold intensity. For each data point, the maximum charge state is considered 1. Anything below X% (relative to the max charge state) will be set to 0. This operates only on each data point, ignoring neighboring data points. The idea is that minor charge state assignments for each peak will be zeroed out, which helps remove artifacts. Smaller values will give a more gentle effect. For example, 0.1 will zero out any charge states below 10% relative intensity, but it allows multiple charge state assignment for each data point, provided they have a reasonable intensity. 0.5 is more extreme, and 0.9 will be similar to the Top N of 1, where usually only 1 charge state is allowed per data point. The exception would be for two peaks cases where a data point is almost equally split between two charge states.
-
Suppression Satellites: The default is 0, which is off. Adding an integer greater than 0 will turn this on. Basically, it scans across a window of +/- n for each data point charge column. If it sees a local maximum, it leaves that alone and sets all others to 0. This allows local suppression of satellite artifacts that are within +/- n from a true charge state assignment. Startit applies but not the suppression percent, since it didn't seem helpful.
-
Suppression Harmonics: The default is 0, which is off. This is a boolean parameter, so setting to 1 will turn it on. The idea here is that for every charge state it can, it will look for it's double. If it finds a double, it will compare the two. If the double is higher, it will assume the lower one is a harmonic. If the double is lower, it will assume the higher one is the harmonic. Whichever is the harmonic is zeroed out. Startit applies but not the suppression percent, since it didn't seem helpful.
-
Suppression Percent: The default here is 0, which says to actually zero out the intensity for each m/z and z pair that is removed. I put in a softer cutoff option so that you could downplay it by multiplying by this value. In other words, you could set this to 0.1, in which case each intensity flagged by either the top X or top N above would be multiplied by 0.1 instead of multiplied by 0. Thus, it still allows intensity for those spots, just tries to downplay them and watch how the rest of the algorithm factors that in. It's a more Bayesian way of thinking. However, it doesn't seem to work. I'm not sure if it's technical or something in the algorithm, but even very low values that should be effectively the same as 0 give messy results. For now, leave it at 0, but let me know if you find any value in playing with this.
-
Suppression Start Iteration: The default here is 3, which says to only start applying this filtering after the algorithm has run 3 iterations. It gives the algorithm time to stabilize a bit, finding which charge states should be the highest naturally, before it starts to remove minor assignments. As opposed to the Cut Percent above, this seems to be pretty flexible. Lower values like 0 give the algorithm less time to settle before applying the cutoffs and can be more extreme. Larger values do less and may waste some computational time as it initially settles into one state only for you to then disrupt it and force it resettle. 3 seems to be a good balance in early testing, but there is a lot of leeway here.
This seems pretty simple and it is. A few people have suggested things like this to me over the years. You may ask why I never tried this. Actually, I did in the very beginning (2015 to be specific). However, it didn't work very well. Here are the key differences.
First, it was placed at the end of the iteration loop rather than the beginning. See the original loop below:
for(m=0;m<abs(config.numit);m++)
{
//if (config.suppression!=0 &&m>0&& m%2==1&&off==0){
// Suppression(newblur, blur, barr,nztab,starttab,endtab,isotopepos,dataMZ,lengthmz,isolength,config,0);
//}
if (config.zsig >= 0){
blur_it_mean(lengthmz,
config.numz,
numclose,
closemzind,
closezzind,
newblur,
blur,
barr);}
else{
blur_it(lengthmz,
config.numz,
numclose,
closemzind,
closezzind,
closeval,
newblur,
blur,
barr);
}
deconvolve_iteration_speedy(lengthmz, config.numz,maxlength,
newblur, blur,barr, config.aggressiveflag, dataInt,
config.damp,isolength, isotopepos, isotopeval, starttab, endtab, mzdist, config.speedyflag);
//Determine the metrics for conversion. Only do this every 10% to speed up.
if ((config.numit<10||m%10==0||m%10==1||m>0.9*config.numit)){
double diff=0;
double tot=0;
//#pragma omp parallel for private(i), schedule(dynamic), reduction(+:diff,tot)
for(i=0;i<lengthmz*config.numz;i++)
{
if(barr[i]==1)
{
diff += pow((blur[i]-oldblur[i]), 2);
tot += blur[i];
}
}
conv=(diff/tot);
printf("Iteration: %d Convergence: %f\n",m,conv);
if(conv<0.000001){
if(off==1&&config.numit>0){
if (config.suppression != 0){
Reconvolve(lengthmz, config.numz, maxlength, starttab, endtab, mzdist, blur, newblur, config.speedyflag);
memcpy(blur, newblur, lengthmz*config.numz*sizeof(double));
Suppression(newblur, blur, barr, nztab, starttab, endtab, isotopepos, dataMZ, lengthmz, isolength, config, 1);
config.suppression = 0;
m = 0;
off = 0;
}
else{
printf("\nConverged in %d iterations.\n\n",m);
break;
}
}
off=1;
}
memcpy(oldblur, blur, lengthmz*config.numz*sizeof(double));
if (m==config.numit-1&&config.suppression != 0){
Reconvolve(lengthmz, config.numz, maxlength, starttab, endtab, mzdist, blur, newblur, config.speedyflag);
memcpy(blur, newblur, lengthmz*config.numz*sizeof(double));
Suppression(newblur, blur, barr, nztab, starttab, endtab, isotopepos, dataMZ, lengthmz, isolength, config, 1);
config.suppression = 0;
m = 0;
off = 0;
}
}
}
Interestingly, you'll note that I tried it at the beginning of the loop, and it didn't seem to help. Also, when I tried moving the new code to the bottom of the loop to match this more closely, it gave fairly similar (only slightly worse) results. So, I don't think placement matters.
Another potential difference was that I had used much harsher smashing. If you notice below, the code for Suppression_Total is nearly identical to a Top N of 1. However, this not only zeroed out the intensity of that point, but it also flipped the boolean array (barr) to 0 for that m/z and z pair, which prevented that assignment from ever coming back again. Our new code doesn't do this, so the rest of the algorithm can override the suggestions of the artifact suppression if it wants.
Here, I think it was a substantial difference. When I switched the code to also flip barr, it made things a lot worse for the new code. Lesson learned: don't kill potential assignments altogether. Zeroing them out helps (somehow more than just lowering their intensity), but make sure they can come back. Interestingly, it looks like I also tried this back in the day. There was a killflag to allow the barr part to be disengaged. For some reason, that's hard coded to 1 above, so it must have been that we tried it. I know for a fact we tried it because it's 0 in the hidden code at the beginning of the loop. So, why didn't that work??? Read on below to find out...
void Suppression_Harmonic(double *blur, double *newblur, char *barr, int *ztab, int lengthmz, int numz,int killflag)
{
deepcopy(newblur, blur, lengthmz*numz);
int zmax = ztab[numz - 1];
int zmin = ztab[0];
int end = (zmax / 2)-zmin;
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < lengthmz; i++)
{
for (int j = 0; j < end; j++)
{
if(barr[index2D(numz, i, j)] == 1)
{
int j2 = 2 * j + zmin;
double val = blur[index2D(numz, i, j)];
double val2 = blur[index2D(numz, i, j2)];
if (val > val2){ newblur[index2D(numz, i, j2)] = 0; if (killflag == 1){ barr[index2D(numz, i, j2)] = 0.0; } }
if (val < val2){ newblur[index2D(numz, i, j)] = 0; if (killflag == 1){ barr[index2D(numz, i, j)] = 0.0; } }
}
}
}
memcpy(blur, newblur, lengthmz*numz*sizeof(double));
}
void Suppression_Satellite(double *blur, double *newblur, char *barr, int lengthmz, int numz,int killflag)
{
deepcopy(newblur, blur, lengthmz*numz);
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < lengthmz; i++)
{
for (int j = 2; j < numz; j++)
{
if (barr[index2D(numz, i, j)] == 1)
{
int j2 = j-1;
int j3 = j - 2;
double val = blur[index2D(numz, i, j)];
double val2 = blur[index2D(numz, i, j2)];
double val3 = blur[index2D(numz, i, j3)];
if (val > val2 && val>val3){
newblur[index2D(numz, i, j2)] = 0;
newblur[index2D(numz, i, j3)] = 0;
if (killflag == 1){ barr[index2D(numz, i, j2)] = 0.0; barr[index2D(numz, i, j3)] = 0.0; }
}
if (val2 > val && val2>val3){
newblur[index2D(numz, i, j)] = 0;
newblur[index2D(numz, i, j3)] = 0;
if (killflag == 1){ barr[index2D(numz, i, j)] = 0.0; barr[index2D(numz, i, j3)] = 0.0; }
}
if (val3 > val && val3>val2){
newblur[index2D(numz, i, j)] = 0;
newblur[index2D(numz, i, j2)] = 0;
if (killflag == 1){ barr[index2D(numz, i, j2)] = 0.0; barr[index2D(numz, i, j)] = 0.0; }
}
}
}
}
memcpy(blur, newblur, lengthmz*numz*sizeof(double));
}
void Suppression_Total(double *blur, double *newblur, char *barr,int lengthmz, int numz,int killflag)
{
deepcopy(newblur, blur, lengthmz*numz);
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < lengthmz; i++)
{
double maxval = 0;
for (int j = 0; j < numz; j++)
{
double val = blur[index2D(numz, i, j)];
if (val > maxval){ maxval = val; }
}
for (int j = 0; j < numz; j++)
{
if (blur[index2D(numz, i, j)] <maxval)
{
newblur[index2D(numz, i, j)] = 0;
if (killflag == 1){ barr[index2D(numz, i, j)] = 0.0; }
}
}
}
memcpy(blur, newblur, lengthmz*numz*sizeof(double));
}
void Suppression_Ultimate(double *blur, double *newblur,char *barr, int * starttab,int * endtab, double * dataMZ, int lengthmz, int numz,double fwhm,int killflag)
{
deepcopy(newblur, blur, lengthmz*numz);
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < lengthmz; i++)
{
double maxval = 0;
double mztop = dataMZ[i];
for (int k = starttab[i]; k <= endtab[i]; k++)
{
if (fabs(dataMZ[k] - mztop) <= fwhm)
{
for (int j = 0; j < numz; j++)
{
double val = blur[index2D(numz, k, j)];
if (val > maxval){ maxval = val; }
}
}
}
for (int k = starttab[i]; k <= endtab[i]; k++)
{
if (fabs(dataMZ[k] - mztop) <=fwhm)
{
for (int j = 0; j < numz; j++)
{
if (blur[index2D(numz, k, j)] < maxval)
{
newblur[index2D(numz, k, j)] = 0;
if (killflag == 1){ barr[index2D(numz, k, j)] = 0.0; }
}
}
}
}
}
memcpy(blur, newblur, lengthmz*numz*sizeof(double));
}
void Suppression_Isotope(double *blur, double *newblur, char * barr, int *isotopepos, int lengthmz, int numz, int isolength,int num,int killflag)
{
deepcopy(newblur, blur, lengthmz*numz);
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < lengthmz; i++)
{
double maxval = 0;
for (int j = 0; j < numz; j++)
{
if (barr[index2D(numz, i, j)] == 1)
{
int pos = isotopepos[index3D(numz, isolength, i, j, num)];
for (int k = i; k < pos; k++)
{
double val = blur[index2D(numz, k, j)];
if (val > maxval){ maxval = val; }
}
}
}
for (int j = 0; j < numz; j++)
{
if (barr[index2D(numz, i, j)] == 1)
{
int pos = isotopepos[index3D(numz, isolength, i, j, num)];
for (int k = i; k < pos; k++)
{
if (blur[index2D(numz, k, j)] < maxval){
newblur[index2D(numz, k, j)] = 0.0;
if (killflag == 1){ barr[index2D(numz, k, j)] = 0.0;}
}
}
}
}
}
memcpy(blur, newblur, lengthmz*numz*sizeof(double));
}
void Suppression(double *newblur, double *blur, char *barr, int * nztab, int * starttab, int * endtab, int * isotopepos, double * dataMZ, int lengthmz, int isolength, Config config,int killflag)
{
double width;
int number = 5;
if (killflag == 1){ printf("Suppressing\n"); }
switch (config.suppression){
case 1:
Suppression_Harmonic(newblur, blur, barr, nztab, lengthmz, config.numz,killflag);
break;
case 2:
Suppression_Satellite(newblur, blur, barr, lengthmz, config.numz,killflag);
break;
case 3:
Suppression_Harmonic(newblur, blur, barr, nztab, lengthmz, config.numz,killflag);
Suppression_Satellite(newblur, blur, barr, lengthmz, config.numz,killflag);
break;
case 4:
Suppression_Total(newblur, blur, barr,lengthmz, config.numz,killflag);
break;
case 5:
if (config.psfun == 0){ width = config.mzsig*2.35482/2.; }
else{ width = config.mzsig/2.; }
Suppression_Ultimate(newblur, blur, barr,starttab, endtab, dataMZ, lengthmz, config.numz, width,killflag);
break;
case 6:
if (isolength > number)
{
Suppression_Isotope(newblur, blur, barr, isotopepos, lengthmz, config.numz, isolength, number,killflag);
}
break;
}
}
In fact, that Suppression_Total code works great and is identical to the Top N when set to 1. What changed in the last 11 years to make this approach viable? The answer is in the Point Smoothing option that was added back in early 2019 with version 3.1. It turns out that the combination of peak smoothing with artifact suppression is the magic bullet to smooth out issues (pun intended) and make things look nice. This is definitely a case where the idea first released in 2015 was sound, but it needed some other pieces to fall into place before it worked well.
Interestingly, the point smoothing idea is pretty much directly taken from the original Bayesian algorithm, just out of order a bit, and is related to the Suppression Ultimate feature above, which tried to account for neighboring m/z points. However, the Suppression Ultimate tried to zero out neighboring points, leaving only one in the middle, whereas the point smoothing tries to make them all agree.
Let me know how it works for you. I'll probably streamline the parameters a bit, but I wanted to get something out there with maximal control for people to play with. If you've read this far, drop me an email because I want to I'm not sure anyone will ever do that :)