Maize landraces genetic map from Seeds data
Luis Avila
1/30/2017
Loading libraries
library(tidyverse)
library(cowplot)
library(qtl)I have recombination events for the individuals from Maize populations
From Peter Bradbury’s imputation and phasing
xo<-read.table(file="xo_Combined_LR13_LR14_parents2_AGPv4_filter1222.txt",
header=TRUE,stringsAsFactors=FALSE)
head(xo)## taxon chr start end
## 1 13FL990708:250295003 1 262738 1023028
## 2 13FL991001:250307264 1 560842 1211062
## 3 13FL991580:250306177 1 560842 1211062
## 4 13FL992495:250308398 1 560842 1211062
## 5 13FL992927:250426784 1 560842 1211062
## 6 14FL991065:250410291 1 560842 1211062dim(xo)## [1] 193444 4num.individuals<-length(unique(xo$taxon))
num.individuals## [1] 5130working with all 10 chromosomes
window.size<-2e7
chromo.keep<-c()
position.keep<-c()
cM.keep<-c()
for (current.chr in 1:10) {
max.end<-max(filter(xo,chr==current.chr)$end)
my.breaks<-seq(from=0,to=max.end,by=window.size)
if(max(my.breaks)<max.end){ my.breaks<-c(my.breaks,max.end)}
current.chr.table<-hist(xo[xo$chr==current.chr,]$end,breaks = my.breaks,plot=FALSE)
#prob.recombination<-current.chr.table$counts/sum(current.chr.table$counts)
prob.recombination<-
current.chr.table$counts/(2*num.individuals)
#pos<- 50 * log(1/(1-2*prob.recombination))
pos<-prob.recombination*100
for (i in 2:length(pos)) {
pos[i]=pos[i]+pos[i-1]
}
pos<-c(0,pos) #adding an start of the chromosome
chromo.keep<-c(chromo.keep,rep(current.chr,length(pos)))
cM.keep<-c(cM.keep,pos)
position.keep<-c(position.keep,my.breaks)
}
roughMaizeLandRaceData<-data.frame(chr=chromo.keep,pos=cM.keep,phys_pos=position.keep)
rownames(roughMaizeLandRaceData) <- paste0(as.character(roughMaizeLandRaceData$chr),"_",round(roughMaizeLandRaceData$phys_pos/1e6,digits=0),"Mb")
head(roughMaizeLandRaceData)## chr pos phys_pos
## 1_0Mb 1 0.00000 0e+00
## 1_20Mb 1 41.96881 2e+07
## 1_40Mb 1 64.64912 4e+07
## 1_60Mb 1 85.11696 6e+07
## 1_80Mb 1 96.56920 8e+07
## 1_100Mb 1 106.09162 1e+08 tail(roughMaizeLandRaceData)## chr pos phys_pos
## 10_60Mb 10 54.98051 60000000
## 10_80Mb 10 61.35478 80000000
## 10_100Mb 10 69.19103 100000000
## 10_120Mb 10 75.90643 120000000
## 10_140Mb 10 100.55556 140000000
## 10_151Mb 10 143.41131 150827693And plotting the map for the ten maize chromosomes
maizeLandracesMap <- table2map(roughMaizeLandRaceData)
plotMap(maizeLandracesMap,show.marker.names=FALSE,main="Maize landraces all")
## Loading maize genetic map from Ogut et al. 2015
OgutMap<-read.table(file="ogut_fifthcM_map_agpv4.txt",
header=FALSE,stringsAsFactors=FALSE)
names(OgutMap)<-c("marker.name.1","marker.name.2","cM","chromosome","position")
head(OgutMap)## marker.name.1 marker.name.2 cM chromosome position
## 1 S_88897 M1 -4.8 1 134914
## 2 S_170386 M2 -4.6 1 223935
## 3 S_251876 M3 -4.4 1 350753
## 4 S_333365 M4 -4.2 1 380021
## 5 S_496344 M6 -3.8 1 541304
## 6 S_577834 M7 -3.6 1 612144dim(OgutMap)## [1] 6740 5dim(OgutMap[OgutMap$chromosome==1,])## [1] 985 5finding markers similar to our simulated markers
window.size<-2e7
chromo.keep<-c()
position.keep<-c()
cM.keep<-c()
for(chr in 1:10) {
current.chromosome<-OgutMap[OgutMap$chromosome==chr,]
current.chromosome.min.cM<-min(current.chromosome$cM)
current.max=max(current.chromosome$position)
max.limit<-(abs(current.max/window.size)+1)*window.size
current.breaks=seq(from=0,to=max.limit,by=window.size)
for (j in 1:length(current.breaks)) {
best.index<-which.min(abs(current.chromosome$position - current.breaks[j]))
#print(best.index)
chromo.keep<-c(chromo.keep,chr)
position.keep<-c(position.keep,current.chromosome[best.index,"position"])
cM.keep<-c( cM.keep,current.chromosome[best.index,"cM"]+abs(current.chromosome.min.cM))
}
}
roughOgutData<-data.frame(chr=chromo.keep,pos=cM.keep,phys_pos=position.keep)
head(roughOgutData)## chr pos phys_pos
## 1 1 0.0 134914
## 2 1 40.6 19912274
## 3 1 62.2 39965698
## 4 1 76.2 60072613
## 5 1 84.8 79738562
## 6 1 91.8 99413055 tail(roughOgutData)## chr pos phys_pos
## 116 10 42.8 63639269
## 117 10 45.6 80508347
## 118 10 49.6 100056068
## 119 10 53.4 120024729
## 120 10 72.4 139943108
## 121 10 112.8 150873537 rounded.phys.pos<-round(roughOgutData$phys_pos/1e6,digits=0)
rownames(roughOgutData) <- paste0(as.character(chromo.keep),"_", rounded.phys.pos,"Mb")
rough.ogut.map <- table2map(roughOgutData)Plotting a maize genetic map for Chr1 from Ogut et al. data
plotMap(rough.ogut.map,chr=c(1),show.marker.names=TRUE,main="Chromosome 1 maize")
One besides the other
Comparing Maps with R/Qtl function
left map is Maize (from Ogut et. al), right map is from the maize landraces
plotMap(rough.ogut.map,maizeLandracesMap,chr=c(1),show.marker.names=TRUE,shift=TRUE, main="Ogut's maize (left) and Seeds landraces (right)")
roughMaizeLandRaceData[roughMaizeLandRaceData$chr==1,]## chr pos phys_pos
## 1_0Mb 1 0.00000 0
## 1_20Mb 1 41.96881 20000000
## 1_40Mb 1 64.64912 40000000
## 1_60Mb 1 85.11696 60000000
## 1_80Mb 1 96.56920 80000000
## 1_100Mb 1 106.09162 100000000
## 1_120Mb 1 110.45809 120000000
## 1_140Mb 1 112.04678 140000000
## 1_160Mb 1 115.96491 160000000
## 1_180Mb 1 124.38596 180000000
## 1_200Mb 1 140.66277 200000000
## 1_220Mb 1 159.91228 220000000
## 1_240Mb 1 175.17544 240000000
## 1_260Mb 1 190.94542 260000000
## 1_280Mb 1 210.21442 280000000
## 1_300Mb 1 244.24951 300000000
## 1_307Mb 1 257.29045 306776194 roughOgutData[roughOgutData$chr==1,]## chr pos phys_pos
## 1_0Mb 1 0.0 134914
## 1_20Mb 1 40.6 19912274
## 1_40Mb 1 62.2 39965698
## 1_60Mb 1 76.2 60072613
## 1_80Mb 1 84.8 79738562
## 1_99Mb 1 91.8 99413055
## 1_120Mb 1 93.8 119583370
## 1_132Mb 1 94.0 131929986
## 1_161Mb 1 95.8 160996844
## 1_180Mb 1 100.4 180044375
## 1_200Mb 1 115.0 199956280
## 1_220Mb 1 131.8 220202926
## 1_240Mb 1 144.0 239783178
## 1_260Mb 1 155.6 259906710
## 1_280Mb 1 171.2 280012317
## 1_300Mb 1 198.2 300092579
## 1_307Mb 1 210.2 306879533Individual comparisons
for(my.chr in 1:10) {
par(mfrow=c(1,2))
plotMap(maizeLandracesMap,chr=c(my.chr), show.marker.names=TRUE, main=paste0("Chromosome ",my.chr ," Seeds maize landraces"))
plotMap(rough.ogut.map,chr=c(my.chr),show.marker.names=TRUE,main=paste0("Chromosome ",my.chr ," maize"))
}
Total distance in cM
max.per.chromosome.teo<-tapply(roughMaizeLandRaceData$pos,roughMaizeLandRaceData$chr,max)
max.per.chromosome.teo## 1 2 3 4 5 6 7 8
## 257.2904 220.3119 199.5517 187.3782 200.6530 168.2261 187.3099 166.5497
## 9 10
## 154.7368 143.4113 total.cM.SeedsLandraces<-sum(max.per.chromosome.teo)
total.cM.SeedsLandraces## [1] 1885.419 ## former negative cM in Ogut's data have been adjusted to 0
## hence the final cM is larger here that in Ogut's
max.per.chromosome.maize<-tapply(roughOgutData$pos,roughOgutData$chr,max)
max.per.chromosome.maize## 1 2 3 4 5 6 7 8 9 10
## 210.2 160.0 161.4 151.8 157.0 111.8 138.4 137.4 130.2 112.8 total.cM.Maize<-sum(max.per.chromosome.maize)
total.cM.Maize## [1] 1471