-
Notifications
You must be signed in to change notification settings - Fork 1
/
Figure2A_correlation.Rmd
85 lines (65 loc) · 3.18 KB
/
Figure2A_correlation.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
---
title: "Correlation between Fold change of B6 vs NZO"
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("ggpubr")
library(VennDiagram)
library(gplots)
library(dplyr)
```
### Function to convert ensembl IDs to gene symbols
```{r convert_gene, 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)
sym.genes <- mapping$symbol[m]
names(sym.genes) <- genes
return(sym.genes)
}
```
```{r main_function, include=TRUE, echo=TRUE}
Celltypes <- c("MEMORY", "NAIVE", "PBL", "SPLEEN")
ap1_genes = c("Jdp2", "Atf6", "Maf", "Batf3", "Mafb", "Mafa", "Mafg", "Atf5", "Atf3", "Batf2", "Atf2", "Atf7", "Atf6b", "Batf", "Junb", "Fos", "Jun", "Fosb", "Jund", "Atf4", "Maff", "Mafk")
for (i in Celltypes) {
fileB6 <- paste0("./data/", i, "_Age18vs3_B6_RNAseq.csv")
fileNZO <- paste0("./data/", i, "_Age18vs3_NZO_RNAseq.csv")
B6_file_orig <- read.delim(fileB6, sep=",", quote = "\"", header=TRUE)
NZO_file_orig <- read.delim(fileNZO, sep=",", quote = "\"", header=TRUE)
B6_NZO_file_orig <- merge(B6_file_orig, NZO_file_orig, by = "Gene.Ensembl")
B6_NZO_file <- B6_NZO_file_orig[(((B6_NZO_file_orig$adj.P.Val.x<0.05) & (abs(B6_NZO_file_orig$logFC.x) > 1)) | ((B6_NZO_file_orig$adj.P.Val.y<0.05) & (abs(B6_NZO_file_orig$logFC.y)>1))),]
print(head(B6_NZO_file))
B6_NZO_file$gene <- ens2gene(B6_NZO_file$Gene.Ensembl)
B6_NZO_file$color_val <- "grey"
B6_NZO_file <- B6_NZO_file%>%mutate(color_val = case_when(B6_NZO_file$gene %in% ap1_genes == "TRUE" ~ "red", TRUE~color_val))
print(i)
print(dim(B6_NZO_file))
print(head(B6_NZO_file))
filename = paste0("./output/",i, "_B6_NZO_merged.txt")
write.table(B6_NZO_file, filename, sep="\t", quote=FALSE)
cor_calc <- cor.test(B6_NZO_file$logFC.x, B6_NZO_file$logFC.y)
print("Correlation test between B6 vs NZO for memory cells")
print(cor_calc)
plot_filename <- paste0("./output/Plot_FC_", i, "_B6_NZO_RNASeq.pdf")
print(plot_filename)
pdf(plot_filename)
print(ggscatter(B6_NZO_file, x = "logFC.x", y = "logFC.y", color = "color_val", label = "gene", repel = TRUE, label.select = ap1_genes, palette = c(grey = "grey", red = "red", blue = "blue"), add = "reg.line", conf.int = TRUE, cor.coef = TRUE, cor.method = "pearson", title = plot_filename, xlab = "FC (RNASeq) for B6", ylab = "FC (RNAseq) for NZO"))
#print(ggscatter(B6_NZO_file, x = "logFC.x", y = "logFC.y", label = B6_NZO_file$gene, repel = TRUE, label.select = ap1_genes, color = "color_val", palette = c("red", "grey"), font.label = "black", add = "reg.line", conf.int = TRUE, cor.coef = TRUE, cor.method = "pearson", xlab = "FC (RNASeq) for B6", ylab = "FC (RNAseq) for NZO"))
dev.off()
}
```