-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathFigure4C_vlnplot.Rmd
125 lines (103 loc) · 5.46 KB
/
Figure4C_vlnplot.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
---
title: "Script to generate expression plots using seurat object"
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(dplyr)
library(Seurat)
library(cowplot)
library(ggplot2)
library(dplyr)
library(forcats)
library(hrbrthemes)
library(viridis)
library(ggpubr)
library(gridBase)
library(gridExtra)
library(tidyr)
library(cowplot)
```
###Generate expression plots.
```{r func_generate_violinplot_piechart, include=TRUE, echo=TRUE}
gene_Vln_Pie <- function(mydata, gene_name, Celltype) {
subset_celltype <- subset(mydata, idents = Celltype)
subset_celltype$celltype.stim <- paste(Idents(subset_celltype), subset_celltype$age, sep = "_")
subset_celltype$celltype <- Idents(subset_celltype)
Idents(subset_celltype) <- "celltype.stim"
df = as.data.frame(subset_celltype[["RNA"]]@counts)
new_df <- df[, df[gene_name,] > 0]
subset_celltype_gene <- CreateSeuratObject(counts = new_df, assay = "RNA")
subset_celltype_gene <- AddMetaData(subset_celltype_gene, [email protected][colnames(subset_celltype_gene), , drop=F])
exprn_val <- subset(subset_celltype_gene[["RNA"]]@counts, rownames(subset_celltype_gene[["RNA"]]@counts) %in% gene_name)
exprn_val_df <- as.data.frame(exprn_val)
rownames(exprn_val_df) <- colnames(subset_celltype_gene[["RNA"]]@counts)
celltype_age <- [email protected]$celltype.stim
celltype <- [email protected]$celltype
df_final_AV <- cbind(exprn_val_df, celltype_age, celltype)
names(df_final_AV) <- c("exprn", "celltype_age", "celltype")
#print(celltype)
#print(table(df_final_AV$celltype_age))
a <- as.data.frame(table(subset_celltype_gene$celltype.stim))
b <- as.data.frame(table(subset_celltype$celltype.stim))
total = merge(b,a, by="Var1", all=TRUE)
total[is.na(total)] <- 0
names(total) <- c("Sample", "Total", "Expressed")
total$Not_expressed <- total$Total - total$Expressed
cells_pie <- total[,c("Expressed", "Not_expressed")]
cells_pie$Celltype <- total$Sample
cells_pie <- cells_pie[c("Celltype", "Expressed", "Not_expressed")]
cells_pie_long <- gather(cells_pie, condition, cells, Expressed:Not_expressed, factor_key=TRUE)
level_order <- c(paste0(Celltype,"_3m"), paste0(Celltype,"_18m"))
#########Violin Plot ######################
vlnplot_n <- df_final_AV %>%
ggplot(aes(y=exprn, x = factor(celltype_age, level = level_order), fill = factor(celltype_age, level = level_order))) +
geom_violin(position = position_dodge(0.8), alpha=0.5) +
scale_fill_manual(values=c("springgreen", "deepskyblue")) + geom_jitter(shape=16, position=position_jitter(0.2),cex=5) +
xlab("Sample")+ylab("Expression")+
theme(plot.margin = unit(c(4,4,20,10), "cm"))+ggtitle(gene_name)+theme(plot.title = element_text(hjust = 0.5, size=80))+ theme(axis.text.x=element_text(size=60, angle=90), axis.text.y = element_text(size = 60), axis.title.x = element_text(size = 60), axis.title.y = element_text(size = 60), legend.title=element_text(size=60), legend.text=element_text(size=60))+theme(legend.key.size = unit(2, "cm"))
###### Pie chart #########################
cells_pie_long$Celltype = factor(cells_pie_long$Celltype, levels=level_order)
pie_all <- ggplot(data = cells_pie_long, aes(x="", y=cells, fill=condition)) + geom_bar(stat = "identity", position = "fill") + theme_void() +
coord_polar(theta = "y", start=0) + facet_wrap(~Celltype, nrow=1) + scale_fill_grey() + theme(legend.position = "right", legend.title = element_text(size = 10), legend.text = element_text(size = 10), legend.key.size = unit(1, 'cm'))
pie_only <- ggplot(data = cells_pie_long, aes(x="", y=cells, fill=condition)) + geom_bar(stat = "identity", position = "fill") + theme_void() +
coord_polar(theta = "y", start=0) + facet_wrap(~Celltype, nrow=1) + scale_fill_grey() + theme(legend.position = "right", legend.title = element_text(size = 10), legend.text = element_text(size = 10), legend.key.size = unit(1, 'cm'))
Celltype = gsub("/", "_or_", Celltype)
outfile = paste0("./output/","VlnPlot_Pie_chart_",Celltype, "_", gene_name, "_new", ".pdf")
outfile1 = paste0("./output/","VlnPlot_",Celltype, "_", gene_name, "_new", ".pdf")
outfile2 = paste0("./output/","Pie_chart_",Celltype, "_", gene_name, "_new", ".pdf")
pdf(file = outfile, width=40, height=40)
print(ggdraw()+draw_plot(vlnplot_n, x=0, y=0, width=1, height=1)+draw_plot(pie_all, x=0.2, y=-0.4, width=0.35))
dev.off()
pdf(file = outfile1, width=40, height=40)
print(ggdraw()+draw_plot(vlnplot_n, x=0, y=0, width=1, height=1))
dev.off()
pdf(file = outfile2, width=30, height=10)
print(ggdraw()+draw_plot(pie_only))
dev.off()
}
```
############################################
```{r , include=TRUE, echo=TRUE}
args <- commandArgs(TRUE)
if (length(args) < 3) {
stop("Seurat object, gene name and Celltype must be supplied (input file).n", call.=FALSE)
} else
{
print("Generating Violinplots and Piecharts for gene expression for young and old ...")
}
mydata = readRDS(args[1]) #spleen_seurat_subset_rename_ident_Oct24.rds
gene_name = args[2] #Jun
Celltype = args[3] #Tcells
#Call function to generate plots.
gene_Vln_Pie(mydata, gene_name, Celltype)
```