/
Fig1.Rmd
143 lines (111 loc) · 4.48 KB
/
Fig1.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
---
title: "Figure 1 for Deleterious Alleles paper"
output: pdf_document
---
```{r setup, include=TRUE}
#setwd("~/Documents/Github/pvpDiallel/")
knitr::opts_knit$set(root.dir=normalizePath('../'))
```
To reproduce the figures, we should set `pvpDiallel` as the root path, i.e. `setwd("~/Documents/Github/pvpDiallel/")`. And then use `knitr` package to knit a pdf file. Or simplely click `Knit PDF` icon in `RStudio`. Note, to generate each panel into seperate pdf files, we should turn `getpdf` into `TRUE` (i.e. `getpdf=TRUE`) when calling the plotting functions.
First of all, determine fond size and set the getpdf option:
```{r}
#par(mar=c(5,4,4,2))
par(font= 2, font.lab= 2, font.axis= 2)
fs <- 1.6 # times bigger than default
ht= 6; wt= 6 #figure height and weight
getpdf <- TRUE # get figures in seperated pdf [TRUE] or not [FALSE]
```
# Figure 1a
```{r, eval=TRUE, fig.width=ht, fig.height=wt}
plotloh <- function(getpdf, outfile, ...){
trait <- read.csv("data/hyb_heterosis.csv")
trait$pMPH <- abs(trait$pMPH*100)
bymed2 <- with(trait, reorder(trait, pMPH, median))
boxplot(pMPH ~ bymed2, data=trait,
xlab = "", ylab= "MPH (100%)", col="antiquewhite3",
...)
if(getpdf == TRUE){
pdf(outfile, width=ht, height=wt)
par(mar=c(5,5,4,2))
boxplot(pMPH ~ bymed2, data=trait,
xlab = "", ylab= "Mid-Parent Heterosis (100%)", col="antiquewhite3",
...)
dev.off()
}
}
plotloh(getpdf, outfile="graphs/Fig1a.pdf",
main="", cex.axis=fs, cex.lab=fs, las=2)
```
# Figure 1b
```{r, eval=TRUE, fig.width=ht, fig.height=wt}
plot_load <- function(getpdf, outfile, ...){
dres <- read.table("data/sup_deleterious_hmp3.txt", header=T)
boxplot(DR ~ geno*ordered, data=dres, border=c("darkgreen","darkblue"), col="grey", ...)
#text(cex=fs, x=x-.25, y=-1.25, H2$Traits, xpd=TRUE, srt=45, pos=2)
if(getpdf == TRUE){
pdf(outfile, width=ht, height=wt)
par(mar=c(5,5,4,2))
boxplot(DR ~ geno*ordered, data=dres, border=c("darkgreen","darkblue"), col="grey", ...)
box()
dev.off()
}
}
plot_load(getpdf, outfile="graphs/Fig1d_v3.pdf",
names=c("LR", "MZ", "LR", "MZ", "LR", "MZ"),
xlab="", ylab="Deleterious Load per bp",
main="", cex.axis=fs, cex.lab=fs)
```
# Figure 1c
```{r, eval=TRUE, fig.width=ht, fig.height=wt}
plotReg <- function(getpdf, outfile, ...){
snptab <- read.csv("cache/daf_gerp2.csv")
snptab <- snptab[order(snptab$GERP2),]
plx <- predict(loess(snptab$meandaf ~ snptab$GERP2), se=T)
x <- snptab$GERP2
y <- snptab$meandaf
plot(x, y, ...)
lines(snptab$GERP2, plx$fit, col="cornflowerblue", lwd=2)
lines(snptab$GERP2, plx$fit - qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
lines(snptab$GERP2, plx$fit + qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
if(getpdf == TRUE){
pdf(outfile, width=wt, height=ht)
par(mar=c(5,5,4,2))
plot(x, y, ...)
lines(snptab$GERP2, plx$fit, col="cornflowerblue", lwd=2)
lines(snptab$GERP2, plx$fit - qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
lines(snptab$GERP2, plx$fit + qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
dev.off()
}
}
plotReg(getpdf, outfile="graphs/Fig1b_v3.pdf",
pch=16, col="antiquewhite3", xlab="GERP Score", ylab="Allele Frequency",
main="", cex.axis=fs, cex.lab=fs)
```
# Figure 1d
```{r, eval=TRUE, fig.width=wt, fig.height=ht}
##########################
plotmgerp <- function(mgerp, getpdf, outfile, ...){
#mgerp <- read.csv("cache/mgerp_cm.csv")
mgerp <- mgerp[order(mgerp$gen),]
plx <- predict(loess(mgerp$mgerp ~ mgerp$gen), se=T)
x <- mgerp$gen
y <- mgerp$mgerp
plot(x, y, ...)
lines(mgerp$gen, plx$fit, col="cornflowerblue", lwd=2)
lines(mgerp$gen, plx$fit - qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
lines(mgerp$gen, plx$fit + qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
if(getpdf == TRUE){
pdf(outfile, width=wt, height=ht)
par(mar=c(5,5,4,2))
plot(x, y, ...)
lines(mgerp$gen, plx$fit, col="cornflowerblue", lwd=2)
lines(mgerp$gen, plx$fit - qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
lines(mgerp$gen, plx$fit + qt(0.975,plx$df)*plx$se, lty=2, lwd=2, col="black")
dev.off()
}
}
mgerp <- read.csv("cache/mgerp_cm.csv")
plotmgerp(mgerp, getpdf=TRUE, outfile="graphs/Fig1d.pdf",
pch=16, col="antiquewhite3", xlab="Recombination Rate (cM/Mb)", ylab="GERP Score",
main="", cex.axis=fs, cex.lab=fs)
```