Another brilliant Super-Elf, with a blast from the past - Numerical Recipes FFT rules 😊
Many thanks for the excellent explanations, Cheers, BR From: James Holton <jmhol...@lbl.gov> Sent: Saturday, June 13, 2020 17:34 To: b...@hofkristallamt.org; CCP4BB@JISCMAIL.AC.UK Subject: Re: [ccp4bb] visual mask editor - why Bernhard, Sounds like you are plotting something similar to what I was tinkering with once. A script you may find useful is this one: https://bl831.als.lbl.gov/~jamesh/bin_stuff/map_func.com I wrote this because although mapmask, maprot, etc have very useful functionalities I found I wanted additional features, such as dividing one map by another, or taking a square root. These are important if you are trying to derive the "signal-to-noise ratio", for example. Once you have a map of rho/sigma(rho) you can convert that into a "pobability" by passing it through the "erf()" and "pow()" functions. This can be a good way of estimating the "probability something is there" or P(rho) for a given map voxel. Specifically: P(rho) = 1-pow(erf(abs(rho/sigma(rho))/sqrt(2)),V/2/d^3) where: rho is the electron density map value (preferably from a Fo-Fc map) sigma(rho) is the error on that voxel V is the unit cell volume (A^3) d is the resolution in A erf() is called the "error function" pow(m,e) is the raise-to-a-power function: m^e The erf() function by itself turns a rho/sigma(rho)=3 peak into 0.997, and a 1-sigma peak into 0.683. What that means is: assuming the noise is Gaussian, you expect voxels in the range -1-sigma to +1-sigma to be ~68.3% of the total. The "probability it isn't noise" (sometimes called a "p-value") is then 1-0.683 = 0.32. Seems like a pretty high probability to give to a 1-sigma peak, but now remember that the map is not just a single observation but thousands. So, the question you really want to ask is: if I generate 100x100x100 = 1e6 Gaussian-random numbers, what are the odds that a 4-sigma peak occurred at random? The answer is: pretty much garanteed. In any collection of 1 million Gaussian-random numbers with rms=1 it is virtually impossible to not have at least one of them > 4. Trust me, I have tried. This is where the "pow()" function comes in. You need to multiply all the individual voxel probabilities together to get the probability of at least one >4-sigma peak happening at random. But, then again, map voxels are hardly independent observations. Finite resolution means that neighboring pixels are highly correlated. So, rather than map grid points, we should be considering "Shannon voxels". All this is is the number of blobs of diameter d, where d is the resolution, that can fit into the volume of the map. For example, if we have a 100 A edge on a cubic cell and 3 A resolution, then we have about 3.7e4 independent "observations" of density, so the probability of a random 4-sigma peak is: P(rho) = erf(4/sqrt(2))**(((100./3)**3)/2) = 0.31 That is, if you make a zillion maps of random data using different seeds each time, 3.7e4 voxels each, and draw from a Gaussian distribution with rms=1, you expect 31% of these maps to have a 4-sigma peak. Randomly. The other 69% will not have anything > 4. So, if you see a 4-sigma peak in a map with 3.7e4 Shannon voxels I'd say it is real about 69% of the time. You might consider 0.69 to be a decent "weight" you should give such a 4-sigma peak. A 5-sigma peak under the same circumstances gets P(rho) = 0.989, and a 3-sigma peak gets P(rho) = 1e-22. aka: probaby noise. It is perhaps worth remembering that at a given resolution large-sigma noise peaks are more common in bigger cells than small ones. So, how do you get sigma(rho)? To be honest, the rms value of the mFo-DFc map is a pretty good estimate. The rms value of the 2mFo-DFc map is not (Lang et al. PNAS 2013). I usually get sigma(rho) empirically. That is, by making a stack of maps: start with your favorite refined structure and introduce random noise from whatever source you want to test. I.E. Gaussian error proprotional to SigI is an obvious one. A more realistic on is the Fo-Fc difference itself. After adding this "extra" noise to the data, re-refine the structure and generate a 2mFo-DFc map. Do this ~50 times with different random number seeds. Then take those 50 maps and compute the mean and standard deviation of "rho" at every voxel. You can do this with mapmask's ADD, MULT and SCALE features, but you can't do the last step, which is taking the square root of the variance. Hence: map_func.com There are lots of other functions supported, including random number genration, etc. Run the script with no arguments to get a list. Oh, but don't try it on an mtz file! mtz files are not maps. -James Holton MAD Scientist On 5/28/2020 12:11 PM, Bernhard Rupp wrote: Yes I have already pilfered useful parts of it in the scripts… Thx, BR From: Boaz Shaanan <mailto:bshaa...@bgu.ac.il> <bshaa...@bgu.ac.il> Sent: Thursday, May 28, 2020 11:59 To: b...@hofkristallamt.org <mailto:b...@hofkristallamt.org> Cc: CCP4BB@JISCMAIL.AC.UK <mailto:CCP4BB@JISCMAIL.AC.UK> Subject: Re: [ccp4bb] visual mask editor - why Hi Bernhard, Did you consider trying 'polder' in the phenix package? Boaz Boaz Shaanan, Ph.D. Department of Life Sciences Ben Gurion University of the Negev Beer Sheva Israel On May 28, 2020 21:17, Bernhard Rupp <hofkristall...@gmail.com <mailto:hofkristall...@gmail.com> > wrote: Maybe I should explain an example: Say coot detects an unmodelled blob (maybe a ligand). Now, I would like to do a number of things without biasing towards a model. Like comparing these map regions, excluding intrusion of a solvent mask, etc. Now could coot for example just generate a mask around what it already knows are blobs? Possible useful items could be a solvent mask not including that regions, or a density map that includes only features with a certain boundary around that blob. I pilfered some kludges together from different sources, but let’s just say inelegant would be a compliment. Best, BR ------------------------ Brief question: Does something like a visual density mask editor exist? Thx, BR ------------------------------------------------------ Bernhard Rupp http://www.hofkristallamt.org/ b...@hofkristallamt.org <mailto:b...@hofkristallamt.org> +1 925 209 7429 +43 676 571 0536 ------------------------------------------------------ Many plausible ideas vanish at the presence of thought ------------------------------------------------------ _____ To unsubscribe from the CCP4BB list, click the following link: https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB <https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1> &A=1 _____ To unsubscribe from the CCP4BB list, click the following link: https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB <https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1> &A=1 ######################################################################## To unsubscribe from the CCP4BB list, click the following link: https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1 This message was issued to members of www.jiscmail.ac.uk/CCP4BB, a mailing list hosted by www.jiscmail.ac.uk, terms & conditions are available at https://www.jiscmail.ac.uk/policyandsecurity/
