R/read_zdata.R

In RZooRoH: Partitioning of Individual Autozygosity into Multiple Homozygous-by-Descent Classes

setClass(Class = "zdata",
         representation(genos = "matrix", bp = "vector", chrnames = "vector", chrbound = "matrix",
                        nind = "numeric", nsnps = "numeric", freqs ="vector", nchr ="numeric",
                        zformat ="character", sample_ids = "vector", map_unit = "character",
                        RecRate = "numeric")
)

is.zdata <- function (x)
{
  res <- (is(x,"zdata") & validObject(x))
  return(res)
}

#### function to get allele frequencies from genotype data

getfreq1 <- function(genoline) {
  n1=length(genoline)
  n2=length(genoline[genoline<3])
  if((n2/n1) > 0.00){
    f0 <- sum(genoline[genoline<3]) / (2*length(genoline[genoline<3]))
  }else{
    f0 <- 0.00
  }
  return(f0)
}

getfreq2 <- function(genoline, genformat){
  if(genformat == 'gp' | genformat== 'gl'){
    i1 <- seq(from=1,to=(length(genoline)),by=3)
    i2 <- seq(from=2,to=(length(genoline)),by=3)
    i3 <- seq(from=3,to=(length(genoline)),by=3)
    g11 <- genoline[i1]
    g12 <- genoline[i2]
    g22 <- genoline[i3]
    g0 <- g11 + g12 + g22
    if(genformat == 'gp'){g0[g0 < 0.05]=0} #### considered missing
    if(genformat == 'gl'){
      g11 <- 10**(-g11/10)
      g12 <- 10**(-g12/10)
      g22 <- 10**(-g22/10)
    }
    gT <- g11 + g12 + g22
    if(genformat == 'gl'){
      g11[g0 > 0] <- g11[g0 > 0]/gT[g0 > 0]
      g12[g0 > 0] <- g12[g0 > 0]/gT[g0 > 0]
      g22[g0 > 0] <- g22[g0 > 0]/gT[g0 > 0]
    }
    g0[g11 >= 0.33 & g12 >= 0.33 & g22 >= 0.33]=0 #### when coding missing as 0.33 / 0.33 / 0.33
    n0 <- length(g0[g0 >0])
    f0 <- sum(2*g11[g0 > 0] + 1*g12[g0 > 0])
    if(n0 > 0)f0 <- f0/(2*n0)
    if(n0 == 0)f0 <- 0.00
  }

  if(genformat == 'ad'){
    i1 <- seq(from=1,to=(length(genoline)),by=2)
    i2 <- seq(from=2,to=(length(genoline)),by=2)
    ad1 <- genoline[i1]
    ad2 <- genoline[i2]
    ad <- ad1 + ad2
    f0 <- sum(ad1[ad > 0]/ad[ad > 0])
    n0 <- length(ad[ad >0])
    if(n0 > 0)f0 <- f0/(n0)
    if(n0 == 0)f0 <- 0.00
  }
  return(f0)
}

#### function to estimate the allele frequencies from haplotype data

getfreq3 <- function(genoline) {
  genoline <- genoline[which(!is.na(genoline))]
  n1=length(genoline)
  n2=length(genoline[genoline<2])
  if((n2/n1) > 0.00){
    f0 <- sum(genoline[genoline<2]) / (length(genoline[genoline<2]))
  }else{
    f0 <- 0.00
  }
  return(f0)
}

#### functions to estimate the allele frequencies with the EM algorithm

getfreqem1 <- function(genoline){
  n <- length(genoline)/3
  f_est <- .Fortran("freqem1",as.double(genoline),as.integer(n),as.double(0))
  return(f_est[3][[1]][1])
}

getfreqem2 <- function(genoline){
  n <- length(genoline)/3
  f_est <- .Fortran("freqem2",as.double(genoline),as.integer(n),as.double(0))
  return(f_est[3][[1]][1])
}

getfreqem3 <- function(genoline){
  n <- length(genoline)/2
  f_est <- .Fortran("freqem3",as.integer(genoline),as.integer(n),as.double(0))
  return(f_est[3][[1]][1])
}

######### wrappers for the C versions of the functions #########

cfreq1 <- function(genos) {
  .Call("getfreq_gt", genos)
}

cfreq2 <- function(genos, genformat) {
  if (genformat %in% c("gp","gl")) return(.Call("getfreq_gp_gl", genos, as.integer(genformat=="gl")))
  if (genformat == "ad") return(.Call("getfreq_ad", genos))
}

cfreq3 <- function(genos) {
  .Call("getfreq_haps", genos)
}

######## wrappers for the C version of the EM versions #######

cfreqem1 <- function(genos) {
  .Call("getfreqem_gp", genos)
}

cfreqem2 <- function(genos) {
  .Call("getfreqem_gl", genos)
}

cfreqem3 <- function(genos) {
  .Call("getfreqem_ad", genos)
}

#'Read the genotype data file
#'
#'Read a data file and convert it to the RZooRoH format in a 'zooin' object
#'required for further analysis.
#'
#'@param genofile The name of the input data file. Note that the model is
#'  designed for autosomes. Other chromosomes and additional filtering (e.g.
#'  call rate, missing, HWE, etc.) should be performed prior to run RZooRoH with
#'  tools such as PLINK or bcftools for instance. The model works on an ordered
#'  map and ignores SNPs with a null position.
#'
#'@param min_maf The minimum allele frequency to keep variants in the analysis
#'  (optional / set to 0.00 by default to keep all markers). Values such as 0.01
#'  allows to exclude monomorphic markers that are not informative and to reduce
#'  the size of the data and computational costs. There is no marker exclusion
#'  on call rate. However, we expect that data filtering is done prior to
#'  RZooRoH with tools such as PLINK or vcftools.
#'
#'@param zformat The code corresponding to the format of the data file ("gt" for
#'  genotypes, "gp" for genotype probabilities, "gl" for genotype likelihoods in
#'  Phred scores, "ad" for allelic depths). For all these formats, markers are
#'  ordered per rows and individuals per columns. Variants should be ordered by
#'  chromosome and position. By default, the format is based on the Oxford/GEN
#'  format, and the first five columns are chromosome identification (e.g, "1",
#'  "chr1"), the name of the marker, the position of the marker in base pairs or
#'  better in cM multiplied when genetic distances are known (see the map_unit
#'  and RecRate options for more details), the first marker allele and the
#'  second marker allele. Information per individual varies according to the
#'  format. With the "gt" format we have one column per individual with 0, 1 and
#'  2 indicating the number of copies of the first allele (and 9 for missing).
#'  With the "gp" format we have three column per individual with the
#'  probabilities of genotype 11 (homozygous for the first allele), genotype 12
#'  and genotype 22 (this corresponds to the oxford GEN format). Similarly, with
#'  the "gl" format, we have three column per individual with the likelihoods
#'  for genotypes 11, 12 and 22 in Phred scores. Finally, with the "ad" format,
#'  we expect two columns per individual: the number of reads for allele 1 and
#'  the number of reads for allele 2. For these three last formats, missing
#'  values must be indicated by setting all elements to 0. If one of the columns
#'  is non-null for one individual, the genotype will be considered non-missing.
#'  Note that the marker alleles specified in columns 4 and 5 are not used.
#'
#'  For applications related to identity-by-descent (IBD) estimation, the ZooRoH
#'  model is applied on two phased chromosomes or two haploid chromosomes. Two
#'  additional formats are therefore available to provide haplotypes (or haploid
#'  data). The first is "vcf" which refers to a phased VCF (for example, the
#'  output from Beagle5). In that case, the VCF must contain only the phased
#'  haplotype information (e.g 0|1, 0|0, 1|1, etc). The phased VCF can only be
#'  used with diploid individuals (after phasing). By default, we assume that
#'  the first and second column indicate chromosome and SNP position and that
#'  the individual haplotypes information starts in column 10. The second format
#'  is "haps" that is similar to the "GEN" format with five columns (chromosome
#'  identification (e.g, "1", "chr1"), the name of the marker, the position of
#'  the marker in base pairs or better in cM multiplied by 1,000,000 when
#'  genetic distances are known, the first marker allele and the second marker
#'  allele) followed by the haplotypes (two columns for diploid individuals and
#'  one in case of haploid chromosome). The alleles are coded as 0 and 1. The
#'  "haps" format is the only format that can be used with haploid data.
#'
#'  It is possible to convert a PLINK PED/BED file or a VCF file to the RZooRoH
#'  genotype format. This can be done using tools such as PLINK or bcftools.
#'  Both of these tools can for example be used to convert to GEN format. With
#'  PLINK2, it is possible to convert the file to a dosage file (0/1/2), which
#'  is similar to the 'GT' format, but with six columns of metadata before the
#'  genotypes. Please note that the GEN format has recently been modified
#'  slightly, with some versions having 5 or 6 additional columns (with or
#'  without chromosome information), and the conversion options have changed
#'  with VCFtools.
#'
#'  For example, conversion to GT format is performed with PLINK2 using the
#'  command plink2 --bfile myinput --export Av --out myoutput. Then, the header
#'  line should be deleted. Ideally, NA genotypes should be replaced with 9,
#'  which can be achieved using a bash command:
#'
#'  sed -e 's/NA/9/g' myoutput > newfile
#'
#'  The 'newfile' can then be read with the zoodata function, but the poscol and
#'  supcol parameters (see below) must be adjusted to poscol=4 and supcol=6.
#'
#'  We also provide examples using previous versions of PLINK (version 1.9) and
#'  VCFtools (conversions are also possible with the new versions). As shown
#'  above, RZooRoH is flexible if there are additional columns, as users can
#'  adjust the reading using options such as 'chrcol', 'poscol' and 'supcol'
#'  (see below).
#'
#'  For ped files, PLINK (v1.9) can recode them to oxford gen format with plink
#'  --file myinput --recode oxford --autosome --out myoutput. The autosome
#'  option keeps only #'  SNPs on autosomes as required by RZooRoH.
#'
#'  For vcf files, bcftools can be used to recode a vcf to the oxford gen format
#'  with the convert option: bcftools convert -t ^chrX,chrY,chrM -g outfile
#'  --chrom --tag GT myfile.vcf. The --chrom option is important to obtain
#'  chromosome number in the first column. The tag option allows to select which
#'  field from the vcf file (GT, PL, GL or GP) is used to generate the genotype
#'  probabilities exported in the oxford gen format. The -t option allows to
#'  exclude chromosomes (this is an example and chromosome names must be adapted
#'  if necessary). The needed output data is then outfile.gen.
#'
#'  If some genotype probabilities are missing, with a value of "-nan", you must
#'  replace them with "0" (triple 0 is considered as missing). This can be done
#'  with this command:
#'
#'  sed -e 's/-nan/0/g' file.gen > newfile.gen
#'
#'  It is important to properly encode missing records, as some missing symbols
#'  are not recognised and will generate errors. Note also that some software
#'  recodes missing values as 0.333, 0.333 and 0.333. This was not previously
#'  considered as missing data. This has no effect on emission probabilities,
#'  but slightly impacts the estimation of allele frequencies when the freqem
#'  option is FALSE. If the three genotype probabilities are now above 0.33, the
#'  genotype is considered missing. Ideally, missing values should be set to 0,
#'  0, 0 (for the 'gp' format).
#'
#'  As for genotype formats, tools such as bcftools can be used to format files
#'  in the 'haps' or phased VCF format.
#'
#'@param chrcol An optional argument that indicates the column number where the
#'  chromosome information is indicated (first column by default for all
#'  formats).
#'
#'@param poscol An optional argument that indicates the column number where the
#'  marker position is indicated (third column by default for all formats except
#'  phased vcf where it is the second column). See also
#'
#'@param supcol An optional argument that indicates the number of additional
#'  columns before the individuals genotypes or haplotypes (five columns are
#'  expected by default as described in the zformat argument description, except
#'  for phased VCF where this value is set to 9). Note that the function
#'  requires at least two information columns: the chromosome number and the
#'  marker position.
#'
#'@param haploid An optional argument that indicates whether we work on an
#'  haploid organism or chromosome (default = FALSE). This is only compatible
#'  with the 'IBD' option from the zoorun function (we don't estimate HBD in
#'  haploid data!). In the case you use haploid data, the only possible format
#'  is "haps".
#'
#'@param allelefreq A vector with allele frequencies for the first marker allele
#'  (optional). By default, the allele frequencies are estimated from the data.
#'  The option allows to skip this computation or to provide external allele
#'  frequencies estimated by another method or on another data set.
#'
#'@param freqem A logical indicating whether allele frequencies should be
#'  estimated with an EM algorithm (TRUE by default). Simpler approaches are
#'  used when set to FALSE. The approach is ignored with the GT, VCF and HAPS
#'  formats. For high confidence genotypes (e.g., genotyping arrays,
#'  high-coverage sequencing data), it is not necessary to use this EM approach
#'  as genotypes or haplotypes are known. The approach is more relevant with
#'  low-fold sequencing for example, and more so with the PL or GP format (the
#'  approximation with the AD format being closer to the EM).
#'
#'@param samplefile A file with names of the samples (optional). It must match
#'  with the number of genotyped individuals. If none is provided, the position
#'  in the genofile is used as ID.
#'
#'@param map_unit Can be "bp" or "cM". By default the variant positions are
#'  assumed to be the physical positions in base pairs (bp). If the user knows
#'  the map in centimorgans (cM), they can provide the positions in cM and set
#'  the map_unit to "cM". Note that the maximum resolution is 1e-6 cM: map
#'  positions will be rounded to this resolution, and two markers separated by
#'  less than 1e-6 cM may therefore end up at the same rounded position, in
#'  which case no recombination will be possible between them. If the positions
#'  are provided in bp but the cM-to-Mb ratio is not equal to 1, the user can
#'  specify the correct ratio using the 'RecRate' option (default: 1 cM/Mb).
#'
#'@param RecRate This gives the number of centimorgans (cM) per megabase (Mb).
#'  By default this value is set to 1 cM per Mb. Therefore, the rates of the
#'  exponential distributions correspond to rates of coacenstry changes per
#'  Morgan, and the expected length of HBD segments is equal to 100 cM / R. This
#'  value will be used to convert the physical position to cM. Note that the
#'  maximum resolution in cM is 1e-6 cM: map positions will be rounded to this
#'  resolution, so two markers separated by less than 1e-6 cM after conversion
#'  (which becomes more likely as the cM-to-Mb ratio decreases below 1) may end
#'  up at the same rounded position, in which case no recombination will be
#'  possible between them. This option is not used when the positions are
#'  provided directly in cM (when "map_unit" is set to "cM").
#'
#'@return The function return a zooin object called containing the following
#'  elements: zooin@@genos a matrix with the genotypes, genotype probabilities
#'  or haplotypes, zooin@@bp an array with marker positions, zooin@@chrbound a
#'  matrix with the first and last marker number for each chromosome,
#'  zooin@@nind the number of individuals, zooin@@nsnps the number of markers
#'  conserved after filtering for minor allele frequency, zooin@@freqs an array
#'  with the marker allele frequencies, zooin@@nchr the number of chromosomes,
#'  zooin@@zformat the format of the data ("gt","gp","gl","ad","vcf","haps") and
#'  zooin@@sample_ids (the names of the samples).
#'
#' @examples
#'
#' # Get the name and location of example files
#'
#' myfile1 <- system.file("exdata","genoex.txt",package="RZooRoH")
#' myfile2 <- system.file("exdata","genosim.txt",package="RZooRoH")
#'
#' # Load your data with default format into a zooin object named "data1":
#'
#' data1 <- zoodata(myfile1)
#'
#' # Load the first data file with default format and filtering out markers with MAF < 0.02
#' # into a zooin object called "data1frq002":
#'
#' data1frq002 <- zoodata(myfile1, min_maf = 0.02)
#'
#' # Load the first data file with default format, with external allele frequencies
#' # (here a random set we create) and filtering out markers with MAF < 0.01:
#'
#' myrandomfreq <- runif(14831)
#' data1c <- zoodata(myfile1, allelefreq = myrandomfreq, min_maf = 0.01)
#'
#' # Load the second data file and indicate your own format (chromosome number in column 1,
#' # map position in column 2, 4 columns before genotypes) and filtering out markers with
#' # MAF < 0.01. The created zooin object is called "Sim5":
#'
#' Sim5 <- zoodata(myfile2, chrcol = 1, poscol =2, supcol = 4, min_maf = 0.01)
#'
#'@export
#'@import data.table

zoodata<-function(genofile, min_maf=0.00, zformat = "gt", chrcol = 1, poscol = 0, supcol = 0,
                  haploid = FALSE, allelefreq = NULL, freqem = TRUE, samplefile = NA,
                  map_unit="bp", RecRate = 1){

  if(supcol==0){
    if(zformat=="gt" | zformat=="gp" | zformat=="gl" | zformat=="ad" | zformat=="haps"){supcol=5}
    else if(zformat=="vcf"){supcol=9}
  }

  if(poscol==0){
    if(zformat=="gt" | zformat=="gp" | zformat=="gl" | zformat=="ad" | zformat=="haps"){poscol=3}
    else if(zformat=="vcf"){poscol=2}
  }

  if(!is.na(samplefile)){
    mysamples <- read.table(samplefile,header=FALSE)
    badnames <- mysamples$V1[grepl("_", mysamples$V1, fixed = TRUE)]
    if(length(badnames) > 0) stop("These individual identifiers contain a '_', which is not allowed: ", paste(badnames, collapse = ", "))
  }

  if(map_unit != "cM" & map_unit != "bp"){stop("The option map_unit must be cM or bp!")}

#  if(freqem){
#    if(zformat=="gt" | zformat=="haps" | zformat=="vcf"){
#      print("Freqem is not usefull with this format, the option is ignored.")
#    }
#  }

  zooin <- new("zdata")
  max_maf <- (1 - min_maf)
  freqs <- allelefreq
#  na_strings <- if (zformat == "vcf") c("NA","na") else c("NA","na",".")
  if(zformat == "vcf"){
#    genfile <- fread(genofile,data.table=FALSE,header=TRUE, na.strings=na_strings)
    genfile <- fread(genofile,data.table=FALSE,header=TRUE)
  }else{
#    genfile <- fread(genofile,data.table=FALSE, na.strings=na_strings)
    genfile <- fread(genofile,data.table=FALSE)
  }

  if(zformat != "gt" && zformat != "gp" && zformat != "gl" && zformat != "ad" && zformat != "vcf" && zformat !="haps"){
    stop(paste("Unknown data format (zformat) ::",zformat,"\n",sep=""))
  }

  if(haploid & zformat != "haps"){
    stop("You must use the haps format with haploid data!")
  }

  print(c('Number of positions in original file ::',length(genfile[,poscol])))

  chr <- as.character(genfile[, chrcol])
  bp <- genfile[, poscol]

  # Test maximum chromosome size
  max_bp_raw <- max(bp, na.rm=TRUE)
  if(max_bp_raw > 1000000000){
    stop("The largest position in your data (", max_bp_raw, " bp) exceeds 1000 Mb, ",
         "which is maximum value allowed.")
  }

  # Test to remain within the representable range of integer positions
  if(map_unit == "bp" & RecRate != 1){
    scaled_max <- max_bp_raw * RecRate
    if(scaled_max > 2000000000){
      stop("The scaled positions (max position x RecRate = ", round(scaled_max), ") exceed the maximum ",
           "allowed value (2e9) by the package. ",
           "As an alternative to scaling the positions, you can keep map_unit = \"bp\" (RecRate = 1, ",
           "default) and instead scale the rates of your model: multiply each rate by RecRate before ",
           "passing it to the model, and interpret the resulting expected HBD segment length as ",
           "100 / rate (in cM), using the original, unscaled rate values.")
    }
  }

  if(map_unit == "bp"){
    if(RecRate != 1) bp <- as.integer(round(bp*RecRate))
  }
  if(map_unit == "cM"){bp <- as.integer(round(bp*1000000))}

  # check for duplicate positions after rounding or conversion - paste FALSE (for first)
  dup_mask <- c(FALSE, diff(bp) == 0 & chr[-1] == chr[-length(chr)])

  if(any(dup_mask)){
    n_dup <- sum(dup_mask)
    dup_idx <- which(dup_mask)
    n_show <- min(10, length(dup_idx))
    example_lines <- paste0("  chr=", chr[dup_idx[1:n_show]], ", pos=", bp[dup_idx[1:n_show]],
                            " (marker index ", dup_idx[1:n_show], ", identical to previous marker)")
    warning(paste0(
      "After rounding to the nearest integer (maximum resolution: 1e-6 cM, and/or after applying RecRate), ",
      n_dup, " marker(s) share an identical position with at least one other marker on the same chromosome. ",
      "Recombination (HBD class change) will not be possible between such markers, effectively treating them as a single block. ",
      "We recommend avoiding this situation (e.g. by removing markers that are too close.\n",
      "First affected marker(s):\n",
      paste(example_lines, collapse = "\n")
    ), call. = FALSE)
  }


  if (zformat == "gt") {
    nind <- length(genfile[1,])-supcol
    genos <- as.matrix(genfile[, (supcol+1):(nind+supcol)])
    rm(genfile)
    gc();malloc_trim_r()
    ### No conversions of 'NA' because properly handled in pemission (and frequency estimation)
#    if(is.null(allelefreq)){freqs <- apply(genos, 1, getfreq1)}
    if(is.null(allelefreq)){freqs <- cfreq1(genos)}
  }

  if (zformat == "gp") {
    nind <- (length(genfile[1,])-supcol)/3
    genos <- as.matrix(genfile[, (supcol+1):(3*nind+supcol)])
    rm(genfile)
    gc();malloc_trim_r()
    invisible(.Call("set_na_to_zero_triplets", genos))
    #    if(is.null(allelefreq) & !freqem){freqs <- apply(genos, 1, getfreq2, genformat = "gp")}
    if(is.null(allelefreq) & !freqem){freqs <- cfreq2(genos, "gp")}
#    if(is.null(allelefreq) & freqem){freqs <- apply(genos,1,getfreqem1)}
    if(is.null(allelefreq) & freqem){freqs <- cfreqem1(genos)}
  }

  if (zformat == "gl") {
    nind <- (length(genfile[1,])-supcol)/3
    genos <- as.matrix(genfile[, (supcol+1):(3*nind+supcol)])
    rm(genfile)
    gc();malloc_trim_r()
    invisible(.Call("set_na_to_zero_triplets", genos))
    #    if(is.null(allelefreq) & !freqem){freqs <- apply(genos, 1, getfreq2, genformat = "gl")}
    if(is.null(allelefreq) & !freqem){freqs <- cfreq2(genos, "gl")}
#    if(is.null(allelefreq) & freqem){freqs <- apply(genos,1,getfreqem2)}
    if(is.null(allelefreq) & freqem){freqs <- cfreqem2(genos)}
  }

  if (zformat == "ad") {
    nind <- (length(genfile[1,])-supcol)/2
    genos <- as.matrix(genfile[, (supcol+1):(2*nind+supcol)])
    rm(genfile)
    gc();malloc_trim_r()
    invisible(.Call("set_na_to_zero_pairs", genos))
    #    if(is.null(allelefreq) & !freqem){freqs <- apply(genos, 1, getfreq2, genformat ="ad")}
    if(is.null(allelefreq) & !freqem){freqs <- cfreq2(genos, "ad")}
#    if(is.null(allelefreq) & freqem){freqs <- apply(genos, 1, getfreqem3)}
    if(is.null(allelefreq) & freqem){freqs <- cfreqem3(genos)}
  }

  if (zformat == "vcf") {
    #genfile[genfile=="."] <- NA
    pvcf <- as.data.frame(genfile[, (supcol+1):ncol(genfile)])
    pvcf[pvcf=="."] <- ".|."
    rm(genfile)
    gc();malloc_trim_r()
    haps1 <- lapply(pvcf, function(col) sapply(strsplit(as.character(col), "\\|"), "[[", 1))
    haps2 <- lapply(pvcf, function(col) sapply(strsplit(as.character(col), "\\|"), "[[", 2))
    haps1[haps1=="."] <- NA
    haps2[haps2=="."] <- NA
    rm(pvcf)
    hap1 <- as.data.frame(lapply(haps1, as.integer))
    hap2 <- as.data.frame(lapply(haps2, as.integer))
    nind <- length(hap1[1,])
    genos <- matrix(9,nrow(hap1),2*nind)
    genos[,seq(1,2*nind-1,by=2)]=as.matrix(hap1)
    genos[,seq(2,2*nind,by=2)]=as.matrix(hap2)
    rm(hap1,hap2)
    gc();malloc_trim_r()
#    if(is.null(allelefreq)){freqs <- apply(genos, 1, getfreq3)}
    if(is.null(allelefreq)){freqs <- cfreq3(genos)}
  }

  if (zformat == "haps") {
    nhap <- length(genfile[1,])-supcol
    if(haploid){nind <- nhap}
    if(!haploid){nind <- nind <- nhap/2}
    genfile[genfile=="."] <- NA
    hapin <- as.data.frame(lapply(genfile[, (supcol+1):(nhap+supcol)], as.integer))
    genos <- as.matrix(hapin)
    rm(genfile)
    gc();malloc_trim_r()
#    if(is.null(allelefreq)){freqs <- apply(genos, 1, getfreq3)}
    if(is.null(allelefreq)){freqs <- cfreq3(genos)}
  }

  gc();malloc_trim_r()

  filterin <- (freqs >= min_maf & freqs <= max_maf & bp > 0)
  chr <- chr[filterin]
  bp <- bp[filterin]
  ff <- freqs[filterin]
  nchr <- length(unique(chr))

  ### solution avoiding loop
  chrbound <- matrix(0, nchr, 2)
  chrnames <- array("",nchr)
  chrbp2 <- c(chr[2:(length(chr))],"chr+")
  chrbound[,2] <- which(chr != chrbp2)
  chrnames <- chr[which(chr != chrbp2)]
  chrbound[1,1] <- 1
  if(nchr > 1){for (i in 2:nchr){chrbound[i,1] <- chrbound[(i-1),2] + 1 }}

  nsnps <- chrbound[nchr,2]

  zooin@nind <- nind
  zooin@genos <- genos[filterin,]
  rm(genos,freqs,filterin)
  gc();malloc_trim_r()
  zooin@nsnps <- nsnps
  zooin@freqs <- ff
  zooin@chrnames <- chrnames
  zooin@chrbound <- chrbound
  zooin@nchr <- nchr
  zooin@bp <- bp
  zooin@zformat <- zformat
  if(!is.na(samplefile)){
    mysamples <- read.table(samplefile,header=FALSE)
    if(length(mysamples$V1) != nind){
      print("The number of sample IDs does not match the number of genotypes !")
    }
    zooin@sample_ids <- as.character(mysamples$V1)
  } else {zooin@sample_ids <- as.character(seq(1:nind))}
  print(c('Number of positions after MAF filtering  ::',zooin@nsnps))

  gc();malloc_trim_r()

  return(zooin)

}