Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Thierry Gosselin

Jérome Goudet

Thibaut Jombart

Klaus Schliep

J Ér

Ôme Goudet

Linking genomics and population genetics with R

Keywords: multivariate analysis, NGS, R, SNP, VCF

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

An exploration into the history of genetics reveals a complex pattern of interconnected concepts, hypotheses, and research programs. A remarkable historical fact is that population geneticists laid the foundations of their discipline several decades before the physical support of heredity was known. [START_REF] Avery | Studies on the chemical nature of the substance inducing transformation of pneumococcal types. Induction of transformation by a desoxyribonucleic acid fraction isolated from Pneumococcus type III[END_REF] demonstrated that DNA was the molecule coding for genes, a fact that was not immediately accepted by their contemporary geneticists for several reasons, perhaps the main one being that, at this time, it was widely accepted that genes most likely consisted of proteins [START_REF] Lederberg | The transformation of genetics by DNA: an anniversary celebration of Avery, MacLeod and McCarty[END_REF][START_REF] Deichmann | Early responses to Avery et al.'s paper on DNA as hereditary material[END_REF]). Population genetics was a mature field at the time of Avery et al.'s finding as witnessed by the long-dated publication of seminal books on the subject (Fisher 1930, e.g.,). Thus, population and molecular genetics have had separate histories at some time and this certainly had an impact on how empirical studies were conducted or motivated. During almost half a century, population genetic studies were based on genotyping individuals on a limited number of loci, especially in conservation-oriented research [START_REF] Ekblom | A field guide to whole-genome sequencing, assembly and annotation[END_REF]. A crucial preliminary of these studies was to find the 'right genetic marker for the study' [START_REF] Sunnucks | Efficient genetic markers for population biology[END_REF].

At the end of the century, genomic data from human populations started to accumulate making possible to address new questions in the emerging field of population genomics, even though at this time we were "a long way from knowing all the SNPs, even in any given population" mainly because "exhaustive typing is currently prohibitive" [START_REF] Goldstein | Population genomics: linkage disequilibrium holds the key[END_REF]. Nevertheless, it appeared clearly that this wealth of human genetic data could not be interpreted without traditional population genetic concepts such as linkage disequilibrium, genetic drift, or coalescent [START_REF] Jorde | Population genomics: a bridge from evolutionary history to genetic medicine[END_REF].

Another result that came out at this time was the peculiarity, from an evolutionary point of view, of human populations compared to other species: low overall genetic diversity, considerable continental-level homogeneity, and small effective population size [START_REF] Jorde | Population genomics: a bridge from evolutionary history to genetic medicine[END_REF]. These features associated with the relatively simple biology of humans (lack of polyploidy, long lifespan) and the considerable interest in medically oriented applications of genomic research resulted in a lot of developments for handling and analysing population genomic data. However, these tools did not appear flexible enough to be applied to organisms with more complex biologies.

In the early twenty-first century, next generation sequencing technologies (NGS) have made possible to access the genotype of an individual at a very large number of loci and even its complete genome sequence [START_REF] Luikart | The power and promise of population genomics: from genotyping to genome typing[END_REF]. Even though there are still few species with their genome sequenced, several approaches have been proposed to apply NGS to any living species [START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF]. Furthermore, the decreasing costs of NGS make increasingly easy to generate data on a large number of loci and to apply these technologies in studies on natural populations [START_REF] Narum | Genotyping-by-sequencing in ecological and conservation genomics[END_REF][START_REF] Andrews | Recent novel approaches for population genomics data analysis[END_REF]. This is a radical change in point of view for population geneticists compared to almost a century of practice of their discipline. In his review, [START_REF] Sunnucks | Efficient genetic markers for population biology[END_REF] listed three components in population genetics studies: genetic markers, statistical methods, and computer programs. We believe that the NGS revolution has solved, at least in principle, the problem of genetic markers. The availability of vast amount of data clearly points to new developments in software. Furthermore, the perspective of merging population and molecular genetics calls for a unified approach to software development.

In this paper, we review the current state of progress in the analysis of population genomic data using R. We illustrate some recent developments with three recently published data sets: the human Y chromosome with 62,042 loci for 1233 individuals (The 1000 Genomes Project Consortium 2015), a set of 1,055,818 phased genotypes for 121 fruit flies (Drosophila melanogaster, Kao et al. 2015ab), and a set of 61,951 SNPs for 113 dolphins (Lagenorhynchus spp., Fernández et al. 2016ab). In the next section, we summarize the challenges ahead and explain how R can help meeting them. The following two sections explore some specific issues related to handling and analyses of large data sets in the context of population genomics, and how the R packages presented in this paper integrate with other genomic bioinformatics software. We argue that R solves the issue of software in population genetics in the current context of NGS and in the perspective of a constant flux of genomic data. In the last section, we discuss the future developments of statistical methods which will, hopefully, complete the merging of population and molecular genetics.

Challenges

NGS contrast sharply with their predecessor (the Sanger method) in two important aspects: they generate far larger quantities of data, which is clearly a source of practical difficulties for data storage and analysis [START_REF] Stephens | Big data: astronomical or genomical?[END_REF], and they evolve much more quickly. The tremendous potential applications of NGS have stimulated a lot of research and development [START_REF] Erlich | A vision for ubiquitous sequencing[END_REF]. A new technology or platform can be created and made available rapidly, but tests and validation for reproducibility take time.

Meanwhile, competing laboratories or companies may initiate another cycle of testing, assessment, and validation for a new method. Furthermore, the commercialisation of a technology and the training of laboratory staff to use this technology on an appropriate scale also take time. All these time-lags result in a form of inertia in the deployment of NGS in research laboratories. A now famous example is provided by the 454 platform which was discontinued in 2013, just six years after being acquired by a large company-some emerging platforms will probably have a shorter lifespan. It appears that we are now in a situation where we can start sampling in the field without knowing exactly what genotyping or sequencing method will be used in the laboratory.

Nowadays, it is possible to genotype a single individual at thousands or even millions of loci, and population genomic studies typically do this for a few individuals (n), at least far less than the number of loci that are analysed (p). This 'small n, large p' problem raises some difficulties, especially in statistical analyses (see section below on multivariate methods). This situation contrasts with traditional population genetic studies which generally used sampling protocols with n > p. With complete genome sequencing the value of p is theoretically bounded by the size of the genome; however, the number of variants (the sites in the genome which are variable among individuals in the sample) depends on the sampling protocol. For instance, the 1000 Genomes Project found almost 38,000,000 variants from 1092 individuals sampled from 14 populations (The 1000 Genomes Project Consortium 2012); this number increased to more than 88,000,000 when 2504 individuals were sampled from 26 populations (The 1000 Genomes Project Consortium 2015). With non-model organisms, we certainly have to expect even more variation since most species are distributed among several subspecies over their range.

The fast pace of NGS development is usually accompanied with a rapid development of bioinformatics tools often resulting in the appearance of a wide variety of data formats. A problem is that very different technologies are used by the different NGS platforms, so that raw data and their assembly differ a lot even though the final goal is the same. For instance, though the FASTQ file format has been adopted as the de facto standard for storing raw DNA reads, the emerging NGS platforms in the late 2000's used subtle differences in their respective FASTQ files making the development of common bioinformatics tools problematic (see: https://en.wikipedia.org/wiki/FASTQ_format for a review of the variation in FASTQ formats). Another difficulty comes from the fact that population genetics software have traditionally been developed independently of each other and often created their own data file formats, resulting in burdensome conversions necessary along bioinformatics workflows (see [START_REF] Lischer | PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs[END_REF], for a tool to work around this problem). Besides, because of the conceptual separation between population and molecular genetics, most software for population genetics has been developed with classical genetic markers and cannot handle large genomic data sets simply because of memory limitations.

Why R

The quest for standardization is a common and recurrent issue in software development, and applications for data analysis are no exception to this. The 1990's were a key period in this respect: efficient hardware (e.g., 32-bit processors) were widely available, Internet became common, and a new generation of software started to spread. R was developed in this exciting context [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF], and started to be progressively adopted by a wide community of scientists for analysing their data, but also for developing new applications [START_REF] Vance | Data analysts captivated by R's power[END_REF]. Several reasons can explain the success of R: we can cite a few that are relevant for population genetics. R is available for all common platforms and operating systems. Its fast graphical capabilities make it ideal for exploratory analyses. Some generic features of R are extremely useful when manipulating data, such as factors or data structures indexed with names. R is modular so that building new packages or pipelines, or conducting simple analyses is straightforward. R can be easily interfaced with code written in C, C++, and some other languages, so that computations which are not efficient in R can be recoded in these languages without losing the advantages of the R environment. R has a standard documentation system which makes very easy to find the relevant information when using a new package or a new function (of course, providing the package developers have written the documentation appropriately). R has been widely adopted in statistics courses in many countries so that most students have a fair knowledge of this language.

Finally, and not the least, several user interfaces (Emacs+ESS, JGR, R-commander, RKWard, RStudio, Tinn-R, among others) are available which considerably ease the work flow while accomodating each researcher's personal taste.

Naturally, all these features (and others) apply to the analysis of population genetic data which has led to the progressive development of a software environment that is attractive to a wide community of population geneticists. Two remarkable features of this software are worth mentioning here. First, these R packages consider all kinds of genetic markers: DNA sequences, SNP, micro-satellites, allozymes. Second, all situations can be handled: genotypes with single or multiple loci, locus with more than two alleles, all levels of ploidy including cases where ploidy level varies within the same data set, phased and unphased genotypes (Table 1). By contrast to some applications which have been developed specifically for diploid organisms with simple life histories and relatively low polymorphism, R packages are developed to meet the needs of a wide community of molecular ecologists. This underlines the philosophy of R to be inclusive in its developments: from the simplest to the most complex or challenging situations.

Data and files

R has a wide range of tools to handle and analyse DNA sequences (mostly in the package ape; [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF], and traditional allelic data (mostly in the package adegenet; [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]. In this section, we detail how some other kinds of data that are more specific to NGS can be handled with R.

SNP Data

Molecular genetic studies have revealed that variation in genetic materials can take many forms. In humans, among 88,332,015 genetic variants identified in a sample of 2504 individuals, 95.53% were biallelic single nucleotide polymorphism (SNP), 4.07% insertions-deletions (indels), 0.33% multiallelic SNPs, and 0.07% structural variants (The 1000 Genomes Project Consortium 2015). Thus, biallelic SNPs are by far the most common form of genetic variation in natural populations (see also the examples below).

Among the technologies developed to acquire this kind of variants, two are most frequently used: sequencing and genotyping arrays, the second one being able to acquire only biallelic SNPs (which are often called 'SNPs' in short). The approach based on genotyping SNPs can be seen as an intermediate stage between the use of traditional population genetic markers and the more modern genomic sequencing-based approaches which consider the whole set of genetic variation in populations.

Typically, SNP data files store genotypes of individuals at many loci in a matrix form with additional information for each locus such as its position along the chromosome. A difficulty with these files is that there is no widely accepted standard format. For instance, a genotype can be coded in the usual form (e.g., A/A, A/T, or T/T), or as the number of minor alleles in the genotype (0, 1, or 2, if T is the minor allele). However, it appears that the VCF format has recently been adopted as a standard for all kinds of SNP data files (see the example with the dolphin data).

The package adegenet provides several tools to read and handle SNP data [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF]. The core of these tools is the class "genlight" designed to store SNP data in a compact way: for example, such an object with one million SNPs for 100 individuals uses about 16 Mb of RAM (or active memory) which is less than 1% of the available RAM on most modern computers. Besides, adegenet has several functions specially designed for the analysis of "genlight" data, such as glPca which performs a principal components analysis (PCA; see below), glPlot which plots a set of SNPs, or glSim which simulates efficiently SNP data. The function read.PLINK in the same package can read some SNP files. More generally, SNPs are stored in plain-text files in tabular form, so standard R functions (e.g., read.table) can be used to parse such files.

VCF files

Over the last few years, the Variant Call Format (VCF) has appeared as a standard for storing population genetic data from NGS sequencing or genotyping [START_REF] Danecek | The variant call format and VCFtools[END_REF]. By contrast to SNP data files, VCF files can contain data on all types of genetic variation. VCF files are text files storing a lot of information, including on the process used to create it from the raw sequencing reads (variant calling). The genetic data are stored in a matrix-like layout where the individuals are the columns and the loci are the rows.

Several R packages have tools to analyse VCF files: two of them, stacksr and vcfR, are described in this issue. Bioconductor has the package VariantAnnotation [START_REF] Obenchain | VariantAnnotation: a Bioconductor package for exploration and annotation of genetic variants[END_REF] for annotating VCF files. We focus here on the functions in pegas because this package is dedicated to population genetics [START_REF] Paradis | pegas: an R package for population genetics with an integrated-modular approach[END_REF]. vcfR provides functions to assess the quality of the VCF file [START_REF] Knaus | VcfR: an R package to manipulate and visualize VCF format data[END_REF] together with functions for data conversion (Table 2).

VCF files can be very big, so a strategy is to first scan the files before reading more detailed information such as the genotypes. pegas follows this strategy by providing several functions: the main ones are VCFloci and read.vcf. The first function scans a VCF file and extracts the information on the loci; the second one reads the individual genotypes with the possibility to select a subset. VCFloci uses a specific mechanism to read big files: the file is scanned sequentially by chunks of one billion bytes (1 Gb; the chunk size can be modified by the user). The function then stores a small table with the name of the file, the size of each chunk adjusted to meet the nearest end of a line, and the number of loci in each chunk. With this mechanism, not more than 1 Gb of RAM is used at the same time whatever the size of the file. The running time increases roughly linearly with the file size: it takes about 20 sec to scan an 8 Gb file (compressed or not) with about 500,000 loci. For smaller files, the running time depends on whether the file is compressed or not: for a 200 Mb file with 62,000 loci, this takes 1 sec or 0.2 sec, respectively.

The output of VCFloci is a data frame (R's usual data table) with nine columns where each row is a locus and the columns are the information for each locus from the VCF file. The information correspond to the mandatory VCF header fields (see [START_REF] Danecek | The variant call format and VCFtools[END_REF]: chromosome (CHROM), position (POS), unique identifier (ID), reference and alternative alleles (REF and ALT), phred-scaled quality score (QUAL), site-filtering information (FILTER), a semicolon-separated list of additional annotation metadata (INFO) and the FORMAT field containing genotype data information (often a semicolon-separated list). The data frame can be manipulated to extract information with standard R functions (Box 1). The function VCFheader prints the header of a VCF file which, usually, contains details on these different fields.

The genotypes are then read by read.vcf. At this stage, it is good to think if it is interesting or feasible to read a complete VCF file. pegas takes roughly four bytes of RAM to store a single genotype; so (using the above notation), about 4np bytes are needed for a given genetic data set. With the rule of thumb that data in R should not use more than 25% of the available RAM to perform "comfortable" analyses, one million loci for 1000 individuals appears as a reasonable upper limit on a computer with 16 Gb of RAM. Furthermore, the user may be interested in only some loci depending on the type of genetic variant, position on the chromosome, or else. A strategy is to locate these loci by finding the rows of the output of VCFloci using standard R data manipulation (Box 1). Similar operations can be done on the field "CHROM" (if loci from several chromosomes are in the same file), "QUAL", or "FILTER". The field "INFO" requires a special attention because it includes detailed information on each locus. To make this information easier to digest, pegas provides the function getINFO to extract a specific information from this character string; by default this is the sequencing depth (DP), but the option what makes possible to extract another information, for instance, the variant type (Box 2).

For small VCF files, read.vcf can be used directly without calling VCFloci first, but the latter is still needed to get the information on the loci.

Analysis of population genomic data with R

Haplotypes and Linkage Disequilibrium

The possibility to investigate the genetics of populations at many loci of individuals brings the perspective to study the dynamics of genomes in populations. An exciting opportunity is to identify parts of species genomes which evolve as a single unit [START_REF] Andrews | Recent novel approaches for population genomics data analysis[END_REF]). Several approaches with NGS allow one to infer genotype phasing for diploid or polyploid organisms, either directly by sequencing long fragments of DNA (e.g., [START_REF] Mccoy | Illumina TruSeq synthetic long-reads empower de novo assembly and resolve complex, highly-repetitive transposable elements[END_REF], or with high coverage sequencing of shorter fragments [START_REF] Nielsen | Genotype and SNP calling from next-generation sequencing data[END_REF]. With phased genotypes, it is straightforward to reconstruct the haplotypes of each individual. pegas has the generic function haplotype which does this operation (it also works with DNA sequences). Linkage disequilibrium (LD) is the most widely used statistical method to identify alleles at different loci that are statistically associated. It is also of central importance in association studies of genetic diseases [START_REF] Goldstein | Population genomics: linkage disequilibrium holds the key[END_REF][START_REF] Jorde | Population genomics: a bridge from evolutionary history to genetic medicine[END_REF]. The function LDscan permits to analyse LD by calculating pairwise correlation coefficients (r 2 ) for a series of biallelic loci. The output of this function can be plotted with LDmap; this function has an option POS to plot the correlation coefficients together with the positions of the loci on the chromosome.

Figure 1 shows the LD at two scales for the chromosome 2L of the fruit fly data: for the first 100 loci and for 100 loci equally spaced on the chromosome (Box 3). Two other functions analyse LD for a pair of loci with any number of alleles: LD and LD2 for phased and unphased genotypes, respectively. These two functions return detailed statistics and tests [START_REF] Schaid | Linkage disequilibrium testing when linkage phase is unknown[END_REF][START_REF] Zaykin | Correlation-based inference for linkage[END_REF].

Multivariate Methods

Multivariate methods are common when analysing big data sets because of their abilities to provide dimension reduction, that is to provide a summary of many variables into a few, usually linear, combinations of them. Thus, population genetics has used these methods for some time [START_REF] Westfall | Allozyme markers in breeding zone designation[END_REF][START_REF] Guinand | Use of a multivariate model using allele frequency distributions to analyse patterns of genetic differentiation among populations[END_REF][START_REF] Moazami-Goudarzi | Is a multivariate consensus representation of genetic relationships among populations always meaningful[END_REF]. The package adegenet implements a lot of tools for the multivariate analysis of genetic data [START_REF] Jombart | Genetic markers in the playground of multivariate analysis[END_REF]. Some of the methods of particular interest here include the spatial principal components analysis (sPCA; [START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF]) and the discriminant analysis of principal components (DAPC; [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF] which is illustrated below.

In spite of their efficiency with large tables, multivariate methods are penalized when analysing very large data sets, for instance, when n and p are larger than 10,000. Some approaches have been developed to decompose very large matrices without using traditional methods such as eigen or singular value decomposition [START_REF] Halko | Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions[END_REF].

These approaches are based on random matrices and aim to approximate the first few principal components. [START_REF] Abraham | Fast principal component analysis of large-scale genome-wide data[END_REF] implemented this approach to perform a PCA on genomic microarray data collected on a large sample of individuals (n > 10,000). Their implementation is available as the R package flashpca distributed on GitHub (https://github.com/gabraham/flashpca). This approach is quite straightforward to code in R and is available in several forms with more or less accuracy [START_REF] Halko | Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions[END_REF].

Discriminant analysis (DA) is a multivariate method which searches for linear combinations of variables resulting in the strongest discrimination of groups identified from the data [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF]. Under some assumptions, it is possible to quantify the reliability of the inferred assignments. Thus, DA offers a powerful approach to assess the presence of structure in genetic data with many loci. However, because it is based on matrix decomposition algorithms, it suffers from the same limitations than PCA with very big data sets, and these methods need to be adapted to be used on genomic data.

The DAPC gives a solution to this problem: it can be used to assess group structure in genetic diversity with large datasets, and is also implemented separately for "genlight" objects [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF]. The main idea behind the DAPC is to first perform a data reduction using a PCA. Typically, less than 100 principal axes are sufficient. If group priors are missing, a sequential k-means clustering with model selection based on the Bayesian information criterion (BIC) is done to infer prior clusters. The next step is to perform a discriminant analysis in the reduced data space using the prior clusters returned by the k-means clustering (or the cluster definition provided by the user if available). The results can be graphically visualized in two ways: by plotting the projections of the individuals on discriminant axes (as usual in multivariate analyses), or by examining the relative posterior probabilities to the different clusters inferred from the discriminant analysis (as usual in stochastic assignment methods).

We consider here the dolphin data and try to show evidence for geographical structure (Fernández et al. 2016a). The data are provided as a VCF file (simply named 'vcf') available from Dryad [START_REF] Fernández | Data from: A genomewide catalogue of single nucleotide polymorphisms in white-beaked and Atlantic white-sided dolphins[END_REF]. Because the multivariate methods are implemented in adegenet, we have to convert the data read by pegas into the appropriate class (Box 4). Remarkably, none of these analyses take significant running time (in fact, the longest operation is to convert from the class "loci" to the class "genind").

Population Genetics and Genomic Bioinformatics

The main strength of R, and maybe the most attractive one for data analysis, is the possibility to integrate a great variety of methods and tools while using the same programming language. In the past few years, the analysis of genetic data has changed quite radically. Until the early twenty-first century, population geneticists used to get their data from the laboratory on their physical support (e.g., electrophoresis gels), then input the data in the computer, and analyse them with specialised software. With the advent of NGS, genetic data are directly acquired on the computer, but an additional step is required: the identification of alleles from the many (hundreds of millions) sequencing reads. This step is called 'variant calling' and requires highly specialised software. Because NGS technologies are evolving very fast (see above), this software is also evolving very fast. Fortunately, the past few years have seen the development of more integrated bioinformatics tools to accomplish such tasks. The R package poRe has been recently released to analyse data from the the mobile nanopore sequencer MinION [START_REF] Watson | poRe: an R package for the visualization and analysis of nanopore sequencing data[END_REF]. This package takes profit of R's tools such as scanning a directory with many files, searching for patterns in file names, and reading files in a specific format (here HDF5), to extract and plot summary statistics from the sequencing runs, as well as export data into FASTQ or FASTA format.

Bioconductor is a long-term effort to produce R packages for the analysis of large-scale genomic data sets (Gentleman et al. 2005;[START_REF] Huber | Orchestrating high-throughput genomic analysis with Bioconductor[END_REF]). Bioconductor's website (www.bioconductor.org) hosts a suite of integrated R packages mostly dedicated to the analysis of expression data. Several packages have been released for the analysis of DNA microarray data (e.g., [START_REF] Dunning | beadarray: R classes and methods for Illumina bead-based data[END_REF][START_REF] Ritchie | R/Bioconductor software for Illumina's Infinium whole-genome genotyping BeadChips[END_REF][START_REF] Morgan | argyle: a R package for analysis of Illumina genotyping arrays[END_REF]. The package ape has functions to convert DNA data from Bioconductor classes into the "DNAbin" class from which population and evolutionary analyses can be performed.

A crucial aspect of R packages which permits the integration of different tasks is the stability of the packages themselves and especially of the data classes defined therein (Tables 1,3). Several R packages for evolutionary genetic analyses are now more than one decade-old and have been widely used and tested, so that their contents are now reliable for a wide range of applications, and package developers have the possibility to work with them to develop new tools and methods (e.g., [START_REF] Paradis | Analysis of Phylogenetics and Evolution with R (Second Edition)[END_REF]Lawrence & Morgan 2014).

Perspectives

Remarkable progress has been accomplished in the development of software for population genetics in the context of analysing data sets from NGS. There is no doubt that this will continue in the years to come where population genetic data will become more and more important to assess biodiversity dynamics and adaptations of species in a changing environment. We discuss below three specific areas where significant progress is likely to be achieved in the coming years.

As apparent from the contributions to this Special Issue, a lot of effort has been devoted to improve the overall quality of the software dedicated to population genetics.

On the other hand, hardware has continued to improve. For instance, multi-core processors are now the rule rather than the exception, so that parallel computing can be accomplished on almost all computers. The adaptation of genetic software to this has been uneven. Lawrence & Morgan (2014) give some examples using Bioconductor, illustrating the issue that implementing a parallel algorithm to accomplish a given task may not always be beneficial. R has a built-in capacity to parallelize computations using the package parallel. The functions in parallel can be used with most R packages described in this article, for instance, by distributing a calculation to several cores of a processor.

However, there is certainly some improvements possible since a lot of the computations used in population genetics implies repetitive and independent calculations (e.g., for allele frequencies by population or for LD).

NGS raw data are more or less noisy because they are generated through some biological and electronical processes which may be affected by random variation. Besides these measurement errors, intrinsic biological processes induce another form of variation, for instance because of cell-specific genetic changes. [START_REF] Lynch | Estimation of nucleotide diversity, disequilibrium coefficients, and mutation rates from high-coverage genome-sequencing projects[END_REF] addressed the issue of inferring population genetic parameters using NGS data (see [START_REF] Korneliussen | Calculation of Tajima's D and other neutrality test statistics from low depth next-generation sequencing data[END_REF], for a more recent contribution). [START_REF] Blischak | Accounting for genotype uncertainty in the estimation of allele frequencies in autopolyploids[END_REF] developed a hierarchical Bayesian model to assess genotype uncertainty in autopolyploids using raw read counts. They implemented their method in the R package polyfreqs (available on CRAN). We take this opportunity to emphasize the importance of open source software for the future of population genomics. O' Rawe et al. (2015) warned about the risk that commercial NGS applications tend to ignore uncertainty in genotyping inference. As a sign of this risk, these authors mention that several applications for haplotype or variant call are proprietary and closed-source. Considering the continuously increasing diversity of NGS technologies and the increasing range of their applications, a special attention to uncertainty in statistical inference using NGS data will be crucial.

Environmental monitoring already allows acquisition of data in real time, and recent developments in NGS (e.g., nanopore sequencing) strongly suggest that similar things will soon be accomplished with population genetic data. This is a crucial challenge for developing efficient software and making sense of the results in real time. Facing this forthcoming challenge will mean having to combine data analysis with theoretical population genetics in order to help us meet issues related to global change, biodiversity conservation, and pathogen dynamics. It will surely be useful to remember the words of caution-if not wisdom-from [START_REF] Gower | The biological stimulus to multidimensional data analysis[END_REF]: "Automatic instrumentation ensures that there is no lack of large data sets. This is something new but I sometimes think that we are getting perilously close to the search for the philosophers stone. If only we could find the right recipe, great truths would be revealed by analysing vast masses of data." disequilibrium with multiple alleles. Genetics, 180, 533-545.

546

Box 1 Scanning and reading VCF files First we scan the human Y chromosome data with VCFloci, print the number of loci and the labels of the columns with names: > info.Y <-VCFloci ("chrY.vcf The number of individuals in the VCF file are obtained with the function VCFlabels which extracts their labels (or identifiers):

> labs <-VCFlabels("chrY.vcf.gz") > length(labs)

[1] 1233 > head(labs)

[1] "HG00096" "HG00101" "HG00103" "HG00105" "HG00107" [6] "HG00108"

Then we read the genotypes with read.vcf. Here, we select the loci in a 1 Mb region of the Y chromosome between positions 5,000,000 and 6,000,000 which can be found with:

> sel <-which(info.Y$POS > 5e6 & info.Y$POS <= 6e6) > length(sel) [1] 5
read.vcf has two ways to specify which loci to read: either using the options from and to (which are set to 1 and 10,000 by default), or using the option which.loci which requires a vector of integers specifying the indices of the loci to read:

> x <-

Box 3 Analysis of linkage disequilibrium

We read the fruitfly data and call the two functions LDscan and LDmap on the first 100 loci (we could use the option to as above to read only 100 loci, but this gives us the opportunity to show how to subset a data frame of genotypes). The positions of the loci are stored in the column POS output by VCFloci (Box 2), so we extract the positions for these 100 loci and pass them as the (optional) second argument of LDmap (Fig. 1a):

> x <-read.vcf(fl) > res <-LDscan(x[, 1:100]) > LDmap(res, info.fly$POS[1:100])

We now select 100 loci regularly spaced from 1 to 224,253 (this value comes from the last tabulation in Box 2) and using ceiling to ensure we have only integers:

> s <-ceiling(seq(1, 224253, length.out = 100))

We then repeat the LD analysis using the appropriate POS information (Fig. 1b):

> xs <-read.vcf(fl, which.loci = s) > res2 <-LDscan(xs) > LDmap(res2, info.fly$POS[s]) 549
Box 4 DAPC with the dolphin data

We read 1000 loci and convert them to the class "genind": > x2 <-read.vcf("vcf", to = 1000) > g2 <-loci2genind(x2)

We can now call the method implemented in the function find.clusters: > o2 <-find.clusters(g2)

When called with no other option, the function displays two graphs (Fig. 2) and asks the user to input the number of principal components (PCs) and the number of clusters to retain. The first graph shows that 40 PCs retain about 80% of the variance, so we select this number. The second graph shows a lowest value of BIC with seven clusters. However, the present data has eight populations (see below) so we select this number of clusters (K) for this analysis (we could use simply the original populations).

The output is a list with several elements, including grp giving the cluster assignment of the individuals. The labels of the individuals include the geographical origin of the samples (as often the case with this kind of data set). So it is simple to extract this information, and cross-tabulate it with the cluster assignment done by find.clusters:

> head(names(o2$grp))

[1] "FAE_2" "FAE_609" "FAE_660" "FAE_667" "FAE_668" "FAE_671" > pop <-gsub("_.*", "", names(o2$grp)) > cluster <-o2$grp > table(pop, cluster) cluster pop 1 2 3 4 5 6 7 8 DEN 0 0 0 0 4 2 0 0 FAE 0 9 1 3 2 0 2 1 FRA 0 0 0 0 1 0 1 0 GER 3 0 0 0 5 3 1 0 ICE 0 0 0 0 15 2 0 0 IRE 0 9 0 0 2 2 3 0 NOR 0 0 0 0 5 2 0 0 SCO 0 5 0 0 24 1 2 0 There does not seem to be a very good match between the originally sampled populations and the cluster assignments. Thus, we repeat the analysis with a smaller number of clusters (K = 4) in a non-interactive way by specifying the required options: > o2b <-find.clusters(g2, n.pca = 40, n.clust = 4)

As above, we extract the 'pop' and 'cluster' information (not shown):

> table(pop, cluster) cluster pop 1 2 3 4 DEN 0 6 0 0 FAE 0 6 4 8 FRA 0 1 0 1 GER 6 5 0 1 ICE 2 15 0 0 IRE 2 4 0 10 NOR 2 5 0 0 SCO 0 26 0 6

We now perform the DAPC with 4 groups asking to output two discriminant axes: > d2b <-dapc(g2, cluster, n.pca = 40, n.d = 2)

The output contains several elements, including the coordinates of the individuals and the group means. The plots can now be done with R's standard functions (see Supplementary Information) or with adegenet's functions (Figs. 3,4): > col <-c("red", "blue", "gold", "darkgreen") > scatter(d2b, posi.da = "topright", col = col) > compoplot(d2b, legend = FALSE, col = col) [START_REF] Goudet | Hierfstat, a package for R to compute and test variance components and F -statistics[END_REF], poppr [START_REF] Kamvar | Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction[END_REF], and mmod [START_REF] Winter | mmod: an R library for the calculation of population differentiation statistics[END_REF].

Package Main methods implemented The linkage coefficients (r 2 ) between each pair of loci are indicated as coloured squares: the squares at the bottom of the triangle are for nearby loci, whereas the square at the top is for the two most distant loci.

adegenet*

  Multivariate methods (sPCA, DAPC) ape Evolutionary distances, distance-based phylogenetics, phylogenetic bootstrap hierfstat* Hierarchical F -statistics pegas Population differentiation (HWE,, AMOVA), haplotype networks, linkage maps poppr* Applications to clonal organisms, distances for microsatellites mmod* Population differentiation (F ST , G ST , Φ ST , D) stacksr Interface to STACKS vcfR Tools for FASTA, GFF, and VCF files *These packages use the class "genind".

Fig. 1

 1 Fig. 1 LD maps at two different genomic scales for the fruit fly data: (a) for the 100 first loci on the chromosome 2L, and (b) for 100 loci regularly spaced along this chromosome. The horizontal axis indicates the position of the loci on the chromosome.

Fig. 2

 2 Fig. 2 Results of finding clusters with the dolphin data. (a) Cumulated variance explained by the principal components (PCs). (b) Bayesian information criterion (BIC) for the different values of K (number of clusters).

Fig. 3

 3 Fig.3Results of the DAPC with the dolphin data.

Fig. 4

 4 Fig. 4 Posterior probabilities of group membership from the DAPC with the dolphin data. The color scheme is the same than in Fig 3.

  Figure 1:

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure 2:

  Box 2 Going further with the INFO field in VCF filesHere we use the information returned by VCFloci, such as the field CHROM (if loci from several chromosomes are in the same file), to select some loci to be read. The field INFO requires a special attention because it includes detailed information on each locus:To make this information easier to digest, the function getINFO helps to extract a specific information from this character string; by default this is the sequencing depth (DP), but the option what makes possible to extract any other information, for instance, the variant type:To illustrate how to use other information returned by VCFloci, we scan the fruitfly data[START_REF] Kao | Data from: Population genomic analysis uncovers African and European admixture in Drosophila melanogaster populations from the south-eastern United States and Caribbean Islands[END_REF]) which include genotypes from all chromosomes in a single file:

	> info.Y$INFO[1]
	[1] "AA=G;AC=22;AF=0.0178427;AN=1233;DP=84761;NS=1233;\
	AMR_AF=0.0000;AFR_AF=0.0000;EUR_AF=0.0000;SAS_AF=0.0000;\
	EAS_AF=0.0451;VT=SNP;EX_TARGET"
	> VT.Y <-getINFO(info.Y, what = "VT")
	> table(VT.Y)
	INDEL	MNP	SNP	SV
	1314	113 60505	110
	The function VCFheader reads the header of the VCF file where details on the abbreviations used
	can be found:
	> cat(VCFheader("chrY.vcf.gz"))
	....
	##INFO=<ID=VT,Number=.,Type=String,Description="indicates what type\
	of variant the line represents"
	....
	Thus, we can now locate which loci are SNPs:
	> sel.snp <-which(VT.Y == "SNP")
	> length(sel.snp)
	[1] 60505
	As usual with R, it is possible to combine these logical comparisons with the AND operator (&),
	for instance, to locate the SNPs within a 2 Mb region:
	> sel <-which(VT.Y == "SNP" & info.Y$POS > 1e6 & info.Y$POS <= 3e6)
	read.vcf("chrY.vcf.gz", which.loci = sel) > length(sel)
	Reading 5 / 5 loci. [1] 1231
	Done.
	> x
	Allelic data frame: 1233 individuals
	5 loci > fl <-"global.pop.GATK.SNP.hard.filters.V3.phased_all.pop.maf.05.recode.vcf.gz"
	> info.fly <-VCFloci(fl) If the user wants to read all loci with a single command, the following can be done: > table(info.fly$CHROM)
	> x <-read.vcf("chrY.vcf.gz", to = nrow(info.Y)) 2L 2R 3L 3R R X
	224253 193675 214235 270619	1 153035

Table 1 :

 1 Major R data classes for population genetics and genomics. The package in parentheses is where the class is defined and which contains the main functions for data manipulation and visualization. The question mark (?) is not part of the IUPAC ambiguity code.†Standard missing value in R.

	Class Data type Main features Missing data	DNAbin (ape) DNA sequences Single or multiples sequences, aligned or not Ambiguity code + ?*	NA †; null allele(s) can genind (adegenet) Allelic (allele-centered) Any level of ploidy; ploidy may vary	be specified at input among loci and among individuals	Not supported genlight (adegenet) biallelic SNP Any level of ploidy	NA †; null allele(s) can loci (pegas) Allelic (locus-centered) Any level of ploidy; ploidy may vary	be specified from R among loci and among individuals;	genotypes may be phased or not

*

Table 2 :

 2 Data inter-operability.

	Class		File	Class conversion
		Input	Output	From	To
	DNAbin FASTA	FASTA	Bioconductor* genlight
		Phylip	Phylip
		Clustal
	genind	Genetix		loci	loci
		Structure
		Fstat
		Genepop
		GenAlex
	genlight PLINK		DNAbin
	loci	Tabular files † Tabular files †	genind	genind
		Genetix
		VCF
	vcfR	VCF	VCF		DNAbin
					genind
					loci

*Several classes are supported.

†Tab-delimined, CSV, and similar formats.

Table 3 :

 3 Main R packages for population genetics and genomics. The packages not cited in the text are hierfstat

Acknowledgements

We are grateful to David Winter and two anonymous reviewers for their constructive comments on a previous version of this paper. The resource reported in this paper started at the Population Genetics in R Hackathon, which was held in March 2015 at the National Evolutionary Synthesis Center (NESCent) in Durham, NC, with the goal of addressing interoperability, scalability, and workflow building challenges for the population genetics package ecosystem in R. The authors were participants in the hackathon, and are indebted to NESCent (NSF #EF-0905606) for hosting and supporting the event. Many thanks to all participants of the hackathon for their feedback and enthusiasm during and after the event. TJ is funded by the Medical Research Council Centre for Outbreak Analysis and the National Institute for Health Research -Health Protection Research Unit (NIHR HPRU) in Modelling Methodology at Imperial College London in partnership with Public Health England (PHE). KS is funded by a grant from the National Science Foundation (DEB 1350474). This is publication ISEM 2016-150.

Membership probability 0.0 0.2 0.4 0.6 0.8 1.0

Supplementary Information

R code for custom plot of the results from the DAPC on the dolphin data. This example shows how to extract the coordinates of the individuals and plot them with a symbol representing the original populations (Fig. S1):

plot(d2b$ind.coord, type = "n") points(d2b$grp.coord, pch = 19, col = "grey50", cex = 4) text(d2b$grp.coord, labels = 1:4, col = "white") points(d2b$ind.coord, pch = (1:8)[pop]) legend("top", legend = levels(pop), pch = 1:8, ncol = 4)


Vincent Angladon
email: vincent.angladon@irit.fr

Simone Gasparini
email: simone.gasparini@irit.fr

Vincent Charvillat
email: vincent.charvillat@irit.fr

Room Floor Plan Generation on a Project Tango Device

Keywords: Floor plan, Mobile device, RGB-D cameras, 2D/3D mapping

This article presents a method to ease the generation of room floor plans with a Project Tango device. Our method takes as input range images as well as camera poses. It is based on the extraction of vertical planar surfaces we label as wall or clutter. While the user is scanning the scene, we estimate the room layout from the labeled planar surfaces by solving a shortest path problem. The user can intervene in this process, and affect the estimated layout by changing the label of the extracted planar surfaces. We compare our approach with other mobile applications and demonstrate its validity.

Introduction

Floor plan generation is the problem of generating a drawing or a digital model to scale of an existing room or building. A common workflow consists in taking individual measurements reported on a freehand sketch of the floor, then in using a CAD software to draw the floor plan. This approach is mostly manual and requires the user to collect all the measurements, typically with a laser range finder. The task is not trivial, as it requires the correct use of the measurement tool and the collection of all the necessary measurements to fully define the floor plan: missing measurements may lead to incomplete plans, thus requiring costly do-overs on site. Moreover, the task can be challenging when dealing with furnished or cluttered environments, preventing, e.g., the direct measurement of certain distances. In order to provide more automatic and efficient solutions, the use of desktop or mobile applications performing the drawing on-site has gained a lot of interest over the last years.

With the latest advancement in 3D scanning technologies, 3D scanners placed on tripods have become an interesting solution for the generation of the floor plan. The generated 3D point cloud enables to perform Building Information Modeling (BIM) and floor plan generation in a semi-automatic way, through the use of a dedicated software. The operating cost is generally higher than the previous method, due to the investment of the scanner and the scanning time which can be quite long when there are several small rooms. The automatic creation of floor plans or BIM from such scanners has been extensively researched for single rooms [START_REF] Budroni | Automated 3d reconstruction of interiors from point clouds[END_REF][START_REF] Xiong | Automatic creation of semantically rich 3d building models from laser scanner data[END_REF] and multiple rooms [START_REF] Xiao | Reconstructing the worlds museums[END_REF][START_REF] Stambler | Building modeling through enclosure reasoning[END_REF][START_REF] Oesau | Indoor scene reconstruction using feature sensitive primitive extraction and graph-cut[END_REF][START_REF] Mura | Automatic room detection and reconstruction in cluttered indoor environments with complex room layouts[END_REF][START_REF] Mura | Piecewise-planar reconstruction of multi-room interiors with arbitrary wall arrangements[END_REF][START_REF] Ochmann | Automatic reconstruction of parametric building models from indoor point clouds[END_REF]. These methods, are called offline because they start the processing after the scan is complete and all the data is available.

At the opposite, for handheld scanners, online methods provide results incrementally during the scanning process. Several online approaches have been proposed [START_REF] Zhang | Walk&sketch: create floor plans with an rgb-d camera[END_REF][START_REF] Zhang | Online structure analysis for real-time indoor scene reconstruction[END_REF][START_REF] Dzitsiuk | De-noising, stabilizing and completing 3d reconstructions on-the-go using plane priors[END_REF] to extract some walls of indoor scenes, but none of them focus on the creation of single or multiple room floor plans.

On common smartphones, user-driven approaches have been proposed [START_REF] Rosser | Modelling of building interiors with mobile phone sensor data[END_REF][START_REF] Pintore | Interactive mapping of indoor building structures through mobile devices[END_REF][START_REF] Sensopia | MagicPlan, Create a floor plan in just a few minutes[END_REF]. Knowing the phone orientation and the vertical distance of the device to the floor (assumed constant and calibrated), they estimate the distances between the user and the room corners. More recently, as mobile devices started to sport depth sensors, userdriven approaches have been proposed for the Project Tango mobile devices: the user is required to manually select the walls with FloorPlanEx [START_REF]Google: java_floor_plan_example create a floor plan by using the depth sensor of a Google Tango device to detect and measure walls in a room[END_REF] or the edge corners of the indoor scene [START_REF]Google: Measure: augmented reality measurement application for Google Tango devices[END_REF].

In this paper, we propose an online approach for the Tango Tablet Development Kit(TDK) Fig. 1a, which is a tablet equipped with a depth sensor and a fish eye camera running on Android. It offers dedicated libraries to perform localization and 3D reconstructions in the form of 3D textured meshes. Our method can generate a room floor plan with optional user interaction. We rely on the ability of this device to localize itself and capture depth maps, in order to have a better understanding of the scene and hence reduce the user efforts. Given the complexity of the task and the diversity of the scenes, the user interaction is yet helpful to improve the robustness of the approach: large occlusions, lack of physical separation between rooms, missing data, incorrect depth perception (e.g. glasses and mirrors), and other issues that could be challenging for a fully automatic system, can be easily solved by the operator. [START_REF] Zlatanova | Problems in indoor mapping and modelling. Acquisition and Modelling of Indoor and Enclosed Environments[END_REF] enumerates and classifies all the issues related to the problem of indoor mapping and modeling. While the Tango TDK addresses most of the acquisitions and sensors problems (variable lighting conditions, sensor fusion, mobility support), our works also tackles acquisition problems (variable occupancy support), data structure and modeling issues (real-time modeling, dynamic abstraction), visualization problems (on mobile visualization, real-time change visualization) and legal issues (user privacy).

Overview of the Method

Hypotheses. We assume the considered rooms are made of a horizontal ceiling and floor, and vertical planar walls, the latters not necessarily orthogonal w.r.t. each other. They can contain clutter (furniture, movable objects, . . . ) occluding the walls. With the Tango TDK, we observed objects located 4.2 m away suffer from a depth uncertainty over 25 ± 9 mm, and lack of texture (e.g. walls and ceiling of uniform color) can lead to incorrect camera poses from the Tango motion tracking. Therefore, our approach is limited to rooms with reasonable texture and where all the walls and the ceiling can be observed by the depth sensor without too many efforts (medium size room with a 5 m ceiling height maximum).

Pipeline. Figure 1b summarizes the proposed pipeline. The Tango Tablet Development Kit middleware provides at each iteration the depth map of the scene, as well as the camera pose. The depth map is processed in order to extract sets of 3D points corresponding to candidate walls, what we call planar patches (detailed in Sect. 3). Thanks to the known camera pose, the planar patches are brought into a global world coordinate system, so that they can be associated and then fused with the existing patches of the model (Sect. 4). Whenever the model is updated, the visibility polygon of the discovered area(s) is also updated (Sect. 5), which can be used as a visual feedback for the user. Furthermore, the visibility polygon is used by the labeling module to automatically classify the planar patches as wall or clutter (Sect. 6), thus allowing the disambiguation between actual walls and other objects that may be lying inside the room. This classifier was trained beforehand against 900 manually labeled vertical planar regions from our training dataset. This last task can take advantage of the interaction of the user, who can correct and change the automatic labeling of the vertical planar patches into wall or clutter (Sect. 7). Finally, the layout generation module computes the room layout from the labeled vertical planes and the boundary given by the visibility polygon (Sect. 8).

Planar Patches Extraction

A planar patch is defined as a list of 3D points associated with the equation of the fitted infinite plane and the boundaries of this patch. We define the planar patches extraction problem as a function which takes as input a range image (also called depth map) and returns a list of planar patches.

The common strategies to extract planar patches are the Hough Transform approaches [START_REF] Borrmann | The 3d hough transform for plane detection in point clouds: A review and a new accumulator design[END_REF], the region growing algorithms [START_REF] Poppinga | Fast plane detection and polygonalization in noisy 3d range images[END_REF], the normal map segmentation [START_REF] Holz | Robot soccer world cup[END_REF] and the split and merge approaches [START_REF] Feng | Fast plane extraction in organized point clouds using agglomerative hierarchical clustering[END_REF].

In our preliminary tests, we noticed that normal map segmentation approaches were too sensitive to the data noise and required higher computational time. We chose instead a region growing approach, similar to the algorithm described by Poppinga et al. [START_REF] Poppinga | Fast plane detection and polygonalization in noisy 3d range images[END_REF]. In this approach, a random seed point and its neighbors are picked and extended by taking into consideration neighboring points. A plane is estimated on this set of points and a new point is considered valid when its distance to the plane is small enough. The planar patch keeps growing iteratively until no valid point can be found in the neighborhood of the patch. A new seed point is picked until all points have been considered. We made some modifications to the original algorithm in order to get better results with the point cloud provided by the Tango TDK. We ignore the 3D points with a distance to the camera farther to 3.7 m as they revealed to be very noisy. To cope with the low density of the depth map, the neighboring selection is not pixel-wise but performed on a 3 × 3 pixel window. In a post-processing step, we remove isolated points of the planar patches with a neighborhood pixel analysis performed in the segmented depth map. This ensure our planar patches are more compact. Finally, we perform a validation of the planar patch with the RANSAC algorithm. Without implementing the optimizations suggested by [START_REF] Poppinga | Fast plane detection and polygonalization in noisy 3d range images[END_REF], we can process 10k points in 47 ms with our C++ implementation.

The ceiling and the floor are incrementally estimated with the minimal and maximal height of the considered vertical planar patches and validated with the fitting of horizontal planes at the considered heights. In the following sections, we consider only the vertical planar patches for which we compute a rectangle boundary.

Planar Patches Model Update

We create a model of planar patches in the world coordinates which is initialized with the patches of the first frame. For each successive frame, we compute planar patches which are associated with the model with the help of an association function. We update the associated planar patches of the model with a fusion function.

Planar Patch Association. Given a list of planar patches from the model and a list of planar patches from the current frame the association function computes a list of planar patches pairs where each patch appears only once. We consider a distance between planar patches which takes into account the normal angle difference of the planes and the mean distance of the points-to-plane distances, using the points of the rectangle boundary of the planar patches. This distance defines the candidate pairs of associated patches, which are validated with an overlapping test between the rectangle boundaries.

Planar Patch Fusion After two planar patches have been associated, we need to create a new planar patch combining the two. We compute and update the covariances matrices of the planes using [START_REF] Lee | Indoor mapping using planes extracted from noisy rgb-d sensors[END_REF]. The plane equations of the fused patch is updated with a Principal Component Analysis, in order to avoid the costly storage of the 3D points associated with the planar patches. We then update the rectangular boundary to include the two planar patches.

Visibility Polygon Computation

During a scan, the user needs to know which parts of the scene he visited or not. We provide this information with a top-down 2D view representation of the area observed by the camera in the form of a visibility polygon, similarly to [START_REF] Zhang | Walk&sketch: create floor plans with an rgb-d camera[END_REF], but computed differently from the vertical planar patches model and the history of the camera positions. We denote s v the line segment obtained from the projection of a vertical planar patch on a horizontal plane. We consider the visibility view polygon P s associated to each line segment s v , which is formed by the union of all the triangles formed by the camera position and the two extremities of s v , see Fig. 2. Each triangle is made of one wall segment and two frustum segments. The geometric union of all the polygons P s forms the visibility polygon. We implemented the computation of the visibility polygon with the help of the geometry engine GEOS [START_REF]OSGeo: GEOS, Geometry Engine[END_REF].

Wall-Clutter Separation

While clutter often consists of irregular shapes, such as plants, sofas, etc., which are eliminated in planar primitive approaches, it can also consist of piecewise planar shapes such as cupboards, radiators, etc. In this section, we propose a method to classify the planar primitive as wall or clutter in order to reduce the number of potential irrelevant planes considered by the room layout estimations component. The segmentation of the primitives can be performed individually [START_REF] Mura | Automatic room detection and reconstruction in cluttered indoor environments with complex room layouts[END_REF][START_REF] Murali | Indoor scan2bim: Building information models of house interiors[END_REF][START_REF] Stambler | Building modeling through enclosure reasoning[END_REF], i.e. using features on each primitive considered independently, or globally [START_REF] Mura | Piecewise-planar reconstruction of multi-room interiors with arbitrary wall arrangements[END_REF][START_REF] Xiong | Automatic creation of semantically rich 3d building models from laser scanner data[END_REF], i.e. considering the adjacent primitives, which provide contextual information.

Our objective is to build the room layout incrementally, during the scan progress, which means all the adjacencies are only known at the end of the scan. Therefore we favored a classification with an individual approach. Our features take the form of a vector (d c , d f , l, d v ). They include the distances d c between the highest point of the vertical planar patch and the estimated ceiling and d f between the lowest point of the vertical planar patch and the estimated floor. Intuitively, a planar patch both close to the ceiling and floor is likely to be a wall. Similarly, a segment s v with a longer length l has a higher probability to correspond to a wall. We also consider the distance d v = max p∈s v d(p, ∂ P v ) between a segment s v and the exterior boundary of the visibility polygon ∂ P v . A high distance d v corresponds to a segment with at least one extremity far from the visibility polygon, which is likely to correspond to clutter.

We compute a wall probability P (s v ) for each segment s v of the model with a Multi-layer Perceptron classifier using one hidden layer and a logistic sigmoid activation, trained on our dataset. When a new frame is ingested, we compute P (s v ) if s v was modified since the previous frame or the ceiling/floor estimation changed. Our implementation based on the Python Scikit-learn module can label 50 planar patches in 29 ms on average.

User Interaction

Automatic approaches for floor plan generation or scan to BIM can be prone to errors. These errors can come from missing data, clutter, sensor noise or some specificities in the scene (e.g. small walls, large openings, cavities in the walls, obstacles fully covering a wall, etc.). Depending on the progress of the scan, it may not be possible for an algorithm to determine whether a planar patch corresponds to a wall or clutter.

The user interaction schemes proposed in the literature are usually corrective [START_REF] Arikan | O-snap: Optimizationbased snapping for modeling architecture[END_REF][START_REF] Mura | Piecewise-planar reconstruction of multi-room interiors with arbitrary wall arrangements[END_REF][START_REF] Pintore | Mobile mapping and visualization of indoor structures to simplify scene understanding and location awareness[END_REF], where the user intervenes at the end to correct a proposed solution (commonly considered for offline approaches), or user-driven [START_REF] Sensopia | MagicPlan, Create a floor plan in just a few minutes[END_REF][START_REF] Rosser | Modelling of building interiors with mobile phone sensor data[END_REF][START_REF] Pintore | Interactive mapping of indoor building structures through mobile devices[END_REF][START_REF]Google: Measure: augmented reality measurement application for Google Tango devices[END_REF][START_REF]Google: java_floor_plan_example create a floor plan by using the depth sensor of a Google Tango device to detect and measure walls in a room[END_REF] where no solution can be computed without user interaction (mostly found in online approaches).

We wanted to propose a collaborative interaction scheme where the interaction would be optional and could be performed at any time. The first way to interact is to move the Tango TDK to make it observe new parts of the scene. A first view displays the camera image augmented with the detected planar patches colored relatively to their wall probability and the estimated layout, while a second view provides a top view of the scene with the visibility polygon and the estimated layout too. The two views shown in Fig. 3 are updated at the frame rate of the depth sensor, giving an immediate feedback to the user who can decide to visit the area(s) with missing data. The user can interchange the wall/clutter labels of the detected patches by touching them in the augmented view, which directly affects the room layout estimation component. At the end of the scan, we let the user perform further editions: suppress forgotten clutter patches and merge planar patches. At the time of writing, we only designed a desktop prototype which can replay a recorded scan or process live data transmitted by the application. The presented results were obtained by performing the interactions on replays of the scans.

The Room Layout Estimation takes as input the model of line segments s v coarsely classified as wall or clutter and generates a simplified floor plan (see Fig. 1b). In order to compute the room layout, we need to retrieve the topological relationships between the walls, i.e. their adjacencies.

Related Work.

A simple approach is to consider the topology of the room can be recovered by connecting adjacent segments [START_REF] Valero | Automatic method for building indoor boundary models from dense point clouds collected by laser scanners[END_REF][START_REF] Xiong | Automatic creation of semantically rich 3d building models from laser scanner data[END_REF][START_REF] Murali | Indoor scan2bim: Building information models of house interiors[END_REF], which works well for rooms with low clutter, and when the walls are well separated from the clutter. The use of cell-complex is very popular for both single rooms [START_REF] Budroni | Automated 3d reconstruction of interiors from point clouds[END_REF][START_REF] Previtali | Towards automatic indoor reconstruction of cluttered building rooms from point clouds[END_REF] and multiple rooms [START_REF] Oesau | Indoor scene reconstruction using feature sensitive primitive extraction and graph-cut[END_REF][START_REF] Mura | Automatic room detection and reconstruction in cluttered indoor environments with complex room layouts[END_REF][START_REF] Mura | Piecewise-planar reconstruction of multi-room interiors with arbitrary wall arrangements[END_REF][START_REF] Ochmann | Automatic reconstruction of parametric building models from indoor point clouds[END_REF] layout estimation. The line segments are replaced by infinite lines which partition the 2D space into polygonal cells. A graph connecting the adjacent cells is defined: instead of considering topological information on the segments, adjacency relationship between the cells is taken into account. The inside/outside label of the cells are computed with a graph-cut algorithm. This approach is robust to missing data because the extension of the segments is automatically considered, but it does not simplifies the layout.

Weighted graph construction. In our approach, the visibility polygon provides the topology, and we try to compute a simple chain of segments s v which explains the visibility polygon. We create a graph from the adjacency relationships between the segments. For simplification purpose, the segments with a low wall probability (inferior to 0.5 in practice) are discarded. The remaining segments S are very likely to be close to the boundary of the visibility polygon, but they might not cover this boundary completely due to missing data. We complete the boundary with segments s c from the relative complement of S with respect to the visibility polygon. A node is created at each segment extremity and each segment intersection as illustrated in Fig. 4a. The nodes belonging to the same segments are linked by an edge. We assign to each edge e a weight proportional to 1 -P (s v ), where s v is the segment associated to e. The edges corresponding to frustum segments and completion segments s c are assigned a penalty weight. With this design, the shortest path, which minimizes the sum of the weights of its edges, gives more importance to segments with a high wall probability and favors solutions with a low number of planar patches.

Solving and discussion. We use the Dijkstra's algorithm to find a cycle in our graph which minimizes the sum of our weights. To avoid trivial solutions, we use a start point on the segment with a high P (s v ), an endpoint which is adjacent, and we remove the edge between them. As shown in Fig. 4a, the graph may contain several cycles which can lead to incorrect layouts. In order to avoid finding an incorrect cycle, we impose the segments along the visibility polygon to correspond to directed edges following a clockwise orientation in the visibility polygon. This method cannot handle the cycles non adjacent to the visibility polygon, which revealed to be non-existent in our dataset. In a post-processing step, we replace, whenever possible, the chains of frustum/completion segments by lines segments extending their adjacent segments, otherwise we simply join them. Figure 4b illustrates the solution computed from a simplified scene. We also apply the Ramer-Douglas-Peucker algorithm [START_REF] Ramer | An iterative procedure for the polygonal approximation of plane curves[END_REF] with a low threshold For a satisfying user experience, the layout of previously seen areas should not change when the user visits a new part of the scene. This behavior cannot be guaranteed with an offline approach. When the ceiling and the floor are detected, our method is suitable for incremental changes of the model: the wall probability P (s v ) of the previously observed segments s v does not change, which means the computed path restricted to the previously seen segments is the same and has the same cost. Fig. 5: Four steps of our room layout estimation : the extracted segments, the generated graph, the shortest path solution, and the layout obtained after post-processing. See Fig. 3 and Fig. 4 for the signification of the colors.

In this section, we compare -the geometry accuracy of-the room floor plans generated with our approach, Magic Plan [START_REF] Sensopia | MagicPlan, Create a floor plan in just a few minutes[END_REF] (run on iPad Air 1, prior to ARKit release) and FloorPlanEx from Google [START_REF]Google: java_floor_plan_example create a floor plan by using the depth sensor of a Google Tango device to detect and measure walls in a room[END_REF].

Evaluation protocol. We considered five indoor scenes Lab1 MW , Lab2, House1 MW , House2 MW and House3 MW , where MW denotes the scenes respecting the Manhattan World assumption. The ground truth room layouts of these scenes were created with a Bosh DLE 50 laser ranger finder. We evaluated the geometry accuracy of the obtained layouts with the ground truth and the reproducibility of the measurements by repeating the measurements five times. Each estimated layout was aligned with the ground truth by computing the transformation which minimizes the distances between their corresponding vertices. The mean of these distances defines our residual error. We also evaluated the user effort during the use of the considered mobile applications. Magic Plan and FloorPlanEx are user-driven applications where the user selects the walls and the corners, respectively. The number of interactions is equal to the number of corners (plus one for Magic Plan). We did not count the interactions required for Magic Plan calibration processes. For our approach, we evaluated the number of labels corrections on the planar patches and the number of post-scan modifications.

Results and analysis. Table 1 and Fig. 6 show the obtained results. Magic Plan estimates the camera-to-corner distances from the device orientation instead of taking advantage of a localization module or a depth sensor. Consequently, even the best results of the application are less accurate than the results obtained with the other approaches. Due to the amount of clutter, most of the corners were captured on the ceiling, which reduces the accuracy of the measurements. Magic Plan assumes the angle between two consecutive walls is 90 • or 45 • , for this reason, the results are unsatisfactory on the scene Lab2 which does not follow the MW assumption. We can also observe the residual increases with the area of the room, which is coherent when there is a small error with the device height estimation.

For selling or renting a property in France, the Alur law defines the maximal error of the measured area to 5%. The area errors from the FloorPlanEx application and our approach are inferior or equal to this threshold, which may not be enough for some official uses. The results show our method is generally more accurate and provides more repeatable results than FloorPlanEx. One explanation is that we consider the 3D points from multiple frames to estimate the planes of the walls when the FloorPlanEx only considers the points from one frame.

The last column of Table 1 describes the degree of interaction. It confirms our approach generally requires fewer screen interactions than user-driven approaches, except for the scene House3 MW which contained a fireplace, high furnitures and many curtains. The Lab1 MW was also quite challenging because of the presence of a high cupboard and pillars which were incorrectly labeled as wall.

Conclusion

We presented an online approach to estimate and measure the room layout of an indoor scene using a Project Tango mobile device. The various components of our pipeline were evaluated on a desktop computer, using scans recorded from a Project Tango mobile device. Compared to existing online works, the proposed method relies on scene understanding to compute the room layout. In order to cope with the complexity of some rooms, the user can interact in a cooperative way to add or remove walls among the observed planes. The comparison with other mobile applications demonstrates that in general, we achieve higher accuracy. In term of efficiency, the number of user interactions depends on the complexity of the scenes, which may contain clutter data incorrectly classified as walls. Some limitations of the proposed method are implicitly associated to the Tango TDK: drift of the device pose and surfaces not handled by the depth sensors, such as mirrors and glasses. Applications of our works include floor plan generation, area estimation and room reconstruction for virtual reality games.

Fig. 1 :

 1 Fig. 1: (a) Tango Tablet Development Kit image sensors: RGB-IR camera (1), Fisheye camera (2), IR pattern emitter (3). (b) Our pipeline.

Fig. 2 :

 2 Fig. 2: Left: the three blue, cyan and green triangles represent the part of the camera frustum viewing s v and associated with the camera poses c 1 , c 2 and c 3 respectively. Their union form P s : the visibility view polygon associated with the line segment s v . Right: the visibility polygon is the union of the visibility view polygons P s , here represented with different colors.

Fig. 3 :

 3 Fig. 3: Visualization of the scan progress of the scene House2. Planar patches classified as wall and clutter are represented in green and blue respectively. Intermediate colors represent intermediate probabilities. Left: augmented view of the device camera with the detected and classified planar patches. The estimated room layout is displayed with black lines. Right: visibility polygon in light gray, camera view polygon in black.

  Creation of a weighted graph (bottom) from the set of segments with their probabilities (top). The visibility polygon is represented by a gray thick line. Some edges are overlapping, for visualization purpose, we represented them with a curved dashed arrow line. Red edges correspond to segments from s c (including frustum edges). (b) Room layout toy example illustrating various cases: open door on the top, missing data in the upper right corner, clutter at the bottom. Orange dashed lines correspond to segments created during the post-processing step.

Fig. 4 :

 4 Fig. 4: Toy examples illustrating the graph creation and the solving. Faded colors represent the segments which are not part of the solution. Red edges correspond to frustum edges, purple for the other edges.

( 1

 1 cm) to simplify near parallel adjacent segments. Contrary to the offline approaches mentioned earlier, ours can deliver an immediate room layout. Our Python implementation of this component takes 200 ms to generate the graph, and 11 ms to compute the shortest path and apply the post-processing steps.

Fig. 6 :

 6 Fig. 6: First row: the point clouds and the labeled planar patches (blue: clutter, green: walls) of the scanned rooms: Lab1 MW , Lab2, House1 MW , House2 MW and House3 MW (from left to right). Second row: comparison of our approach in orange with the ground truth in blue. Third and fourth row : comparison of the room layout obtained with FloorPlanEx and Magic Plan in orange and the ground truth in blue.

Table 1 :

 1 Magic Plan 15.4 % 24.9 % 233 mm 79 mm 163 mm 344 mm 7 Results of the geometry accuracy and reproducibility comparison experiment.

	Scene	Method	Mean	Max	Mean	σ resid.	Min	Max	Number
			area err.	area err.	residual		residual	residual	Interact.
	Lab1 MW (25 m 2 )	Ours FloorPlanEx 2.2 % 2.3 % Magic Plan 12 %	4.2 % 29 mm 14 mm 14 mm 48 mm 4.5% 47 mm 26 mm 18 mm 84 mm 17 % 164 mm 52 mm 106 mm 231 mm	3.25 4 5
	Lab2 (47 m 2 )	Ours FloorPlanEx 3.3 % 1.1 % Magic Plan 15 %	2.4 % 38 mm 3 mm 35 mm 43 mm 4.3 % 73 mm 19 mm 45 mm 100 mm 24% 264 mm 65 mm 199 mm 329 mm	0.75 6 7
	House1 MW (11 m 2 )	Ours FloorPlanEx 2.6 % 1.8 % Magic Plan 4.9 %	2.4 % 30 mm 12 mm 14 mm 44 mm 4.5 % 53 mm 9 mm 38 mm 60 mm 8.8 % 66 mm 18 mm 46 mm 87 mm	0.5 6 7
	House2 MW (13 m 2 )	Ours FloorPlanEx 3.3 % 3.0 % Magic Plan 5.4 %	5.0 % 4.2 % 9.2 %	29 mm 12 mm 17 mm 49 mm 41 mm 11 mm 28 mm 58 mm 50 mm 23 mm 22 mm 89 mm	0 4 5
	House3 MW (48 m 2 )	Ours FloorPlanEx 2.8 % 1.9 %	2.3 % 32 mm 5 mm 27 mm 37 mm 4.0 % 105 mm 23 mm 78 mm 144 mm	7.5 6