Code for VDJ plots


#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)