add fusion summary for LGAT

analyses/fusion-summary/01-fusion-summary.Rmd

@@ -1,8 +1,8 @@
 ---
 title: "Generate Fusion Summary Files"
 output: html_notebook
-author: Daniel Miller (D3b) and Jaclyn Taroni (CCDL)
-date: January 2020
+author: Daniel Miller (D3b), Jaclyn Taroni (CCDL), Jo Lynne Rokita (D3b)
+date: January 2020, November 2020
 ---
 
 Generate fusion files specifically for consumption by molecular subtyping analyses
@@ -44,6 +44,30 @@ filterFusion <- function(df, bioid, fuses, genes) {
   return(df %>% select(Sample, FusionName))
 }
 
+#create function specifically for lgat, since we want more information from the putative oncogenic file
+filterLGATFusion <- function(df, bioid, fuses, genes) {
+  if (!missing(bioid)) {
+    df <- filter(df, Sample %in% bioid)
+  }
+  if (!missing(fuses) & !missing(genes)) {
+    df <- filter(df, FusionName %in% fuses |
+                   Gene1A %in% genes |
+                   Gene2A %in% genes |
+                   Gene1B %in% genes |
+                   Gene2B %in% genes)
+  } else if (!missing(fuses)) {
+    df <- filter(df, FusionName %in% fuses)
+  } else if (!missing(genes)) {
+    df <- filter(df,
+                 Gene1A %in% genes |
+                   Gene2A %in% genes |
+                   Gene1B %in% genes |
+                   Gene2B %in% genes)
+  }
+  return(df %>% 
+           select(Sample, FusionName, Fusion_Type, Gene1A, Gene1B, Gene2A,  Gene2B, Gene1A_anno, Gene1B_anno, LeftBreakpoint, RightBreakpoint,      reciprocal_exists, DomainRetainedGene1A, DomainRetainedGene1B) %>%
+           distinct())
+}
 
 #' Generate matrix with fusion counts
 #' @param fuseDF Filtered fusion data frame
@@ -67,10 +91,14 @@ prepareOutput <- function(fuseDF, bioid) {
 
 ```{r}
 dataDir <- file.path("..", "..", "data")
+fusDir <- file.path("..", "..", "analyses", "fusion_filtering", "results")
+annotDir <- file.path("..", "..", "analyses", "fusion_filtering", "references")
 #' The putative oncogenic fusion file is what we'll use to check for the 
+#' annotation file
+annot <- read.delim2(file.path(annotDir, "genelistreference.txt"), sep = "\t", header = T)
 #' presence or absence of the fusions.
 putativeOncogenicDF <- 
-  read_tsv(file.path(dataDir, "pbta-fusion-putative-oncogenic.tsv"))
+  read_tsv(file.path(fusDir, "pbta-fusion-putative-oncogenic.tsv"))
 #' However, some biospecimens are not represented in this filtered, prioritized
 #' file but *are* present in the original files -- this will cause them to be
 #' "missing" in the final files for consumption which could mislead analysts.
@@ -88,11 +116,12 @@ if (!dir.exists(resultsDir)) {
 ependFile <- file.path(resultsDir, "fusion_summary_ependymoma_foi.tsv")
 embryFile <- file.path(resultsDir, "fusion_summary_embryonal_foi.tsv")
 ewingsFile <- file.path(resultsDir, "fusion_summary_ewings_foi.tsv")
+lgatFile <- file.path(resultsDir, "fusion_summary_lgat_foi.tsv")
 ```
 
 ## Fusions and genes of interest
 
-Taken from [`AlexsLemonade/OpenPBTA-analysis#245`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/245), [`AlexsLemonade/OpenPBTA-analysis#251`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/251), and [`AlexsLemonade/OpenPBTA-analysis#623`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/623) respectively. 
+Taken from [`AlexsLemonade/OpenPBTA-analysis#245`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/245), [`AlexsLemonade/OpenPBTA-analysis#251`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/251), [`AlexsLemonade/OpenPBTA-analysis#623`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/623) respectively, and [`AlexsLemonade/OpenPBTA-analysis#808`](https://github.com/AlexsLemonade/OpenPBTA-analysis/issues/808) 
 
 ```{r}
 #' **Filters**
@@ -132,6 +161,28 @@ ewingsFuses<- c(
   "FUS--ERG"
 )
 
+#' 4: Exact match a list of fusions common in low-grade astrocytic (LGAT) tumors
+#' as well as fusions containing a particular gene with any other gene
+lgatFuses <- c(
+  "KIAA1549--BRAF",
+  "FGFR1--TACC1",
+  "MYB--QKI"
+)
+lgatGenes <- c(
+  "BRAF",
+  "ALK",
+  "ROS1",
+  "NTRK1",
+  "NTRK2",
+  "NTRK3",
+  "PDGFRA",
+  "FGFR2",
+  "FGFR1",
+  "MYB",
+  "MYBL1",
+  "RAF1"
+)
+
 ```
 
 ### Filter putative oncogenic fusions list
@@ -145,6 +196,9 @@ allFuseEmbry <- filterFusion(df = putativeOncogenicDF,
                              genes = embryGenes)
 allFuseEwing <- filterFusion(df = putativeOncogenicDF,
                              fuses = ewingsFuses)
+allFuseLGAT <- filterLGATFusion(df = putativeOncogenicDF,
+                             fuses = lgatFuses,
+                             genes = lgatGenes)
 
 ```
 
@@ -203,9 +257,104 @@ allFuseEwing %>%
   write_tsv(ewingsFile)
 ```
 
+#### Perform selection for LGAT fusions
+
+```{r}
+# Which genes/fusions are not kinases, but in the list?
+# Separate LGAT fusions into genes, combine with gene list, check for not kinase
+lgatFuses_df <- as.data.frame(lgatFuses) %>%
+  separate(lgatFuses, into = c("Gene1A", "Gene1B"), remove = F)
+kinases <- annot %>%
+  filter(type == "Kinase")
+lgatFuses_df$Gene1A_anno <- ifelse(lgatFuses_df$Gene1A %in% kinases$Gene_Symbol, "kinase", "")
+lgatFuses_df$Gene1B_anno <- ifelse(lgatFuses_df$Gene1B %in% kinases$Gene_Symbol, "kinase", "")
+# Only pull fusions that do not contain kinase genes, as ones with kinases will be dealt with sepaately later
+lgatFuses_df <- lgatFuses_df %>%
+  filter(lgatFuses_df$Gene1A_anno != "kinase" & lgatFuses_df$Gene1B_anno != "kinase") %>%
+  select(lgatFuses)
+nonkinase_lgatFuses <- lgatFuses_df$lgatFuses
+
+# Identify non-kinase genes in LGAT goi list
+nonkinase_lgatGenes <- setdiff(lgatGenes, kinases$Gene_Symbol)
+
+# Pull LGAT non-kinase fusions
+nonkinaseLGAT <- filterFusion(df = putativeOncogenicDF,
+                             fuses = nonkinase_lgatFuses,
+                             genes = nonkinase_lgatGenes) %>%
+  distinct()
+
+# First, filter all fusion dataframe for 3' kinases which are in-frame and retain the kinase domain - keep these
+three_prime_kinase_inframe <- allFuseLGAT %>%
+  filter(grepl("Kinase", Gene1B_anno) & Fusion_Type == "in-frame" & DomainRetainedGene1B == "Yes") %>%
+  select(Sample, FusionName, Gene1A, Gene1B) %>%
+  distinct()
+
+# Which fusions are not in-frame?
+three_prime_kinase_outframe <- allFuseLGAT %>%
+  filter(grepl("Kinase", Gene1B_anno) & Fusion_Type != "in-frame") %>%
+  distinct()
+# Check they don't retain the kinase domain. They do not, so leave out.
+table(three_prime_kinase_outframe$DomainRetainedGene1B)
+
+# Filter all fusions for 5' kinase fusions which have reciprocal fusions and have lost the kinase domain. 
+five_prime_domain_lost <- allFuseLGAT %>%
+  filter(grepl("Kinase", Gene1A_anno) & reciprocal_exists == "TRUE") %>%
+  select(Sample, FusionName, Gene1A, Gene1B, Fusion_Type, DomainRetainedGene1A) %>%
+  filter(DomainRetainedGene1A == "No") %>%
+  distinct()
+
+# Keep those with kinase domain retained and fusion in-frame - keep this list
+five_prime_domain_intact <- allFuseLGAT %>%
+  filter(grepl("Kinase", Gene1A_anno) & reciprocal_exists == "TRUE") %>%
+  select(Sample, FusionName, Gene1A, Gene1B, Fusion_Type, DomainRetainedGene1A) %>%
+  filter(DomainRetainedGene1A == "Yes") %>%
+  distinct()
+
+# Identify reciprocal matched fusions which are in-frame
+five_prime_domain_lost$reciprocal <- paste(five_prime_domain_lost$Gene1B, five_prime_domain_lost$Gene1A, sep ="--")
+five_prime_reciprocals <- five_prime_domain_lost[,c("Sample", "reciprocal")] %>%
+  rename(FusionName = reciprocal)
+five_prime_reciprocals_full <- five_prime_reciprocals %>%
+  left_join(allFuseLGAT, by = c("Sample", "FusionName")) %>%
+  filter(Fusion_Type == "in-frame")
+# Keep the list of the original 5' kinase fusions with matched reciprocals
+five_prime_kinase_keep <- five_prime_reciprocals_full[,c("Sample", "Gene1A", "Gene1B")] %>%
+  distinct()
+  five_prime_kinase_keep$FusionName <- paste(five_prime_kinase_keep$Gene1B, five_prime_kinase_keep$Gene1A, sep ="--")
+
+# Which ones are not in-frame? These are mostly genes fused to self
+out_frame_reciprocals <- setdiff(five_prime_reciprocals[,c("Sample", "FusionName")], five_prime_reciprocals_full[,c("Sample", "FusionName")])
+out_frame_reciprocals
+
+#Check these do not retain the kinase domain. They do not, so will leave out
+out_frame_reciprocals_full <- out_frame_reciprocals %>%
+  left_join(allFuseLGAT, by = c("Sample", "FusionName")) %>%
+  select(Sample, FusionName, Gene1A, Gene1B, Fusion_Type, DomainRetainedGene1A, DomainRetainedGene1B)
+
+# Rbind lists for final table of LGAT fusions of interest
+subsetFuseLGAT <- bind_rows(nonkinaseLGAT, three_prime_kinase_inframe[,c("Sample", "FusionName")],
+                         five_prime_domain_intact[,c("Sample", "FusionName")],  
+                         five_prime_kinase_keep[,c("Sample", "FusionName")]) %>%
+  distinct()
+```
+
+#### Write LGAT fusions to file
+
+```{r}
+subFuseLGAT <- subFuseLGAT %>%
+  prepareOutput(specimensUnion)
+
+# Which fusions of interest are not present in any samples?
+missingLgatFusion <- setdiff(lgatFuses, colnames(subFuseLGAT)[-1])
+# For the fusions that are not present, fill those columns with 0
+subFuseLGAT[, missingLgatFusion] <- 0
+
+subFuseLGAT %>%
+  write_tsv(lgatFile)
+```
+
 ## Session Info
 
 ```{r}
 sessionInfo()
 ```
-