Export utils

antigenomics · Jun 3, 2018 · a79dad3 · a79dad3
1 parent e9ca778
commit a79dad3
Show file tree

Hide file tree

Showing 2 changed files with 95 additions and 16 deletions.
diff --git a/misc/misc.Rmd b/misc/misc.Rmd
@@ -19,3 +19,78 @@ fread("../database/vdjdb.slim.txt") %>%
   get_segment_parts %>%
   fwrite("../res/segments.aaparts.txt", sep = "\t")
 ```
+
+Export TCRs by chain & epitope for human
+
+```{r}
+fread("../database/vdjdb.slim.txt") %>%
+  filter(species == "HomoSapiens") %>%
+  proc_slim_vdjdb %>%
+  as.vdjtools.df -> dt.vdjdb
+
+epi.count <- dt.vdjdb %>%
+  group_by(gene, antigen.epitope) %>%
+  summarise(count = n())
+
+epi.good <- epi.count %>%
+  filter(count >= 30)
+
+system("mkdir export/")
+for (gg in c("TRA", "TRB")) {
+  for (ee in epi.good %>% filter(gene == gg) %>% .$antigen.epitope) {
+    dt.vdjdb %>% 
+      filter(gene == gg, antigen.epitope == ee) %>%
+      select(-species, -gene, -antigen.epitope) %>%
+      fwrite(paste0("export/", gg, ".", ee, ".txt"), sep = "\t")
+  }
+}
+```
+```
+dt.vdjdb.hb <- dt.vdjdb %>%
+  filter(species == "HomoSapiens", cdr3.beta != "", 
+         v.beta != "", j.beta != "") %>%
+  select(cdr3.beta, v.beta, j.beta, antigen.epitope) %>%
+  unique
+
+dt.epi.count.b <- dt.vdjdb.hb %>%
+  group_by(antigen.epitope) %>%
+  summarise(count = n())
+
+epi.sel.b <- dt.epi.count.b %>% 
+  filter(count >= 30) %>% 
+  .$antigen.epitope %>%
+  unique
+
+dt.vdjdb.hb.1 <- dt.vdjdb.hb %>%
+  filter(antigen.epitope %in% epi.sel.b)
+
+dt.vdjdb.ha <- dt.vdjdb %>%
+  filter(species == "HomoSapiens", cdr3.alpha != "",
+         v.alpha != "", j.alpha != "") %>%
+  select(cdr3.alpha, v.alpha, j.alpha, antigen.epitope) %>%
+  unique
+
+dt.epi.count.a <- dt.vdjdb.ha %>%
+  group_by(antigen.epitope) %>%
+  summarise(count = n())
+
+epi.sel.a <- dt.epi.count.a %>% 
+  filter(count >= 30) %>% 
+  .$antigen.epitope %>%
+  unique
+
+dt.vdjdb.ha.1 <- dt.vdjdb.ha %>%
+  filter(antigen.epitope %in% epi.sel.a)
+
+system("mkdir export/")
+for (epi in epi.sel.b) {
+  dt.vdjtools.b %>% 
+    filter(antigen.epitope == epi) %>%
+    fwrite(paste0("export/TRB.", epi, ".txt"), sep = "\t")
+}
+for (epi in epi.sel.a) {
+  dt.vdjtools.a %>% 
+    filter(antigen.epitope == epi) %>%
+    fwrite(paste0("export/TRA.", epi, ".txt"), sep = "\t")
+}
+```
diff --git a/misc/utils.R b/misc/utils.R
@@ -2,6 +2,17 @@ library(dplyr)
 library(stringr)
 library(data.table)
 
+### VDJdb processing
+
+proc_slim_vdjdb <- function(.df) {
+  .df %>%
+    filter(v.segm != "", j.segm != "") %>%
+    mutate(v.segm = str_split_fixed(v.segm, "[*,]", 2)[,1],
+           j.segm = str_split_fixed(j.segm, "[*,]", 2)[,1]) %>%
+    select(species, gene, antigen.epitope, cdr3, v.segm, j.segm) %>%
+    unique
+}
+
 ### Segment processing
 get_segment_parts <- function(.df) {
   .df.v <- .df %>%
@@ -26,8 +37,6 @@ get_segment_parts <- function(.df) {
 
   rbind(.df.v, .df.j)
 }
-
-
 
 ### VDJtools export
 
@@ -47,23 +56,18 @@ mock_back_translate <- function(x) {
 
 # "CASS" %>% strsplit('') %>% lapply(mock_back_translate)
 
-as.vdjtools.df <- function(.df, .chain = c("beta", "alpha")) {
-  if (.chain == "beta") {
-    .df$cdr3aa <- .df$cdr3.beta
-    .df$v <- .df$v.beta
-    .df$j <- .df$j.beta
-  } else {
-    .df$cdr3aa <- .df$cdr3.alpha
-    .df$v <- .df$v.alpha
-    .df$j <- .df$j.alpha
-  }
-
-  .df$cdr3nt <- cdr3.beta %>% 
+as.vdjtools.df <- function(.df) {
+  .df$cdr3aa <- .df$cdr3
+  .df$cdr3nt <- .df$cdr3aa %>% 
     strsplit('') %>% 
     lapply(mock_back_translate)
 
   .df %>%
-    mutate(count = 1, freq = 1 / n(), d = "",
+    mutate(count = 1, v = v.segm, d = "", j = j.segm,
            vend = -1, dstart = -1, dend = -1, jstart = -1) %>%
-    select(count, freq, cdr3nt, cdr3aa, v, d, j, vend, dstart, dend, jstart)
+    group_by(species, gene, antigen.epitope) %>%
+    mutate(freq = 1 / n()) %>%
+    ungroup %>%
+    select(count, freq, cdr3nt, cdr3aa, v, d, j, vend, dstart, dend, jstart,
+           species, gene, antigen.epitope)
 }