Commit 20fec586 authored by Sterken, Mark's avatar Sterken, Mark
Browse files

Commit BioRxiv version 4temp paper by M.I. Maualana et al.

parents
.Rproj.user
.Rhistory
.RData
.Ruserdata
###function for broad-sense heritability
#currently 3 methods are supported:
# ANOVA (very basic method) [ref?]
# REML (based on lme4 package) [ref?]
# Keurentjes_2007 (based on Brem, 2005 and Keurentjes, 2007)
# all methods support permutation
###input
#trait.value: a vector with trait values
#strain: a vector with strain names matching the trait values
#trait: a vector with trait names (e.g you can enter multiple traits this way)
#method: one of three methods: ANOVA, REML, or Keurentjes_2007 (standard is ANOVA because it's quick)
#parental.strain: names of parental strains (only for Keurentjes; otherwise parental strains should be removed from the analysis)
#n.perm: number of permutations to be done (standard 1000)
#Vg.factor: fraction of the measured genotypic variance should be counted towards total variance (standard is 1)
###output
#a dataframe with traits per row and columns with the heritability, underlying values and an FDR0.05 threshold
H2.broad <- function(trait.value,strain,trait,method,parental.strain,n.perm,Vg.factor){
if(missing(trait.value)){ stop("requires trait value")}
if(missing(strain)){ stop("requires strain values")}
if(missing(trait)){ trait <- rep("A",length(trait.value))}
if(missing(method)){ method <- "ANOVA"}
if(!method %in% c("ANOVA","REML","Keurentjes_2007")){ stop("method should be either ANOVA, REML, or Keurentjes_2007")}
if(method=="Keurentjes_2007"){
if(missing(parental.strain)){ stop("when using Keurentjes_2007, parental strain names should be provided")}
if(sum(parental.strain %in% strain) != length(parental.strain)){
stop("not all parental strains are in the data")}
}
if(missing(n.perm)){ n.perm <- 1000}
if(missing(Vg.factor)){ Vg.factor <- 1}
###calculates H2 based on REML (part of H2.broad)
H2.REML <- function(input,Vg.factor){
if(missing(Vg.factor)){Vg.factor <- 1}
##test if lme4 package is installed
if(!require("lme4")){
install.packages("lme4")
library("lme4")
}
REML <- lmer(trait.value ~ 1 + (1|strain),data=input)
Variations <- as.data.frame(VarCorr(REML))$vcov
H2 <- c((Variations[1]*Vg.factor)/sum(Variations*c(Vg.factor,1)),(Variations*c(Vg.factor,1)))
names(H2) <- c("H2_REML","Vg","Ve")
return(H2)
}
###Calculates H2 based on ANOVA (part of H2.broad)
H2.ANOVA <- function(input,Vg.factor){
if(missing(Vg.factor)){Vg.factor <- 1}
ANOVA <- anova(lm(trait.value~strain,data=input))
Variations <- ANOVA[[2]]
H2 <- c((Variations[1]*Vg.factor)/sum(Variations*c(Vg.factor,1)),(Variations*c(Vg.factor,1)))
names(H2) <- c("H2_ANOVA","Vg","Ve")
return(H2)
}
###Calculates H2 as described in Keurentjes, 2007, in case error has to be estimated on limited number of strains
###Permutation is based on Brem, 2005
H2.Keurentjes <- function(input,parental.strain,Vg.factor){
if(missing(Vg.factor)){Vg.factor <- 1}
selc.parental <- input$strain %in% parental.strain
###Ve is estimated based on parental replicates
VAR.PL <- tapply(input$trait.value[selc.parental],as.character(unlist(input$strain[selc.parental])),var)
N.PL <- tapply(input$trait.value[selc.parental],as.character(unlist(input$strain[selc.parental])),length)
VAR.E <- sum(VAR.PL*(N.PL-1))/sum((N.PL-1))
##Vt is estimated based on RILs
VAR.T <- var(input$trait.value[!selc.parental])
H2 <- c(((VAR.T-VAR.E)*Vg.factor)/VAR.T,((VAR.T-VAR.E)*Vg.factor),VAR.E)
names(H2) <- c("H2_keurentjes","Vg","Ve")
return(H2)
}
###Prep the data
input <- data.frame(cbind(trait.value=trait.value,strain=strain))
input$trait.value <- as.numeric(as.character(unlist(input$trait.value)))
input$strain <- as.factor(as.character(unlist(input$strain)))
input <- split(input,trait)
if(method == "ANOVA"){
H2 <- t(sapply(input,H2.ANOVA,Vg.factor))
perm <- NULL
for(i in 1:n.perm){
input.perm <- lapply(input,function(x){x[,1] <- x[order(runif(nrow(x))),1];return(x)})
H2.perm <- t(sapply(input.perm,H2.ANOVA,Vg.factor))
perm <- rbind(perm,H2.perm[,1])
}
H2 <- cbind(H2,FDR0.05_ANOVA=apply(perm,2,function(x){sort(x)[ceiling(n.perm*0.95)]}))
}
if(method == "REML"){
H2 <- t(sapply(input,H2.REML,Vg.factor))
perm <- NULL
for(i in 1:n.perm){
input.perm <- lapply(input,function(x){x[,1] <- x[order(runif(nrow(x))),1];return(x)})
H2.perm <- t(sapply(input.perm,H2.REML,Vg.factor))
perm <- rbind(perm,H2.perm[,1])
}
H2 <- cbind(H2,FDR0.05_REML=apply(perm,2,function(x){sort(x)[ceiling(n.perm*0.95)]}))
}
if(method == "Keurentjes_2007"){
H2 <- t(sapply(input,H2.Keurentjes,parental.strain = parental.strain,Vg.factor))
perm <- NULL
for(i in 1:n.perm){
input.perm <- lapply(input,function(x){x[,1] <- x[order(runif(nrow(x))),1];return(x)})
H2.perm <- t(sapply(input.perm,H2.Keurentjes,parental.strain = parental.strain,Vg.factor))
perm <- rbind(perm,H2.perm[,1])
}
H2 <- cbind(H2,FDR0.05_Keurentjes=apply(perm,2,function(x){sort(x)[ceiling(n.perm*0.95)]}))
}
output <- as.data.frame(cbind(trait=names(input),round(H2,digits=3)))
output[,1] <- as.character(unlist(output[,1]))
output[,2] <- as.numeric(as.character(unlist(output[,2])))
output[,3] <- as.numeric(as.character(unlist(output[,3])))
output[,4] <- as.numeric(as.character(unlist(output[,4])))
output[,5] <- as.numeric(as.character(unlist(output[,5])))
return(output)
}
###Data preparation for QTL mapping; start function for most of subsequent functions
###input
#trait.matrix (matrix; traits per row and strains per column)
#strain.trait (vector with strain names per trait. Should match strain.map)
#strain.map (markers in rows, genotypes in columns)
#strain.marker (markers in rows, with columns name, chromosome and position)
###output
#a list witht the entries Trait, Map, and Marker. Where the Traits are aligned with the map
###Description
#Arranges the trait and the map data so they are aligned and runs some basic checks on the input data
###See also
#QTL.map.1
###Data preparation
QTL.data.prep <- function(trait.matrix,strain.trait,strain.map,strain.marker){
if(missing(trait.matrix)){ stop("requires trait data, traits per row and strains per column")}
if(missing(strain.trait)){ stop("requires the strain names matching the trait data")}
if(missing(strain.map)){ stop("requires the genetic map of the population")}
if(missing(strain.marker)){ stop("requires the marker info (name, chromosome, basepair) of the map")}
if(ncol(trait.matrix) != length(strain.trait)){ stop("the number of strains does not correspond with the trait matrix")}
if(sum(tolower(strain.trait) %in% tolower(colnames(strain.map))) < length(strain.trait)){
warning(paste("The strains: ",paste(strain.trait[tolower(strain.trait) %in% tolower(colnames(strain.map))],collapse=", ")," are not in the strain.map file"))
}
if(nrow(strain.map) != nrow(strain.marker)){ stop("the number of markers (strain.map) does not equal the number of marker descriptions (strain.marker)")}
if(ncol(strain.marker) != 3){ stop("make sure the strain.marker file contains a column with the marker name, chromosome, and bp location")}
strain.map <- data.matrix(strain.map)
###Make matrix
map.nu <- matrix(NA,nrow=nrow(strain.map),ncol=length(strain.trait))
for(j in 1:length(strain.trait)){
map.nu[,j] <- strain.map[,tolower(colnames(strain.map)) == tolower(strain.trait[j])]
}
colnames(map.nu) <- strain.trait
rownames(map.nu) <- rownames(strain.map)
###add rownames to the trait matrix
if(length(rownames(trait.matrix)) == 0){
rownames(trait.matrix) <- 1:nrow(trait.matrix)
}
###Make output
output <- NULL
output[[1]] <- trait.matrix
output[[2]] <- map.nu
output[[3]] <- data.frame(strain.marker)
names(output) <- c("Trait","Map","Marker")
###return output
return(output)
}
\ No newline at end of file
###function for single marker mapping
#Snoek & Sterken, 2017; Sterken, 2017
###input
#trait.matrix (traits in rows, genotypes in columns)
#strain.map (markers in rows, genotypes in columns)
#strain.marker (markers in rows, with columns name, chromosome and position)
###output
#a list with names: LOD, Effect, Trait, Map, and Marker.
###Description
# map.1 is optimized for dataset in which many markers are un-informative
# these markers are not mapped but the p-vals and effect are copied from the previous marker.
QTL.map.1 <- function(trait.matrix,strain.map,strain.marker){
if(missing(trait.matrix)){ stop("specify trait.matrix, matrix with traits in rows and genotypes in columns")}
if(missing(strain.map)){ stop("specify strain.map (genetic map)")}
if(ncol(strain.map) != dim(as.matrix(trait.matrix))[2] &
ncol(strain.map) != prod(dim(as.matrix(trait.matrix)))){ stop("number of strains is not equal to number of genotypes in map")}
if(missing(strain.marker)){ strain.marker <- data.frame(cbind(name=1:dim(strain.map)[1],chromosome=NA,bp=1:dim(strain.map)[1]))
rownames(strain.marker) <- 1:dim(strain.map)[1]}
###Check for NAs
NAs <- sum(is.na(trait.matrix))>0 | sum(is.na(strain.map))>0
small <- FALSE
if(dim(as.matrix(trait.matrix))[1] ==1){
traitnames <- rownames(trait.matrix)
trait.matrix <- rbind(trait.matrix,trait.matrix)
rownames(trait.matrix) <- c("A","B")
small <- TRUE
}
eff.out <- matrix(NA,nrow(trait.matrix),nrow(strain.map))
pval.out <- matrix(NA,nrow(trait.matrix),nrow(strain.map))
for(i in 1:nrow(strain.map)){
noseg <- length(unique(strain.map[i,])) ==1
if(noseg){
output.tmp <- matrix(c(0,0),byrow=TRUE,ncol=2,nrow=nrow(trait.matrix))
}
if( i == 1 & !noseg){
if(!NAs){output.tmp <- fast.lm.1(trait.matrix,strain.map[i,])}
if(NAs){output.tmp <- lm.1(trait.matrix,strain.map[i,])}
}
if( i != 1 & sum(abs(strain.map[i-1,]-strain.map[i,]),na.rm=T) != 0 & !noseg){
if(!NAs){output.tmp <- fast.lm.1(trait.matrix,strain.map[i,])}
if(NAs){output.tmp <- lm.1(trait.matrix,strain.map[i,])}
}
if( i != 1 & sum(abs(strain.map[i-1,]-strain.map[i,]),na.rm=T) == 0 & !noseg){
output.tmp <- output.tmp
}
eff.out[,i] <- output.tmp[,2]
pval.out[,i] <- output.tmp[,1]
}
if(!small){
colnames(eff.out) <- rownames(strain.marker)
rownames(eff.out) <- rownames(trait.matrix)
colnames(pval.out) <- rownames(strain.marker)
rownames(pval.out) <- rownames(trait.matrix)
}
if(small){
eff.out <- matrix(eff.out[1,],nrow=1)
rownames(eff.out) <- traitnames
colnames(eff.out) <- rownames(strain.map)
pval.out <- matrix(pval.out[1,],nrow=1)
rownames(pval.out) <- traitnames
colnames(pval.out) <- rownames(strain.map)
trait.matrix <- matrix(trait.matrix[1,],nrow=1)
rownames(trait.matrix) <- traitnames
colnames(trait.matrix) <- colnames(strain.map)
}
output <- NULL; output <- as.list(output)
output[[1]] <- round(pval.out,digits=2)
output[[2]] <- round(eff.out,digits=3)
output[[3]] <- trait.matrix
output[[4]] <- strain.map
output[[5]] <- strain.marker
names(output) <- c("LOD","Effect","Trait","Map","Marker")
return(output)
}
###Permutation function for single marker mapping
#Snoek & Sterken, 2017; Sterken, 2017
###input
#trait.matrix (traits in rows, genotypes in columns)
#strain.map (markers in rows, genotypes in columns)
#strain.marker (markers in rows, with columns name, chromosome and position)
#n.perm (the number of times each trait should be permutated)
###output
#a list with names: LOD, Effect, Trait_perm, Map, and Marker.
###Description
# map.1 is optimized for dataset in which many markers are un-informative
# these markers are not mapped but the p-vals and effect are copied from the previous marker.
QTL.map.1.perm <- function(trait.matrix,strain.map,strain.marker,n.perm){
if(missing(trait.matrix)){ stop("specify trait.matrix, matrix with traits in rows and genotypes in columns")}
if(missing(strain.map)){ stop("specify strain.map (genetic map)")}
if(ncol(strain.map) != dim(as.matrix(trait.matrix))[2] &
ncol(strain.map) != prod(dim(as.matrix(trait.matrix)))){ stop("number of strains is not equal to number of genotypes in map")}
if(missing(strain.marker)){ strain.marker <- data.frame(cbind(name=1:dim(strain.map)[1],chromosome=NA,bp=1:dim(strain.map)[1]))
rownames(strain.marker) <- 1:dim(strain.map)[1]}
if(missing(n.perm)){ n.perm <- 1000}
if(n.perm == 1 & dim(as.matrix(trait.matrix))[2] ==1){ stop("cannot conduct 1 permutation on 1 trait, set n.perm > 1")}
NAs <- sum(is.na(trait.matrix))>0 | sum(is.na(strain.map))>0
small <- FALSE
if(dim(as.matrix(trait.matrix))[1] ==1){
traitnames <- rownames(trait.matrix)
trait.matrix <- rbind(trait.matrix,trait.matrix)
rownames(trait.matrix) <- c(traitnames,"B")
small <- TRUE
}
trait.matrix <- trait.matrix[rep(1:nrow(trait.matrix),each=n.perm),]
trait.matrix <- permutate.traits(trait.matrix)
rownames(trait.matrix) <- paste(rownames(trait.matrix),1:n.perm,sep="_")
if(small){
trait.matrix <- trait.matrix[1:n.perm,]
}
eff.out <- matrix(NA,nrow(trait.matrix),nrow(strain.map))
pval.out <- matrix(NA,nrow(trait.matrix),nrow(strain.map))
for(i in 1:nrow(strain.map)){
noseg <- length(unique(strain.map[i,])) ==1
if(noseg){
output.tmp <- matrix(c(0,0),byrow=TRUE,ncol=2,nrow=nrow(trait.matrix))
}
if( i == 1 & !noseg){
if(!NAs){output.tmp <- fast.lm.1(trait.matrix,strain.map[i,])}
if(NAs){output.tmp <- lm.1(trait.matrix,strain.map[i,])}
}
if( i != 1 & sum(abs(strain.map[i-1,]-strain.map[i,]),na.rm=T) != 0 & !noseg){
if(!NAs){output.tmp <- fast.lm.1(trait.matrix,strain.map[i,])}
if(NAs){output.tmp <- lm.1(trait.matrix,strain.map[i,])}
}
if( i != 1 & sum(abs(strain.map[i-1,]-strain.map[i,]),na.rm=T) == 0 & !noseg){
output.tmp <- output.tmp
}
eff.out[,i] <- output.tmp[,2]
pval.out[,i] <- output.tmp[,1]
}
colnames(eff.out) <- rownames(strain.marker)
rownames(eff.out) <- rownames(trait.matrix)
colnames(pval.out) <- rownames(strain.marker)
rownames(pval.out) <- rownames(trait.matrix)
output <- NULL; output <- as.list(output)
output[[1]] <- round(pval.out,digits=2)
output[[2]] <- round(eff.out,digits=3)
output[[3]] <- trait.matrix
output[[4]] <- strain.map
output[[5]] <- strain.marker
names(output) <- c("LOD","Effect","Trait_perm","Map","Marker")
return(output)
}
###function for transforming output of QTL.map.1 to a dataframe
###input
#QTL.map.1 output
#small (logic) if a single trait was mapped set to TRUE
###output
#dataframe with the columns trait, qtl_chromosome, qtl_bp, qtl_marker, qtl_significance, and qtl_effect
###Description
#This transforms the relevant output to a dataframe, making handling in dplyr and plotting using ggplot easier
#This takes output of map.1 functions (list file with LOD, effect, Trait, and Map and converts it to a list.
#Optimized (a bit) by generating an empty matrix and filling the columns one-by-one.
#added an option for single traits
QTL.map1.dataframe <- function(map1.output,small){
if(missing(map1.output)){ stop("Missing map1.output ([[LOD]],[[Effect]],[[Trait]],[[Map]],[[Marker]])")}
if(missing(small)){ small <- FALSE}
if(class(map1.output[[1]])=="numeric"){
small <- TRUE
}
if(small){
###Split up to long format
output <- matrix(NA,length(map1.output[[1]]),6)
output[,1] <- rep(rownames(map1.output[[3]])[1],each=length(map1.output[[1]]))
output[,2] <- as.character(unlist(map1.output[[5]][,2]))
output[,3] <- as.numeric(as.character(unlist(map1.output[[5]][,3])))
output[,4] <- c(1:length(map1.output[[1]]))
output[,5] <- c(map1.output[[1]])
output[,6] <- c(map1.output[[2]])
output <- as.data.frame(output)
}
if(!small){
###Split up to long format
output <- matrix(NA,ncol(map1.output[[1]])*nrow(map1.output[[1]]),6)
output[,1] <- rep(rownames(map1.output[[1]]),each=ncol(map1.output[[1]]))
output[,2] <- rep(as.character(unlist(map1.output[[5]][,2])),times=nrow(map1.output[[1]]))
output[,3] <- rep(as.numeric(as.character(unlist(map1.output[[5]][,3]))),times=nrow(map1.output[[1]]))
output[,4] <- rep(c(1:ncol(map1.output[[1]])),times=nrow(map1.output[[1]]))
output[,5] <- c(t(map1.output[[1]]))
output[,6] <- c(t(map1.output[[2]]))
output <- as.data.frame(output)
}
###overwrite dataframe formats
output[,1] <- as.character(unlist(output[,1]))
output[,2] <- as.character(unlist(output[,2]))
output[,3] <- as.numeric(as.character(unlist(output[,3])))
output[,4] <- as.numeric(as.character(unlist(output[,4])))
output[,5] <- as.numeric(as.character(unlist(output[,5])))
output[,6] <- as.numeric(as.character(unlist(output[,6])))
colnames(output) <- c("trait","qtl_chromosome","qtl_bp","qtl_marker","qtl_significance","qtl_effect")
###Return
return(output)
}
###function for single marker mapping
#Snoek & Sterken, 2017; Sterken, 2017
###input
#map1.dataframe; output of the QTL.map1.dataframe function
#threshold; the significance treshold for calling QTL (in -log10(p))
#LOD.drop; the decrease in LOD-score compared to the peak to call the boundaries of the QTL. Standard is 1.5
###output
#the original dataframe with the added columns: qtl_peak, qtl_marker_left, qtl_bp_left, qtl_marker_right, qtl_bp_right
#these columns are only non-NA where a QTL peak is located.
###Description
#This takes output of mapping.to.list function and identiefies peaks and confidence intervals v1.3
#Based on a given threshold and a 1.5 LOD-drop.
#v1.2 made the function more efficient, the peak borders are selected in one processing
#v1.3 fixed the problem calling the peak if the peak borders the edge of the chromosome
#v1.4 fixed 'larger then, smaller then' bug
#v1.5 fixed 'incompatible types' bug; NA is not recognized as numeric
QTL.peak.finder <- function(map1.dataframe,threshold,LOD.drop){
if(missing(map1.dataframe)){ stop("Missing map1.dataframe (trait, qtl_chromosome, qtl_bp, qtl_marker, qtl_significance, qtl_effect)")}
if(missing(threshold)){ stop("Missing threshold")}
if(missing(LOD.drop)){ LOD.drop <- 1.5}
###locate the peak
peaks <- dplyr::filter(map1.dataframe, qtl_significance >= threshold-LOD.drop) %>%
dplyr::group_by(trait,qtl_chromosome) %>%
dplyr::filter(qtl_significance == max(qtl_significance),max(qtl_significance)>=threshold) %>%
dplyr::filter(abs(qtl_marker - mean(qtl_marker)) == min(abs(qtl_marker - mean(qtl_marker)))) %>%
dplyr::filter(qtl_marker == min(qtl_marker)) %>%
as.data.frame() %>%
dplyr::mutate(qtl_peak=qtl_marker)
###locate left boundary
peaks_left <- dplyr::filter(map1.dataframe,trait %in% unique(peaks$trait)) %>%
dplyr::group_by(trait,qtl_chromosome) %>%
dplyr::mutate(Border=ifelse((qtl_significance < max(qtl_significance)-LOD.drop & max(qtl_significance)>=threshold),T,
ifelse((qtl_bp == min(qtl_bp) & min(qtl_bp) %in% qtl_bp[qtl_significance >= max(qtl_significance)-LOD.drop & max(qtl_significance)>=threshold]),T,
ifelse((qtl_bp == max(qtl_bp) & max(qtl_bp) %in% qtl_bp[qtl_significance >= max(qtl_significance)-LOD.drop & max(qtl_significance)>=threshold]),T,F))),
peakloc=qtl_marker[which.max(qtl_significance)]) %>%
dplyr::filter(Border) %>%
dplyr::mutate(qtl_marker_left=as.numeric(ifelse(qtl_marker == max(qtl_marker[qtl_marker <= peakloc]),qtl_marker,NA)),
qtl_bp_left=as.numeric(ifelse(qtl_bp == max(qtl_bp[qtl_marker <= peakloc]),qtl_bp,NA)),
qtl_marker_right=as.numeric(ifelse(qtl_marker == min(qtl_marker[qtl_marker >= peakloc]),qtl_marker,NA)),
qtl_bp_right=as.numeric(ifelse(qtl_bp == min(qtl_bp[qtl_marker >= peakloc]),qtl_bp,NA))) %>%
dplyr::filter(!is.na(qtl_bp_left) | !is.na(qtl_bp_right)) %>%
dplyr::select(-Border,-peakloc)
peaks_right <- dplyr::filter(peaks_left,!is.na(qtl_bp_right))
peaks_left <- dplyr::filter(peaks_left,!is.na(qtl_bp_left))
###Merge the findings
peaks <- cbind(peaks,peaks_left$qtl_marker_left,peaks_left$qtl_bp_left,peaks_right$qtl_marker_right,peaks_right$qtl_bp_right)
###Prepare input
output <- dplyr::mutate(map1.dataframe,qtl_peak = NA,qtl_marker_left = NA,qtl_bp_left = NA,qtl_marker_right = NA,qtl_bp_right = NA)
###Merge peaks and input
output[match(paste(peaks[,1],peaks[,2],peaks[,3]),paste(output[,1],output[,2],output[,3])),7:11] <- peaks[,7:11]
output <- dplyr::mutate(output,qtl_bp_left=ifelse(qtl_bp_left > qtl_bp,qtl_bp,qtl_bp_left),qtl_bp_right=ifelse(qtl_bp_right<qtl_bp,qtl_bp,qtl_bp_right))
###return called peaks
return(output)
}
###function for determining variance explained by QTL
###input
#QTL.table.file (should contain marker locations: qtl_marker, qtl_marker_1, qtl_marker_2 and trait names)
#QTL.list.file (list with QTL data, should contain entries Map and Trait)
###output
#the R-squared is added to the QTL.table.file
###Description
# This function will add the sum of squares of the used model to a QTL table output, it also calculates the r-squared of the model
QTL.table.addR2 <- function(QTL.table.file,QTL.list.file){
if(missing(QTL.table.file)){ stop("Give a QTL summary table, should contain marker locations (qtl_marker, qtl_marker_1, qtl_marker_2) and trait names")}
if(missing(QTL.list.file)){ stop("Give the input on which was mapped (output of QTL.data.prep)")}
###R-squared function
R.sq <- function(x,y){return(cor(x,y,use="na.or.complete")^2)}
###True if single marker mapping
if("qtl_marker" %in% colnames(QTL.table.file)){
output <- NULL
for(i in 1:nrow(QTL.table.file)){
###get marker
mrk <- QTL.list.file$Map[as.numeric(as.character(unlist(QTL.table.file$qtl_marker[i]))),]
###get expression
y <- unlist(QTL.list.file$Trait[rownames(QTL.list.file$Trait)==as.character(unlist(QTL.table.file$trait[i])),])
###run model
output <- c(output,R.sq(y,mrk))
}
output <- data.frame(cbind(QTL.table.file,qtl_R2_sm=output))
}
###True if multiple marker mapping
if("marker_1" %in% colnames(QTL.table.file) & "marker_2" %in% colnames(QTL.table.file)){
output <- cbind(QTL.table.file,SS_int_mrk1=NA,SS_int_mrk2=NA,SS_int=NA,SS_int_res=NA,SS_add_mrk1=NA,SS_add_mrk2=NA,SS_add_res=NA)
for(i in 1:nrow(QTL.table.file)){
###get marker
mrk1 <- QTL.list.file$Map[QTL.table.file$marker_1[i],]
mrk2 <- QTL.list.file$Map[QTL.table.file$marker_2[i],]
###get expression
y <- unlist(QTL.list.file$Trait[rownames(QTL.list.file$Trait)==QTL.table.file$trait[i],])
###run model
output[i,colnames(output)%in%c("SS_int_mrk1","SS_int_mrk2","SS_int","SS_int_res")] <- anova(lm(y ~ mrk1*mrk2))$"Sum Sq"
output[i,colnames(output)%in%c("SS_add_mrk1","SS_add_mrk2","SS_add_res")] <- anova(lm(y ~ mrk1+mrk2))$"Sum Sq"
}
output <- dplyr::mutate(output,r2_add=(SS_add_mrk1+SS_add_mrk2)/(SS_add_mrk1+SS_add_mrk2+SS_add_res),r2_int=SS_int/(SS_int_mrk1+SS_int_mrk2+SS_int+SS_int_res))
}