-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathScript_08_Functional annotation of aDMCs and aVMCs.Rmd
193 lines (151 loc) · 8.66 KB
/
Script_08_Functional annotation of aDMCs and aVMCs.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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
---
title: "Functional annotation of aDMCs and aVMCs"
author: "Yunfeng"
date: "1/16/2023"
output: html_document
---
#Script information
This script was written to annotate XCI related features of aDMCs and aVMCs.
```{r}
#set library path.
.libPaths("~/researchdrive/yliu/Rlibs")
```
```{r}
# load data
load(file = "Output_05_Catalogue_aDMCs_validated.RData")
load(file = "Output_06_Catalogue_aVMCs_validated.RData")
load(file = "Output_07_450K_chrX_annotation.RData")
```
```{r}
# Load necessary libraries
library(tidyverse)
```
```{r}
### Prepare aDMCs annotation data.
# aDMCs in females
ann450K_aDMCs_females<-ann450K_chrX[row.names(dat_aDMCs_females_validated),]
table(ann450K_aDMCs_females$CGI_Feature)
# CGI:8 Shore:7 non-CGI:18
# aDMCs in males
ann450K_aDMCs_males<-ann450K_chrX[row.names(dat_aDMCs_males_validated),]
table(ann450K_aDMCs_males$CGI_Feature)
# CGI:83 Shore:67 non-CGI:40
```
```{r}
### Prepare aDMCs annotation data.
# aVMCs in females
ann450K_aVMCs_females<-ann450K_chrX[row.names(dat_aVMCs_females_validated),]
table(ann450K_aVMCs_females$CGI_Feature)
# CGI:473 Shore:299 non-CGI:215
# aVMCs in males
ann450K_aVMCs_males<-ann450K_chrX[row.names(dat_aVMCs_males_validated),]
table(ann450K_aVMCs_males$CGI_Feature)
# CGI:2 Shore:21 non-CGI:14
```
```{r}
# CGI annotation of aDMCs and aVMCs in females and males.
ann450K_aDMCs_females<-as.data.frame(ann450K_aDMCs_females)
ann450K_aDMCs_males<-as.data.frame(ann450K_aDMCs_males)
ann450K_aVMCs_females<-as.data.frame(ann450K_aVMCs_females)
ann450K_aVMCs_males<-as.data.frame(ann450K_aVMCs_males)
ann450K_chrX_reference<-as.data.frame(ann450K_chrX)
ann450K_aDMCs_aVMCs_females_males<-rbind(ann450K_chrX_reference,ann450K_aDMCs_females,ann450K_aVMCs_females,ann450K_chrX_reference,ann450K_aDMCs_males,ann450K_aVMCs_males)
ann450K_aDMCs_aVMCs_females_males$symbol<-c(rep("Control",9777),rep("aDMCs",33),rep("aVMCs",987),rep("Control",9777),rep("aDMCs",316),rep("aVMCs",37))
ann450K_aDMCs_aVMCs_females_males$sex<-c(rep("Females",10797),rep("Males",10130))
ann450K_aDMCs_aVMCs_females_males$symbol<-factor(ann450K_aDMCs_aVMCs_females_males$symbol,levels = unique(ann450K_aDMCs_aVMCs_females_males$symbol))
ann450K_aDMCs_aVMCs_females_males$sex<-factor(ann450K_aDMCs_aVMCs_females_males$sex,levels = c("Males","Females"))
ann450K_aDMCs_aVMCs_females_males$CGI_Feature<-factor(ann450K_aDMCs_aVMCs_females_males$CGI_Feature,levels = c("non-CGI","Shore","CGI"))
save(ann450K_aDMCs_aVMCs_females_males,file = "CGI_XCI_feature_annotation.RData")
tiff("Output_08_CGI annotation of aDMCs and aVMCs.tiff", units="in", width=15, height=9, res=600,compression ='lzw')
ggplot(ann450K_aDMCs_aVMCs_females_males) +
geom_bar(aes(x=symbol,fill=CGI_Feature), position = "fill") +
facet_grid(~sex,switch ="y",scales = "fixed") +
labs(fill = "CGI_Feature") +
scale_fill_manual(values=c("#3CAEA3","#F6D55C","#ED553B"))+
xlab("") +
ylab("Fraction") +
theme_bw()+
theme(panel.spacing.x = unit(.5, "cm"),legend.position = "bottom",axis.text.x = element_text(angle = 45,hjust = 1,size = 34),axis.text.y = element_text(size = 34),axis.title.y = element_text(size = 34),legend.text = element_text(size = 29),legend.title = element_blank(),strip.text.x = element_text(size = 30))
dev.off()
```
```{r}
# XCI annotation based on TSS meta calls.
ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls<-gsub("escapes XCI","Escape XCI",ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls)
ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls<-gsub("subject to XCI","Subject to XCI",ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls)
ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls<-gsub("variably escapes from XCI","Variably escape XCI",ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls)
ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls<-factor(ann450K_aDMCs_aVMCs_females_males$XCI.Status.meta.calls,levels = c("Escape XCI","Variably escape XCI","Subject to XCI"))
ann450K_aDMCs_aVMCs_females_males<-filter(ann450K_aDMCs_aVMCs_females_males,Distance.meta.calls<=2000)
xtabs(~XCI.Status.meta.calls+symbol+sex,ann450K_aDMCs_aVMCs_females_males)
tiff("Output_08_XCI annotation_meta_calls.tiff", units="in", width=15, height=9, res=600,compression ='lzw')
ggplot(ann450K_aDMCs_aVMCs_females_males) +
geom_bar(aes(x=symbol,fill=XCI.Status.meta.calls), position = "fill") +
facet_grid(~sex,switch ="y",scales = "fixed") +
labs(fill = "XCI status") +
scale_fill_manual(values=c("#3CAEA3","#F6D55C","#ED553B"))+
xlab("") +
ylab("Fraction") +
theme_bw()+
theme(panel.spacing.x = unit(.5, "cm"),legend.position = "bottom",axis.text.x = element_text(angle = 45,hjust = 1,size = 34),axis.text.y = element_text(size = 34),axis.title.y = element_text(size = 34),legend.text = element_text(size = 29),legend.title = element_blank(),strip.text.x = element_text(size = 30))
dev.off()
```
```{r}
# Calculate the mean methylation of aDMCs in females and males.
load(file = "Output_01_chrX_betas.RData")
aDMCs_females<-t(assay(betas_Xfemales))[,rownames(dat_aDMCs_females_validated)]
aDMCs_females<-as.data.frame(aDMCs_females)
dat_mean_aDMCs_females<-as.data.frame(rowMeans(t(aDMCs_females),na.rm = TRUE))
colnames(dat_mean_aDMCs_females)<-"Beta value"
table(dat_mean_aDMCs_females$`Beta value`<0.25) # 13
table(dat_mean_aDMCs_females$`Beta value`>0.7) # 1
aDMCs_males<-t(assay(betas_Xmales))[,rownames(dat_aDMCs_males_validated)]
aDMCs_males<-as.data.frame(aDMCs_males)
dat_mean_aDMCs_males<-as.data.frame(rowMeans(t(aDMCs_males),na.rm = TRUE))
colnames(dat_mean_aDMCs_males)<-"Beta value"
table(dat_mean_aDMCs_males$`Beta value`<0.25) # 161
table(dat_mean_aDMCs_males$`Beta value` >0.7) # 63
```
```{r}
# Calculate the mean methylation of aVMCs in females and males.
aVMCs_females<-t(assay(betas_Xfemales))[,rownames(dat_aVMCs_females_validated)]
aVMCs_females<-as.data.frame(aVMCs_females)
dat_mean_aVMCs_females<-as.data.frame(rowMeans(t(aVMCs_females),na.rm = TRUE))
colnames(dat_mean_aVMCs_females)<-"Beta value"
table(dat_mean_aVMCs_females$`Beta value`<0.25) # 112
table(dat_mean_aVMCs_females$`Beta value`>0.7) # 117
aVMCs_males<-t(assay(betas_Xmales))[,rownames(dat_aVMCs_males_validated)]
aVMCs_males<-as.data.frame(aVMCs_males)
dat_mean_aVMCs_males<-as.data.frame(rowMeans(t(aVMCs_males),na.rm = TRUE))
colnames(dat_mean_aVMCs_males)<-"Beta value"
table(dat_mean_aVMCs_males$`Beta value`<0.25) # 2
table(dat_mean_aVMCs_males$`Beta value`>0.7) # 4
```
```{r}
# Calculate the mean methylation of 9777 CpGs in females and males.
# Xfemales
dat_mean_Xfemales<-as.data.frame(rowMeans(assay(betas_Xfemales),na.rm = TRUE))
colnames(dat_mean_Xfemales)<-"Beta value"
# Xmales
dat_mean_Xmales<-as.data.frame(rowMeans(assay(betas_Xmales),na.rm = TRUE))
colnames(dat_mean_Xmales)<-"Beta value"
```
```{r}
dat_mean_aDMCs_aVMCs_females_males<-rbind(dat_mean_Xfemales,dat_mean_aDMCs_females,dat_mean_aVMCs_females,dat_mean_Xmales,dat_mean_aDMCs_males,dat_mean_aVMCs_males)
dat_mean_aDMCs_aVMCs_females_males$symbol<-c(rep("Control",9777),rep("aDMCs",33),rep("aVMCs",987),rep("Control",9777),rep("aDMCs",316),rep("aVMCs",37))
dat_mean_aDMCs_aVMCs_females_males$sex<-c(rep("Females",10797),rep("Males",10130))
dat_mean_aDMCs_aVMCs_females_males$Methylation.level<-with(dat_mean_aDMCs_aVMCs_females_males,ifelse(`Beta value` < 0.25, "Hypomethylated", ifelse(`Beta value` >= 0.25 & `Beta value` < 0.7,"Intermediate methylated","Hypermethylated")))
dat_mean_aDMCs_aVMCs_females_males$symbol<-factor(dat_mean_aDMCs_aVMCs_females_males$symbol,levels = c("Control","aDMCs","aVMCs"))
dat_mean_aDMCs_aVMCs_females_males$Methylation.level<-factor(dat_mean_aDMCs_aVMCs_females_males$Methylation.level,levels = c("Hypomethylated","Hypermethylated","Intermediate methylated"))
dat_mean_aDMCs_aVMCs_females_males$sex<-factor(dat_mean_aDMCs_aVMCs_females_males$sex,levels = c("Males","Females"))
save(dat_mean_aDMCs_aVMCs_females_males,file = "dat_mean_aDMCs_aVMCs_females_males.RData")
tiff("Output_08_Methylation level of aDMCs and aVMCs.tiff", units="in", width=15, height=9, res=600,compression ='lzw')
ggplot(dat_mean_aDMCs_aVMCs_females_males) +
geom_bar(aes(x=symbol,fill=Methylation.level),position = "fill") +
labs(fill = "Methylation level") +
facet_grid(~sex,switch ="y",scales = "fixed") +
scale_fill_manual(values=c("#3CAEA3","#F6D55C","#ED553B"))+
xlab("") +
ylab("Fraction") +
theme_bw()+
theme(panel.spacing.x = unit(.5, "cm"),legend.position = "bottom",axis.text.x = element_text(angle = 45,hjust = 1,size = 34),axis.text.y = element_text(size = 34),axis.title.y = element_text(size = 34),legend.text = element_text(size = 29),legend.title = element_blank(),strip.text.x = element_text(size = 30))
dev.off()
```