-
Notifications
You must be signed in to change notification settings - Fork 0
/
heterogeneity_bulkRNA.Rmd
67 lines (50 loc) · 2.13 KB
/
heterogeneity_bulkRNA.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
---
title: "Selected examples of inter-individual heterogeneity"
author: "Anthony Hung"
date: "2021-01-21"
output: workflowr::wflow_html
editor_options:
chunk_output_type: console
---
# Introduction
Here, we plot the expression values of three genes for which we identified inter-individual hetereogeneity in the response to mechanical strain.
# Load data and libraries
```{r load data}
library(ggplot2)
#load in V2 results from voom in DE analysis
v2 <- readRDS("output/voom_results.rds")
```
# Plot a dot plot for each of the three genes
```{r plotting function}
# this function takes in the v2 object and the name of a gene (Ensmbl ID) and outputs a plot
plot_gene <- function(v, g) {
# v - An EList object containing log2 counts per million
# g - character vector of a single gene
stopifnot(class(v) == "EList",
is.character(g), length(g) == 1)
library("tidyr")
single_gene <- v$E[g, ]
single_gene_long <- as.data.frame(single_gene)
colnames(single_gene_long) <- "log2cpm"
single_gene_long$sample <- rownames(single_gene_long)
single_gene_long <- separate(single_gene_long, col = "sample", sep = "_",
into = c("Individual", "Replicate", "Treatment"))
single_gene_long$Treatment <- gsub('S', 'Strain', single_gene_long$Treatment)
single_gene_long$Treatment <- gsub('U', 'Control', single_gene_long$Treatment)
single_gene_long$Treatment <- factor(single_gene_long$Treatment, levels(factor(single_gene_long$Treatment))[c(1,2)])
return(ggplot(single_gene_long, aes(x = Individual, y = log2cpm, fill = Individual)) +
labs(title = g, y = expression("Expression level (" * log[2] * " cpm)")) +
geom_dotplot(binaxis = "y", stackdir = "center", key_glyph = 'rect', dotsize = 0.75) +
facet_wrap(~Treatment, strip.position = 'bottom') +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult=1)) +
theme(axis.title.x=element_blank())
)
}
```
# Make plots
```{r}
gene_list <- c("ENSG00000157227", "ENSG00000181789", "ENSG00000120699")
for (gene in gene_list) {
print(plot_gene(v2, as.character(gene)))
}
```