Hi James, I like your approach with dummy atoms and occupancy refinement. Dealing with actual maps sounds like hell to me indeed (especially given that we deal with weighted Fourier maps!). Reading this as someone who immediately translates this into a computer code (in my mind), a few things that aren't clear to me: - Where exactly inside the blob of density do you place these dummy atoms? - How many do you place? - Is there a risk of placing them too close to the boundary of the blob (in which case the question remains: what is the boundary?)? - I guess O or C as an atom type should do it, but what about B factor (would you refine B as well?)? - if you refine B, how do you deconvolute occupancy from refined B values (and eventually from effects of positional errors of your DAs)? -.... - How all these choices are going to affect the result?
All the best! Pavel On Mon, Aug 15, 2022 at 4:38 PM James Holton <jmhol...@lbl.gov> wrote: > There are several programs for integrating electron density, but please > let me assure you that it is almost always the wrong thing to do. > > A much better strategy is occupancy refinement. Throw in dummy atoms, > turn off non-bonded interactions to them, and refine their occupancy until > it a) stops changing (may be more than one round), and b) there are no > Fo-Fc differences left in the region of interest. Then all you do is add > up the occupancies, multiply by the relevant atomic number (usually 8), and > voila! you get the best-fit number of electrons in your blob. You may want > to try re-running with random starting points to get an idea of the error > bars. > > What is wrong with integrating density? Well, for one, it is hard to know > where to set the boundaries. Integrated density can be VERY sensitive to > the choice of radius, making your arbitrary decision of which radius to use > a source of error. Too small and you miss stuff. Too big and you add > unnecessary noise. Also, neighboring atoms have tails, and if you don't > subtract them properly, that is another source of error. Also, because of > the missing F000 term, there is an offset, which adds a term proportional > to the integration volume. For example, an integral resulting in zero > "electrons" does NOT mean you have vacuum. It just means that the area you > integrated has the same average density as the entire map. This may not be > the number you want. > > The beauty of occupancy refinement is that it automatically handles all > these problems. The "vacuum level" and F000 are known quantities in the > calculated map. The B factors given to the dummy atoms als o allow the > borders of your integration region to be "soft": down-weighting the > contribution of map noise far from your region of interest. And, finally, > by including atoms in the green density, neighboring atoms won't be sucked > into it. > > Think of it as fitting a smooth curve to noisy data and the number of > electrons is just a parameter in that fit, rather than trying to integrate > the noisy data itself. This is not an analogy. Refinement programs are > really just very sophisticated curve-fitting programs. And if you have a > forest of overlapping peaks and you are trying to de-convolute the > area/volume of one peak, it is best to include that peak in the fit, rather > than leave it out. Shoulder peaks especially tend to get "eaten" by large > neighboring peaks. > > How do you turn off non-bonds? Well, there is documentation for refmac: > http://www.ysbl.york.ac.uk/refmac/data/refmac_keywords.html > and phenix: > https://phenix-online.org/documentation/reference/refinement.html > > All that said, to answer the original question: > One very easy thing to do within the CCP4 suite is to use "mapmask" to > make a mask corresponding to your "sphere", or other region of interest. > Perhaps place a water at the center of your peak, and either use the > "border" feature of mapmask, or use "sfall" to compute a calculated map and > convert that to a mask using "threshold" in mapmask. This mask should have > values of 0 or 1 at every voxel. (or, if you feel like being clever, > something between 0 and 1 to reflect how much you want to weight a given > voxel). You can check it in coot. If you then multiply this mask by your > mFo-DFc map the result will have a non-zero average value. This will be > printed out in the log file. Multiply this average value by the volume of > the unit cell and you have your integrated number of electrons. Yes, its > that simple. > One issue you may have is map parameter compatibility (grid spacing, axis > order, xyz limits, etc.). You get around these by using the same grid in > all your fft or sfall runs, and then use mapmask to make the axis and > limits match before you multiply the map and mask. The only other issue > here might be the average value being a very small number and rounded off > by the default print precision. You can fix this by multiplying the map by > a large constant (again, using mapmask), then the printed value will have > lots of digits. > > This may seem complicated, but the use of masks can be a very valuable > skill to develop. In fact, one way to simplify, stabilize and accelerate > the occupancy refinement described above is to use a mask to isolate the > region of interest. That is, take the mFo-DFc map, zero out everything far > away from your peak, and convert the result to structure factors. You can > then call these structure factors "Fobs" (alongside the original > sigma(Fobs)) in a new refinement. The Rwork/Rfree then becomes a local > statistic, indicative of the % error in your refined total occupancy. One > caveat is that if every atom in the new refinement is having its occupancy > refined you will lose the absolute scale. To fix this, you need to add back > at least one well-ordered atom into "Fobs", and also include it in the > model. For example, take a well-ordered helix, extract those atoms, > calculate a map using "sfall", and add it to the masked-off difference map > before doing the "new Fobs" structure factor calculation. Include these > same atoms in the new refinement. They won't move, but they will keep the > scale stable. Oh, and don't forget to turn off the bulk solvent > correction! The bulk solvent has already been subtracted in the mFo-DFc > map. > > Hope this all makes sense and feel free to ask follow-up questions, > > -James Holton > MAD Scientist > > > On 8/10/2022 9:59 AM, Neno Vuksanovic wrote: > > Dear All, > > I would like to quantify electron density inside of positive Fo-Fc blobs > in active sites of multiple protomers in the map and compare them. I am > aware that I can interpolate maps and obtain density values at coordinate > points using either MapMan, Chimera or Coot, but I would like to know if > there is a tool that would let me designate a sphere of a certain volume > and calculate total electron density inside it? > > Best Regards, > Neno > > ------------------------------ > > To unsubscribe from the CCP4BB list, click the following link: > https://www.jiscmail.ac.uk/cgi-bin/WA-JISC.exe?SUBED1=CCP4BB&A=1 > > > > ------------------------------ > > To unsubscribe from the CCP4BB list, click the following link: > https://www.jiscmail.ac.uk/cgi-bin/WA-JISC.exe?SUBED1=CCP4BB&A=1 > ######################################################################## To unsubscribe from the CCP4BB list, click the following link: https://www.jiscmail.ac.uk/cgi-bin/WA-JISC.exe?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/
