Import 1segSite haplo

[david20011999]

I have tried to import an haplotype of just one segSite but an error occurs. It is due when each haplotype for each chromosome is selected (file R/importData.R, line 279) you are substraccting elements from a matrix. Nevertheless, by default R convert it to a simpler object. It just happens when it is substracting 1 segSite (1 column) and it becomes a vector. Being a vector instead a matrix newMapPop function aborts. You can solve it easlily with the function as.matrix or adding haplo[,take, drop = F] on the mentioned line.

Thanks ;)

If you want to test it, here you have an example:

nsegsites <- 1 

#Creating the map
genMap = data.frame(
  marker     = paste0("x", 1:nsegsites),
  chromosome = 1,
  position   = seq(from = 0, to = (nsegsites / 2) - 0.5, by = 0.5)
)

#number of individuals in the founder population
nind <- 1000

#allele frequencies
p = 0.3
q = 1 - p

#genotype frequencies (in this case HWE)
geno <- c(rep("11", times = round(p * p * nind)),
          rep("10", times = round(p * q * nind)),
          rep("01", times = round(q * p * nind)),
          rep("00", times = round(q * q * nind)))
#if there are multiple segsites, paste whole individual genotypes

#sampling genos
geno <- data.frame(V1 = sample(geno, size = nsegsites * nind, replace = F))

# Load the haplotypes and formatting
aux <- 1:(nsegsites * 2)
library(stringr)
haplo12 <- data.frame(str_split_fixed(geno$V1, "", 2 * nsegsites))
haplo1 <- data.frame(V1 = as.numeric(haplo12[,aux[aux %% 2 != 0]]))
haplo2 <- data.frame(V1 = as.numeric(haplo12[,aux[aux %% 2 == 0]]))
aux <- c((1:nind)*2-1,(1:nind)*2)
haplo <- rbind(haplo1,haplo2)
haplo <- data.frame(haplo [aux,])
haplo <- as.matrix(haplo)

# Assign marker names to the haplotypes
colnames(haplo) = genMap$marker

# Load a pedigree 
ped = data.frame(id = 1:nind,
                 mother = rep(0, nind),
                 father = rep(0, nind))

# Create the founder population
founderPop = importHaplo(haplo  = haplo,
                         genMap = genMap,
                         ploidy = 2,
                         ped    = ped)

# Set simulation parameters
SP = SimParam$new(founderPop)

# Load QTL effects 
qtlEffects = data.frame(marker = c("x1"),
                        aditiveEffect = c(1),
                        dominanceEffect = c(0.5))

# Import in SimParam
SP$importTrait(markerNames = qtlEffects$marker,
               addEff = qtlEffects$aditiveEffect,
               domEff = qtlEffects$dominanceEffect,
               name = "Trait1")
SP$setSexes("yes_sys")

[david20011999]

I've not sent a Pull Request due it is a small thing

[gregorgorjanc]

@david20011999 cool catch! Do provide a PR so you simplify the work for @gaynorr ;)

[gaynorr]

@david20011999, thanks for reporting this. Sadly, it's not the first time I've been hit with this bug.