-
Notifications
You must be signed in to change notification settings - Fork 1
/
Figure2G_boxplot_pbl.Rmd
144 lines (118 loc) · 6.37 KB
/
Figure2G_boxplot_pbl.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
---
title: "Boxplots of fold change based on AP1 subfamily"
author: "Neerja Katiyar"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: workflowr::wflow_html
editor_options:
chunk_output_type: console
---
### Prepare the environment
```{r setup, include=FALSE}
require(knitr)
knitr::opts_chunk$set(echo = TRUE)
opts_knit$set(root.dir = "/Users/katiyn/Dropbox (JAX)/MouseAging_clean/Mice_aging_NK_resubmission/code/") #set root dir!
```
```{r libraries, include=TRUE, echo=TRUE}
library(ggplot2)
library(ggpubr)
library(rstatix)
```
### Function for mapping Ensembl IDs to gene symbols
```{r convert_geneID, include=TRUE, echo=TRUE}
ens2gene <- function(genes){
genome <- annotables::grcm38
# ens to gene symbol mapping
mapping <-
base::subset(genome,
genome$ensgene %in% genes,
select = c('ensgene', 'symbol'))
m <- match(genes, mapping$ensgene)
ens.genes <- mapping$symbol[m]
names(ens.genes) <- genes
return(ens.genes)
}
```
### Classify the data into AP1 subfamilies.
```{r AP1_family, include=TRUE, echo=TRUE}
data_AP1_all <- read.table("./data/Table_B6_NZO_FC_RNA_AP1_all_genes.txt", sep="\t", header=T, row.names=1)
data_AP1_family <- data_AP1_all[data_AP1_all$tissue == "PBL",]
data_AP1_family_Jun <- data_AP1_family[data_AP1_family$gene == "Jun",]
data_AP1_family_Junb <- data_AP1_family[data_AP1_family$gene == "Junb",]
data_AP1_family_Jund <- data_AP1_family[data_AP1_family$gene == "Jund",]
data_AP1_family_Fos <- data_AP1_family[data_AP1_family$gene == "Fos",]
data_AP1_family_Fosb <- data_AP1_family[data_AP1_family$gene == "Fosb",]
data_AP1_family_Fra1 <- data_AP1_family[data_AP1_family$gene == "Fra1",]
data_AP1_family_Fra2 <- data_AP1_family[data_AP1_family$gene == "Fra2",]
data_AP1_family_Atf <- data_AP1_family[data_AP1_family$gene == "Atf",]
data_AP1_family_Atf2 <- data_AP1_family[data_AP1_family$gene == "Atf2",]
data_AP1_family_Atf3 <- data_AP1_family[data_AP1_family$gene == "Atf3",]
data_AP1_family_Atf4 <- data_AP1_family[data_AP1_family$gene == "Atf4",]
data_AP1_family_Atf5 <- data_AP1_family[data_AP1_family$gene == "Atf5",]
data_AP1_family_Atf6 <- data_AP1_family[data_AP1_family$gene == "Atf6",]
data_AP1_family_Atf6b <- data_AP1_family[data_AP1_family$gene == "Atf6b",]
data_AP1_family_Atf7 <- data_AP1_family[data_AP1_family$gene == "Atf7",]
data_AP1_family_Batf <- data_AP1_family[data_AP1_family$gene == "Batf",]
data_AP1_family_Batf2 <- data_AP1_family[data_AP1_family$gene == "Batf2",]
data_AP1_family_Batf3 <- data_AP1_family[data_AP1_family$gene == "Batf3",]
data_AP1_family_Maf <- data_AP1_family[data_AP1_family$gene == "Maf",]
data_AP1_family_Mafa <- data_AP1_family[data_AP1_family$gene == "Mafa",]
data_AP1_family_Mafb <- data_AP1_family[data_AP1_family$gene == "Mafb",]
data_AP1_family_Maff <- data_AP1_family[data_AP1_family$gene == "Maff",]
data_AP1_family_Mafg <- data_AP1_family[data_AP1_family$gene == "Mafg",]
data_AP1_family_Mafk <- data_AP1_family[data_AP1_family$gene == "Mafk",]
######################################################################
data_AP1_family_Jun$AP1_subfamily <- c(rep("Jun"))
data_AP1_family_Junb$AP1_subfamily <- c(rep("Jun"))
data_AP1_family_Jund$AP1_subfamily <- c(rep("Jun"))
data_AP1_family_Fos$AP1_subfamily <- c(rep("Fos"))
data_AP1_family_Fosb$AP1_subfamily <- c(rep("Fos"))
#data_AP1_family_Fra1$AP1_subfamily <- c(rep("Fos"))
#data_AP1_family_Fra2$AP1_subfamily <- c(rep("Fos"))
#data_AP1_family_Atf$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf2$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf3$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf4$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf5$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf6$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf6b$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Atf7$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Batf$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Batf2$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Batf3$AP1_subfamily <- c(rep("Atf"))
data_AP1_family_Maf$AP1_subfamily <- c(rep("Maf"))
data_AP1_family_Mafa$AP1_subfamily <- c(rep("Maf"))
data_AP1_family_Mafb$AP1_subfamily <- c(rep("Maf"))
data_AP1_family_Maff$AP1_subfamily <- c(rep("Maf"))
data_AP1_family_Mafg$AP1_subfamily <- c(rep("Maf"))
data_AP1_family_Mafk$AP1_subfamily <- c(rep("Maf"))
df_AP1_family_all <- rbind(data_AP1_family_Jun, data_AP1_family_Junb, data_AP1_family_Jund, data_AP1_family_Fos, data_AP1_family_Fosb, data_AP1_family_Atf2, data_AP1_family_Atf3, data_AP1_family_Atf4, data_AP1_family_Atf5, data_AP1_family_Atf6, data_AP1_family_Atf6b, data_AP1_family_Atf7, data_AP1_family_Batf, data_AP1_family_Batf2, data_AP1_family_Batf3, data_AP1_family_Maf, data_AP1_family_Mafa, data_AP1_family_Mafb, data_AP1_family_Maff, data_AP1_family_Mafb, data_AP1_family_Mafg, data_AP1_family_Mafk)
write.table(df_AP1_family_all, "./output/PBL_FC_AP1_all.txt", sep="\t", quote=FALSE)
```
#########################################################
### Generate boxplots to compare fold change with age in individual AP1 subfamily.
```{r generate_AP1_family_boxplot, include=TRUE, echo=TRUE}
df <- df_AP1_family_all
df_B6 <- df[df$strain == "B6",]
df_NZO <- df[df$strain == "NZO",]
data_Jun_B6 <- df_B6[df_B6$AP1_subfamily == "Jun",]
data_Jun_NZO <- df_NZO[df_NZO$AP1_subfamily == "Jun",]
data_Fos_B6 <- df_B6[df_B6$AP1_subfamily == "Fos",]
data_Fos_NZO <- df_NZO[df_NZO$AP1_subfamily == "Fos",]
data_Atf_B6 <- df_B6[df_B6$AP1_subfamily == "Atf",]
data_Atf_NZO <- df_NZO[df_NZO$AP1_subfamily == "Atf",]
data_Maf_B6 <- df_B6[df_B6$AP1_subfamily == "Maf",]
data_Maf_NZO <- df_NZO[df_NZO$AP1_subfamily == "Maf",]
#One sample t-test
Jun_B6_test <- t.test(data_Jun_B6$FC, mu=0, alternative = "greater")
Jun_NZO_test <- t.test(data_Jun_NZO$FC, mu=0, alternative = "greater")
Fos_B6_test <- t.test(data_Fos_B6$FC, mu=0, alternative = "greater")
Fos_NZO_test <- t.test(data_Fos_NZO$FC, mu=0, alternative = "greater")
Atf_B6_test <- t.test(data_Atf_B6$FC, mu=0, alternative = "greater")
Atf_NZO_test <- t.test(data_Atf_NZO$FC, mu=0, alternative = "greater")
Maf_B6_test <- t.test(data_Maf_B6$FC, mu=0, alternative = "greater")
Maf_NZO_test <- t.test(data_Maf_NZO$FC, mu=0, alternative = "greater")
p <- ggboxplot(df_AP1_family_all, x="AP1_subfamily", y="FC", color = "strain", palette = "jco", add="jitter", facet.by = "tissue", nrow=1, short.panel.labs = FALSE)
pdf("./output/Boxplot_PBL_FC_RNA_DE_not_filtered_all_24_AP1_genes.pdf", height=10, width=10)
p
dev.off()
```