forked from KatiePelletier/WingShapeBSA
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathds_wild_fst_plottingByPop.R
More file actions
164 lines (136 loc) · 7.74 KB
/
ds_wild_fst_plottingByPop.R
File metadata and controls
164 lines (136 loc) · 7.74 KB
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
#plotting by population fst plots for wild pools (left vs right)
library(data.table)
library(tidyverse)
#my plotting functions
source("KP_genomescan_source.R")
ds.wild <- fread("../Data/wild_lvr_nosync_100window.fst")
ds.wild <- populaionFst_cleanup(ds.wild, x = c('CMO', 'FVW13', 'FVW14', 'PHO'))
#Reordering and numbering the chr for plotting.
fixed_ds.wild <- chrNumbering(ds.wild)
chrlabel.ds <- middleChr(fixed_ds.wild)
#Quick Plots
#ID suggested adding either the mean fst to each chr or a smoothing spline and increasing alpha
cmo.xMean <- mean(fixed_ds.wild[chr == "X", CMO])
cmo.2LMean <- mean(fixed_ds.wild[chr == "2L", CMO])
cmo.2RMean <- mean(fixed_ds.wild[chr == "2R", CMO])
cmo.3LMean <- mean(fixed_ds.wild[chr == "3L", CMO])
cmo.3RMean <- mean(fixed_ds.wild[chr == "3R", CMO])
cmo.4Mean <- mean(fixed_ds.wild[chr == "4", CMO])
#can resuse these all for all of them (because the numbering is actually common)
cmo.xstart <- min(fixed_ds.wild[chr == "X", number])
cmo.2Lstart <- min(fixed_ds.wild[chr == "2L", number])
cmo.2Rstart <- min(fixed_ds.wild[chr == "2R", number])
cmo.3Lstart <- min(fixed_ds.wild[chr == "3L", number])
cmo.3Rstart <- min(fixed_ds.wild[chr == "3R", number])
cmo.4start <- min(fixed_ds.wild[chr == "4", number])
CMO <- ggplot(data = fixed_ds.wild, aes(x=number, y=CMO, color=chr)) +
geom_point(size=1, show.legend = F, alpha = 0.2) +
theme(panel.background = element_blank()) +
xlab("Chromosome") +
ylab("meanFst") +
scale_x_discrete(limits=c(chrlabel.ds),
labels = c("X","2L", "2R", '3L', '3R', '4')) +
scale_colour_manual(values=c('black', 'grey46', 'black', 'grey46', 'black','grey46')) +
geom_segment(x = cmo.xstart, xend = cmo.2Lstart, y = cmo.xMean, yend = cmo.xMean , col = "red") +
geom_segment(x = cmo.2Lstart, xend = cmo.2Rstart, y = cmo.2LMean, yend = cmo.2LMean , col = "red") +
geom_segment(x = cmo.2Rstart, xend = cmo.3Lstart, y = cmo.2RMean, yend = cmo.2RMean , col = "red") +
geom_segment(x = cmo.3Lstart, xend = cmo.3Rstart, y = cmo.3LMean, yend = cmo.3LMean , col = "red") +
geom_segment(x = cmo.3Rstart, xend = cmo.4start, y = cmo.3RMean, yend = cmo.3RMean , col = "red") +
geom_segment(x = cmo.4start, xend = max(fixed_ds.wild$number), y = cmo.4Mean, yend = cmo.4Mean , col = "red") +
ylim(0,0.75) +
theme(text = element_text(size=15),
axis.text.x= element_text(size=12),
axis.text.y= element_text(size=12),
panel.border = element_rect(colour = "black",
fill=NA, size=0.5))
f13.xMean <- mean(fixed_ds.wild[chr == "X", FVW13])
f13.2LMean <- mean(fixed_ds.wild[chr == "2L", FVW13])
f13.2RMean <- mean(fixed_ds.wild[chr == "2R", FVW13])
f13.3LMean <- mean(fixed_ds.wild[chr == "3L", FVW13])
f13.3RMean <- mean(fixed_ds.wild[chr == "3R", FVW13])
f13.4Mean <- mean(fixed_ds.wild[chr == "4", FVW13])
FVW13 <- ggplot(data = fixed_ds.wild, aes(x=number, y=FVW13, color=chr)) +
geom_point(size=1, show.legend = F, alpha = 0.2) +
theme(panel.background = element_blank()) +
xlab("Chromosome") +
ylab("meanFst") +
scale_x_discrete(limits=c(chrlabel.ds),
labels = c("X","2L", "2R", '3L', '3R', '4')) +
scale_colour_manual(values=c('black', 'grey46', 'black', 'grey46', 'black','grey46')) +
ylim(0,0.75) +
geom_segment(x = cmo.xstart, xend = cmo.2Lstart, y = f13.xMean, yend = f13.xMean , col = "red") +
geom_segment(x = cmo.2Lstart, xend = cmo.2Rstart, y = f13.2LMean, yend = f13.2LMean , col = "red") +
geom_segment(x = cmo.2Rstart, xend = cmo.3Lstart, y = f13.2RMean, yend = f13.2RMean , col = "red") +
geom_segment(x = cmo.3Lstart, xend = cmo.3Rstart, y = f13.3LMean, yend = f13.3LMean , col = "red") +
geom_segment(x = cmo.3Rstart, xend = cmo.4start, y = f13.3RMean, yend = f13.3RMean , col = "red") +
geom_segment(x = cmo.4start, xend = max(fixed_ds.wild$number), y = f13.4Mean, yend = f13.4Mean , col = "red") +
theme(text = element_text(size=15),
axis.text.x= element_text(size=12),
axis.text.y= element_text(size=12),
panel.border = element_rect(colour = "black",
fill=NA, size=0.5))
f14.xMean <- mean(fixed_ds.wild[chr == "X", FVW14])
f14.2LMean <- mean(fixed_ds.wild[chr == "2L", FVW14])
f14.2RMean <- mean(fixed_ds.wild[chr == "2R", FVW14])
f14.3LMean <- mean(fixed_ds.wild[chr == "3L", FVW14])
f14.3RMean <- mean(fixed_ds.wild[chr == "3R", FVW14])
f14.4Mean <- mean(fixed_ds.wild[chr == "4", FVW14])
FVW14 <- ggplot(data = fixed_ds.wild, aes(x=number, y=FVW14, color=chr)) +
geom_point(size=1, show.legend = F, alpha = 0.2) +
theme(panel.background = element_blank()) +
xlab("Chromosome") +
ylab("meanFst") +
scale_x_discrete(limits=c(chrlabel.ds),
labels = c("X","2L", "2R", '3L', '3R', '4')) +
scale_colour_manual(values=c('black', 'grey46', 'black', 'grey46', 'black','grey46')) +
geom_segment(x = cmo.xstart, xend = cmo.2Lstart, y = f14.xMean, yend = f14.xMean , col = "red") +
geom_segment(x = cmo.2Lstart, xend = cmo.2Rstart, y = f14.2LMean, yend = f14.2LMean , col = "red") +
geom_segment(x = cmo.2Rstart, xend = cmo.3Lstart, y = f14.2RMean, yend = f14.2RMean , col = "red") +
geom_segment(x = cmo.3Lstart, xend = cmo.3Rstart, y = f14.3LMean, yend = f14.3LMean , col = "red") +
geom_segment(x = cmo.3Rstart, xend = cmo.4start, y = f14.3RMean, yend = f14.3RMean , col = "red") +
geom_segment(x = cmo.4start, xend = max(fixed_ds.wild$number), y = f14.4Mean, yend = f14.4Mean , col = "red") +
ylim(0,0.75) +
theme(text = element_text(size=15),
axis.text.x= element_text(size=12),
axis.text.y= element_text(size=12),
panel.border = element_rect(colour = "black",
fill=NA, size=0.5))
pho.xMean <- mean(fixed_ds.wild[chr == "X", PHO])
pho.2LMean <- mean(fixed_ds.wild[chr == "2L", PHO])
pho.2RMean <- mean(fixed_ds.wild[chr == "2R", PHO])
pho.3LMean <- mean(fixed_ds.wild[chr == "3L", PHO])
pho.3RMean <- mean(fixed_ds.wild[chr == "3R", PHO])
pho.4Mean <- mean(fixed_ds.wild[chr == "4", PHO])
PHO <- ggplot(data = fixed_ds.wild, aes(x=number, y=PHO, color=chr)) +
geom_point(size=1, show.legend = F, alpha = 0.2) +
theme(panel.background = element_blank()) +
xlab("Chromosome") +
ylab("meanFst") +
scale_x_discrete(limits=c(chrlabel.ds),
labels = c("X","2L", "2R", '3L', '3R', '4')) +
scale_colour_manual(values=c('black', 'grey46', 'black', 'grey46', 'black','grey46')) +
ylim(0,0.75) +
geom_segment(x = cmo.xstart, xend = cmo.2Lstart, y = pho.xMean, yend = pho.xMean , col = "red") +
geom_segment(x = cmo.2Lstart, xend = cmo.2Rstart, y = pho.2LMean, yend = pho.2LMean , col = "red") +
geom_segment(x = cmo.2Rstart, xend = cmo.3Lstart, y = pho.2RMean, yend = pho.2RMean , col = "red") +
geom_segment(x = cmo.3Lstart, xend = cmo.3Rstart, y = pho.3LMean, yend = pho.3LMean , col = "red") +
geom_segment(x = cmo.3Rstart, xend = cmo.4start, y = pho.3RMean, yend = pho.3RMean , col = "red") +
geom_segment(x = cmo.4start, xend = max(fixed_ds.wild$number), y = pho.4Mean, yend = pho.4Mean , col = "red") +
theme(text = element_text(size=15),
axis.text.x= element_text(size=12),
axis.text.y= element_text(size=12),
panel.border = element_rect(colour = "black",
fill=NA, size=0.5))
library(cowplot)
png("../Figures/wild_ds_byPoolFst.png")
plot_grid(FVW13, FVW14, CMO, PHO, nrow = 2, labels = c("FVW13", "FVW14", "CMO", "PHO"))
dev.off()
mean(fixed_ds.wild$CMO)
mean(fixed_ds.wild$PHO)
mean(fixed_ds.wild$FVW13)
mean(fixed_ds.wild$FVW14)
###Maybe I want to look at some of those higher windows? I dont know? The couple hits in fvw13 look sketch.
#taking some of the high hits in fvw13 quickly.
#these all have pretty low coverage.
f13.top <- filter(fixed_ds.wild, FVW13 >= 0.45)
dim(f13.top)