-
Notifications
You must be signed in to change notification settings - Fork 0
/
Sex_bias_plotsetc_withSB_from_Vi_and_Md_females.R
398 lines (304 loc) · 23.2 KB
/
Sex_bias_plotsetc_withSB_from_Vi_and_Md_females.R
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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
library(ggplot2)
library(pheatmap)
library(cowplot)
library(grid)
library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)
library(cowplot)
library(gtable)
library(RColorBrewer)
library(pvclust)
library("VennDiagram")
library("SuperExactTest")
library(raster)
theme_update(plot.title = element_text(hjust = 0.5))
print (sessionInfo())
# R version 3.5.1 (2018-07-02)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS 10.14.3
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# attached base packages:
# [1] grid stats graphics grDevices utils datasets methods base
# other attached packages:
# [1] raster_2.8-4 sp_1.3-1 SuperExactTest_1.0.4 VennDiagram_1.6.20 futile.logger_1.4.3 pvclust_2.0-0 RColorBrewer_1.1-2
# [8] gtable_0.2.0 lattice_0.20-38 gridExtra_2.3 cowplot_0.9.3 pheatmap_1.0.10 ggplot2_3.1.0
# loaded via a namespace (and not attached):
# [1] Rcpp_1.0.0 pillar_1.3.0 compiler_3.5.1 formatR_1.5 plyr_1.8.4 bindr_0.1.1 futile.options_1.0.1
# [8] tibble_1.4.2 pkgconfig_2.0.2 rlang_0.3.0.1 bindrcpp_0.2.2 withr_2.1.2 dplyr_0.7.8 tidyselect_0.2.5
# [15] glue_1.3.0 R6_2.3.0 purrr_0.2.5 lambda.r_1.2.3 magrittr_1.5 scales_1.0.0 codetools_0.2-15
# [22] assertthat_0.2.0 colorspace_1.3-2 lazyeval_0.2.1 munsell_0.5.0 crayon_1.3.4
# >
################################################################################################################################################
##### read in data (virgins)
dat_Tbi_WB_RBBH_SB_asex_VI = read.csv("Output/DE_joined_Virgin/Tbi_WB_RBBH_disp_allsepar_Vi_SB_asex.csv")
dat_Tce_WB_RBBH_SB_asex_VI = read.csv("Output/DE_joined_Virgin/Tce_WB_RBBH_disp_allsepar_Vi_SB_asex.csv")
dat_Tcm_WB_RBBH_SB_asex_VI = read.csv("Output/DE_joined_Virgin/Tcm_WB_RBBH_disp_allsepar_Vi_SB_asex.csv")
dat_Tpa_WB_RBBH_SB_asex_VI = read.csv("Output/DE_joined_Virgin/Tpa_WB_RBBH_disp_allsepar_Vi_SB_asex.csv")
dat_Tps_WB_RBBH_SB_asex_VI = read.csv("Output/DE_joined_Virgin/Tps_WB_RBBH_disp_allsepar_Vi_SB_asex.csv")
head(dat_Tbi_WB_RBBH_SB_asex_VI)
################################################################################################################################################
##### read in data (mated)
dat_Tbi_WB_RBBH_SB_asex_MT = read.csv("Output/DE_joined/Tbi_WB_RBBH_disp_allsepar_SB_asex.csv")
dat_Tce_WB_RBBH_SB_asex_MT = read.csv("Output/DE_joined/Tce_WB_RBBH_disp_allsepar_SB_asex.csv")
dat_Tcm_WB_RBBH_SB_asex_MT = read.csv("Output/DE_joined/Tcm_WB_RBBH_disp_allsepar_SB_asex.csv")
dat_Tpa_WB_RBBH_SB_asex_MT = read.csv("Output/DE_joined/Tpa_WB_RBBH_disp_allsepar_SB_asex.csv")
dat_Tps_WB_RBBH_SB_asex_MT = read.csv("Output/DE_joined/Tps_WB_RBBH_disp_allsepar_SB_asex.csv")
head(dat_Tbi_WB_RBBH_SB_asex_VI, n = 20)
head(dat_Tbi_WB_RBBH_SB_asex_MT, n = 20)
#################################################################################################################################################
###### refine sex bias by FC (not log FC)
FC_cutoff = 2
dat_Tbi_WB_RBBH_SB_asex_MT$Tbi_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tbi_WB_RBBH_SB_asex_MT$Tbi_WB_log2FC_SB * dat_Tbi_WB_RBBH_SB_asex_MT$Tbi_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tbi_WB_RBBH_SB_asex_MT$Tbi_WB_sexbias), "Unbiased")
dat_Tce_WB_RBBH_SB_asex_MT$Tce_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tce_WB_RBBH_SB_asex_MT$Tce_WB_log2FC_SB * dat_Tce_WB_RBBH_SB_asex_MT$Tce_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tce_WB_RBBH_SB_asex_MT$Tce_WB_sexbias), "Unbiased")
dat_Tcm_WB_RBBH_SB_asex_MT$Tcm_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tcm_WB_RBBH_SB_asex_MT$Tcm_WB_log2FC_SB * dat_Tcm_WB_RBBH_SB_asex_MT$Tcm_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tcm_WB_RBBH_SB_asex_MT$Tcm_WB_sexbias), "Unbiased")
dat_Tpa_WB_RBBH_SB_asex_MT$Tpa_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tpa_WB_RBBH_SB_asex_MT$Tpa_WB_log2FC_SB * dat_Tpa_WB_RBBH_SB_asex_MT$Tpa_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tpa_WB_RBBH_SB_asex_MT$Tpa_WB_sexbias), "Unbiased")
dat_Tps_WB_RBBH_SB_asex_MT$Tps_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tps_WB_RBBH_SB_asex_MT$Tps_WB_log2FC_SB * dat_Tps_WB_RBBH_SB_asex_MT$Tps_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tps_WB_RBBH_SB_asex_MT$Tps_WB_sexbias), "Unbiased")
#################################################################################################################################################
###### refine sex bias by FC (not log FC)
FC_cutoff = 2
dat_Tbi_WB_RBBH_SB_asex_VI$Tbi_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tbi_WB_RBBH_SB_asex_VI$Tbi_WB_log2FC_SB * dat_Tbi_WB_RBBH_SB_asex_VI$Tbi_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tbi_WB_RBBH_SB_asex_VI$Tbi_WB_sexbias), "Unbiased")
dat_Tce_WB_RBBH_SB_asex_VI$Tce_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tce_WB_RBBH_SB_asex_VI$Tce_WB_log2FC_SB * dat_Tce_WB_RBBH_SB_asex_VI$Tce_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tce_WB_RBBH_SB_asex_VI$Tce_WB_sexbias), "Unbiased")
dat_Tcm_WB_RBBH_SB_asex_VI$Tcm_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tcm_WB_RBBH_SB_asex_VI$Tcm_WB_log2FC_SB * dat_Tcm_WB_RBBH_SB_asex_VI$Tcm_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tcm_WB_RBBH_SB_asex_VI$Tcm_WB_sexbias), "Unbiased")
dat_Tpa_WB_RBBH_SB_asex_VI$Tpa_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tpa_WB_RBBH_SB_asex_VI$Tpa_WB_log2FC_SB * dat_Tpa_WB_RBBH_SB_asex_VI$Tpa_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tpa_WB_RBBH_SB_asex_VI$Tpa_WB_sexbias), "Unbiased")
dat_Tps_WB_RBBH_SB_asex_VI$Tps_WB_sexbias2 <- ifelse(2^(sqrt(dat_Tps_WB_RBBH_SB_asex_VI$Tps_WB_log2FC_SB * dat_Tps_WB_RBBH_SB_asex_VI$Tps_WB_log2FC_SB)) > FC_cutoff , as.character(dat_Tps_WB_RBBH_SB_asex_VI$Tps_WB_sexbias), "Unbiased")
########## Merge together - make new datasets
dat_Tbi_WB_RBBH_SB_asex_VIMT <- merge(dat_Tbi_WB_RBBH_SB_asex_VI, dat_Tbi_WB_RBBH_SB_asex_MT, by="genename",sort = TRUE)
dat_Tce_WB_RBBH_SB_asex_VIMT <- merge(dat_Tce_WB_RBBH_SB_asex_VI, dat_Tce_WB_RBBH_SB_asex_MT, by="genename",sort = TRUE)
dat_Tcm_WB_RBBH_SB_asex_VIMT <- merge(dat_Tcm_WB_RBBH_SB_asex_VI, dat_Tcm_WB_RBBH_SB_asex_MT, by="genename",sort = TRUE)
dat_Tpa_WB_RBBH_SB_asex_VIMT <- merge(dat_Tpa_WB_RBBH_SB_asex_VI, dat_Tpa_WB_RBBH_SB_asex_MT, by="genename",sort = TRUE)
dat_Tps_WB_RBBH_SB_asex_VIMT <- merge(dat_Tps_WB_RBBH_SB_asex_VI, dat_Tps_WB_RBBH_SB_asex_MT, by="genename",sort = TRUE)
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT <- paste(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2.x, dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2.y, sep = "_" )
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT <- paste(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2.x, dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2.y, sep = "_" )
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT <- paste(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2.x, dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2.y, sep = "_" )
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT <- paste(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2.x, dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2.y, sep = "_" )
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT <- paste(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2.x, dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2.y, sep = "_" )
levels(as.factor(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT ))
levels(as.factor(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT ))
levels(as.factor(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT ))
levels(as.factor(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT ))
levels(as.factor(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT ))
###### class genes as sex-biased if they are SB in BOTH datasets
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_both <-
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "female_biased_female_biased", "female_biased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "female_biased_Unbiased" , "Unbiased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "male_biased_male_biased" , "male_biased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "male_biased_Unbiased" , "Unbiased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "Unbiased_female_biased" , "Unbiased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "Unbiased_male_biased" , "Unbiased",
ifelse(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_VIMT == "Unbiased_Unbiased", "Unbiased", "ERROR"
)))))))
dat_Tbi_WB_RBBH_SB_asex_MT_both <-
as.data.frame(cbind(
as.character(dat_Tbi_WB_RBBH_SB_asex_VIMT$genename),
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_log2FC_SB.y,
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_FDR_SB.y,
as.character(dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_sexbias2_both),
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_log2FC_SA.y,
dat_Tbi_WB_RBBH_SB_asex_VIMT$Tbi_WB_FDR_SA.y
))
colnames(dat_Tbi_WB_RBBH_SB_asex_MT_both) <- c("genename", "Tbi_WB_log2FC_SB", "Tbi_WB_FDR_SB", "Tbi_WB_sexbias2", "Tbi_WB_log2FC_SA", "Tbi_WB_FDR_SA")
dat_Tbi_WB_RBBH_SB_asex_MT_both$genename <- as.character(dat_Tbi_WB_RBBH_SB_asex_MT_both$genename)
dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_log2FC_SB <- as.numeric(as.character(dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_log2FC_SB))
dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_FDR_SB <- as.numeric(as.character(dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_FDR_SB ))
dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_log2FC_SA <- as.numeric(as.character(dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_log2FC_SA))
dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_FDR_SA <- as.numeric(as.character(dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_FDR_SA))
###### class genes as sex-biased if they are SB in BOTH datasets
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_both <-
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "female_biased_female_biased", "female_biased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "female_biased_Unbiased" , "Unbiased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "male_biased_male_biased" , "male_biased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "male_biased_Unbiased" , "Unbiased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "Unbiased_female_biased" , "Unbiased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "Unbiased_male_biased" , "Unbiased",
ifelse(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_VIMT == "Unbiased_Unbiased", "Unbiased", "ERROR"
)))))))
dat_Tce_WB_RBBH_SB_asex_MT_both <-
as.data.frame(cbind(
as.character(dat_Tce_WB_RBBH_SB_asex_VIMT$genename),
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_log2FC_SB.y,
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_FDR_SB.y,
as.character(dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_sexbias2_both),
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_log2FC_SA.y,
dat_Tce_WB_RBBH_SB_asex_VIMT$Tce_WB_FDR_SA.y
))
colnames(dat_Tce_WB_RBBH_SB_asex_MT_both) <- c("genename", "Tce_WB_log2FC_SB", "Tce_WB_FDR_SB", "Tce_WB_sexbias2", "Tce_WB_log2FC_SA", "Tce_WB_FDR_SA")
dat_Tce_WB_RBBH_SB_asex_MT_both$genename <- as.character(dat_Tce_WB_RBBH_SB_asex_MT_both$genename)
dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_log2FC_SB <- as.numeric(as.character(dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_log2FC_SB))
dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_FDR_SB <- as.numeric(as.character(dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_FDR_SB ))
dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_log2FC_SA <- as.numeric(as.character(dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_log2FC_SA))
dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_FDR_SA <- as.numeric(as.character(dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_FDR_SA))
###### class genes as sex-biased if they are SB in BOTH datasets
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_both <-
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "female_biased_female_biased", "female_biased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "female_biased_Unbiased" , "Unbiased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "male_biased_male_biased" , "male_biased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "male_biased_Unbiased" , "Unbiased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "Unbiased_female_biased" , "Unbiased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "Unbiased_male_biased" , "Unbiased",
ifelse(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_VIMT == "Unbiased_Unbiased", "Unbiased", "ERROR"
)))))))
dat_Tcm_WB_RBBH_SB_asex_MT_both <-
as.data.frame(cbind(
as.character(dat_Tcm_WB_RBBH_SB_asex_VIMT$genename),
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_log2FC_SB.y,
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_FDR_SB.y,
as.character(dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_sexbias2_both),
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_log2FC_SA.y,
dat_Tcm_WB_RBBH_SB_asex_VIMT$Tcm_WB_FDR_SA.y
))
colnames(dat_Tcm_WB_RBBH_SB_asex_MT_both) <- c("genename", "Tcm_WB_log2FC_SB", "Tcm_WB_FDR_SB", "Tcm_WB_sexbias2", "Tcm_WB_log2FC_SA", "Tcm_WB_FDR_SA")
dat_Tcm_WB_RBBH_SB_asex_MT_both$genename <- as.character(dat_Tcm_WB_RBBH_SB_asex_MT_both$genename)
dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_log2FC_SB <- as.numeric(as.character(dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_log2FC_SB))
dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_FDR_SB <- as.numeric(as.character(dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_FDR_SB ))
dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_log2FC_SA <- as.numeric(as.character(dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_log2FC_SA))
dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_FDR_SA <- as.numeric(as.character(dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_FDR_SA))
###### class genes as sex-biased if they are SB in BOTH datasets
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_both <-
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "female_biased_female_biased", "female_biased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "female_biased_Unbiased" , "Unbiased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "male_biased_male_biased" , "male_biased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "male_biased_Unbiased" , "Unbiased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "Unbiased_female_biased" , "Unbiased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "Unbiased_male_biased" , "Unbiased",
ifelse(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_VIMT == "Unbiased_Unbiased", "Unbiased", "ERROR"
)))))))
dat_Tpa_WB_RBBH_SB_asex_MT_both <-
as.data.frame(cbind(
as.character(dat_Tpa_WB_RBBH_SB_asex_VIMT$genename),
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_log2FC_SB.y,
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_FDR_SB.y,
as.character(dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_sexbias2_both),
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_log2FC_SA.y,
dat_Tpa_WB_RBBH_SB_asex_VIMT$Tpa_WB_FDR_SA.y
))
colnames(dat_Tpa_WB_RBBH_SB_asex_MT_both) <- c("genename", "Tpa_WB_log2FC_SB", "Tpa_WB_FDR_SB", "Tpa_WB_sexbias2", "Tpa_WB_log2FC_SA", "Tpa_WB_FDR_SA")
dat_Tpa_WB_RBBH_SB_asex_MT_both$genename <- as.character(dat_Tpa_WB_RBBH_SB_asex_MT_both$genename)
dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_log2FC_SB <- as.numeric(as.character(dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_log2FC_SB))
dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_FDR_SB <- as.numeric(as.character(dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_FDR_SB ))
dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_log2FC_SA <- as.numeric(as.character(dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_log2FC_SA))
dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_FDR_SA <- as.numeric(as.character(dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_FDR_SA))
###### class genes as sex-biased if they are SB in BOTH datasets
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_both <-
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "female_biased_female_biased", "female_biased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "female_biased_Unbiased" , "Unbiased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "male_biased_male_biased" , "male_biased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "male_biased_Unbiased" , "Unbiased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "Unbiased_female_biased" , "Unbiased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "Unbiased_male_biased" , "Unbiased",
ifelse(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_VIMT == "Unbiased_Unbiased", "Unbiased", "ERROR"
)))))))
dat_Tps_WB_RBBH_SB_asex_MT_both <-
as.data.frame(cbind(
as.character(dat_Tps_WB_RBBH_SB_asex_VIMT$genename),
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_log2FC_SB.y,
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_FDR_SB.y,
as.character(dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_sexbias2_both),
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_log2FC_SA.y,
dat_Tps_WB_RBBH_SB_asex_VIMT$Tps_WB_FDR_SA.y
))
colnames(dat_Tps_WB_RBBH_SB_asex_MT_both) <- c("genename", "Tps_WB_log2FC_SB", "Tps_WB_FDR_SB", "Tps_WB_sexbias2", "Tps_WB_log2FC_SA", "Tps_WB_FDR_SA")
dat_Tps_WB_RBBH_SB_asex_MT_both$genename <- as.character(dat_Tps_WB_RBBH_SB_asex_MT_both$genename)
dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_log2FC_SB <- as.numeric(as.character(dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_log2FC_SB))
dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_FDR_SB <- as.numeric(as.character(dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_FDR_SB ))
dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_log2FC_SA <- as.numeric(as.character(dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_log2FC_SA))
dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_FDR_SA <- as.numeric(as.character(dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_FDR_SA))
dat_Tbi_WB_RBBH_SB_asex_MT_both$Tbi_WB_sexbias <- rep("NA", length(dat_Tbi_WB_RBBH_SB_asex_MT_both[,1]))
dat_Tce_WB_RBBH_SB_asex_MT_both$Tce_WB_sexbias <- rep("NA", length(dat_Tce_WB_RBBH_SB_asex_MT_both[,1]))
dat_Tcm_WB_RBBH_SB_asex_MT_both$Tcm_WB_sexbias <- rep("NA", length(dat_Tcm_WB_RBBH_SB_asex_MT_both[,1]))
dat_Tpa_WB_RBBH_SB_asex_MT_both$Tpa_WB_sexbias <- rep("NA", length(dat_Tpa_WB_RBBH_SB_asex_MT_both[,1]))
dat_Tps_WB_RBBH_SB_asex_MT_both$Tps_WB_sexbias <- rep("NA", length(dat_Tps_WB_RBBH_SB_asex_MT_both[,1]))
################################################################################################################################################
############### boxplot SB
################################################################################################################################################
make_long_table = function(dat_Tbi,dat_Tce,dat_Tcm,dat_Tpa,dat_Tps,tiss){
sp_u = "Tbi"
df_t1 = as.data.frame(cbind(
eval(parse(text=paste(dat_Tbi,'$', sp_u, '_' ,tiss,'_log2FC_SA',sep=''))),
as.character(eval(parse(text=paste(dat_Tbi,'$', sp_u, '_' ,tiss,'_sexbias',sep='')))),
as.character(eval(parse(text=paste(dat_Tbi,'$', sp_u, '_' ,tiss,'_sexbias2',sep=''))))))
df_t1$sp = rep(sp_u, length(df_t1[,1]))
df_t1$group = paste(df_t1$sp, df_t1$V2, sep = "_")
df_t1$group2 = paste(df_t1$sp, df_t1$V3, sep = "_")
sp_u = "Tce"
df_t2 = as.data.frame(cbind(
eval(parse(text=paste(dat_Tce,'$', sp_u, '_' ,tiss,'_log2FC_SA',sep=''))),
as.character(eval(parse(text=paste(dat_Tce,'$', sp_u, '_' ,tiss,'_sexbias',sep='')))),
as.character(eval(parse(text=paste(dat_Tce,'$', sp_u, '_' ,tiss,'_sexbias2',sep=''))))))
df_t2$sp = rep(sp_u, length(df_t2[,1]))
df_t2$group = paste(df_t2$sp, df_t2$V2, sep = "_")
df_t2$group2 = paste(df_t2$sp, df_t2$V3, sep = "_")
sp_u = "Tcm"
df_t3 = as.data.frame(cbind(
eval(parse(text=paste(dat_Tcm,'$', sp_u, '_' ,tiss,'_log2FC_SA',sep=''))),
as.character(eval(parse(text=paste(dat_Tcm,'$', sp_u, '_' ,tiss,'_sexbias',sep='')))),
as.character(eval(parse(text=paste(dat_Tcm,'$', sp_u, '_' ,tiss,'_sexbias2',sep=''))))))
df_t3$sp = rep(sp_u, length(df_t3[,1]))
df_t3$group = paste(df_t3$sp, df_t3$V2, sep = "_")
df_t3$group2 = paste(df_t3$sp, df_t3$V3, sep = "_")
sp_u = "Tpa"
df_t4 = as.data.frame(cbind(
eval(parse(text=paste(dat_Tpa,'$', sp_u, '_' ,tiss,'_log2FC_SA',sep=''))),
as.character(eval(parse(text=paste(dat_Tpa,'$', sp_u, '_' ,tiss,'_sexbias',sep='')))),
as.character(eval(parse(text=paste(dat_Tpa,'$', sp_u, '_' ,tiss,'_sexbias2',sep=''))))))
df_t4$sp = rep(sp_u, length(df_t4[,1]))
df_t4$group = paste(df_t4$sp, df_t4$V2, sep = "_")
df_t4$group2 = paste(df_t4$sp, df_t4$V3, sep = "_")
sp_u = "Tps"
df_t5 = as.data.frame(cbind(
eval(parse(text=paste(dat_Tps,'$', sp_u, '_' ,tiss,'_log2FC_SA',sep=''))),
as.character(eval(parse(text=paste(dat_Tps,'$', sp_u, '_' ,tiss,'_sexbias',sep='')))),
as.character(eval(parse(text=paste(dat_Tps,'$', sp_u, '_' ,tiss,'_sexbias2',sep=''))))))
df_t5$sp = rep(sp_u, length(df_t5[,1]))
df_t5$group = paste(df_t5$sp, df_t5$V2, sep = "_")
df_t5$group2 = paste(df_t5$sp, df_t5$V3, sep = "_")
df_t = rbind(df_t1,df_t2,df_t3,df_t4,df_t5)
colnames(df_t) <- c("log2FC_SA", "sexbias", "sexbias2", "sp", "SBgroup", "SBgroup2")
df_t$log2FC_SA <- as.numeric(as.character(df_t$log2FC_SA))
df_t$SBgroup <- as.factor(df_t$SBgroup)
df_t$SBgroup2 <- as.factor(df_t$SBgroup2)
df_t$SBgroup_ord <- ordered(df_t$SBgroup, levels = c(
"Tbi_female_biased", "Tbi_male_biased", "Tbi_Unbiased",
"Tce_female_biased", "Tce_male_biased", "Tce_Unbiased",
"Tps_female_biased", "Tps_male_biased", "Tps_Unbiased",
"Tcm_female_biased", "Tcm_male_biased", "Tcm_Unbiased",
"Tpa_female_biased", "Tpa_male_biased", "Tpa_Unbiased"
))
df_t$SBgroup2_ord <- ordered(df_t$SBgroup2, levels = c(
"Tbi_female_biased", "Tbi_male_biased", "Tbi_Unbiased",
"Tce_female_biased", "Tce_male_biased", "Tce_Unbiased",
"Tps_female_biased", "Tps_male_biased", "Tps_Unbiased",
"Tcm_female_biased", "Tcm_male_biased", "Tcm_Unbiased",
"Tpa_female_biased", "Tpa_male_biased", "Tpa_Unbiased"
))
df_t$sp_ord <- ordered(df_t$sp, levels = c("Tbi", "Tce", "Tps", "Tcm", "Tpa" ))
return(df_t)
}
dat_ALL_WB_RBBH_SB_asex_long_MT_both <- make_long_table("dat_Tbi_WB_RBBH_SB_asex_MT_both", "dat_Tce_WB_RBBH_SB_asex_MT_both", "dat_Tcm_WB_RBBH_SB_asex_MT_both", "dat_Tpa_WB_RBBH_SB_asex_MT_both", "dat_Tps_WB_RBBH_SB_asex_MT_both", "WB")
dat_ALL_WB_RBBH_SB_asex_long_MT_both$sexbias2_ord <- ordered(dat_ALL_WB_RBBH_SB_asex_long_MT_both$sexbias2, levels = c("female_biased", "Unbiased", "male_biased"))
head(dat_ALL_WB_RBBH_SB_asex_long_MT_both)
##### plot boxplot
WB_SA_box_s2_noout_both <- ggplot(dat_ALL_WB_RBBH_SB_asex_long_MT_both, aes(sp_ord, log2FC_SA)) +
theme_classic() +
geom_boxplot(aes(fill = factor(sexbias2_ord)),position=position_dodge(0.6), width = 0.5, outlier.size = 0, lwd = 0.5, fatten = 1, outlier.shape = NA) +
coord_cartesian(ylim=c(-3.5,3.5)) +
ylab ("log2 fold change in asexual females") +
xlab ("Gene-class by species") +
scale_fill_manual(values=c("firebrick2", "grey", "royalblue2")) + geom_hline(yintercept = 0) + labs(fill='') +
ggtitle("Whole-body | SB class = FDR <0.05 FC > 2\nSB in Vi AND Md")
## output
dir.create("Output/SB_in_Vi_and_Md")
## pdf
pdf(paste("Output/SB_in_Vi_and_Md/Sex_asex_WB_SBinViANDMd_boxplot_s2_FDR005FC2.pdf",sep = ""), width = 6, height = 8)
WB_SA_box_s2_noout_both
dev.off()
getwd() ## wh
########################################################################################################################################################################
####### output session info
print (sessionInfo())
writeLines(capture.output(sessionInfo()), "sex_bias_plotsetc_withVIfemales_sessionInfo.txt")