Figures for Contributions of RLRs to WNV innate immune signaling in BMMs

0a. Intro
1. Set Locations
2. Load Libraries and supplemental files
3. Create data structures needed for figure preparation
4. Heatmap of differentially expressed genes in WNV infection
5. Create Heatmaps of “wild-type” DE genes from Immune and Innate Immune Signaling GOBP categories
7. RLR-dependent WNV-induced transcriptional response
8. Identify RLR-specific gene signatures
9. Identify pathway and polarization profiles among gene sets

0a. Intro

Macrophages were derived from freshly isolated bone marrow from mice of various genotypes. These cells were then infected with WNV-Tx i.c. p3 at an MOI of 2.5 for 24 or 48 hours (plus a 48h mock). RNA was then isolated and sent to Expression Analysis, Inc. The targets file of the analysis script describes the replication strategy and layout of all samples. The next steps are differential expression analysis to determine gene signatures specifically induced (or lacking) in each knockout. This script creates figures of interest from the differential gene expression results.

1. Set Locations

# Input files 

in.expression <- "results/ExpressMatrix.txt"

in.dataMask <- "results/sig_dataMask.txt"



# Input files for pathway annotation

in.immune <- "raw_data/GO_categories/GO0006955_immune_response.csv"

in.innate <- "raw_data/GO_categories/GO0045087_innate_immune_response.csv"

in.isgs <- "raw_data/ISGs_BMM_AH.csv"

in.th1 <- "raw_data/GO_categories/GO0042088_Th1_type_immune_response.txt"

in.th2 <- "raw_data/GO_categories/GO0042092_Th2_type_immune_response.txt"

in.m1 <- "raw_data/GO_categories/M1_Becker.txt"

in.m2 <- "raw_data/GO_categories/M2_Becker.txt"

2. Load Libraries and supplemental files

# Install libraries if needed

#source("http://www.bioconductor.org/biocLite.R")

#biocLite(pkgs = c("biomaRt", "xlsx", "limma", "plotrix", "reshape2", "ggplot2"))



# Load libraries

library(biomaRt)

library(xlsx)

## Loading required package: rJava

## Loading required package: xlsxjars

library(limma)

library(plotrix)

library(reshape2)

library(ggplot2)



# Source in heatmap.F.4 function

source("../../MiscScripts/heatmap.F.R")

## 

## Attaching package: 'gplots'

## The following object is masked from 'package:plotrix':

## 

##     plotCI

## The following object is masked from 'package:stats':

## 

##     lowess

## Loading required package: fastcluster

## 

## Attaching package: 'fastcluster'

## The following object is masked from 'package:stats':

## 

##     hclust

##

## ==========================================================================

## *

## *  Package WGCNA 1.51 loaded.

## *

## *    Important note: It appears that your system supports multi-threading,

## *    but it is not enabled within WGCNA in R. 

## *    To allow multi-threading within WGCNA with all available cores, use 

## *

## *          allowWGCNAThreads()

## *

## *    within R. Use disableWGCNAThreads() to disable threading if necessary.

## *    Alternatively, set the following environment variable on your system:

## *

## *          ALLOW_WGCNA_THREADS=<number_of_processors>

## *

## *    for example 

## *

## *          ALLOW_WGCNA_THREADS=8

## *

## *    To set the environment variable in linux bash shell, type 

## *

## *           export ALLOW_WGCNA_THREADS=8

## *

## *     before running R. Other operating systems or shells will

## *     have a similar command to achieve the same aim.

## *

## ==========================================================================

## 

## Attaching package: 'WGCNA'

## The following object is masked from 'package:stats':

## 

##     cor

# Print out session data

sessionInfo()

## R version 3.3.2 (2016-10-31)

## Platform: x86_64-w64-mingw32/x64 (64-bit)

## Running under: Windows 7 x64 (build 7601) Service Pack 1

## 

## locale:

## [1] LC_COLLATE=English_United States.1252 

## [2] LC_CTYPE=English_United States.1252   

## [3] LC_MONETARY=English_United States.1252

## [4] LC_NUMERIC=C                          

## [5] LC_TIME=English_United States.1252    

## 

## attached base packages:

## [1] stats     graphics  grDevices utils     datasets  methods   base     

## 

## other attached packages:

##  [1] WGCNA_1.51            fastcluster_1.1.22    dynamicTreeCut_1.63-1

##  [4] gplots_3.0.1          ggplot2_2.2.1         reshape2_1.4.2       

##  [7] plotrix_3.6-4         limma_3.30.10         xlsx_0.5.7           

## [10] xlsxjars_0.6.1        rJava_0.9-8           biomaRt_2.30.0       

## 

## loaded via a namespace (and not attached):

##  [1] Rcpp_0.12.9           lattice_0.20-34       GO.db_3.4.0          

##  [4] gtools_3.5.0          assertthat_0.1        rprojroot_1.2        

##  [7] digest_0.6.12         foreach_1.4.3         plyr_1.8.4           

## [10] backports_1.0.5       acepack_1.4.1         stats4_3.3.2         

## [13] RSQLite_1.1-2         evaluate_0.10         lazyeval_0.2.0       

## [16] data.table_1.10.4     gdata_2.17.0          S4Vectors_0.12.1     

## [19] rpart_4.1-10          Matrix_1.2-8          preprocessCore_1.36.0

## [22] checkmate_1.8.2       rmarkdown_1.3         splines_3.3.2        

## [25] stringr_1.1.0         foreign_0.8-67        htmlwidgets_0.8      

## [28] RCurl_1.95-4.8        munsell_0.4.3         BiocGenerics_0.20.0  

## [31] base64enc_0.1-3       htmltools_0.3.5       nnet_7.3-12          

## [34] htmlTable_1.9         tibble_1.2            gridExtra_2.2.1      

## [37] Hmisc_4.0-2           IRanges_2.8.1         codetools_0.2-15     

## [40] matrixStats_0.51.0    XML_3.98-1.5          bitops_1.0-6         

## [43] grid_3.3.2            gtable_0.2.0          DBI_0.5-1            

## [46] magrittr_1.5          scales_0.4.1          KernSmooth_2.23-15   

## [49] impute_1.48.0         stringi_1.1.2         doParallel_1.0.10    

## [52] latticeExtra_0.6-28   Formula_1.2-1         RColorBrewer_1.1-2   

## [55] iterators_1.0.8       tools_3.3.2           Biobase_2.34.0       

## [58] parallel_3.3.2        survival_2.40-1       yaml_2.1.14          

## [61] AnnotationDbi_1.36.2  colorspace_1.3-2      cluster_2.0.5        

## [64] caTools_1.17.1        memoise_1.0.0         knitr_1.15.1

3. Create data structures needed for figure preparation

Read in differential expression matrix for union of differentially expressed genes

# DE = > 2 fold change up or down over mock, < 0.05 BH-adjusted p-val

expression <- as.matrix(read.table(in.expression, header = TRUE, row.names = 1))



# "masked" expression matrix: all non-significant genes set to 0

sigMask <- as.matrix(read.table(in.dataMask, header = TRUE, row.names = 1))



# Create DE results matrix from sigMask

DE_results <- apply(sigMask, MARGIN = c(1, 2), FUN = function(x) {

  if (x >= 1) {

    x = 1

  } else if (x <= -1) {

    x = -1

  } else {

    x = 0

  }

})

Annotate genes according to ENSEMBL ID using biomaRt package

# Extract annotation from ENSEMBL database

ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl", host = "www.ensembl.org")

gene_map <- getBM(attributes = c("external_gene_name", "ensembl_gene_id"), filters = "external_gene_name",

                  values = row.names(expression), mart = ensembl)

Create function to collect heatmap modules

# gene_list is a named vector of genes mapped to respective modules

# exprs_matrix is matrix used to create heatmap

collect_modules <- function(gene_list, exprs_matrix) {

  # Create filename for results (file named after expression matrix used to create heatmap)

  filename <- paste("results/modules/", deparse(substitute(exprs_matrix)), ".xlsx", sep = "")

  

  # Split vector of genes by module color

  module_list <- split(gene_list, as.factor(gene_list))

  

  # Collect the gene names mapping to each module color (vector names)

  modules <- lapply(module_list, function(x) { names(x) })



  # Collect relevant data for each module

  module_expression <- lapply(modules, function(x) {

    # collect annotation and expression data

    annotation_data <- gene_map[(match(x, gene_map$external_gene_name)), ]

    express_data <- exprs_matrix[row.names(exprs_matrix) %in% x, ]

    

    # reorder expression data to match annotation

    express_data <- express_data[(match(x, row.names(express_data))), ]

    cbind.data.frame(geneName = x, annotation_data$ensembl_gene_id, express_data)

  })

  

  # Write results to file

  for (i in 1:length(modules)) {

    toAppend <- ifelse (i == 1, FALSE, TRUE)

    module_color <- names(modules[i])

    write.xlsx(module_expression[[i]], filename, sheetName = module_color, showNA = TRUE, col.names = TRUE,

               row.names = FALSE, append = toAppend)

  }

}

Create function to annotate gene list

# gene_list is a named vector of genes

# exprs_matrix is matrix of fold changes

annotate_gene_list <- function(gene_list, exprs_matrix) {

  # collect annotation and expression data

  annotation_data <- gene_map[(match(gene_list, gene_map$external_gene_name)), ]

  express_data <- exprs_matrix[row.names(exprs_matrix) %in% gene_list, ]

    

  # reorder expression data to match annotation

  express_data <- express_data[(match(gene_list, row.names(express_data))), ]

  annotated <- cbind.data.frame(geneName = gene_list, ensembl = annotation_data$ensembl_gene_id, express_data)

}

4. Heatmap of differentially expressed genes in WNV infection

Visualize union of WT differentially expressed genes

# Number of differentially expressed genes in at least one sample

length(row.names(expression))

## [1] 6823

# Create heatmap of expression profile

express_map <- heatmap.F.4(expression, cutoff = 3, distmethod = "euclidean",clustermethod = "ward.D", 

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Write copy of heatmap to file

pdf("figures/Heatmap_DE_genes.pdf")

express_map <- heatmap.F.4(expression, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(express_map, expression)

Plot numbers of upregulated and downregulated genes

# Calculate numbers of upregulated and downregulated genes

updown <- t(apply(DE_results, MARGIN = 2, FUN = table))

updown <- as.data.frame(updown[, colnames(updown) != 0])

updown <- cbind(updown, sample = row.names(updown), ord = seq(1, length(row.names(updown))))

updown

##                   -1   1          sample ord

## BMM_WT_24h      3511 861      BMM_WT_24h   1

## BMM_LGP2_24h    3757 795    BMM_LGP2_24h   2

## BMM_MDA5_24h    3160 801    BMM_MDA5_24h   3

## BMM_DKO_24h      238 321     BMM_DKO_24h   4

## BMM_RIGI.WT_24h 3562 922 BMM_RIGI.WT_24h   5

## BMM_RIGI_24h    3350 840    BMM_RIGI_24h   6

## BMM_WT_48h      2783 759      BMM_WT_48h   7

## BMM_LGP2_48h    2694 762    BMM_LGP2_48h   8

## BMM_MDA5_48h    2886 810    BMM_MDA5_48h   9

## BMM_RIGI.WT_48h 2160 731 BMM_RIGI.WT_48h  10

## BMM_RIGI_48h    2369 746    BMM_RIGI_48h  11

# Reshape data frame for ggplot2 requirements

updown <- melt(updown, id.vars = c("sample", "ord"), variable.name = "direction", value.name = "count")

for (i in 1:(length(row.names(updown))/2)) {

  updown[i, 4] <- updown[i, 4] * -1

}

updown$pos <- updown$count >= 0



# Create Plot

DE_plot <- ggplot(updown, aes(x = reorder(sample, ord), y = count, fill = pos)) +

  geom_bar(stat = "identity", pos = "identity") +

  scale_fill_manual(values = c("#0000FF", "#FF0000"), guide = FALSE) +

  ylab("Number of Differentially Expressed Genes") + 

  theme_bw() +

  theme(panel.grid.major.x = element_blank(),

        panel.grid.minor.x = element_blank(),

        panel.grid.major.y = element_blank(),

        panel.grid.minor.y = element_blank(),

        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),

        axis.title.x = element_blank())

print(DE_plot)

# Create PDF of figure

pdf("figures/Barplot_DE_numbers.pdf")

print(DE_plot)

dev.off()

## png 

##   2

5. Create Heatmaps of “wild-type” DE genes from Immune and Innate Immune Signaling GOBP categories

Immune response (GO:0006955)

# Read in gene list (single column)

immune <- read.csv(in.immune, header = FALSE)$V1



# Extract DEs mapping to each category

immune_DEs <- expression[which(toupper(row.names(expression)) %in% immune), ]



# Number of differentially expressed immune genes

length(row.names(immune_DEs))

## [1] 532

# Create heatmap of expression profile

immune_map <- heatmap.F.4(immune_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                          symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from immune modules

pdf("figures/Heatmap_Immune_GO0006955_DE.pdf")

immune_map <- heatmap.F.4(immune_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                          symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(immune_map, immune_DEs)

Innate immune response (GO:0045087)

# Read in gene list (single column)

innate <- read.csv(in.innate, header = FALSE)$V1



# Extract DEs mapping to each category

innate_DEs <- expression[which(toupper(row.names(expression)) %in% innate), ]



# Number of differentially expressed innate immune genes

length(row.names(innate_DEs))

## [1] 281

# Create heatmap of expression profile

innate_map <- heatmap.F.4(innate_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                          symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from innate modules

pdf("figures/Heatmap_Innate_GO0045087_DE.pdf")

innate_map <- heatmap.F.4(innate_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                          symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(innate_map, innate_DEs)

Create heatmap of all WNV-induced ISGs among differentially expressed genes ISGs defined per IFN-beta treatment of BMM for 6 or 24 hours inducing > 2-fold gene induction (Al Huang data set)

# Read in measurement of IFN-beta treatment of BMM: 6 and 24 hours

ISG_expression <- read.csv(in.isgs, header = TRUE, row.names = 1)



# Which genes are induced > 2-fold at either 6 or 24hr post-IFN-beta treatment in BMM?

ISGs <- row.names(ISG_expression[which(ISG_expression["mac_beta_vs_mock_6h"] > 1 |

                                       ISG_expression["mac_beta_vs_mock_24h"] > 1), ])



# Extract WNV DEs mapping to ISG list

ISG_DEs <- expression[which(toupper(row.names(expression)) %in% toupper(ISGs)), ]



# Number of differentially expressed ISGs

length(row.names(ISG_DEs))

## [1] 341

# Create heatmap of expression profile

ISG_map <- heatmap.F.4(ISG_DEs, cutoff = 4, distmethod = "euclidean", clustermethod = "ward.D",

                       symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from ISG lists

pdf("figures/Heatmap_ISG_DE.pdf")

ISG_map <- heatmap.F.4(ISG_DEs, cutoff = 4, distmethod = "euclidean", clustermethod = "ward.D",

                       symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(ISG_map, ISG_DEs)

Non-ISGs differentially expressed following WNV infection (all genes NOT in above list)

# Extract WNV DEs NOT mapping to ISG list

non_ISG_DEs <- expression[-which(toupper(row.names(expression)) %in% toupper(ISGs)), ]



# Number of differentially expressed non-ISGs

length(row.names(non_ISG_DEs))

## [1] 6482

# Create heatmap of expression profile

non_ISG_map <- heatmap.F.4(non_ISG_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from ISG lists

pdf("figures/Heatmap_non_ISG_DE.pdf")

non_ISG_map <- heatmap.F.4(non_ISG_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(non_ISG_map, non_ISG_DEs)

7. RLR-dependent WNV-induced transcriptional response

Visualize those genes that are differentially expressed in B6 WT but not RIG-I/MDA5 DKO BMM (24h only)

# Use sigMask to define filter for genes that are DE in B6 WT but not DKO at 24h

WT_24h_up <- row.names(sigMask[which(sigMask[, "BMM_WT_24h"] > 0), ]) 

WT_24h_down <- row.names(sigMask[which(sigMask[, "BMM_WT_24h"] < 0), ])

WT_48h_up <- row.names(sigMask[which(sigMask[, "BMM_WT_48h"] > 0), ]) 

WT_48h_down <- row.names(sigMask[which(sigMask[, "BMM_WT_48h"] < 0), ])

DKO_24h_up <- row.names(sigMask[which(sigMask[, "BMM_DKO_24h"] > 0), ])

DKO_24h_down <- row.names(sigMask[which(sigMask[, "BMM_DKO_24h"] < 0), ])



RLR_dependent_up <- WT_24h_up[-(WT_24h_up %in% DKO_24h_up)]

RLR_dependent_down <- WT_24h_down[-(WT_24h_down %in% DKO_24h_down)]

RLR_dependent_genes <- c(RLR_dependent_up, RLR_dependent_down)



# Extract expression matrix mapping to RLR-dependent genes (only show 24h samples)

RLR_dependent_DEs <- expression[which(row.names(expression) %in% RLR_dependent_genes), 

                                c("BMM_WT_24h", "BMM_LGP2_24h", "BMM_MDA5_24h", "BMM_DKO_24h",

                                  "BMM_RIGI.WT_24h", "BMM_RIGI_24h")]



# Number of RLR-dependent DEs

length(row.names(RLR_dependent_DEs))

## [1] 4370

# Create heatmap of expression profile

RLR_map <- heatmap.F.4(RLR_dependent_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                       symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from RLR-dependent gene lists

pdf("figures/Heatmap_RLR_dependent_DE.pdf")

RLR_map <- heatmap.F.4(RLR_dependent_DEs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                       symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(RLR_map, RLR_dependent_DEs)

RLR-dependent ISGs (does the pattern of RLR-dependent ISGs look like that of non-ISGs?)

# Extract RLR-dependent DEs mapping to ISG list

RLR_ISGs <- RLR_dependent_DEs[which(toupper(row.names(RLR_dependent_DEs)) %in% toupper(ISGs)), ]



# Number of RLR-dependent ISGs

length(row.names(RLR_ISGs))

## [1] 239

# Create heatmap of expression profile

RLR_ISG_map <- heatmap.F.4(RLR_ISGs, cutoff = 4, distmethod = "euclidean", clustermethod = "ward.D",

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from ISG lists

pdf("figures/Heatmap_RLR_dependent_ISGs.pdf")

RLR_ISG_map <- heatmap.F.4(RLR_ISGs, cutoff = 4, distmethod = "euclidean", clustermethod = "ward.D",

                           symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(RLR_ISG_map, RLR_ISGs)

RLR-dependent non-ISGs (does the pattern of RLR-dependent ISGs look like that of non-ISGs?)

# Extract RLR-dependent DEs mapping to ISG list

RLR_non_ISGs <- RLR_dependent_DEs[-which(toupper(row.names(RLR_dependent_DEs)) %in% toupper(ISGs)), ]



# Number of RLR-dependent ISGs

length(row.names(RLR_non_ISGs))

## [1] 4131

# Create heatmap of expression profile

RLR_non_ISG_map <- heatmap.F.4(RLR_non_ISGs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                               symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

# Create heatmaps from ISG lists

pdf("figures/Heatmap_RLR_dependent_non_ISGs.pdf")

RLR_non_ISG_map <- heatmap.F.4(RLR_non_ISGs, cutoff = 3, distmethod = "euclidean", clustermethod = "ward.D",

                               symbreaks = TRUE, margins = c(10, 1), keytitle = "log2 Fold Change")

dev.off()

## png 

##   2

# Write results of each module to file

collect_modules(RLR_non_ISG_map, RLR_non_ISGs)

8. Identify RLR-specific gene signatures

24h time point

# Use sigMask to define filter for genes that are DE in B6 WT but not KO at 24h

LGP2_up_24 <- sigMask[, "BMM_WT_24h"] > 0 & !sigMask[, "BMM_LGP2_24h"] > 0

LGP2_down_24 <- sigMask[, "BMM_WT_24h"] < 0 & !sigMask[, "BMM_LGP2_24h"] < 0

MDA5_up_24 <- sigMask[, "BMM_WT_24h"] > 0 & !sigMask[, "BMM_MDA5_24h"] > 0

MDA5_down_24 <- sigMask[, "BMM_WT_24h"] < 0 & !sigMask[, "BMM_MDA5_24h"] < 0

RIGI_up_24 <- sigMask[, "BMM_RIGI.WT_24h"] > 0 & !sigMask[, "BMM_RIGI_24h"] > 0

RIGI_down_24 <- sigMask[, "BMM_RIGI.WT_24h"] < 0 & !sigMask[, "BMM_RIGI_24h"] < 0



# Venn diagrams

KO_genes_up_24 <- cbind(LGP2_up_24, MDA5_up_24, RIGI_up_24)

KO_counts_up_24 <- vennCounts(KO_genes_up_24)

test <- vennDiagram(KO_counts_up_24)

pdf("figures/Venn_KO_differential_24h_up.pdf")

vennDiagram(KO_counts_up_24)

dev.off()

## png 

##   2

KO_genes_down_24 <- cbind(LGP2_down_24, MDA5_down_24, RIGI_down_24)

KO_counts_down_24 <- vennCounts(KO_genes_down_24)

vennDiagram(KO_counts_down_24)

pdf("figures/Venn_KO_differential_24h_down.pdf")

vennDiagram(KO_counts_down_24)

dev.off()

## png 

##   2

# Get gene names for each RLR differential list

LGP2_24h_up <- names(LGP2_up_24[which(LGP2_up_24 == TRUE)])

LGP2_24h_down <- names(LGP2_down_24[which(LGP2_down_24 == TRUE)])

MDA5_24h_up <- names(MDA5_up_24[which(MDA5_up_24 == TRUE)])

MDA5_24h_down <- names(MDA5_down_24[which(MDA5_down_24 == TRUE)])

RIGI_24h_up <- names(RIGI_up_24[which(RIGI_up_24 == TRUE)])

RIGI_24h_down <- names(RIGI_down_24[which(RIGI_down_24 == TRUE)])



# Write tables of venn diagram data

write.table(KO_genes_up_24, "results/KO_genes_uo_24.txt", sep = "\t", quote = FALSE, col.names = NA)

write.table(KO_genes_down_24, "results/KO_genes_down_24.txt", sep = "\t", quote = FALSE, col.names = NA)

48h time point

# Use sigMask to define filter for genes that are DE in B6 WT but not KO at 48h

LGP2_up_48 <- sigMask[, "BMM_WT_48h"] > 0 & !sigMask[, "BMM_LGP2_48h"] > 0

LGP2_down_48 <- sigMask[, "BMM_WT_48h"] < 0 & !sigMask[, "BMM_LGP2_48h"] < 0

MDA5_up_48 <- sigMask[, "BMM_WT_48h"] > 0 & !sigMask[, "BMM_MDA5_48h"] > 0

MDA5_down_48 <- sigMask[, "BMM_WT_48h"] < 0 & !sigMask[, "BMM_MDA5_48h"] < 0

RIGI_up_48 <- sigMask[, "BMM_RIGI.WT_48h"] > 0 & !sigMask[, "BMM_RIGI_48h"] > 0

RIGI_down_48 <- sigMask[, "BMM_RIGI.WT_48h"] < 0 & !sigMask[, "BMM_RIGI_48h"] < 0



# Venn diagrams

KO_genes_up_48 <- cbind(LGP2_up_48, MDA5_up_48, RIGI_up_48)

KO_counts_up_48 <- vennCounts(KO_genes_up_48)

vennDiagram(KO_counts_up_48)

pdf("figures/Venn_KO_differential_48h_up.pdf")

vennDiagram(KO_counts_up_48)

dev.off()

## png 

##   2

KO_genes_down_48 <- cbind(LGP2_down_48, MDA5_down_48, RIGI_down_48)

KO_counts_down_48 <- vennCounts(KO_genes_down_48)

vennDiagram(KO_counts_down_48)

pdf("figures/Venn_KO_differential_48h_down.pdf")

vennDiagram(KO_counts_down_48)

dev.off()

## png 

##   2

# Get gene names for each RLR differential list

LGP2_48h_up <- names(LGP2_up_48[which(LGP2_up_48 == TRUE)])

LGP2_48h_down <- names(LGP2_down_48[which(LGP2_down_48 == TRUE)])

MDA5_48h_up <- names(MDA5_up_48[which(MDA5_up_48 == TRUE)])

MDA5_48h_down <- names(MDA5_down_48[which(MDA5_down_48 == TRUE)])

RIGI_48h_up <- names(RIGI_up_48[which(RIGI_up_48 == TRUE)])

RIGI_48h_down <- names(RIGI_down_48[which(RIGI_down_48 == TRUE)])



# Write tables of venn diagram data

write.table(KO_genes_up_48, "results/KO_genes_uo_48.txt", sep = "\t", quote = FALSE, col.names = NA)

write.table(KO_genes_down_48, "results/KO_genes_down_48.txt", sep = "\t", quote = FALSE, col.names = NA)

9. Identify pathway and polarization profiles among gene sets

Create radar plot showing patterns of genes in each list of KO vs WT that fall into particular categories

# Read in additional pathway lists (all single columns)

Th1_genes <- read.table(in.th1, sep = "\t", header = FALSE)$V1

Th2_genes <- read.table(in.th2, sep = "\t", header = FALSE)$V1

M1_genes <- read.table(in.m1, sep = "\t", header = FALSE)$V1

M2_genes <- read.table(in.m2, sep = "\t", header = FALSE)$V1



# Get percentage of Th1, Th2, M1, M2, or ISG for each gene list

getProportion <- function(geneList, category) {

  sum(toupper(geneList) %in% toupper(category))/length(category)

}



# Identify gene lists to compare

geneLists_KO_24h_up <- c("LGP2_24h_up", "MDA5_24h_up", "RIGI_24h_up")

geneLists_KO_24h_down <- c("LGP2_24h_down", "MDA5_24h_down", "RIGI_24h_down")

geneLists_KO_48h_up <- c("LGP2_48h_up", "MDA5_48h_up", "RIGI_48h_up")

geneLists_KO_48h_down <- c( "LGP2_48h_down", "MDA5_48h_down", "RIGI_48h_down")



# Identify colors for plots

KO_cols <- c("red", "blue", "orange")



createRadarData <- function(geneList) {

  Th1 <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), Th1_genes)})

  Th2 <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), Th2_genes)})

  M1 <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), M1_genes)})

  M2 <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), M2_genes)})

  ISG <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), ISGs)})

  Innate <- sapply(geneList, function(x) { getProportion(eval(parse(text = x)), innate)})

  byCategory <- cbind.data.frame(Th1, Th2, M1, M2, ISG, Innate)

  byCategory # return data frame ready for radar plot

}



# KO samples 24h up 

radarPlot_KO_24h_up <- createRadarData(geneLists_KO_24h_up)

radial.plot(radarPlot_KO_24h_up, labels = colnames(radarPlot_KO_24h_up),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_24h_up)))

legend(-0.2, 0.15, row.names(radarPlot_KO_24h_up), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO samples 24h down

radarPlot_KO_24h_down <- createRadarData(geneLists_KO_24h_down)

radial.plot(radarPlot_KO_24h_down, labels = colnames(radarPlot_KO_24h_down),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_24h_down)))

legend(-0.2, 0.15, row.names(radarPlot_KO_24h_down), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO samples 48h up 

radarPlot_KO_48h_up <- createRadarData(geneLists_KO_48h_up)

radial.plot(radarPlot_KO_48h_up, labels = colnames(radarPlot_KO_48h_up),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_48h_up)))

legend(-0.06, 0.045, row.names(radarPlot_KO_48h_up), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO samples 48h down

radarPlot_KO_48h_down <- createRadarData(geneLists_KO_48h_down)

radial.plot(radarPlot_KO_48h_down, labels = colnames(radarPlot_KO_48h_down),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_48h_down)))

legend(-0.07, 0.05, row.names(radarPlot_KO_48h_down), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

Radar plots of genes independent of each individual RLR (present in each knockout)

# Use sigMask to define filter for genes that are DE in B6 WT but not KO at 24h

LGP2_present_up_24 <- names(which(sigMask[, "BMM_LGP2_24h"] > 0) == TRUE)

LGP2_present_down_24 <- names(which(sigMask[, "BMM_LGP2_24h"] < 0) == TRUE)

MDA5_present_up_24 <- names(which(sigMask[, "BMM_MDA5_24h"] > 0) == TRUE)

MDA5_present_down_24 <- names(which(sigMask[, "BMM_MDA5_24h"] < 0) == TRUE)

RIGI_present_up_24 <- names(which(sigMask[, "BMM_RIGI_24h"] > 0) == TRUE)

RIGI_present_down_24 <- names(which(sigMask[, "BMM_RIGI_24h"] < 0) == TRUE)

KO_present_24h_up <- c("LGP2_present_up_24", "MDA5_present_up_24", "RIGI_present_up_24")

KO_present_24h_down <- c("LGP2_present_down_24", "MDA5_present_down_24", "RIGI_present_down_24")



LGP2_present_up_48 <- names(which(sigMask[, "BMM_LGP2_48h"] > 0) == TRUE)

LGP2_present_down_48 <- names(which(sigMask[, "BMM_LGP2_48h"] < 0) == TRUE)

MDA5_present_up_48 <- names(which(sigMask[, "BMM_MDA5_48h"] > 0) == TRUE)

MDA5_present_down_48 <- names(which(sigMask[, "BMM_MDA5_48h"] < 0) == TRUE)

RIGI_present_up_48 <- names(which(sigMask[, "BMM_RIGI_48h"] > 0) == TRUE)

RIGI_present_down_48 <- names(which(sigMask[, "BMM_RIGI_48h"] < 0) == TRUE)

KO_present_48h_up <- c("LGP2_present_up_48", "MDA5_present_up_48", "RIGI_present_up_48")

KO_present_48h_down <- c("LGP2_present_down_48", "MDA5_present_down_48", "RIGI_present_down_48")



# KO genes present 24h up 

radarPlot_KO_present_24h_up <- createRadarData(KO_present_24h_up)

radial.plot(radarPlot_KO_present_24h_up, labels = colnames(radarPlot_KO_present_24h_up),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_present_24h_up)))

legend(-1, 0.8, row.names(radarPlot_KO_present_24h_up), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO genes present 24h down

radarPlot_KO_present_24h_down <- createRadarData(KO_present_24h_down)

radial.plot(radarPlot_KO_present_24h_down, labels = colnames(radarPlot_KO_present_24h_down),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_present_24h_down)))

legend(-0.7, 0.6, row.names(radarPlot_KO_present_24h_down), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO genes present 48h up

radarPlot_KO_present_48h_up <- createRadarData(KO_present_48h_up)

radial.plot(radarPlot_KO_present_24h_up, labels = colnames(radarPlot_KO_present_48h_up),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_present_48h_up)))

legend(-1, 0.8, row.names(radarPlot_KO_present_48h_up), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

# KO genes present 48h down

radarPlot_KO_present_48h_down <- createRadarData(KO_present_48h_down)

radial.plot(radarPlot_KO_present_48h_down, labels = colnames(radarPlot_KO_present_48h_down),

            line.col = KO_cols, lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_KO_present_48h_down)))

legend(-0.7, 0.6, row.names(radarPlot_KO_present_48h_down), col = KO_cols, lty = 1, cex = 0.8, bty = "n")

WT (B6) radar plots

# WT 24h up

radarPlot_WT_24h_up <- createRadarData(c("WT_24h_up"))

radial.plot(radarPlot_WT_24h_up, labels = colnames(radarPlot_WT_24h_up),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_WT_24h_up)))

legend(-0.9, 0.6, row.names(radarPlot_WT_24h_up), col = "black", lty = 1, cex = 0.8, bty = "n")

# WT 24h down

radarPlot_WT_24h_down <- createRadarData(c("WT_24h_down"))

radial.plot(radarPlot_WT_24h_down, labels = colnames(radarPlot_WT_24h_down),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_WT_24h_down)))

legend(-0.45, 0.35, row.names(radarPlot_WT_24h_down), col = "black", lty = 1, cex = 0.8, bty = "n")

# WT 48h up

radarPlot_WT_48h_up <- createRadarData(c("WT_48h_up"))

radial.plot(radarPlot_WT_48h_up, labels = colnames(radarPlot_WT_48h_up),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_WT_48h_up)))

legend(-0.9, 0.7, row.names(radarPlot_WT_48h_up), col = "black", lty = 1, cex = 0.8, bty = "n")

# WT 48h down

radarPlot_WT_48h_down <- createRadarData(c("WT_48h_down"))

radial.plot(radarPlot_WT_48h_down, labels = colnames(radarPlot_WT_48h_down),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_WT_48h_down)))

legend(-0.55, 0.4, row.names(radarPlot_WT_48h_down), col = "black", lty = 1, cex = 0.8, bty = "n")

RLR-dependent/independent(present in DKO cells) gene radar plots

# RLR-dependent 24h up

radarPlot_RLR_dep_24h_up <- createRadarData(c("RLR_dependent_up"))

radial.plot(radarPlot_RLR_dep_24h_up, labels = colnames(radarPlot_RLR_dep_24h_up),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_RLR_dep_24h_up)))

legend(-0.9, 0.7, row.names(radarPlot_RLR_dep_24h_up), col = "black", lty = 1, cex = 0.8, bty = "n")

# RLR-dependent 24h down

radarPlot_RLR_dep_24h_down <- createRadarData(c("RLR_dependent_down"))

radial.plot(radarPlot_RLR_dep_24h_down, labels = colnames(radarPlot_RLR_dep_24h_down),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_RLR_dep_24h_down)))

legend(-0.5, 0.35, row.names(radarPlot_RLR_dep_24h_down), col = "black", lty = 1, cex = 0.8, bty = "n")

# DKO 24h up (RLR-independent)

radarPlot_DKO_24h_up <- createRadarData(c("DKO_24h_up"))

radial.plot(radarPlot_DKO_24h_up, labels = colnames(radarPlot_DKO_24h_up),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_DKO_24h_up)))

legend(-0.35, 0.3, row.names(radarPlot_DKO_24h_up), col = "black", lty = 1, cex = 0.8, bty = "n")

# DKO 24h down (RLR-independent)

radarPlot_DKO_24h_down <- createRadarData(c("DKO_24h_down"))

radial.plot(radarPlot_DKO_24h_down, labels = colnames(radarPlot_DKO_24h_down),

            line.col = "black", lwd = 2, label.prop = 1.15, mar = c(3, 8, 3, 2),

            boxed.radial = FALSE, rp.type = "p", show.grid.labels = 3,

            radial.lim = c(0, max(radarPlot_DKO_24h_down)))

legend(-0.025, 0.02, row.names(radarPlot_DKO_24h_down), col = "black", lty = 1, cex = 0.8, bty = "n")