swsoyee / r3dmol

🧬 An R package for visualizing molecular data in 3D
https://swsoyee.github.io/r3dmol/
Other
88 stars 4 forks source link

Color by b-factor (spectrum) #45

Closed dschust-r closed 2 years ago

dschust-r commented 2 years ago

Thank you for the great package, it's really nice to use and makes protein vis in R easy!

I was wondering if maybe you know a way to color structures by b-factor using r3dmol? I've been checking 3dmol but I don't think they have implemented that functionality yet. Unfortunately, I don't know any javascript and I am guessing this would have to be used as a colorfunction in the m_style_cartoon() function.

Do you maybe have any suggestions for this?

swsoyee commented 2 years ago

@dschust-r Thank you for using {r3dmol}. You can using some color function like this example:

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = "
        function(atom) {
          return atom.resi % 2 == 0 ? 'white' : 'green';
        }"
  )) %>%
  m_zoom_to()

something like:

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = "
        function(atom) {
          if (<THIS atom IS A B-FACTOR>) {
            return '<SOME COLOR NAME>';
          }
        }"
  )) %>%
  m_zoom_to()

You need to find a way to distinguish which atom belongs to the b-factor (e.g. know the positions of all the b-factors), then just set a color name or code (hex) for the position.

I am very sorry I don't know what b-factor is, so I don't know if this answer will help you to reach your goal.

@BradyAJohnston Do you have any better advice?

dschust-r commented 2 years ago

Thank you @swsoyee for the super fast reply!

The b-factor is contained as a column in the structure file (pdb, cif, etc) and it can more or less tell you how flexible these regions are. AlphaFold models contain the confidence scores in the b-factor column, so you could color them by that in e.g. pyMOL or ChimeraX.

In a PDB file it would be this column: image

swsoyee commented 2 years ago

@dschust-r Thanks for the explanation. Here is a quick solution.

b_factor_position <- 1:30 # Find the location of all your b-factors and save them as a variable

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = paste0("
        function(atom) {
          return [", paste0(b_factor_position, collapse = ","), "].includes(atom.resi) ? 'red' : 'white';
        }"
  ))) %>%
  m_zoom_to() %>%
  m_png()

image

dschust-r commented 2 years ago

Thank you for the reply, that's not quite what I meant. Sorry, I should have clarified.

The b-factor is a continuous number, assigned to each position. So you would have different colors on different regions throughout all of the protein sequence.

For example an Alphafold prediction would look like this (low b-factor orange with a gradient to high b-factor in blue) image

BradyAJohnston commented 2 years ago

I hadn't thought to look into colouring by b-factor, but you are right that the support for it isn't obvious in 3Dmol.js.

I have had a play around, and it's possible! We don't have a dedicated function for it (but can certainly make one) but you can pass a list to anything that has the colorScheme argument (such as m_style_surface() and m_style_ballnstick()).

r3dmol() %>% 
  m_add_model(data = pdb_1j72) %>% 
  m_add_surface(
    type = "SAS",
    style = m_style_surface(
      colorScheme = list(prop = "b", gradient = "roygb", min = 10, max = 200)
    )
  ) %>%
  m_zoom_to()

image

r3dmol() %>% 
  m_add_model(data = pdb_1j72) %>% 
  m_add_style(
    m_style_ballnstick(
      colorScheme = list(prop = "b", gradient = "roygb", min = 10, max = 200)
    )
  ) %>% 
  m_zoom_to()

image

In these examples, I am using one of the pre-defined colorschemes that already exist within 3Dmol.js, though they aren't well documented. image

I'm sure we could probably implement user-defined gradients a bit better, but @swsoyee would know more about the javascript side of that than I would (I'm the structural biologist, he's the proper javascript developer).

I don't quite know how to implement the gradients for the cartoon representations, I am showing an example below, but @swsoyee will need to help us out with being able to define a gradient inside of the function.

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = "
        function(atom) {
          return atom.b < 80 ? 'red' : 'white';
        }"
  )) %>%
  m_zoom_to()

image

BradyAJohnston commented 2 years ago

For clarification as well @swsoyee, for the cartoon javascript function, the b-factors would range from 0 to 100 and we would like to select a colour from a vector based on that. I know how I would do it in R,

colour_gradient_vector[round(protein$b)]

But I don't know an easy way to define a vector of hex codes or something similar in javascript.

dschust-r commented 2 years ago

@BradyAJohnston thank you so much! That already helps a lot.

I think the b-factors can even go higher than 100 for "real" structures. For AF predictions they go from 0-100 - as far as I know.

swsoyee commented 2 years ago

@BradyAJohnston Thank you for your detailed additions!

@dschust-r You could define the color for each atom according their score or something, just like this:

# Set color for each atom
color_vector <- paste0(c(rep('"blue"', 30), rep('"red"', 20), rep('"white"', 40), rep('"yellow"', 10)), collapse = ",")

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = paste0("
        function(atom) {
          const color = [", color_vector, "];
          return color[atom.resi];
        }"
  ))) %>%
  m_zoom_to() %>%
  m_png()

image

Or, If you know the what props are in atom, you could do this directly through Javascript functions Tips: all props of atom are list here: https://3dmol.csb.pitt.edu/doc/types.html#AtomSpec

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = "
        function(atom) {
          if (atom.b < 50) {return 'red'};
          if (atom.b < 80 && atom.b >= 50) {return 'green'};
          return 'white';
        }"
  )) %>%
  m_zoom_to() %>%
  m_png()

image

BradyAJohnston commented 2 years ago

Okay so it turns out that colorscheme is also supported for cartoons, but it isn't documented at all inside of 3Dmol.js so we haven't explicitly added support for it inside of m_style_cartoon() function (yet!).

A workaround in the meantime however is as follows (and this shouldn't break with future updates / fixes either):

The result of m_style_cartoon() however is just a list, so you can add in the colorscheme (all lowercase for the javascript-facing stuff) and define the gradient there as in the other examples:

library(r3dmol)
cartoon_styles <- m_style_cartoon()
cartoon_styles$cartoon$colorscheme <- list(prop = "b", gradient = "roygb", min = 10, max = 200)

r3dmol() %>% 
  m_add_model(pdb_1j72) %>% 
  m_add_style(style = cartoon_styles) %>% 
  m_zoom_to()

image

dschust-r commented 2 years ago

Thank you both @swsoyee and @BradyAJohnston that is really helpful and works like a charm for my data!

BradyAJohnston commented 2 years ago

Based on @swsoyee's response as well, you can define your own colour gradient in R as below, and use that in the colorfunc argument.

make_colours <- function(x, start = "red", end = "green") {
  colorfun <- colorRampPalette(c(start, end))

  paste0('"', paste(colorfun(x), collapse = '", "'), '"')
}

r3dmol() %>%
  m_add_model(data = pdb_1j72) %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = paste0("
        function(atom) {
          const color = [", make_colours(120, "blue", "tomato"),"]
          return color[Math.round(atom.b)]
        }")
  )) %>%
  m_zoom_to()

image

Thanks for opening the issue, we'll be sure to add in proper support for the colour gradients :)

dschust-r commented 2 years ago

That's amazing! Thank you so much, works great 🥳

So happy that it works!

swsoyee commented 2 years ago

It looks like this problem has been solved, so I'll close it for now. Feel free to reopen it again if needed.