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 <bshaa...@bgu.ac.il>
*Sent:* Thursday, May 28, 2020 11:59
*To:* b...@hofkristallamt.org
*Cc:* 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&A=1
------------------------------------------------------------------------
To unsubscribe from the CCP4BB list, click the following link:
https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&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/
