-
Notifications
You must be signed in to change notification settings - Fork 0
/
application.r
83 lines (71 loc) · 2.51 KB
/
application.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
library(PopGenome)
genome =
readVCF("AGC_refHC_bialSNP_AC2_2DPGQ.3L_V2.CHRcode2.DRYAD.vcf.gz",
10000,"4",1,45000000, include.unknown=TRUE)
# Define the populations
Aquad = c("SRS408143", "SRS408145", "SRS408151", "SRS408155",
"SRS408966", "SRS408969", "SRS408972", "SRS408973", "SRS408983",
"SRS420578")
Amela = c("SRS408142", "SRS408185", "SRS408994")
Ameru = c("SRS408186", "SRS408187", "SRS408967", "SRS408974",
"SRS408992", "SRS410266","SRS410284", "SRS410286", "SRS410290",
"SRS420577")
# Define the outgroup
Chris = c("CHRISTYI")
# Set the populations
genome = set.populations(genome, list(Aquad,Amela,Ameru),diploid=TRUE)
# Set the outgroup
genome = set.outgroup(genome, Chris, diploid=TRUE)
#
inversionStart=1.45e07
inversionEnd=3.575e07
# Split the data into 50kb windows
slide =
sliding.window.transform(genome,jump=5000,width=5000,type=2)
# Perform the introgression statistics
slide = introgression.stats(slide, lambda=5)
df = slide@df
df_theta = slide@df_theta
df_bayes = slide@df_bayes
pos = seq(1, 45000000, by= 10000)-1
RES = cbind(pos, df, df_theta, df_bayes)
colnames(RES) = c("position","df","df_theta","df_BF")
#plots
library(ggplot2)
library(reshape)
RES = as.data.frame(RES)
p1 = ggplot(RES, aes(x=position, y=df, color=df_BF)) +
geom_point() +
geom_smooth(method = "loess") +
scale_color_gradient(low="black", high="orange")+
ylim(c(-1,1)) +
ylab("Estimate of introgression (df)") +
xlab("Genomic position (kB)")
p2 = ggplot(RES, aes(x=position, y=df_theta, color=df_BF)) +
geom_point() +
geom_smooth(method = "loess") +
scale_color_gradient(low="black", high="orange")+
ylim(c(0,1)) +
ylab("Estimate of introgression (df_theta)") +
xlab("Genomic position (kB)")
p3 = ggplot(RES, aes(x=position, y=df_BF, color=df_BF)) +
geom_point() +
geom_smooth(method = "loess") +
scale_color_gradient(low="black", high="orange")+
#ylim(c(0,1)) +
ylab("Evidence of introgression (df_BF)") +
xlab("Genomic position (kB)")
# ---------------- #
p4 = ggplot(RES, aes(y=df, x=df_BF, color=df_BF)) +
geom_point() +
scale_color_gradient(low="black", high="orange")+
#ylim(c(0,200)) +
ylab("Estimate of introgression (df)") +
xlab("Evidence of introgression (df_BF)")
library(gridExtra)
library(grid)
library(lattice)
library(ggpubr)
#grid.arrange(grobs=list(p1,p2,p3,p4), ncol=2)
ggarrange(p1, p2, p3, p4, ncol=2, nrow=2,
common.legend = TRUE, legend="bottom")