-
Notifications
You must be signed in to change notification settings - Fork 0
/
Overall_CpG_analysis.Rmd
96 lines (72 loc) · 3.17 KB
/
Overall_CpG_analysis.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
---
title: "Untitled"
author: "Shelly Trigg"
date: "12/22/2019"
output: html_document
---
```{r}
library(ggplot2)
library(colorRamps)
```
read in dataa
```{r}
CpGs <- read.table("allc_5x_CpG.txt", header = FALSE, sep = "\t",stringsAsFactors = FALSE)
```
add column names
```{r}
colnames(CpGs) <- c("Sample","mCpG_5x", "CpG_5x", "perc_meth" )
```
add column for lane info
```{r}
CpGs$Lane <- gsub(".*_L","L", CpGs$Sample)
```
remove lane info from sample name
```{r}
CpGs$Sample <- gsub("_L.*","", CpGs$Sample)
```
read in meta data
```{r}
meta_data <- read.csv("Sample.Info.csv", stringsAsFactors = FALSE)
```
create sample.ID column in CpG data frame
```{r}
CpGs$Sample.ID <- gsub("EPI-","EPI_", CpGs$Sample)
CpGs$Sample.ID <- gsub("_S.*", "", CpGs$Sample.ID)
```
merge meta_data with CpG data
```{r}
CpGs <- merge(CpGs, meta_data, by = "Sample.ID")
```
create column with primary and secondary treatment info
```{r}
CpGs$Treatment <- paste(CpGs$Initial.Treatment, CpGs$Secondary.Treatment, sep = "_")
```
plot CpGs
```{r}
#plot mCpGs
jpeg("5x_num_mCpG_boxplot.jpg", width = 6, height = 4, units = "in", res = 300)
ggplot(CpGs, aes(TimePoint,mCpG_5x/1000000))+geom_boxplot() + theme_bw() + theme(axis.text.x = element_text(size = 10, hjust= 1,angle = 45)) + ylab("millions of mCpGs")
dev.off()
#plot percent methylation
jpeg("5xCovPercMeth_boxplot.jpg", width = 6, height = 4, units = "in", res = 300)
ggplot(CpGs, aes(x=CpGs$TimePoint,CpGs$perc_meth))+geom_boxplot()+geom_point(aes(color=Treatment),position = position_jitter(0.25)) + theme_bw() + theme(axis.text.x = element_text(size = 10, hjust= 1,angle = 60)) + ylab("CpG % methylation") + scale_color_manual(values = colorRamps::rgb.tables(length(unique(CpGs$Treatment))))
jpeg("5xCovPercMeth_facet_boxplots.jpg", width = 6, height = 4, units = "in", res = 300)
ggplot(CpGs, aes(x=CpGs$Treatment,CpGs$perc_meth))+geom_boxplot()+geom_point(aes(color=Treatment),position = position_jitter(0.25)) + theme_bw() + theme(axis.text.x = element_text(size = 10, hjust= 1,angle = 60)) + ylab("CpG % methylation") + scale_color_manual(values = colorRamps::rgb.tables(length(unique(CpGs$Treatment)))) + facet_wrap(~TimePoint)
dev.off()
```
```{r}
jpeg("5xCovPercMeth_boxplots.jpg", width = 6, height = 4, units = "in", res = 300)
ggplot(CpGs, aes(x=paste(CpGs$TimePoint,CpGs$Treatment),CpGs$perc_meth))+geom_boxplot()+geom_point(aes(color=Treatment),position = position_jitter(0.25)) + theme_bw() + theme(axis.text.x = element_text(size = 10, hjust= 1,angle = 60)) + ylab("CpG % methylation") + scale_color_manual(values = colorRamps::rgb.tables(length(unique(CpGs$Treatment))))
dev.off()
```
run a test for differences between groups
```{r}
CpGaov <- aov(mCpG_5x~group, data = CpGs)
CpGaov2w <- aov(mCpG_5x~temperature*salinity, data = CpGs)
CpGaov <- aov(Tperc_meth~group, data = CpGs0.75x)
CpGglm <- glm(matrix(c(mCpG_5x0.75X, CpG_5x0.75X), ncol=2) ~ group, family = binomial, data = CpGs0.75x[-grep("CTRL",CpGs0.75x$group),] )
CpGglm0 <- glm(matrix(c(mCpG_5x0.75X, CpG_5x0.75X), ncol=2) ~ 1, family = binomial, data = CpGs0.75x[-grep("CTRL",CpGs0.75x$group),] )
p.val<- pchisq(deviance(CpGglm0)-deviance(CpGglm),
df.residual(CpGglm0)-df.residual(CpGglm),
lower.tail=FALSE)
```