Sponsored by the AIRR Community

Code for VDJ plots

Protocols
#1

As requested by some folks at the AIRR meeting last week I am happy to provide the R code for the VDJ plots I have presented in my talk. A fancier way of making these plots (such as Fig. S2 in our J Immunol paper: http://www.jimmunol.org/content/jimmunol/194/1/252.full.pdf?with-ds=yes) needs some longer code which I can also provide if people are interested.

Load libraries

x <- c(“tidyverse”, “alakazam”, “ggthemes”)
sapply(x, require, character.only = T)
rm(x)

Download example files from below URL:

http://changeo.readthedocs.io/en/version-0.3.9---defineclones-brownbag/examples/igblast.html#example-data

Load .tab data into R

out <- readChangeoDb("~/Downloads/Changeo_Example/S43_db-pass_parse-select.tab")

Modify V/D/J_CALL

out$V_GENE <- as.character(substr(lapply(strsplit(as.character(out$V_CALL), “\"), “[”, 1), 8, 20))
out$D_GENE <- as.character(substr(lapply(strsplit(as.character(out$D_CALL), "\”), “[”, 1), 8, 20))
out$J_GENE <- as.character(substr(lapply(strsplit(as.character(out$J_CALL), “\*”), “[”, 1), 8, 20))

Add Isotype information from CPRIMER and make a factor

out$ISOTYPE <- out$CPRIMER %>% gsub("-PCR", “”, .) %>% factor(levels = c(“IgD”, “IgM”, “IgA”, “IgG”, “IgE”))

VDJ plot

a <- out %>% group_by(ISOTYPE, MID, V_GENE, D_GENE, J_GENE) %>% summarise(Frequency = sum(DUPCOUNT))

a[0, ] %>%
ggplot(aes(V_GENE, D_GENE, group = J_GENE, color = J_GENE, size = Frequency)) +
theme_few(8) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5), strip.text = element_text(size = 12, face = “bold”)) +
geom_point(data = a %>% filter(J_GENE == “IGHJ1”), position = position_nudge(x = -.2, y = .1), alpha = .5, shape = 16) +
geom_point(data = a %>% filter(J_GENE == “IGHJ2”), position = position_nudge(x = 0, y = .1), alpha = .5, shape = 16) +
geom_point(data = a %>% filter(J_GENE == “IGHJ3”), position = position_nudge(x = .2, y = .1), alpha = .5, shape = 16) +
geom_point(data = a %>% filter(J_GENE == “IGHJ4”), position = position_nudge(x = -.2, y = -.1), alpha = .5, shape = 16) +
geom_point(data = a %>% filter(J_GENE == “IGHJ5”), position = position_nudge(x = 0, y = -.1), alpha = .5, shape = 16) +
geom_point(data = a %>% filter(J_GENE == “IGHJ6”), position = position_nudge(x = .2, y = -.1), alpha = .5, shape = 16) +
facet_grid(MID ~ ISOTYPE)

7 Likes
#2

Hi jtrueck，
The download URL that provided by you is out of date.And I’m interested in the plots in [quote=“jtrueck, post:1, topic:494”]
A fancier way of making these plots (such as Fig. S2 in our J Immunol paper: http://www.jimmunol.org/content/jimmunol/194/1/252.full.pdf?with-ds=yes ) needs some longer code which I can also provide if people are interested.
[/quote]

Could you please help me provide the plot code?

Many thanks,
Arida

#3

Thanks for your interest. Here is the updated URL for the example files: http://clip.med.yale.edu/immcantation/examples/Changeo_Example.tar.gz

I’ll also provide the code for the “fancier” plot in due course.

#4

Thank you so much！！！

#5

Ok here is the complete updated code for both plots as there were some changes made to the example data and the original code no longer worked. To achieve the “glowing” effect of the more fancy plot, I am using a loop to graph the same points 30 times with differing sizes and transparency values. So if you have a lot of points, this can take a while or even crash the system. There may be a better / easier way to achieve the same thing - if anyone discovers this, please let me know :slight_smile:

Load libraries

x <- c("tidyverse", "alakazam", "ggthemes")
sapply(x, require, character.only = T)

Download example files from below URL, unzip and load .tab data into R

http://clip.med.yale.edu/immcantation/examples/Changeo_Example.tar.gz

out <- readChangeoDb("~/Downloads/Changeo_Example/S43_db-pass_parse-select.tab")

Modify V/D/J_CALL & add Isotype information from CPRIMER as factor

out$V_GENE <- as.character(substr(lapply(strsplit(as.character(out$V_CALL), "\\*"), "[", 1), 1, 20))
out$D_GENE <- as.character(substr(lapply(strsplit(as.character(out$D_CALL), "\\*"), "[", 1), 1, 20))
out$J_GENE <- as.character(substr(lapply(strsplit(as.character(out$J_CALL), "\\*"), "[", 1), 1, 20))
out$ISOTYPE <- out$CPRIMER %>% gsub("-PCR", "", .) %>% factor(levels = c("IgD", "IgM", "IgA", "IgG", "IgE"))

VDJ plot

a <- out %>% group_by(ISOTYPE, MID, V_GENE, D_GENE, J_GENE) %>% summarise(Frequency = sum(DUPCOUNT))

a[0, ] %>% 
  ggplot(aes(V_GENE, D_GENE, group = J_GENE, color = J_GENE, size = Frequency)) + 
  theme_few(8) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5), strip.text = element_text(size = 12, face = "bold")) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ1"), position = position_nudge(x = -.2, y = .1), alpha = .5, shape = 16) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ2"), position = position_nudge(x = 0, y = .1), alpha = .5, shape = 16) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ3"), position = position_nudge(x = .2, y = .1), alpha = .5, shape = 16) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ4"), position = position_nudge(x = -.2, y = -.1), alpha = .5, shape = 16) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ5"), position = position_nudge(x = 0, y = -.1), alpha = .5, shape = 16) + 
  geom_point(data = a %>% filter(J_GENE == "IGHJ6"), position = position_nudge(x = .2, y = -.1), alpha = .5, shape = 16) +
  facet_grid(MID ~ ISOTYPE)

More fancy version of VDJ plot with glowing appearance of points

create J nudging data.frame

J.nudge <- data.frame(J_GENE = paste0("IGHJ", 1:6), J.x = rep(c(-.2, 0, .2), 2), J.y = rep(c(.1, -.1), each = 3), color = c("#F8766D", "#B79F00", "#00BA38", "#00BFC4", "#619CFF", "#F564E3"), stringsAsFactors = F)

setting “frame” for the VDJ plot

VDJplot <- a[0, ] %>% 
  ggplot(aes(V_GENE, D_GENE, group = J_GENE, color = J_GENE)) + theme_few(8) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5), strip.text = element_text(size = 12, face = "bold")) + 
  geom_point(data = a, alpha = 0) + 
  facet_grid(MID ~ ISOTYPE) + guides(size = F, color = guide_legend(override.aes = list(alpha = 1))) + scale_size(range = c(0, 10)) # scale_size adjusts the size of the points

add the points to the VDJ plot (loop)

for (i in J.nudge$J_GENE){ 
  VDJplot <- VDJplot + 
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*30), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/30, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*29), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/29, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*28), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/28, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*27), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/27, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*26), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/26, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*25), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/25, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*24), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/24, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*23), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/23, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*22), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/22, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*21), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/21, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*20), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/20, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*19), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/19, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*18), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/18, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*17), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/17, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*16), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/16, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*15), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/15, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*14), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/14, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*13), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/13, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*12), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/12, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*11), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/11, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*10), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/10, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*9), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/9, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*8), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/8, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*7), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/7, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*6), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/6, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*5), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/5, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*4), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/4, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*3), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/3, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) +
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency*2), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1/2, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i]) + 
    geom_point(data = a %>% filter(J_GENE == i), aes(size = Frequency), position = position_nudge(x = J.nudge$J.x[J.nudge$J_GENE == i], y = J.nudge$J.y[J.nudge$J_GENE == i]), alpha = 1, shape = 16, color = J.nudge$color[J.nudge$J_GENE == i])
}
print(VDJplot)
rm(x, a, out, VDJplot, J.nudge, i)