9 IAV Kegg Pathway Candidate Gene Analysis

9.1 Libraries

library(limma)
library(org.Hs.eg.db)
library(AnnotationDbi) 
library(tidyverse)
library(ggkegg)
library(tidygraph)

9.2 Data

t_postfilter <- read.table("Output Files/txome_female_edgeR_postfilter.txt")
  
t_acute <- read.csv("Output Files/txome_female_edgeR_DEG_acute.csv", row.names=1)
t_peak <- read.csv("Output Files/txome_female_edgeR_DEG_peak.csv", row.names=1)
t_late <- read.csv("Output Files/txome_female_edgeR_DEG_late.csv", row.names=1)

g_postfilter <- read.table("Output Files/gnome_female_edgeR_postfilter.txt")
g_acute <- read.csv("Output Files/gnome_female_edgeR_DEG_acute.csv", row.names=1)
g_peak <- read.csv("Output Files/gnome_female_edgeR_DEG_peak.csv", row.names=1)
g_late <- read.csv("Output Files/gnome_female_edgeR_DEG_late.csv", row.names=1)

t_virus_postfilter <- read.table("Output Files/txome_female_edgeR_viralRNA_postfilter.txt")
t_virus_degs <- read.csv("Output Files/txome_female_edgeR_viralRNA_DEG.csv", row.names=1)

g_virus_postfilter <- read.table("Output Files/gnome_female_edgeR_viralRNA_postfilter.txt")
g_virus_degs <- read.csv("Output Files/gnome_female_edgeR_viralRNA_DEG.csv", row.names=1)

9.3 Extract genes from KEGG Pathway

iav <- 
  getGeneKEGGLinks(species="hsa") %>% 
  filter(PathwayID == "hsa05164")

iav$Symbol <- mapIds(org.Hs.eg.db, iav$GeneID,
                     column="SYMBOL", keytype="ENTREZID")

9.4 Overlap with Gene Lists

9.4.1 Postfilter

postfilter <- 
  rbind(t_postfilter, g_postfilter, t_virus_postfilter) %>% 
  filter(!duplicated(V1))

iav_postfilter <-
  intersect(iav$Symbol, postfilter$V1)

9.4.2 Differentially Expressed Genes

acute <- 
  data.frame(
    genes=c(intersect(iav$Symbol, t_acute$gene), 
            intersect(iav$Symbol, g_acute$gene)), 
    category="acute")

peak <- 
  data.frame(
    genes=c(intersect(iav$Symbol, t_peak$gene), 
            intersect(iav$Symbol, g_peak$gene)), 
    category="peak")

late <- 
  data.frame(
    genes=c(intersect(iav$Symbol, t_late$gene),
            intersect(iav$Symbol, g_late$gene)),
    category="late") %>% 
  filter(!duplicated(genes))

virus <- 
  data.frame(
    genes=c(intersect(iav$Symbol, t_virus_degs$gene),
            intersect(iav$Symbol, g_virus_degs$gene)),
    category="virus") %>% 
  filter(!duplicated(genes))

9.4.3 Compile gene IDs

deg <-
  rbind(acute, peak, late, virus) %>% 
  mutate(id = paste0("hsa:", ifelse(genes %in% iav$Symbol, 
                     iav[match(genes, iav$Symbol),1],
                     genes)))

expressed <-
  data.frame(genes=iav_postfilter) %>%  
  filter(!genes %in% deg$genes) %>% 
  mutate(id = paste0("hsa:", ifelse(genes %in% iav$Symbol, 
                     iav[match(genes, iav$Symbol),1],
                     genes)))

notexpressed <- 
  data.frame(iav) %>% 
  filter(!Symbol %in% expressed$genes) %>% 
  filter(!Symbol %in% deg$genes) %>% 
  mutate(id = paste0("hsa:", ifelse(Symbol %in% iav$Symbol, 
                     iav[match(Symbol, iav$Symbol),1],
                     Symbol)))

9.5 Plot pathway map

** no longer used - cannot figure out coloring for some genes. Created map online.

graph_info <-
  pathway("hsa05164") %>% 
  activate(nodes) %>% 
  mutate(convert_hsa=convert_id("hsa"), 
         convert_map=convert_id("pathway"),
         bgcolor=ifelse(convert_hsa %in% expressed$id, "#ffffe0", 
                        ifelse(convert_hsa %in% deg$id, "#fbc69d", "white")))

graph <-
  ggraph(graph_info, x=x, y=y) +
   geom_node_rect(aes(filter=type=="gene", 
                      fill=I(bgcolor)),
                 color="black") + 
  overlay_raw_map() +
  theme_void() 
graph 

#ggkeggsave(filename="Figures/iavkegg.jpg", graph, dpi=300)

9.5.1 Create Legend

png(filename="Figures/iavkeggleg.png", width = 6, height=7, units="in", res=300)
  plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', , xlim=c(0,100), ylim=c(0,10))
  legend("topleft", 
         legend =c("Acute", "Peak", "Late", "Virus", "Expressed", "Not Expressed"), 
         pch=15, pt.cex=7, cex=3, bty='o',
         col = c("darkslategray1", "darkslategray3", "darkslategray4", "gray80", "coral", "wheat2"), 
         horiz=F)
dev.off()