QTL Mapping of ERA8

Title: Exome sequencing of bulked segregants identified a novel TaMKK3-A allele linked to the wheat ERA8 ABA-hypersensitive germination phenotype

Authors: Shantel A. Martinez, Oluwayesi Shorinola, Samantha Conselman, Deven See, Daniel Z. Skinner, Cristobal Uauy, and Camille M. Steber

This repository contains all datasets for the paper, Martinez et al. 2020, and can be used to enable reproducibility of the QTL Mapping and figures created.

The datasets generated during and/or analysed during the current study are available in the GitHub repository, github.com/shantel-martinez/ERA8-Mapping.

Load All Data

To download all of the datasets used in this analysis:

Once downloaded, the environment can be loaded:

load("ERA8-Mapping-Data.RData")

Load Libraries

library(qtl)
library(tidyverse)
font<-element_text(face = "bold",  size = 16) #ggplot aesthetic
# library(kableExtra)
# library(knitr)
library(DT) 
sketch = htmltools::withTags(table(    #table aesthetic
  class = 'display',
  thead(
    tr(
      th(rowspan = 2, 'Species'),
      th(colspan = 2, 'Sepal'),
      th(colspan = 2, 'Petal')
    ),
    tr(
      lapply(rep(c('Length', 'Width'), 2), th)
    )
  )
)) 

Zak/ZakERA8 BC3F2:3

After-ripening timecourse

ZParABA <- read.csv("./ZZ8ARTC_E1.csv", header = TRUE, na.string="NA") 

ZParABA$ABA_uM <- as.character(ZParABA$ABA_uM)
ZParABA$Geno <- as.character(ZParABA$Geno)
ParABAf <-
  ZParABA %>%
    filter(Geno != "Zak..2.2.6", Geno != "95.3.91.4" , Geno != "Zak..2.2.8")

ParABA2 <- subset(ParABAf, ParABAf$ABA_uM == "0" )
ggplot(data = ParABA2 , aes(x=AR_wk, y=PG_d5, color = Geno)) +
  geom_line(size=1.5) +
  geom_point(aes(size = 3))+
  geom_errorbar(aes(ymin=PG_d5-PG_Std, ymax=PG_d5+PG_Std), width=1,
                 position=position_dodge(0.05))+
  ylim(0,100)+ theme_bw()+ theme(axis.text = font, axis.title = font)+
  # labs(title = "Figure S2a", subtitle = "0uM ABA (No Hormone)", y = "Percent Germination", x = "Weeks After-ripened")+
  labs(title = "0uM ABA (No Hormone)", y = "Percent Germination", x = "Weeks After-ripened")+ theme(legend.position="none") +
  scale_color_manual(values=c("#440154","#A2CD5A"))

ggsave("20190923_ZZ8_Parents_MES.png", width = 5, height = 5, units = "in")

ParABA2 <- subset(ParABAf, ParABAf$ABA_uM == "2")
ggplot(data = ParABA2 , aes(x=AR_wk, y=PG_d5, color = Geno)) +
  geom_line(size=1.5) +
  geom_point(aes(size = 3))+
  geom_errorbar(aes(ymin=PG_d5-PG_Std, ymax=PG_d5+PG_Std), width=1,
                 position=position_dodge(0.05))+
  ylim(0,100)+ theme_bw()+ theme(axis.text = font, axis.title = font)+
  # labs(title = "Figure S2b", subtitle = "2uM ABA", y = "Percent Germination", x = "Weeks After-ripened")+
  labs(title = "2uM ABA", y = "Percent Germination", x = "Weeks After-ripened")+ theme(legend.position="none") +
  scale_color_manual(values=c("#440154","#A2CD5A"))

ggsave("20190923_ZZ8_Parents_ABA.png", width = 5, height = 5, units = "in")

The ZZ8 dataset contains Zak/ZakERA8 BC3 F2:3 genotypic data, phenotypic data, and marker positions on chromosome 4A.

ZZ8 <- read.cross("csv",".","./ZZ8_SNP1-SNP_30_all.csv",na.strings=c("-","NA","."),BC.gen=3, F.gen=2)

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(3)F(2) cross
## 
##     No. individuals:    424 
## 
##     No. phenotypes:     3 
##     Percent phenotyped: 100 97.9 97.9 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      11 
##     No. markers:        11 
##     Percent genotyped:  92.6 
##     Genotypes (%):      AA:27.9  AB:51.9  BB:20.2  not BB:0.0  not AA:0.0

Phenotypes include: Percent Germination (PG) and Germination Index (GI)

Composite Interval Mapping

data_2 <- calc.genoprob(ZZ8, step=0)
out.cimPGavg <- cim(data_2, pheno.col="PG",map.function=c("kosambi"),method=c("hk"))
out.cimPGavg$Trait <- ("PG") 
out.cimGIavg <- cim(data_2, pheno.col="GI",map.function=c("kosambi"),method=c("hk"))
out.cimGIavg$Trait <- ("GI")

Determine Significance Threshhold

a <- scanone(ZZ8, method="hk", pheno.col = "PG", n.perm=1000)
b <- scanone(ZZ8, method="hk", pheno.col = "GI", n.perm=1000)
summary(a) 
summary(b)

Wrangle Results

ZZ8_all <- rbind(out.cimPGavg, out.cimGIavg)
ZZ8_all$Population <- ("all")

---

Repeat CIM for the X5.1 subset only

data5.1<-read.cross("csv",".","./ZZ8_SNP1-SNP_30_5.1.csv",na.strings=c("-","NA","."),BC.gen=3, F.gen=2)

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(3)F(2) cross
## 
##     No. individuals:    122 
## 
##     No. phenotypes:     3 
##     Percent phenotyped: 100 94.3 94.3 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      11 
##     No. markers:        11 
##     Percent genotyped:  91.7 
##     Genotypes (%):      AA:28.1  AB:53.7  BB:18.2  not BB:0.0  not AA:0.0

data_2 <- calc.genoprob(data5.1, step=0)
out.cimPGavg <- cim(data_2, pheno.col="PG",map.function=c("kosambi"),method=c("hk"))
out.cimPGavg$Trait <- ("PG") 
out.cimGIavg <- cim(data_2, pheno.col="GI",map.function=c("kosambi"),method=c("hk"))
a <- scanone(data5.1, method="hk", pheno.col = "PG", n.perm=1000)
b <- scanone(data5.1, method="hk", pheno.col = "GI", n.perm=1000)
summary(a) 
summary(b) 
ZZ8_5.1 <- rbind(out.cimPGavg, out.cimGIavg)
ZZ8_5.1$Population <- ("5.1")

---

Repeat CIM for the X5.2a subset only

data5.2a<-read.cross("csv",".","./ZZ8_SNP1-SNP_30_5.2a.csv",na.strings=c("-","NA","."),BC.gen=3, F.gen=2)

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(3)F(2) cross
## 
##     No. individuals:    242 
## 
##     No. phenotypes:     3 
##     Percent phenotyped: 100 99.2 99.2 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      11 
##     No. markers:        11 
##     Percent genotyped:  93.2 
##     Genotypes (%):      AA:27.6  AB:50.8  BB:21.5  not BB:0.0  not AA:0.0

data_2 <- calc.genoprob(data5.2a, step=0)
out.cimPGavg <- cim(data_2, pheno.col="PG",map.function=c("kosambi"),method=c("hk"))
out.cimPGavg$Trait <- ("PG") 
out.cimGIavg <- cim(data_2, pheno.col="GI",map.function=c("kosambi"),method=c("hk"))
out.cimGIavg$Trait <- ("GI")
a <- scanone(data5.2a, method="hk", pheno.col = "PG", n.perm=1000)
b <- scanone(data5.2a, method="hk", pheno.col = "GI", n.perm=1000)
summary(a) 
summary(b) 
ZZ8_5.2a <- rbind(out.cimPGavg, out.cimGIavg)
ZZ8_5.2a$Population <- ("5.2a")

Repeat CIM for the X5.2b subset only

data5.2b <- read.cross("csv",".","./ZZ8_SNP1-SNP_30_5.2b.csv",na.strings=c("-","NA","."),BC.gen=3, F.gen=2)

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(3)F(2) cross
## 
##     No. individuals:    60 
## 
##     No. phenotypes:     3 
##     Percent phenotyped: 100 100 100 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      11 
##     No. markers:        11 
##     Percent genotyped:  92 
##     Genotypes (%):      AA:28.8  AB:52.6  BB:18.6  not BB:0.0  not AA:0.0

data_2 <- calc.genoprob(data5.2b, step=0)
out.cimPGavg <- cim(data_2, pheno.col="PG",map.function=c("kosambi"),method=c("hk"))
out.cimPGavg$Trait <- ("PG") 
out.cimGIavg <- cim(data_2, pheno.col="GI",map.function=c("kosambi"),method=c("hk"))
out.cimGIavg$Trait <- ("GI")
a <- scanone(data5.2b, method="hk", pheno.col = "PG", n.perm=1000)
b <- scanone(data5.2b, method="hk", pheno.col = "GI", n.perm=1000)
summary(a) 
summary(b) 
ZZ8_5.2b <- rbind(out.cimPGavg, out.cimGIavg)
ZZ8_5.2b$Population <- ("5.2b")

Graph Results

QTL Graph for Germination Index (GI)

ZZ8_QTL <- rbind(ZZ8_all, ZZ8_5.1, ZZ8_5.2b, ZZ8_5.2a)
GI <- ZZ8_QTL[(ZZ8_QTL$Trait %in% "GI"),]

ggplot(data = GI , aes(x=GI$pos, y=GI$lod)) + 
  geom_line(aes(colour=GI$Population), size=1)+geom_hline(aes(yintercept = 1.85))+
  theme_bw()+ theme(axis.text = font,  axis.title = font)+ 
  labs(x = "Chromosome 4A ERA8 Region (cM)", y = "LOD")+
  scale_color_manual(values=c("#78c679", "#006837", "#31a354","#A2CD5A")) + theme(legend.position="none") +
  ggsave("ZZ8_4A_QTL.png", width = 3, height = 3, units = "in")

Louise/ZakERA8 RIL F6:7

After-ripening timecourse

ParABA <- read.csv("LZ8ARTC_2014.csv", header = TRUE, na.string="NA") 

ParABA$ABA_uM <- as.character(ParABA$ABA_uM)
ParABA$Geno <- as.character(ParABA$Geno)
ParABAf <-
  ParABA %>%
    filter(Geno != "ERA8.1", Geno != "ERA8.3" , Geno != "ERA8.4")

ParABA2 <- subset(ParABAf, ParABAf$ABA_uM == "2" & ParABAf$AR_wk != 7)
ggplot(data = ParABA2 , aes(x=AR_wk, y=PG_d5, color = Geno)) +
  geom_line(size=1.5) +
  geom_point(aes(size = 3))+
  geom_errorbar(aes(ymin=PG_d5-PG_Std, ymax=PG_d5+PG_Std), width=1,
                 position=position_dodge(0.05))+
  ylim(0,100)+ theme_bw()+ theme(axis.text = font, axis.title = font)+
  # labs(title = "Figure 5a")+
  scale_color_manual(values=c("#440154","#008B8B", "#A2CD5A"))

ggsave("20190313_1-6ARWhole_Parents.png", width = 7, height = 7, units = "in")


ParABAc <- subset(ParABA, ParABA$AR_wk == 7)
target <- c("0", "2", "5", "10")
require(gdata)
ParABAc$ABA_uM <- reorder.factor(ParABAc$ABA_uM, new.order=target)
ggplot(data = ParABAc , aes(x=ABA_uM, y=PG_d5, color = Geno, shape = ABA_uM)) +
  geom_point(aes(size = 4))+
  geom_errorbar(aes(ymin=PG_d5-PG_Std, ymax=PG_d5+PG_Std), width=.2,
                 position=position_dodge(0.05))+
  ylim(0,100)+ theme_bw()+ theme(axis.text = font, axis.title = font)+
  # labs(title = "Figure 5b")+
  scale_color_manual(values=c("#440154", "#008B8B", "#A2CD5A"))+
  scale_shape_manual(values=c(1,16,17,15))  

ggsave("20190313_7ARCut_Parents.png", width = 3, height = 5, units = "in")

EMS-derived SNPs

# LZ8 <- read.cross("csv",".","./LZ8_SNP1-SNP30_MKK3.csv",na.strings=c("-","NA","."),BC.gen=0, F.gen=5)
LZ8

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(0)F(5) cross
## 
##     No. individuals:    207 
## 
##     No. phenotypes:     21 
##     Percent phenotyped: 100 100 100 100 100 100 100 88.9 88.9 88.4 88.4 
##                         82.1 82.1 84.5 84.1 84.5 84.5 84.5 83.6 99.5 99.5 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      13 
##     No. markers:        13 
##     Percent genotyped:  92.7 
##     Genotypes (%):      AA:50.1  AB:5.5  BB:44.4  not BB:0.0  not AA:0.0

data_2 <- calc.genoprob(data5.2b, step=0)
data_2 <- calc.genoprob(LZ8, step=0)

out.cimE1PGD1 <- cim(data_2, pheno.col="E1D1",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD1$Trait <- ("E1PGD1") 
out.cimE1PGD2 <- cim(data_2, pheno.col="E1D2",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD2$Trait <- ("E1PGD2")
out.cimE1PGD3 <- cim(data_2, pheno.col="E1D3",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD3$Trait <- ("E1PGD3")
out.cimE1PGD4 <- cim(data_2, pheno.col="E1D4",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD4$Trait <- ("E1PGD4")
out.cimE1PGD5 <- cim(data_2, pheno.col="E1D5",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD5$Trait <- ("E1PGD5")
out.cimE1GIavg <- cim(data_2, pheno.col="E1GI",map.function=c("kosambi"),method=c("hk"))
out.cimE1GIavg$Trait <- ("E1GIavg")

out.cimE2PGD1 <- cim(data_2, pheno.col="E2D1",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD1$Trait <- ("E2PGD1")
out.cimE2PGD2 <- cim(data_2, pheno.col="E2D2",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD2$Trait <- ("E2PGD2")
out.cimE2PGD3 <- cim(data_2, pheno.col="E2D3",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD3$Trait <- ("E2PGD3")
out.cimE2PGD4 <- cim(data_2, pheno.col="E2D4",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD4$Trait <- ("E2PGD4")
out.cimE2PGD5 <- cim(data_2, pheno.col="E2D5",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD5$Trait <- ("E2PGD5")
out.cimE2GIavg <- cim(data_2, pheno.col="E2GI",map.function=c("kosambi"),method=c("hk"))
out.cimE2GIavg$Trait <- ("E2GIavg")

out.cimE3PGD1 <- cim(data_2, pheno.col="E3D1",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD1$Trait <- ("E3PGD1")
out.cimE3PGD2 <- cim(data_2, pheno.col="E3D2",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD2$Trait <- ("E3PGD2")
out.cimE3PGD3 <- cim(data_2, pheno.col="E3D3",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD3$Trait <- ("E3PGD3")
out.cimE3PGD4 <- cim(data_2, pheno.col="E3D4",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD4$Trait <- ("E3PGD4")
out.cimE3PGD5 <- cim(data_2, pheno.col="E3D5",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD5$Trait <- ("E3PGD5")
out.cimE3GIavg <- cim(data_2, pheno.col="E3GI",map.function=c("kosambi"),method=c("hk"))
out.cimE3GIavg$Trait <- ("E3GIavg")

QTL.LZ8.ERA8 <- rbind(out.cimE1PGD1, out.cimE1PGD2, out.cimE1PGD3, out.cimE1PGD4, out.cimE1PGD5, 
                      out.cimE1GIavg, out.cimE2PGD1, out.cimE2PGD2, out.cimE2PGD3, out.cimE2PGD4, 
                      out.cimE2PGD5, out.cimE2GIavg,out.cimE3PGD1, out.cimE3PGD2, out.cimE3PGD3, 
                      out.cimE3PGD4, out.cimE3PGD5, out.cimE3GIavg)
QTL.LZ8.Day1 <- rbind(out.cimE1PGD1, out.cimE2PGD1,out.cimE3PGD1)
QTL.LZ8.GI <- rbind(out.cimE1GIavg, out.cimE2GIavg,out.cimE3GIavg)

Calculating the significant LOD score threshold for each trait

a<-scanone(LZ8, method="hk", pheno.col = "E1D1", n.perm=1000)
c<-scanone(LZ8, method="hk", pheno.col = "E1D2", n.perm=1000)
d<-scanone(LZ8, method="hk", pheno.col = "E1D3", n.perm=1000)
e<-scanone(LZ8, method="hk", pheno.col = "E1D4", n.perm=1000)
f<-scanone(LZ8, method="hk", pheno.col = "E1D5", n.perm=1000)
g<-scanone(LZ8, method="hk", pheno.col = "E1GI", n.perm=1000)
b1<-scanone(LZ8, method="hk", pheno.col = "E2D1", n.perm=1000)
c1<-scanone(LZ8, method="hk", pheno.col = "E2D2", n.perm=1000)
d1<-scanone(LZ8, method="hk", pheno.col = "E2D3", n.perm=1000)
e1<-scanone(LZ8, method="hk", pheno.col = "E2D4", n.perm=1000)
f1<-scanone(LZ8, method="hk", pheno.col = "E2D5", n.perm=1000)
g1<-scanone(LZ8, method="hk", pheno.col = "E2GI", n.perm=1000)
h<-scanone(LZ8, method="hk", pheno.col = "E3D1", n.perm=1000)
i<-scanone(LZ8, method="hk", pheno.col = "E3D2", n.perm=1000)
j<-scanone(LZ8, method="hk", pheno.col = "E3D3", n.perm=1000)
k<-scanone(LZ8, method="hk", pheno.col = "E3D4", n.perm=1000)
l<-scanone(LZ8, method="hk", pheno.col = "E3D5", n.perm=1000)
m<-scanone(LZ8, method="hk", pheno.col = "E3GI", n.perm=1000)
summary(a) #2.12 E1D1 
summary(c) #2.10
summary(d) #2.20
summary(e) #2.15
summary(f) #2.02
summary(g) #1.98
summary(b1) #2.06 E2D1
summary(c1) #2.17
summary(d1) #2.03
summary(e1) #1.96
summary(f1) #2.17
summary(g1) #2.10
summary(h) #2.01 E3D1
summary(i) #2.01
summary(j) #2.21
summary(k) #2.12
summary(l) #2.11
summary(m) #2.19

QTL Graph for GI

ggplot(QTL.LZ8.GI, aes(x=QTL.LZ8.GI$pos, y=QTL.LZ8.GI$lod))+
  geom_line(aes(colour = factor(Trait)),size = 1)+
  theme_bw()+ theme(axis.text = font,  axis.title = font)+
  geom_hline(aes(yintercept = 2.10))+ 
  labs(x = "Chromosome 4A (cM)", y = NULL) +
  ylim (0,10) +
  scale_color_manual(values=c("#96CDCD", "#20B2AA", "#00868B")) + theme(legend.position="none") + 
  ggsave("LZ8_4A_QTL.png", width = 3, height = 3, units = "in")

Genotyping-by-sequencing SNPs

# LZ8GBS <- read.cross("csv",".","./LZ8_GBS_all.csv",na.strings=c("-","NA","."),BC.gen=0, F.gen=5)
LZ8GBS

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(0)F(5) cross
## 
##     No. individuals:    209 
## 
##     No. phenotypes:     21 
##     Percent phenotyped: 99 99 99 99 99 99 99 88 88 87.6 87.6 81.3 81.3 83.7 
##                         83.3 83.7 83.7 83.7 82.8 98.6 98.6 
## 
##     No. chromosomes:    45 
##         Autosomes:      1A.1 1B.1 1B.2 1B.3 1D.1 1D.2 2A.1 2A.2 2B.1 2D.1 
##                         3A.1 3A.2 3A.3 3A.4 3B.1 3D.1 3D.2 4A.1 4A.2 4A.3 
##                         4B.1 4B.2 4D.1 5A.1 5A.2 5B.1 5B.2 5B.3 5D.1 5D.2 
##                         6A.1 6A.2 6B.1 6B.2 6D.1 6D.2 6D.3 7A.1 7A.2 7B.1 
##                         7B.2 7D.1 7D.2 7D.3 UNK.1 
## 
##     Total markers:      2234 
##     No. markers:        231 132 6 9 8 49 14 84 49 57 80 41 5 54 64 19 27 
##                         127 92 5 124 6 11 54 8 66 13 5 18 9 129 40 19 20 16 
##                         13 5 106 132 201 5 20 8 33 20 
##     Percent genotyped:  99 
##     Genotypes (%):      AA:45.2  AB:3.5  BB:51.4  not BB:0.0  not AA:0.0

View genotypic data

plot.map(LZ8GBS)

geno.image(LZ8GBS)

Omit Skewed Markers with potential segregation distortion

gt<-geno.table(LZ8GBS)
g <- gt[gt$P.value<0.05/totmar(LZ8GBS),]
datatable(g,rownames = FALSE, container = sketch, caption = htmltools::tags$caption(style = 'text-align: left;', 'Table A - Skewed GBS markers '))

drop <- c( "A15328_3A_7A_70" ,"A23596_unk_70" ,"A1484_unk_70" ,"A13016_1B_70" ,"A15940_1B_100" ,"A9591_1B_100" ,"A576_1B_70" ,"A18104_1D_100" ,"A27014_1B_2B_100","A17967_2D_100" ,"A13978_2D_70" ,"A25742_3A_100" ,"A18332_3A_70" ,"A12777_3D_100" ,"A273_4A_100" ,"A21245_4A_70" ,"A10491_4A_70" ,"A18377_4A_4B_70" ,"A25000_4B_70" ,"A25441_4B_100" ,"A4035_4D_100" ,"A18379_4D_100" ,"A18009_4D_70" ,"A19099_4D_7D_100","A16328_4D_70" ,"A11848_4D_70" ,"A12209_4D_100" ,"A27652_4D_100" ,"A21275_5A_100" ,"A25877_5A_100" ,"A25439_5A_70" ,"A25438_unk_70" ,"A20969_5A_70" ,"A8097_5A_100" ,"A19723_5A_100" ,"A11602_5A_100" ,"A7781_5A_100" ,"A11185_5D_100" ,"A27769_5D_100" ,"A16263_3D_5D_100","A18607_5D_100" ,"A25725_5D_100" ,"A1900_6A_70" ,"A11910_3A_6A_70" ,"A25264_6A_70" ,"A20866_6A_100" ,"A11585_6B_70" ,"A11967_6B_100" ,"A20255_6B_100" ,"A19959_6B_70" ,"A27878_6B_70" ,"A28179_6B_100" ,"A28536_6B_70" ,"A1460_6B_100" ,"A9811_6B_100" ,"A16328_6B_100" ,"A13442_6B_70" ,"A26650_3A_4B_70" ,"A8101_6B_100" ,"A7718_6B_6D_70" ,"A10138_6B_70" ,"A27566_unk_70" ,"A15195_unk_70" ,"A15639_unk_100" ,"A10596_6B_100" ,"A27867_unk_100" ,"A10599_unk_100" ,"A20619_6B_100" ,"A7070_6B_100" ,"A6713_1B_6B_70" ,"A22215_6B_100" ,"A499_5D_100" ,"A21788_6D_70" ,"A21995_6D_100" ,"A18389_6D_100" ,"A11691_6D_100" ,"A11316_6D_70" ,"A28649_6D_70" ,"A28881_6D_100" ,"A18071_6D_100" ,"A8016_7A_7D_100" ,"A6637_unk_70" ,"A7471_7A_100" ,"A12324_4B_7B_70" ,"A13388_7B_100" ,"A28147_7B_100" ,"A27851_7B_70" ,"A11842_7D_70" ,"A28228_7D_100" ,"A21486_7D_100" ,"A27928_7D_70" ,"A18989_unk_100","A7100_7A_70","A23150_unk_70")
LZ8GBS <- drop.markers(LZ8GBS, drop)
gt<-geno.table(LZ8GBS)
gt[gt$P.value<0.05/totmar(LZ8GBS),]

## [1] chr     missing AA      AB      BB      not.BB  not.AA  P.value
## <0 rows> (or 0-length row.names)

Composite Interval Mapping

Across all traits

data_2 <- calc.genoprob(LZ8GBS, step=0)
out.cimE1PGD1 <- cim(data_2, pheno.col="E1D1",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD1$Trait <- ("E1PGD1") 
out.cimE1PGD2 <- cim(data_2, pheno.col="E1D2",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD2$Trait <- ("E1PGD2")
out.cimE1PGD3 <- cim(data_2, pheno.col="E1D3",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD3$Trait <- ("E1PGD3")
out.cimE1PGD4 <- cim(data_2, pheno.col="E1D4",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD4$Trait <- ("E1PGD4")
out.cimE1PGD5 <- cim(data_2, pheno.col="E1D5",map.function=c("kosambi"),method=c("hk"))
out.cimE1PGD5$Trait <- ("E1PGD5")
out.cimE1GIavg <- cim(data_2, pheno.col="E1GI",map.function=c("kosambi"),method=c("hk"))
out.cimE1GIavg$Trait <- ("E1GIavg")
out.cimE2PGD1 <- cim(data_2, pheno.col="E2D1",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD1$Trait <- ("E2PGD1")
out.cimE2PGD2 <- cim(data_2, pheno.col="E2D2",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD2$Trait <- ("E2PGD2")
out.cimE2PGD3 <- cim(data_2, pheno.col="E2D3",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD3$Trait <- ("E2PGD3")
out.cimE2PGD4 <- cim(data_2, pheno.col="E2D4",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD4$Trait <- ("E2PGD4")
out.cimE2PGD5 <- cim(data_2, pheno.col="E2D5",map.function=c("kosambi"),method=c("hk"))
out.cimE2PGD5$Trait <- ("E2PGD5")
out.cimE2GIavg <- cim(data_2, pheno.col="E2GI",map.function=c("kosambi"),method=c("hk"))
out.cimE2GIavg$Trait <- ("E2GIavg")
out.cimE3PGD1 <- cim(data_2, pheno.col="E3D1",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD1$Trait <- ("E3PGD1")
out.cimE3PGD2 <- cim(data_2, pheno.col="E3D2",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD2$Trait <- ("E3PGD2")
out.cimE3PGD3 <- cim(data_2, pheno.col="E3D3",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD3$Trait <- ("E3PGD3")
out.cimE3PGD4 <- cim(data_2, pheno.col="E3D4",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD4$Trait <- ("E3PGD4")
out.cimE3PGD5 <- cim(data_2, pheno.col="E3D5",map.function=c("kosambi"),method=c("hk"))
out.cimE3PGD5$Trait <- ("E3PGD5")
out.cimE3GIavg <- cim(data_2, pheno.col="E3GI",map.function=c("kosambi"),method=c("hk"))
out.cimE3GIavg$Trait <- ("E3GIavg")

QTL.LZ8.ERA8 <- rbind(out.cimE1PGD1, out.cimE1PGD2, out.cimE1PGD3, out.cimE1PGD4, out.cimE1PGD5, out.cimE1GIavg, out.cimE2PGD1, out.cimE2PGD2, out.cimE2PGD3, out.cimE2PGD4, out.cimE2PGD5, out.cimE2GIavg,out.cimE3PGD1, out.cimE3PGD2, out.cimE3PGD3, out.cimE3PGD4, out.cimE3PGD5, out.cimE3GIavg)
# write.table(QTL.LZ8.ERA8, "./QTLLZ8All.txt", sep="\t")

a<-scanone(LZ8GBS, method="hk", pheno.col = "E1D1", n.perm=1000)
c<-scanone(LZ8GBS, method="hk", pheno.col = "E1D2", n.perm=1000)
d<-scanone(LZ8GBS, method="hk", pheno.col = "E1D3", n.perm=1000)
e<-scanone(LZ8GBS, method="hk", pheno.col = "E1D4", n.perm=1000)
f<-scanone(LZ8GBS, method="hk", pheno.col = "E1D5", n.perm=1000)
g<-scanone(LZ8GBS, method="hk", pheno.col = "E1GI", n.perm=1000)
b1<-scanone(LZ8GBS, method="hk", pheno.col = "E2D1", n.perm=1000)
c1<-scanone(LZ8GBS, method="hk", pheno.col = "E2D2", n.perm=1000)
d1<-scanone(LZ8GBS, method="hk", pheno.col = "E2D3", n.perm=1000)
e1<-scanone(LZ8GBS, method="hk", pheno.col = "E2D4", n.perm=1000)
f1<-scanone(LZ8GBS, method="hk", pheno.col = "E2D5", n.perm=1000)
g1<-scanone(LZ8GBS, method="hk", pheno.col = "E2GI", n.perm=1000)
h<-scanone(LZ8GBS, method="hk", pheno.col = "E3D1", n.perm=1000)
i<-scanone(LZ8GBS, method="hk", pheno.col = "E3D2", n.perm=1000)
j<-scanone(LZ8GBS, method="hk", pheno.col = "E3D3", n.perm=1000)
k<-scanone(LZ8GBS, method="hk", pheno.col = "E3D4", n.perm=1000)
l<-scanone(LZ8GBS, method="hk", pheno.col = "E3D5", n.perm=1000)
m<-scanone(LZ8GBS, method="hk", pheno.col = "E3GI", n.perm=1000)
summary(a) #
summary(c) #
summary(d) #
summary(e) #
summary(f) #
summary(g) #
summary(b1) #
summary(c1) #
summary(d1) #
summary(e1) #
summary(f1) #
summary(g1) #
summary(h) #
summary(i) #
summary(j) #
summary(k) #
summary(l) #
summary(m) #

Heading and Height

QTL Analysis

data_2 <- calc.genoprob(LZ8GBS, step=0)
out.cimHeight <- cim(data_2, pheno.col="Height",map.function=c("kosambi"),method=c("hk"))
out.cimHeight$Trait <- ("Height") 
out.cimHeading <- cim(data_2, pheno.col="Heading",map.function=c("kosambi"),method=c("hk"))
out.cimHeading$Trait <- ("Heading")

Calculate Significance Threshold

a<-scanone(LZ8GBS, method="hk", pheno.col = "Heading", n.perm=1000)
summary(a) #4.43 5%
b<-scanone(LZ8GBS, method="hk", pheno.col = "Height", n.perm=1000)
summary(b) #4.38 5%

Graph the Manhatten Plot

For graphing purposes only, add a column to give unique positions for each marker per chromosome.

# uniqpos <- read.table("./GBS_UniquePos.txt",header=TRUE,na.string=c("", " ", "NA", "na"))
uniqpos$chr <- NULL
uniqpos$pos <- NULL
QTL.LZ8.ERA8 <- as.data.frame(QTL.LZ8.ERA8)
LZ8  <- cbind(QTL.LZ8.ERA8, uniqpos)

ggplot(LZ8, aes(x=LZ8$order, y=LZ8$lod))+
  geom_point(aes(colour = factor(Genome)),size = 2)+ylim(0,6)+theme_bw()+ theme(axis.text = font,  axis.title = font, axis.text.x=element_blank())+ scale_color_manual(values=c("#0A0A0A", "#8A8A8A", "#D1D1D1","#A2CD5A"))+ 
  labs(x = "Chromosome 1A - 7D and Unknown" , y = "LOD",subtitle="ABA Sensitivity QTL")

# uniqposAHH <- read.table("./GBS_UniquePos_ABAHeadHei.txt",header=TRUE,na.string=c("", " ", "NA", "na"))
uniqposAHH$chr <- NULL
uniqposAHH$pos <- NULL
QTL.LZ8.ERA8 <- rbind(out.cimE1PGD1, out.cimE1PGD2, out.cimE1PGD3, out.cimE1PGD4, 
                      out.cimE1PGD5, out.cimE1GIavg, out.cimE2PGD1, out.cimE2PGD2, 
                      out.cimE2PGD3, out.cimE2PGD4, out.cimE2PGD5, out.cimE2GIavg,
                      out.cimE3PGD1, out.cimE3PGD2, out.cimE3PGD3, out.cimE3PGD4, 
                      out.cimE3PGD5, out.cimE3GIavg, out.cimHeight, out.cimHeading )
QTL.LZ8.ERA8 <- as.data.frame(QTL.LZ8.ERA8)
QTL.LZ8  <- cbind(QTL.LZ8.ERA8, uniqposAHH)
ggplot(QTL.LZ8, aes(x=QTL.LZ8$order, y=QTL.LZ8$lod))+
  geom_point(aes(colour = factor(Genome)),size = 2)+ylim(0,6)+theme_bw()+ 
  theme(axis.text = font,  axis.title = font, axis.text.x=element_blank())+ 
  scale_color_manual(values=c("#0A0A0A", "#8A8A8A", "#D1D1D1","#A2CD5A")) + geom_hline(aes(yintercept = 4.4))+ 
  labs(x = "Chromosome 1A - 7D and Unknown" , y = "LOD",  title = "Figure S5b", subtitle="ABA, Heading Date, and Height QTL")

QTL.LZ8.4A <-  subset(QTL.LZ8 ,  QTL.LZ8$chr=="4A.1" | QTL.LZ8$chr=="4A.2" | QTL.LZ8$chr=="4A.3" )
ggplot(QTL.LZ8.4A, aes(x=QTL.LZ8.4A$order, y=QTL.LZ8.4A$lod))+
  geom_point(aes(colour = factor(chr)),size = 2)+ylim(0,6)+theme_bw()+ 
  theme(axis.text = font,  axis.title = font, axis.text.x=element_blank())+ 
  scale_color_manual(values=c("#0A0A0A", "#8A8A8A", "#D1D1D1")) + geom_hline(aes(yintercept = 4.4))+ 
  labs(x = "Chromosome 4A" , y = "LOD",  title = "Figure S5c", subtitle="ABA, Heading Date, and Height QTL")

List of Significant QTL on Chromosome 4A

QTL.LZ8.4A %>% 
  mutate(lod = round(lod,2)) %>%
  filter(lod > 4.4) %>%
  arrange(str_extract(chr, '^.'), pos) %>%
  dplyr::select(Marker, chr, pos,  Trait, lod)  %>%
  datatable()#rownames = FALSE, container = sketch, caption = htmltools::tags$caption(style = 'text-align: left;', 'Table S11 - Significant QTL on chromosome 4A'))

Significant Alleles in LZ8 population

ATol<-read.table("AllTolQTLAlleles.txt", head = T,na.string=c("."," ","NA","na","NaN"))

ggplot(ATol, aes(x = OverallTolAlleles, y = E1D2)) +
  geom_point(aes(colour = ATol$Color), size = 2)+ theme_bw() + 
  theme(panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "#7D7D7D"), panel.background = element_blank(),
        axis.line = element_line(size=0.5))   +
  geom_smooth(method=lm,colour = "#7A7A7A")+ ylim(0, 100) +
  scale_colour_manual(values = c("#008B8B", "#0A0A0A" ,"#7A7A7A","#551A8B")) + 
  # theme(legend.position="none")+
  labs(y = "Percent Germination", title = "E1D2", subtitle = "Figure S4a", x = "")

ggsave("LZ8_4A_HaplotypeE1.png", width = 3, height = 3, units = "in")

ggplot(ATol, aes(x = OverallTolAlleles, y = E2D3)) +
  geom_point(aes(colour = ATol$Color), size = 2)+ theme_bw() + 
  theme(panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "#7D7D7D"), panel.background = element_blank(),
        axis.line = element_line(size=0.5))   +
  geom_smooth(method=lm,colour = "#7A7A7A")+ ylim(0, 100) +
  scale_colour_manual(values = c("#008B8B", "#0A0A0A" ,"#7A7A7A","#551A8B")) + 
  # theme(legend.position="none")+
  labs(y = "Percent Germination", title = "E2D3", subtitle = "Figure S4b", x = "Number of Favorable Alleles")

ggsave("LZ8_4A_HaplotypeE2.png", width = 3, height = 3, units = "in")

ggplot(ATol, aes(x = OverallTolAlleles, y = E3D1)) +
  geom_point(aes(colour = ATol$Color), size = 2)+ theme_bw() + 
  theme(panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "#7D7D7D"), panel.background = element_blank(),
        axis.line = element_line(size=0.5))   +
  geom_smooth(method=lm,colour = "#7A7A7A")+ ylim(0, 100) +
  scale_colour_manual(values = c("#008B8B", "#0A0A0A" ,"#7A7A7A","#551A8B"))+ 
  # theme(legend.position="none")+
  labs(y = "Percent Germination", title = "E3D1", subtitle = "Figure S4c", x = "")

ggsave("LZ8_4A_HaplotypeE3.png", width = 3, height = 3, units = "in")

Otis/ZakERA8 F2:3

The OZ8 dataset contains the Otis/ZakERA8 F2: F3 genotypic data, phenotypic data, and marker positions on chromosome 4A.

OZ8<-read.cross("csv",".","./OZ8_SNP6-SNP29.csv",na.strings=c("-","NA","."),BC.gen=0, F.gen=2)

OZ8

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     BC(0)F(2) cross
## 
##     No. individuals:    108 
## 
##     No. phenotypes:     7 
##     Percent phenotyped: 100 98.1 98.1 98.1 98.1 98.1 98.1 
## 
##     No. chromosomes:    1 
##         Autosomes:      4A 
## 
##     Total markers:      4 
##     No. markers:        4 
##     Percent genotyped:  99.5 
##     Genotypes (%):      AA:24.0  AB:49.1  BB:27.0  not BB:0.0  not AA:0.0

# geno.image(OZ8)
# summary(OZ8$pheno)

Check for skewed markers

gt<-geno.table(OZ8)
gt[gt$P.value<0.05/totmar(OZ8),]  

## [1] chr     missing AA      AB      BB      not.BB  not.AA  P.value
## <0 rows> (or 0-length row.names)

No markers are skewed.

The phenotypes include percent germination (PG) is separated by day (d1…d5) and Germination Index (GI) is calculated over all 5 days:
> PG_d1, PG_d2, PG_d3, PG_d4, PG_d5, GI

Compite Interval Mapping

data_2 <- calc.genoprob(OZ8, step=0)
out.cimGI <- cim(data_2, pheno.col="GI",map.function=c("kosambi"),method=c("hk"))
a<-scanone(OZ8, method="hk", pheno.col = "GI", n.perm=1000)
summary(a) 

out.cimd5 <- cim(data_2, pheno.col="PG_d5",map.function=c("kosambi"),method=c("hk"))
out.cimd4 <- cim(data_2, pheno.col="PG_d4",map.function=c("kosambi"),method=c("hk"))
out.cimd3 <- cim(data_2, pheno.col="PG_d3",map.function=c("kosambi"),method=c("hk"))
out.cimd2 <- cim(data_2, pheno.col="PG_d2",map.function=c("kosambi"),method=c("hk"))
out.cimd1 <- cim(data_2, pheno.col="PG_d1",map.function=c("kosambi"),method=c("hk"))

Wrangle Results

out.cimd5$Pheno <- "d5"
out.cimd4$Pheno <- "d4"
out.cimd3$Pheno <- "d3"
out.cimd2$Pheno <- "d2"
out.cimd1$Pheno <- "d1"
out.cimGI$Pheno <- "GI"
Otera8 <- rbind(out.cimGI ,out.cimd5,out.cimd4,out.cimd3,out.cimd2,out.cimd1 )

Graph Results

ggplot(Otera8, aes(x=Otera8$pos, y=Otera8$lod, colour = factor(Otera8$Pheno)))+
  geom_line(size = 1.5)+
  theme_bw()+ theme(axis.text = font,  axis.title = font)+
  geom_hline(aes(yintercept = 1.73, colour = ""))+ 
  labs(x = "Chromosome 4A (cM)", y = "LOD") 

out.cimGI2 <- rbind(out.cimGI, out.cimd5)
ggplot(out.cimGI2, aes(x=out.cimGI2$pos, y=out.cimGI2$lod))+
  geom_line(aes(colour = factor(Pheno)),size = 1)+
  theme_bw()+ theme(axis.text = font,  axis.title = font)+
  geom_hline(aes(yintercept = 1.78))+ 
  labs(x = "Chromosome 4A (cM)", y = "LOD")  +
  ylim (0,10) +xlim (0,11)+
  scale_color_manual(values=c("#CD6600","#FFB90F")) + theme(legend.position="none")+
  ggsave("OZ8_4A_QTL.png", width = 3, height = 3, units = "in")

Interval Map of 4A Region

imap <- read.csv("D:/Shantel/Research Projects/3. Zak ERA8/1_Mapping Zak ERA8 Data/ERA8 Interval Mapping/ZZ8 5.1 5.2 Group Data for R.csv", na.string="-")

imap$Group <-  as.character(imap$Group)
imapf <-  filter(imap, Group == "ERA8-g" | Group == "I") 
t.test(PG_avg ~ Group, data=imapf)

## 
##  Welch Two Sample t-test
## 
## data:  PG_avg by Group
## t = -4.5576, df = 41.888, p-value = 4.435e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -41.06224 -15.85678
## sample estimates:
## mean in group ERA8-g      mean in group I 
##             24.73684             53.19635

Repeat for other groups

imapf <-  filter(imap, Group == "Zak" | Group == "I") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "Zak-g" | Group == "I") 
t.test(PG_avg ~ Group, data=imapf)

imapf <-  filter(imap, Group == "Zak-g" | Group == "Zak") 
t.test(PG_avg ~ Group, data=imapf)

imapf <-  filter(imap, Group == "ERA8-g" | Group == "II") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "Zak" | Group == "II") 
t.test(PG_avg ~ Group, data=imapf)



imapf <-  filter(imap, Group == "Zak-g" | Group == "III") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "ERA8-g" | Group == "III") 
t.test(PG_avg ~ Group, data=imapf)

imapf <-  filter(imap, Group == "Zak-g" | Group == "IV") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "ERA8-g" | Group == "IV") 
t.test(PG_avg ~ Group, data=imapf)

imapf <-  filter(imap, Group == "ERA8-g" | Group == "VI") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "Zak-g" | Group == "VI") 
t.test(PG_avg ~ Group, data=imapf)
imapf <-  filter(imap, Group == "ERA8" | Group == "VI") 
t.test(PG_avg ~ Group, data=imapf)

Plot Interval Map

imap$PG_avg <- as.integer(imap$PG_avg)
target <- c( "VII","VI","V","IV", "III", "II", "I", "Zak-g", "ERA8-g", "Zak", "ERA8")
require(gdata)
imap$Group<- reorder.factor( imap$Group, new.order=target)
library(ggplot2)
ggplot(imap, aes(imap$Group, imap$PG_avg)) +
  geom_boxplot(aes(fill = factor(color)))+ theme_bw() +
  scale_fill_manual(values=c( "#7A7A7A","#A2CD5A","#440154" ))+
  geom_jitter(width = 0.05)+
  scale_y_continuous(breaks = seq(0, 100,10))+
  coord_flip() + labs(title = "Figure 4e")

ggsave("20190430_ZZ8_IntervalMap.PNG", width = 10, height = 8, units = "in")

Session Information

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gdata_2.18.0         DT_0.5               forcats_0.4.0       
##  [4] stringr_1.4.0        dplyr_0.8.0.1        purrr_0.3.2         
##  [7] readr_1.3.1          tidyr_0.8.3          tibble_2.1.1        
## [10] ggplot2_3.1.1        tidyverse_1.2.1      qtl_1.44-9          
## [13] base64enc_0.1-3      RevoUtils_11.0.3     RevoUtilsMath_11.0.0
## 
## loaded via a namespace (and not attached):
##  [1] gtools_3.8.1     tidyselect_0.2.5 xfun_0.6         haven_2.1.0     
##  [5] lattice_0.20-38  colorspace_1.4-1 generics_0.0.2   htmltools_0.3.6 
##  [9] yaml_2.2.0       rlang_0.3.4      later_0.8.0      pillar_1.3.1    
## [13] glue_1.3.1       withr_2.1.2      modelr_0.1.4     readxl_1.3.1    
## [17] plyr_1.8.4       munsell_0.5.0    gtable_0.3.0     cellranger_1.1.0
## [21] rvest_0.3.3      htmlwidgets_1.3  evaluate_0.13    labeling_0.3    
## [25] knitr_1.22       httpuv_1.5.1     crosstalk_1.0.0  parallel_3.5.3  
## [29] broom_0.5.2      Rcpp_1.0.1       xtable_1.8-3     promises_1.0.1  
## [33] scales_1.0.0     backports_1.1.4  jsonlite_1.5     mime_0.6        
## [37] hms_0.4.2        digest_0.6.18    stringi_1.4.3    shiny_1.3.1     
## [41] grid_3.5.3       cli_1.1.0        tools_3.5.3      magrittr_1.5    
## [45] lazyeval_0.2.2   crayon_1.3.4     pkgconfig_2.0.2  xml2_1.2.0      
## [49] lubridate_1.7.4  assertthat_0.2.1 rmarkdown_1.12   httr_1.4.0      
## [53] rstudioapi_0.10  R6_2.3.0         nlme_3.1-137     compiler_3.5.3