README.md

+# PopSimul
+
+## Authors
+Cécile Triay (IRD)
+Alice Boizet (CIRAD)
+Mathias Lorieux (IRD)
+Mathieu Triay (Atelier Triay)
+
+## Description
+This tool enables to simulate a VCF of a bi-parental population. It can be used to produce a VCF under controlled conditions with setting parameters such as average depth, marker density and expected error rate. It produces 2 VCF files with the simulated genotypes, one with simulated data, and one with the same simulated data but with simulated genotyping errors. 
+The Breakpoints csv file contains the positions of exact genotypes transitions based on the simulated data. This tool has been used to test NOISYmputer (https://gitlab.cirad.fr/noisymputer/noisymputerstandalone) to compare exact Breakpoints position to Breakpoints positions determined on imputed data.
+
+## Requirements
+python 3
+
+## Running popsimul
+
+``` python popsimul.py -o test1511.vcf -ni 10 -s 44000000 -cM 180 -nm 220000 -a 3 -m 10 -eA 0.03 -eB 0.05```
+
+### Parameters description
+
+You can see all parameters by running this command  
+``` python popsimul.py -h```
+
+| Short  | Long | Description | Default value |
+|---|---|---|---|  
+| -o  | --output | Name for the 3 generated files :  [output].vcf, [output]_errors.vcf, [output]_Breakpoints.csv | popsimul |
+| -ni  | --nb_individuals | number of individuals | 10 |
+| -bp  | --chr_bp | chromosome size | 44000000|
+| -cM | --cM_size | size of genetic map | 180 |
+| -nm  | --nb_markers | number of markers| 220000 |
+| -a  | --average_depth | average depth | 3 |
+| -m  | --max_depth | max depth | 10 |
+| -eA | --error_rate_A | error rate for A | 0.05 |
+| -eB | --error_rate_B | error rate for B | 0.05 |
+
+

popsimul.py

+import random
+import math
+import numpy as np
+import argparse
+
+#def generate_chr_for_one_ind (chromosome_size, marker_density, err_rate, mean_depth, max_depth, markers_positions, conversion_factor):
+def generate_chr_for_one_ind (mean_depth, markers_positions, conversion_factor, errA, errB):
+    '''
+    Generate a chromosome for a given size, density of markers, and depth (follwing a gaussian distrib with a mean_depth and sd_depth).
+    
+    Args:
+        mean_depth (float): The mean sequencing depth of the chromosome
+        markers_positions (lst) : liste of marker positions (picked on a grid)
+        conversion_factor (float) : Value of the bp per cM conversion (bp chromosome size / cM chromosome size)
+        errA & errB (float): Respectively the error rate of observing a B whereas the genotype is truely a A and a A whereas the genotype is truely a B
+    
+    Returns:
+        segment: An array containing for each snp of the chromosome, this array : [initial_gt, genotype, genotype_error, site_depth, x_error, y_error]        
+        - intial_gt : simulated genotype with no noise and no sequencing error
+        - genotype : simulated genotype with noise 
+        - genotype_error : simulated genotype with noise and sequencing error
+        - site_depth : simulated read depth 
+        - x_error : simulated allelic depth for A
+        - y_error : simulated allelic depth for B
+    '''
+    
+    # Create an empty list to store the chromosome
+    segment = []
+    breakpoints = []
+
+    # Iterate through the list of markers positions (fixed)
+    previous_marker_gt = ""
+    for i in range(0, len(markers_positions)):
+        recomb1 = False
+        recomb2 = False
+        # If the current marker is not the first one
+        if i > 0:
+            # Calculate the interval between the current marker and the previous one
+            #using the Kosambi inverse function to estimate the chance of recombination
+            IntervalWithPreviousMarker = markers_positions[i] - markers_positions[i-1]
+            tempcM = 2 * (IntervalWithPreviousMarker / 100) / conversion_factor         ## Conversion factor is used to have the bp per cM 
+            IntervalWithPreviousMarkerInRF = 0.5 * ((math.exp(tempcM) - math.exp(-tempcM)) / (math.exp(tempcM) + math.exp(-tempcM)))
+            rnd = random.random()
+
+            # If the random number is greater than the probability of a recombination a recombination occurs
+            # Done for genotype 1
+            if rnd > 1 - IntervalWithPreviousMarkerInRF: 
+                recomb1 = True
+                # If the previous marker is A, change the genotype to B
+                if genotype1 == "A":
+                    genotype1 = "B"
+                else :
+                    genotype1 = "A"
+                    
+            # Done for genotype 2
+            rnd = random.random()
+            if rnd > 1 - IntervalWithPreviousMarkerInRF:
+                recomb2 = True
+                # If the previous marker is B, change the genotype to A
+                if genotype2 == "A":
+                    genotype2 = "B"
+                else :
+                    genotype2 = "A"
+                    
+        # If the current marker is the first one, Set the two first genotypes (random).
+        else:
+            rnd = random.random()
+            # If the random number is greater than 0.5
+            if rnd > 0.5:
+                genotype1 = "A"
+            else : 
+                genotype1 = "B"
+
+            rnd = random.random()
+            if rnd > 0.5:
+                genotype2 = "B"
+            else : 
+                genotype2 = "A"
+            
+            # initialize previous_marker_gt
+            if genotype1 == "A" and genotype2 == "A":
+                previous_marker_gt = "A"
+            elif genotype1 == "B" and genotype2 == "B":
+                previous_marker_gt = "B"
+            else :
+                previous_marker_gt = "H"
+
+        # Breakpoints for this individual
+        if recomb1 or recomb2:
+            #this is a breakpoint [beforeBkpAfter, transition]
+            transition = previous_marker_gt + " => "
+
+            if genotype1 == "A" and genotype2 == "A":
+                previous_marker_gt = "A"
+            elif genotype1 == "B" and genotype2 == "B":
+                previous_marker_gt = "B"
+            else :
+                previous_marker_gt = "H"
+            
+            transition = transition + previous_marker_gt
+        
+            breakpoints.append([i,transition]) 
+
+        # Calculate the site depth of the current marker    
+        g = np.random.poisson(mean_depth)
+        # If the depth of the current marker is less than 0 (safeguard)
+        if g < 0:
+            g = 0
+        # Round the depth of the current marker
+        site_depth = round(g)
+        # Initialize x and y (with and without error)
+        x = 0
+        y = 0
+        x_error = 0
+        y_error = 0
+
+        # If the site is homozygous A (ref, 0/0)
+        if genotype1 == "A" and genotype2 == "A":
+            x = site_depth
+            y = 0
+            initial_gt  = "0/0"
+            genotype = "0/0"
+            x_error = 0
+            y_error = 0
+            # Considering the sequencing and mapping error at each site for a given site depth
+            for j in range(0, site_depth):
+                rnd = random.random()
+                # If the random number is smaller the error rate, we increase the wronge allele depth
+                if rnd < errA:
+                    x_error = x_error + 0
+                    y_error = y_error + 1
+                else : 
+                    x_error = x_error + 1
+                    y_error = y_error + 0
+            if y_error == site_depth:
+                genotype_error = "1/1"
+            elif y_error > 0:
+                genotype_error = "0/1"
+            else :
+                genotype_error = "0/0"
+
+        # If the current marker is B and the previous marker is B, genotype if homozygote B (alt, 1/1)
+        elif genotype1 == "B" and genotype2 == "B":
+            x = 0
+            y = site_depth
+            initial_gt  = "1/1"
+            genotype = "1/1"
+            x_error = 0
+            y_error = 0
+            # Considering the sequencing and mapping error at each site for a given site depth
+            for j in range(0,site_depth):
+                rnd = random.random()
+                # If the random number is smaller the error rate, we increase the wronge allele depth
+                if rnd < errB:
+                    x_error = x_error + 1
+                    y_error = y_error + 0
+                else : 
+                    x_error = x_error + 0
+                    y_error = y_error + 1
+            if x_error == site_depth:
+                genotype_error = "0/0"
+            elif x_error > 0:
+                genotype_error = "0/1"
+            else :
+                genotype_error = "1/1"
+                
+        # If the current marker is neither A nor B, genotype is heterozygous H (0/1)
+        else:
+            initial_gt  = "0/1"
+            # Generate a random number between 0 and the depth of the current marker
+            x = random.randint(0,site_depth)
+            y = site_depth - x
+            # Error of A and B compensate themselves in heterozygous site so no need to change the x_error and y_error
+            x_error = x
+            y_error = y
+            # If the depth of x of the current marker is 0, it's seen as homozygous site (alt, 1/1)
+            if x == 0:
+                genotype = "1/1"
+                genotype_error = "1/1"
+            # If the depth of x of the current marker is 0, it's seen as homozygous site (ref, 0/0)
+            elif y == 0:
+                genotype = "0/0"
+                genotype_error = "0/0"
+            # If the depth x and y of the current marker is not 0, it's seen as heterozygous (0/1)
+            else :
+                genotype = "0/1"
+                genotype_error = "0/1"  
+
+        if site_depth == 0:
+            genotype = "./."
+            genotype_error = "./."
+        
+        # data is an array that stores genotypes and depth for one snp and one individual
+        data = [initial_gt , genotype, genotype_error, site_depth, x_error, y_error]
+
+        segment.append(data)
+
+    return segment, breakpoints
+
+# Generate a list of individuals
+def generate_individuals (nb_individuals, chromosome_size, nb_markers, mean_depth, conversion_factor, errA, errB):
+    matrix = [] #matrix is an array of genotypes for each individual
+    #matrix_error = []
+    
+    # Generate a sorted list of markers positions
+    #markers_positions = sorted(random.sample(range(chromosome_size),size))
+    markers_positions = np.linspace(1, chromosome_size, nb_markers, dtype="int")
+
+    with open(args.output + '_Breakpoints.csv', 'w') as breakpointFile:
+        breakpointFile.write(",".join(["sample","average_bkp_position", "bkp_start_position", "bkp_stop_position", "transitionType"]) + "\n");
+        for i in range(nb_individuals):
+            print("individual " + str(i+1) + " / " + str(nb_individuals), end="\r")
+            bkps_nb = 0
+            result = []
+
+            while bkps_nb < 1 : #resimulate until we get at least one crossing over for this individual
+                result = generate_chr_for_one_ind(mean_depth, markers_positions, conversion_factor, errA, errB)
+                bkps_nb = len(result[1])
+
+            matrix.append(result[0])
+
+            for bkp in result[1]:
+                startBkp = markers_positions[bkp[0] - 1]
+                stopBkp = markers_positions[bkp[0]]
+                averageBkpPosition = startBkp + round((stopBkp - startBkp)/2)
+                bkp_row = [str(i), str(averageBkpPosition), str(startBkp), str(stopBkp), bkp[1]]
+                breakpointFile.write(",".join(bkp_row) + "\n");
+    print("The file", args.output + "_Breakpoint.csv" , "has been generated")
+    
+    return matrix, markers_positions
+
+# Write VCF file
+def write_vcf(output_file_name, matrix, gt_data_position, with_depth):
+    # Header to add to the VCF for it to be recognized as such file.    
+    header = [
+    "##fileformat=VCFv4.1\n",
+    "##fileDate=20090805\n",
+    "##source=myImputationProgramV3.1\n",
+    "##reference=/shared/projects/recombinationlandscape/REFERENCE_GENOME/Osat_Azucena_AGI_chrOK_uniline_WithNIPBARorganelles.fasta\n",
+    "##contig=<ID=chr01,length=44011168>\n",
+    "##phasing=partial\n",
+    "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n",
+    "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Read Depth\">\n",
+    "##FORMAT=<ID=AD,Number=2,Type=Integer,Description=\"Allelic Depth\">\n"
+    ]
+    nb_individuals = len(matrix)
+    nb_snp = len(matrix[0])
+
+    with open(output_file_name + '.vcf', 'w') as file:
+        # Write the header of the file
+        file.writelines(header)
+        
+        # Write the header of the file
+        file.write("\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT"] + [str(n) for n in range(nb_individuals)] + ["Parent1", "Parent2"]) + "\n")
+        
+        # Iterate over each row in the table
+        for m in range(0, nb_snp):            
+            if with_depth:
+                format_info = "GT:DP:AD"
+            else :
+                format_info = "GT"
+            row = ["chr01", str(markers_positions[m]), ".", "A", "T", ".", ".", ".", format_info]
+
+            # add all individual genotypes
+            for i in range(0, nb_individuals):
+                vcf_data = str(matrix[i][m][gt_data_position])
+                if with_depth:
+                    vcf_data = vcf_data + ":" + str(matrix[i][m][3]) + ":" + str(matrix[i][m][4]) + "," + str(matrix[i][m][5])
+                row.append(vcf_data)
+            
+            # add parent genotypes
+            if with_depth:
+                p1_data =  "0/0:"+ str(round(mean_depth)) + ":" + str(round(mean_depth)) + ",0"
+                p2_data =  "1/1:"+ str(round(mean_depth)) + ":0," + str(round(mean_depth))
+            else :
+                p1_data =  "0/0"
+                p2_data =  "1/1"
+
+            row.append(p1_data)
+            row.append(p2_data)
+            file.write("\t".join(row) + "\n")
+
+
+# # Check if a row is correct
+# def is_correct (row):
+#     l = len(row)
+
+#     count_h = row.count("H") / l
+#     count_a = row.count("A") / l
+#     count_b = row.count("B") / l
+    
+#     return [count_a, count_b, count_h, 0.45 <= count_h <= 0.55 and 0.2 <= count_b <= 0.3 and 0.2 <= count_a <= 0.3]
+
+
+# define commandline arguments
+parser = argparse.ArgumentParser(description='PopSimul')
+
+parser.add_argument('-o', '--output', help='Output files name', default="popsimul")
+parser.add_argument('-ni', '--nb_individuals', help='Number of individuals', default=10)
+parser.add_argument('-bp', '--chr_bp', help='Chromosome size in bp', default=44000000)
+parser.add_argument('-cM', '--cM_size', help='Size of genetic map', default=180)
+parser.add_argument('-nm', '--nb_markers', help='Number of markers', default=220000)
+parser.add_argument('-a', '--average_depth', help='Depth average', default=3)
+parser.add_argument('-m', '--max_depth', help='Maximum possible depth', default=10)
+parser.add_argument('-eA', '--error_rate_A', help='Error rate', default=0.05)
+parser.add_argument('-eB', '--error_rate_B', help='Error rate', default=0.05)
+
+args = parser.parse_args()
+
+# Set the parameters for the script
+nb_individuals = int(args.nb_individuals)
+chromosome_size = int(args.chr_bp)
+cMsize = int(args.cM_size)    ## Size of genetic map
+conversion_factor = chromosome_size/cMsize   ## Corresponds to a bpPercM conversion ! Needs to be fixed... Does not produce the correct division!
+#conversion_factor = 233000
+#marker_density = 255000/44000000
+nb_markers = int(args.nb_markers)
+mean_depth = float(args.average_depth)
+max_depth = float(args.max_depth) #TODO better estimate for max ?
+errA = float(args.error_rate_A)
+errB = float(args.error_rate_B)
+
+print("-----USED PARAMETERS-----")
+print("nb_individuals =",str(nb_individuals))
+print("chr_bp =",str(chromosome_size))
+print("cMsize =",str(cMsize))
+print("nb_markers =",str(nb_markers))
+print("average_depth =",str(mean_depth))
+print("max_depth =",str(max_depth))
+print("error_rate_A =",str(errA))
+print("error_rate_B =",str(errB))
+print("------------------------")
+
+## Generate the pop and associate the results to the matrixes with and without errors
+pop = generate_individuals(nb_individuals, chromosome_size, nb_markers, mean_depth, conversion_factor, errA, errB)
+matrix = pop[0]
+markers_positions = pop[1]
+
+# Header to add to the VCF for it to be recognized as such file.    
+header = [
+"##fileformat=VCFv4.1\n",
+"##fileDate=20090805\n",
+"##source=myImputationProgramV3.1\n",
+"##reference=/shared/projects/recombinationlandscape/REFERENCE_GENOME/Osat_Azucena_AGI_chrOK_uniline_WithNIPBARorganelles.fasta\n",
+"##contig=<ID=chr01,length=44011168>\n",
+"##phasing=partial\n",
+"##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n",
+"##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Read Depth\">\n",
+"##FORMAT=<ID=AD,Number=2,Type=Integer,Description=\"Allelic Depth\">\n"
+]
+
+## Write simulated VCF file with no noise and no sequencing error
+write_vcf(args.output, matrix, 0, True)
+print("The file", args.output + ".vcf has been generated")
+
+## Write simulated VCF file with noise but no sequencing error
+filename = args.output + "_noise"
+write_vcf(filename, matrix, 1, True)
+print("The file", filename + ".vcf has been generated")
+
+## Write simulated VCF file with noise and sequencing error
+filename = args.output + "_noise_errors"
+write_vcf(filename, matrix, 2, True)
+print("The file", filename + ".vcf has been generated")

sumulation.vba

+'Simulation module for genetic maps and populations
+
+Function gauss()
+'simulated normally-distributed data
+'use: myvalue = gauss * sd + mean
+
+Dim fac As Double, r As Double, V1 As Double, V2 As Double
+    
+10  V1 = 2 * Rnd - 1
+    V2 = 2 * Rnd - 1
+    r = V1 ^ 2 + V2 ^ 2
+    If (r >= 1) Then GoTo 10
+    fac = Sqr(-2 * Log(r) / r)
+    gauss = V2 * fac
+    
+End Function
+
+Sub testgauss()
+    Dim mean As Variant
+    Dim sd As Variant
+
+    mean = 0
+    sd = 1
+    myvalue = gauss * sd + mean
+    MsgBox myvalue
+End Sub
+
+Public Sub Simulate_Map()
+'For i = 1 To 1000
+Application.StatusBar = "Simulating map..."
+Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+ThisWorkbook.Sheets("SimulMap").Activate
+
+Dim SimulatedChromosome As Integer
+Dim NumberOfSimulatedChromosomes As Integer
+Dim SimulatedMarker As Long
+Dim NumberOfSimulatedMarkers As Long
+Dim NumberOfMarkersInThisChromosome As Long
+Dim SimulRow As Long
+Dim SimulatedMapSize As Double
+Dim ChromosomeSize As Double
+Dim MaxVariationPositionRate As Double
+Dim VariationPosition As Double
+Dim SimulatedMarkerPosition As Double
+Dim AverageChromosomeSize As Double
+Dim SimulatedMarkerDensity As Double
+Dim SimulatedChromosomeSizeVariation As Double
+Dim TotalMapSize As Double
+Dim StartChrSize As Double
+Dim GenomeSize As LongLong 'long is max 2,147,483,647
+Dim nbOfMb As LongLong
+Dim bp_cM_ratio As Double
+
+'Map parameters
+NumberOfSimulatedChromosomes = [N9].Value
+SimulatedMapSize = [N10].Value
+SimulatedMarkerDensity = [N11].Value
+SimulatedChromosomeSizeVariation = [N12].Value / 100
+nbOfMb = [j22].Value
+GenomeSize = nbOfMb * 1000000
+bp_cM_ratio = GenomeSize / SimulatedMapSize
+
+If SimulatedMarkerDensity < 0 Then EquallyDispersed = 1
+If [j11] = "yes" Then EquallyDispersed = 2
+
+If EquallyDispersed = 1 Then
+    SimulatedMarkerDensity = 0 - SimulatedMarkerDensity
+ElseIf EquallyDispersed = 0 Then
+    SimulatedMarkerDensity = SimulatedMarkerDensity + SimulatedMarkerDensity * (NumberOfSimulatedChromosomes / 1000)
+    'gives better results than SimulatedMarkerDensity alone
+ElseIf EquallyDispersed = 2 Then
+    If SimulatedMarkerDensity < 0 Then SimulatedMarkerDensity = (0 - SimulatedMarkerDensity)
+End If
+
+MaxVariationPositionRate = [J12]
+
+'Map Simulation
+'Define Chromosome sizes
+AverageChromosomeSize = SimulatedMapSize / NumberOfSimulatedChromosomes
+'ExpectedMarkerNumber = Int(SimulatedMapSize / SimulatedMarkerDensity) + 1
+'ExpectedMarkerNumberPerChromosome = Int(ExpectedMarkerNumber / NumberOfSimulatedChromosomes) + 1
+
+Dim SimulatedMap() As Variant
+ReDim SimulatedMap(1000000, 4)
+
+If [b6] <> "" Then
+    msg = "Are you sure you want to clear this simulated map?"
+    Style = vbYesNo + vbCritical + vbDefaultButton1
+    Title = "Clear results"
+    Ctxt = 1000
+    response = MsgBox(msg, Style, Title)
+    If response = vbYes Then
+        ThisWorkbook.Sheets("SimulMap").Range("B6", "F100005").ClearContents
+        ThisWorkbook.Sheets("SimulMap").Range("N22:N24").ClearContents
+        ThisWorkbook.Sheets("SimulMap").Range("J27:O1000").ClearContents
+        ThisWorkbook.Sheets("SimulMap").Range("R22:R25").ClearContents
+    Else
+        GoTo 1000
+    End If
+End If
+    
+SimulRow = 0
+NumberOfSimulatedMarkers = 0
+TotalMapSize = 0
+
+[K6] = "Simulating chromosomes... "
+Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+
+For SimulatedChromosome = 1 To NumberOfSimulatedChromosomes
+
+    StartChrSize = 0
+    Randomize (Rnd)
+    If Rnd > 0.5 Then SimulatedChromosomeSizeVariation = 0 - SimulatedChromosomeSizeVariation
+'    Randomize (Rnd)
+    ChromosomeSize = AverageChromosomeSize + AverageChromosomeSize * (Rnd * SimulatedChromosomeSizeVariation)
+    NumberOfMarkersInThisChromosome = Int(ChromosomeSize / SimulatedMarkerDensity) + 1
+    For SimulatedMarker = 1 To NumberOfMarkersInThisChromosome
+        NumberOfSimulatedMarkers = NumberOfSimulatedMarkers + 1
+        If EquallyDispersed = 1 Then 'fixed positions
+            SimulatedMarkerPosition = (ChromosomeSize / (NumberOfMarkersInThisChromosome - 1)) * (SimulatedMarker - 1)
+        ElseIf EquallyDispersed = 0 Then ' random positions
+            Randomize (Rnd)
+            SimulatedMarkerPosition = Rnd * ChromosomeSize
+        ElseIf EquallyDispersed = 2 Then ' semi-random positions
+            Randomize (Rnd)
+            SimulatedMarkerPosition = (ChromosomeSize / (NumberOfMarkersInThisChromosome - 1)) * (SimulatedMarker - 1)
+            If MaxVariationPositionRate > 0 Then
+                If Rnd > 0.5 Then
+                    VariationPosition = (ChromosomeSize / (NumberOfMarkersInThisChromosome - 1)) * Rnd * MaxVariationPositionRate
+                Else
+                    VariationPosition = 0 - ((ChromosomeSize / (NumberOfMarkersInThisChromosome - 1)) * Rnd * MaxVariationPositionRate)
+                End If
+            Else
+                VariationPosition = 0
+            End If
+            SimulatedMarkerPosition = SimulatedMarkerPosition + VariationPosition
+        End If
+        ''Debug.Print  SimulatedChromosome
+        SimulatedMap(SimulRow, 0) = SimulatedChromosome
+        SimulatedMap(SimulRow, 1) = NumberOfSimulatedMarkers
+        fakeMkphysicalPosition = Int(SimulatedMarkerPosition * bp_cM_ratio)
+'        SimulatedMap(SimulRow, 2) = "Chr" & SimulatedChromosome & "_M" & SimulatedMarker
+        SimulatedMap(SimulRow, 2) = "Chr" & SimulatedChromosome & "_" & fakeMkphysicalPosition
+
+        If SimulatedMarker = 1 Then
+                SimulatedMap(SimulRow, 3) = 0
+        ElseIf SimulatedMarker = NumberOfMarkersInThisChromosome Then
+            SimulatedMap(SimulRow, 3) = ChromosomeSize
+'            SimulatedMap(SimulRow, 3) = Format(ChromosomeSize, "####.00")
+        Else
+            SimulatedMap(SimulRow, 3) = SimulatedMarkerPosition
+'            SimulatedMap(SimulRow, 3) = Format(SimulatedMarkerPosition, "####.00") 'very slow
+        End If
+'        'If SimulatedMarkerPosition > StartChrSize Then StartChrSize = SimulatedMarkerPosition
+        SimulRow = SimulRow + 1
+    Next SimulatedMarker
+    'FinalChromosomeSize = StartChrSize
+    TotalMapSize = TotalMapSize + ChromosomeSize
+    Cells(26 + SimulatedChromosome, 10) = SimulatedChromosome
+    Cells(26 + SimulatedChromosome, 11) = Format(ChromosomeSize, "####.00")
+    Cells(26 + SimulatedChromosome, 12) = NumberOfMarkersInThisChromosome
+    Cells(26 + SimulatedChromosome, 14) = Format(ChromosomeSize / NumberOfMarkersInThisChromosome, "####.000")
+'    chrcM_Mbp = ChromosomeSize / (fakeMkphysicalPosition / 1000000)
+'    Cells(26 + SimulatedChromosome, 15) = Format(chrcM_Mbp, "####.00")
+Next SimulatedChromosome
+
+[K6] = ""
+
+[N22] = TotalMapSize
+[N23] = NumberOfSimulatedMarkers
+[N24] = Format(TotalMapSize / NumberOfSimulatedMarkers, "####.000")
+
+[R22] = [N9]
+[r23] = [N10]
+[R24] = [N11]
+[R25] = [N12]
+
+    Cells(26, 10) = "Chr"
+    Cells(26, 11) = "cM"
+    Cells(26, 12) = "Mk"
+    Cells(26, 14) = "density"
+
+'Application.StatusBar = "Writing map..."
+'Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+'ThisWorkbook.Sheets("SimulMap").Activate
+'Range("B6:E100005").Select
+'Selection.ClearContents
+
+ThisWorkbook.Sheets("SimulMap").Range("B6").Resize(NumberOfSimulatedMarkers, 4) = SimulatedMap
+Erase SimulatedMap
+
+'Application.StatusBar = "Sorting map by position..."
+'Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+
+Range("B6:E100005").Select
+Selection.Sort Key1:=Range("B6"), Order1:=xlAscending, Key2:=Range("E6") _
+, Order2:=xlAscending, Header:=xlGuess, OrderCustom:=1, MatchCase:= _
+False, Orientation:=xlTopToBottom
+
+[b6].Select
+
+'Application.StatusBar = "Sorting map by locus number..."
+'Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+
+Range("C6:D100000").Select
+Selection.Sort Key1:=Range("C6"), Order1:=xlAscending, Header:=xlGuess, _
+    OrderCustom:=1, MatchCase:=False, Orientation:=xlTopToBottom
+    
+Range("N22").Select
+Application.ScreenUpdating = True
+'Next i
+1000
+Application.StatusBar = False
+End Sub
+
+
+Public Sub Chr_Start_at_0()
+Application.ScreenUpdating = False
+
+If [b6] <> "" Then
+    msg = "Are you sure you want to do this? All marker positions will be recalculated in order to let chromosomes start at position 0"
+    Style = vbYesNo + vbCritical + vbDefaultButton1
+    Title = "Start at 0"
+    Ctxt = 1000
+    response = MsgBox(msg, Style, Title)
+    If response = vbNo Then
+        GoTo 1000
+    End If
+End If
+
+ThisWorkbook.Sheets("SimulMap").Activate
+
+'first, counts the chr parameters (same than Read_Chr_Simulated)
+Range("B6", Range("B1048576").End(xlUp)).Select
+NbOfRows = Selection.Rows.Count
+
+NumbOfMarkInChrTotal = 0
+TotalMapSize = 0
+NumberOfChr = 0
+For i = 6 To NbOfRows + 6
+    If Cells(i + 1, 2) <> Cells(i, 2) Then 'new chr is detected
+        NumberOfChr = NumberOfChr + 1
+        Cells(26 + Cells(i, 2), 10) = Cells(i, 2) 'chr number
+        ChrSize = Cells(i, 5)
+        Cells(26 + Cells(i, 2), 11) = ChrSize 'chr size
+        TotalMapSize = TotalMapSize + Cells(i, 5)
+        NumbOfMarkInChr = i - 5 - NumbOfMarkInChrTotal
+        Cells(26 + Cells(i, 2), 12) = NumbOfMarkInChr
+        MarkerDensity = ChrSize / NumbOfMarkInChr
+        Cells(26 + Cells(i, 2), 14) = MarkerDensity
+        NumbOfMarkInChrTotal = NumbOfMarkInChrTotal + NumbOfMarkInChr
+    End If
+Next
+[N22] = TotalMapSize
+[N23] = NumbOfMarkInChrTotal
+[N24] = TotalMapSize / NumbOfMarkInChrTotal
+[R22] = NumberOfChr
+
+'second, reformats the marker positions to start at 0
+firstRow = 6
+NbOfProcessedRows = 0
+For chromosome = 1 To NumberOfChr
+    NumbOfMarkInChr = Cells(26 + chromosome, 12)
+    FirstPosition = Cells(NbOfProcessedRows + 6, 5)
+    For mark = 1 To NumbOfMarkInChr
+        Cells(NbOfProcessedRows + mark + 5, 5) = Cells(NbOfProcessedRows + mark + 5, 5) - FirstPosition
+    Next
+    NbOfProcessedRows = NbOfProcessedRows + NumbOfMarkInChr
+Next chromosome
+
+'third, re-reads and writes the chr parameters
+NumbOfMarkInChrTotal = 0
+TotalMapSize = 0
+NumberOfChr = 0
+For i = 6 To NbOfRows + 6
+    If Cells(i + 1, 2) <> Cells(i, 2) Then 'new chr is detected
+        NumberOfChr = NumberOfChr + 1
+        Cells(26 + Cells(i, 2), 10) = Cells(i, 2) 'chr number
+        ChrSize = Cells(i, 5)
+        Cells(26 + Cells(i, 2), 11) = ChrSize 'chr size
+        TotalMapSize = TotalMapSize + Cells(i, 5)
+        NumbOfMarkInChr = i - 5 - NumbOfMarkInChrTotal
+        Cells(26 + Cells(i, 2), 12) = NumbOfMarkInChr
+        MarkerDensity = ChrSize / NumbOfMarkInChr
+        Cells(26 + Cells(i, 2), 14) = MarkerDensity
+        NumbOfMarkInChrTotal = NumbOfMarkInChrTotal + NumbOfMarkInChr
+    End If
+Next
+[N22] = TotalMapSize
+[N23] = NumbOfMarkInChrTotal
+[N24] = TotalMapSize / NumbOfMarkInChrTotal
+[R22] = NumberOfChr
+
+1000 [N22].Select
+Application.CutCopyMode = False
+Application.ScreenUpdating = True
+End Sub
+
+
+Sub Add_Trait_To_Population_In_Data()
+
+Application.ScreenUpdating = False
+Init
+
+popType = Sheets("Data").[D4]
+popsize = Sheets("Data").[D5]
+nbloci = Sheets("Data").[D6]
+NbTraits = Sheets("Data").[D7]
+
+nbHeaders = 15
+
+outcellrow = nbloci + NbTraits + nbHeaders + 1
+
+mkQTL1 = ThisWorkbook.Sheets("SimulMap").[R39]
+NbOfTraitsToAdd = ThisWorkbook.Sheets("SimulMap").[R40]
+
+Dim SimulatedTrait() As Variant
+ReDim SimulatedTrait(1 To NbOfTraitsToAdd, 1 To popsize + 2)
+
+MeanAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R31]
+MeanAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R32]
+MeanBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R33]
+VarCoeffAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R34]
+VarCoeffAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R35]
+VarCoeffBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R36]
+SD_AA_QTL1 = MeanAA_QTL1 * VarCoeffAA_QTL1
+SD_AB_QTL1 = MeanAB_QTL1 * VarCoeffAB_QTL1
+SD_BB_QTL1 = MeanBB_QTL1 * VarCoeffBB_QTL1
+
+Randomize (Rnd)
+
+If popType = "SSD" Or popType = "DH" Then
+    For j = 1 To NbOfTraitsToAdd
+        SimulatedTrait(j, 1) = NbTraits + j
+        SimulatedTrait(j, 2) = "SimulatedTrait_" & NbTraits + j
+        For i = 1 To popsize
+            Genotype = Data(i - 1, mkQTL1 - 1)
+            If Genotype = "A" Then
+                myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                SimulatedTrait(j, i + 2) = Format(myAAvalue, "###.00")
+            ElseIf Genotype = "B" Then
+                myBBvalue = gauss * SD_BB_QTL1 + MeanBB_QTL1
+                SimulatedTrait(j, i + 2) = Format(myBBvalue, "###.00")
+            Else
+                SimulatedTrait(j, i + 2) = "-"
+            End If
+        Next
+    Next
+ElseIf popType = "BC1" Then
+    For j = 1 To NbOfTraitsToAdd
+        SimulatedTrait(j, 1) = NbTraits + j
+        SimulatedTrait(j, 2) = "SimulatedTrait_" & NbTraits + j
+        For i = 1 To popsize
+            Genotype = Data(i - 1, mkQTL1 - 1)
+            If Genotype = "A" Then
+                myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                SimulatedTrait(j, i + 2) = Format(myAAvalue, "###.00")
+            ElseIf Genotype = "H" Then
+                myABvalue = gauss * SD_AB_QTL1 + MeanAB_QTL1
+                SimulatedTrait(j, i + 2) = Format(myABvalue, "###.00")
+            Else
+                SimulatedTrait(j, i + 2) = "-"
+            End If
+        Next
+    Next
+ElseIf popType = "F2" Then
+    For j = 1 To NbOfTraitsToAdd
+        SimulatedTrait(j, 1) = NbTraits + j
+        SimulatedTrait(j, 2) = "SimulatedTrait_" & NbTraits + j
+        For i = 1 To popsize
+            Genotype = Data(i - 1, mkQTL1 - 1)
+            If Genotype = "H" Then
+                myABvalue = gauss * SD_AB_QTL1 + MeanAB_QTL1
+                SimulatedTrait(j, i + 2) = Format(myABvalue, "###.00")
+            ElseIf Genotype = "A" Then
+                myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                SimulatedTrait(j, i + 2) = Format(myAAvalue, "###.00")
+            ElseIf Genotype = "B" Then
+                myBBvalue = gauss * SD_BB_QTL1 + MeanBB_QTL1
+                SimulatedTrait(j, i + 2) = Format(myBBvalue, "###.00")
+            Else
+                SimulatedTrait(j, i + 2) = "-"
+            End If
+        Next
+    Next
+End If
+
+
+Sheets("Data").Activate
+Cells(outcellrow, 1).Resize(NbOfTraitsToAdd, popsize + 2) = SimulatedTrait
+
+Sheets("Data").[D7] = NbTraits + NbOfTraitsToAdd
+End Sub
+
+Sub Add_Random_Trait_To_Population_In_Data()
+
+Application.ScreenUpdating = False
+Init
+
+popType = Sheets("Data").[D4]
+popsize = Sheets("Data").[D5]
+nbloci = Sheets("Data").[D6]
+NbTraits = Sheets("Data").[D7]
+
+nbHeaders = 15
+
+outcellrow = nbloci + NbTraits + nbHeaders + 1
+
+NbOfTraitsToAdd = ThisWorkbook.Sheets("SimulMap").[R43]
+
+Dim SimulatedTrait() As Variant
+ReDim SimulatedTrait(1 To NbOfTraitsToAdd, 1 To popsize + 2)
+
+MeanAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R31]
+MeanAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R32]
+MeanBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R33]
+VarCoeffAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R34]
+VarCoeffAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R35]
+VarCoeffBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R36]
+SD_AA_QTL1 = MeanAA_QTL1 * VarCoeffAA_QTL1
+SD_AB_QTL1 = MeanAB_QTL1 * VarCoeffAB_QTL1
+SD_BB_QTL1 = MeanBB_QTL1 * VarCoeffBB_QTL1
+
+Randomize (Rnd)
+
+For j = 1 To NbOfTraitsToAdd
+    SimulatedTrait(j, 1) = NbTraits + j
+    SimulatedTrait(j, 2) = "RandomTrait_" & NbTraits + j
+    For i = 1 To popsize
+        myvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+        SimulatedTrait(j, i + 2) = Format(myvalue, "###.00")
+    Next
+Next
+
+Sheets("Data").Activate
+Cells(outcellrow, 1).Resize(NbOfTraitsToAdd, popsize + 2) = SimulatedTrait
+
+Sheets("Data").[D7] = NbTraits + NbOfTraitsToAdd
+End Sub
+
+Public Sub Read_Simulation_Parameters()
+'only for autobatchmap
+SimulatedPopulationType = VirtualMap(16, 14) 'ONLY BC1 & F2 pops at this stage
+
+SimulatedPopulationSize = VirtualMap(17, 14)
+
+MapFunctionSimul = VirtualMap(18, 14)
+
+ChromosomeSDL1 = VirtualMap(9, 18)
+If ChromosomeSDL1 > 0 Then
+    PositionSDL1 = VirtualMap(10, 18)
+    ViabU = VirtualMap(15, 18) 'differential viability fA/fB for SDL1
+    SDL1AbsoluteNumber = 0
+    For i = 1 To ChromosomeSDL1 - 1
+        SDL1AbsoluteNumber = SDL1AbsoluteNumber + VirtualMap(26 + i, 12)
+    Next
+    SDL1AbsoluteNumber = SDL1AbsoluteNumber + PositionSDL1
+End If
+
+ChromosomeSDL2 = VirtualMap(11, 18)
+If ChromosomeSDL2 > 0 Then
+    PositionSDL2 = VirtualMap(12, 18)
+    ViabV = VirtualMap(17, 18) 'differential viability fA/fB for SDL2
+    SDL2AbsoluteNumber = 0
+    For i = 1 To ChromosomeSDL2 - 1
+        SDL2AbsoluteNumber = SDL2AbsoluteNumber + VirtualMap(26 + i, 12)
+    Next
+    SDL2AbsoluteNumber = SDL2AbsoluteNumber + PositionSDL2
+End If
+
+NumberOfSimulatedTraits = VirtualMap(19, 14)
+NumberOfSimulatedMarkers = VirtualMap(23, 14)
+NumberOfSimulatedChromosomes = VirtualMap(22, 18)
+RandomErrorRate = VirtualMap(18, 10)
+RandomMissingDataRate = VirtualMap(19, 10)
+
+End Sub
+
+Public Sub Simulate_Population_From_Map()
+
+'version for excel 2007+ ONLY
+'On Error GoTo 900
+Dim SegDistortersAreIndependant As Boolean, SDL2ConditionalToSDL1 As Boolean
+Dim genotypeNames() As Variant
+
+Application.ScreenUpdating = False
+
+Read_Parameters
+
+SegDistortersAreIndependant = True
+SDL2ConditionalToSDL1 = False
+
+'QTL1 parameters
+ChrQTL1 = ThisWorkbook.Sheets("SimulMap").[R29]
+MarkerNumberQTL1 = ThisWorkbook.Sheets("SimulMap").[R30]
+MeanAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R31]
+MeanAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R32]
+MeanBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R33]
+VarCoeffAA_QTL1 = ThisWorkbook.Sheets("SimulMap").[R34]
+VarCoeffAB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R35]
+VarCoeffBB_QTL1 = ThisWorkbook.Sheets("SimulMap").[R36]
+SD_AA_QTL1 = MeanAA_QTL1 * VarCoeffAA_QTL1
+SD_AB_QTL1 = MeanAB_QTL1 * VarCoeffAB_QTL1
+SD_BB_QTL1 = MeanBB_QTL1 * VarCoeffBB_QTL1
+Mean_QTL1_for_Random_Backcross = 0.5 * (MeanAA_QTL1 + MeanAB_QTL1)
+SD_AA_QTL1_for_Random_Backcross = Mean_QTL1_for_Random_Backcross * 0.5 * (VarCoeffAA_QTL1 + VarCoeffAB_QTL1)
+Mean_QTL1_for_Random_DH = 0.5 * (MeanAA_QTL1 + MeanBB_QTL1)
+Mean_QTL1_for_Random_F2 = 0.5 * (MeanAA_QTL1 + MeanBB_QTL1)
+
+If AutoBatchMapOn = 0 Then Application.StatusBar = "Starting..."
+
+'Population Parameters
+Sheets("SimulMap").[K6] = "Simulating population... ": Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+
+If AutoBatchMapOn = 1 Then
+    NumberOfSimulatedPopulations = 1
+    Read_Simulation_Parameters
+Else
+    ThisWorkbook.Sheets("SimulMap").Activate
+    NumberOfSimulatedPopulations = ThisWorkbook.Sheets("SimulMap").[N15]
+    saveInFiles = ThisWorkbook.Sheets("SimulMap").[O15]
+    VirtualMap = ThisWorkbook.Sheets("SimulMap").UsedRange.Value
+    
+    SimulatedPopulationType = ThisWorkbook.Sheets("SimulMap").[N16] 'ONLY BC1 & F2 pops at this stage
+    If SimulatedPopulationType <> "BC1" And SimulatedPopulationType <> "F2" _
+        And SimulatedPopulationType <> "BC2F1" And SimulatedPopulationType <> "SSD" _
+            Then MsgBox "BC1 or F2 or SSD populations only, sorry.": GoTo 1000
+    SimulatedPopulationSize = ThisWorkbook.Sheets("SimulMap").[N17] 'NOT more than 508 individuals at this stage
+'    If MaxNumberOfColumns = 256 Then
+        If SimulatedPopulationSize > 16382 And NumberOfSimulatedPopulations = 1 Then
+            SimulatedPopulationSize = 1
+            MsgBox "I can't simulate more than 16382 individuals. If you need more individuals, please simulate 2 or more populations; this will generate Mapmaker/EXP compatible files that you can analyze using the BatchMap command."
+        End If
+'    End If
+    
+    ChromosomeSDL1 = ThisWorkbook.Sheets("SimulMap").[R9]
+    If ChromosomeSDL1 > 0 Then
+        PositionSDL1 = ThisWorkbook.Sheets("SimulMap").[R10]
+        ViabU = ThisWorkbook.Sheets("SimulMap").[R15]  'differential viability fA/fB for SDL1 (BC1 or F2) - gametic selection
+        ViabU1 = ThisWorkbook.Sheets("SimulMap").[R16] 'differential viability fAA/fBB for SDL1 (F2) - zygotic selection
+        ViabU2 = ThisWorkbook.Sheets("SimulMap").[R17] 'differential viability fAB/fBB for SDL1 (F2) - zygotic selection
+        ViabU3 = ViabU2 / ViabU1                       'differential viability fAB/fAA for SDL1 (F2) - zygotic selection
+        SDL1AbsoluteNumber = 0
+        For i = 1 To ChromosomeSDL1 - 1
+            SDL1AbsoluteNumber = SDL1AbsoluteNumber + ThisWorkbook.Sheets("SimulMap").Cells(26 + i, 12)
+        Next
+        SDL1AbsoluteNumber = SDL1AbsoluteNumber + PositionSDL1
+    End If
+    
+    ChromosomeSDL2 = ThisWorkbook.Sheets("SimulMap").[R11]
+    If ChromosomeSDL2 > 0 Then
+        PositionSDL2 = ThisWorkbook.Sheets("SimulMap").[R12]
+        ViabV = ThisWorkbook.Sheets("SimulMap").[R18]
+        ViabV1 = ThisWorkbook.Sheets("SimulMap").[R19]
+        ViabV2 = ThisWorkbook.Sheets("SimulMap").[R20]
+        ViabV3 = ViabV2 / ViabV1                        'differential viability fAB/fAA for SDL2 (F2) - zygotic selection
+        SDL2AbsoluteNumber = 0
+        For i = 1 To ChromosomeSDL2 - 1
+            SDL2AbsoluteNumber = SDL2AbsoluteNumber + Cells(26 + i, 12)
+        Next
+        SDL2AbsoluteNumber = SDL2AbsoluteNumber + PositionSDL2
+    End If
+    
+    'F2 options for segregation distortion
+    SelTypeGameticMaleSDL1 = ThisWorkbook.Sheets("SimulMap").[X15]
+    SelTypeGameticFemaleSDL1 = ThisWorkbook.Sheets("SimulMap").[X16]
+    SelTypeZygoticSDL1 = ThisWorkbook.Sheets("SimulMap").[X17]
+    SelTypeGameticMaleSDL2 = ThisWorkbook.Sheets("SimulMap").[X18]
+    SelTypeGameticFemaleSDL2 = ThisWorkbook.Sheets("SimulMap").[X19]
+    SelTypeZygoticSDL2 = ThisWorkbook.Sheets("SimulMap").[X20]
+    
+    NumberOfSimulatedTraits = 1
+    NumberOfSimulatedRandomTraits = ThisWorkbook.Sheets("SimulMap").[N19]
+    
+    nbOfSupPops = 1
+    If NumberOfSimulatedPopulations > 1 And saveInFiles = True Then
+        OperatingSystem = ThisWorkbook.Sheets("hide").Cells(8, 1)
+        OperatingSystem = UCase(Left(OperatingSystem, 3))
+        If OperatingSystem = "WIN" Then
+            FileSaveName = Application.GetSaveAsFilename(Title:="Locate the folder to export the data (Don't rename the MyPop.txt file please): ", InitialFileName:="MyPop.txt", fileFilter:="Text Files (*.txt), *.txt")
+        ElseIf OperatingSystem = "MAC" Then
+            FileSaveName = Application.GetSaveAsFilename(Title:="Locate the folder to export the data (Don't rename the MyPop.txt file please): ", InitialFileName:="MyPop.txt", fileFilter:="TEXT")
+        End If
+        If FileSaveName = False Then GoTo 1000
+        
+        LengthOfPath = Len(FileSaveName)
+        OutputPath = Left(FileSaveName, LengthOfPath - 9)
+        'MsgBox "The output folder will be " & Chr(13) & OutputPath
+    ElseIf NumberOfSimulatedPopulations > 1 And saveInFiles = False Then
+        nbOfSupPops = NumberOfSimulatedPopulations
+        subPopSize = Int(SimulatedPopulationSize / nbOfSupPops)
+        remainderSubPopSize = SimulatedPopulationSize - nbOfSupPops * subPopSize
+        NumberOfSimulatedPopulations = 1
+        
+        ReDim genotypeNames(1 To 1, 1 To SimulatedPopulationSize)
+        outcolsub = 1
+        For subPop = 1 To nbOfSupPops
+            For mygenot = 1 To subPopSize
+                genotypeNames(1, outcolsub) = "pop" & subPop & "_genot" & mygenot
+                outcolsub = outcolsub + 1
+            Next
+        Next
+        If remainderSubPopSize > 0 Then
+            genotnumber = subPopSize + 1
+            For mygenot = nbOfSupPops * subPopSize + 1 To SimulatedPopulationSize
+                genotypeNames(1, mygenot) = "pop" & nbOfSupPops & "_genot" & genotnumber
+                genotnumber = genotnumber + 1
+            Next
+        End If
+        MapFunctionSimul = ThisWorkbook.Sheets("SimulMap").[N18]
+    ElseIf NumberOfSimulatedPopulations = 1 And saveInFiles = False Then
+        ReDim genotypeNames(1 To 1, 1 To SimulatedPopulationSize)
+        For mygenot = 1 To SimulatedPopulationSize
+            genotypeNames(1, mygenot) = "pop1" & "_genot" & mygenot
+        Next
+
+        MapFunctionSimul = ThisWorkbook.Sheets("SimulMap").[N18]
+        
+    End If
+
+    
+    'Reading map parameters
+    Application.StatusBar = "Reading map parameters"
+    
+    'TotalMapSize = [N22]
+    NumberOfSimulatedMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+    NumberOfSimulatedChromosomes = ThisWorkbook.Sheets("SimulMap").[R22]
+    RandomErrorRate = [j18]
+    RandomMissingDataRate = [j19]
+    
+End If
+'If NumberOfSimulatedPopulations > 1 And AutoBatchMapOn = 0 Then  ' simulate genotypes and write in txt files
+
+Dim SimulatedGenotypes() As Variant
+
+Randomize (Rnd)
+
+For SimulatedPopulation = 1 To NumberOfSimulatedPopulations
+    
+        If AutoBatchMapOn = 0 Then Application.StatusBar = "Simulating population " & SimulatedPopulation & " of " & NumberOfSimulatedPopulations
+        
+        'ReDim SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedPopulationSize + 1) As Variant
+        ReDim SimulatedGenotypes(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits + 1, SimulatedPopulationSize + 1)
+        
+'BACKCROSS
+        If SimulatedPopulationType = "BC1" Then 'fake, because does not simulate properly the 4 possible gametes. However gives correct RFs (+ or -)
+
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+10
+                If AutoBatchMapOn = 0 Then Application.StatusBar = "Pop " & SimulatedPopulation & " of " & NumberOfSimulatedPopulations & ": Genot " & SimulatedGenotype & " of " & SimulatedPopulationSize
+                StartRow = 1
+                
+                For chromosome = 1 To NumberOfSimulatedChromosomes
+                
+                    ChrSize = VirtualMap(26 + chromosome, 11)
+                    NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        MarkerPositionOld = 0
+                        If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                        SimulatedGenotypes(StartRow, SimulatedGenotype) = Genotype
+                        For Marker = 2 To NbOfMrksInChr
+                        'Application.StatusBar = "Pop " & SimulatedPopulation & " of " & NumberOfSimulatedPopulations & " - Chr " & chromosome & " - Mk " & marker
+        
+                            MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                            IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                            If MapFunctionSimul = 1 Then
+                                IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                            ElseIf MapFunctionSimul = 2 Then
+                                temp = 2 * (IntervalWithPreviousMarker / 100)
+                                IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                            End If
+                            'IntervalWithPreviousMarkerInRF = IntervalWithPreviousMarker / 100
+                            
+                            If Rnd <= IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                                If Genotype = "A" Then
+                                    Genotype = "H"
+                                ElseIf Genotype = "H" Then
+                                    Genotype = "A"
+                                End If
+                            End If
+                            
+                            SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                            
+                            'simulate non-random trait
+                            If ChrQTL1 = chromosome Then
+                                If Marker = MarkerNumberQTL1 Then
+                                    If Genotype = "A" Then
+                                        myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                                        SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myAAvalue, "###.00")
+                                    ElseIf Genotype = "H" Then
+                                        myABvalue = gauss * SD_AB_QTL1 + MeanAB_QTL1
+                                        SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myABvalue, "###.00")
+                                    End If
+                                End If
+                            End If
+                            MarkerPositionOld = MarkerPosition
+                        Next 'Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                Next 'chromosome
+                
+                GameteIsDead = 0
+                GameteIsDeadSDL1 = 0
+                GameteIsDeadSDL2 = 0
+                
+                If ChromosomeSDL1 > 0 Then ' there is a segregation distorter
+                    FreqA = ViabU / (1 + ViabU) 'final frequency, 0 to 0.5
+                    'SurvivalA = FreqA * 2 '0 to 1
+                    FreqH = 1 / (1 + ViabU) 'final frequency
+                    'SurvivalH = FreqH * 2
+                    myRandkillSDL1 = Rnd
+                    If ViabU < 1 Then 'selection against AA genotypes or A gametes
+'                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > FreqA Then 'you're dead, gamete A!
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "A" And myRandkillSDL1 > ViabU Then  'you're dead, gamete A!
+                            GameteIsDeadSDL1 = 1
+                        End If
+                    ElseIf ViabU > 1 Then 'selection against H genotypes or B gametes
+'                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill > FreqH Then  ' you're dead, gamete B!
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "H" And myRandkillSDL1 > 1 / ViabU Then ' you're dead, gamete B!
+                            GameteIsDeadSDL1 = 1
+                        End If
+                    End If
+                End If
+                
+                If ChromosomeSDL2 > 0 Then ' there is another segregation distorter
+                    FreqA = ViabV / (1 + ViabV) 'final frequency, 0 to 0.5
+                    'SurvivalA = FreqA * 2 '0 to 1
+                    FreqH = 1 / (1 + ViabV) 'final frequency
+                    'SurvivalH = FreqH * 2
+                    myRandkillSDL2 = Rnd
+                    If ViabV < 1 Then 'selection against AA genotypes or A gametes
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "A" And myRandkillSDL2 > ViabV Then  '... you're dead, gamete A!
+                            GameteIsDeadSDL2 = 1
+                        End If
+                    ElseIf ViabV > 1 Then 'selection against H genotypes or B gametes
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "H" And myRandkillSDL2 > 1 / ViabV Then   '... you're dead, gamete B!
+                            GameteIsDeadSDL2 = 1
+                        End If
+                    End If
+                End If
+                
+                If SegDistortersAreIndependant = True Then 'no epistasis
+                    If GameteIsDeadSDL1 = 1 Or GameteIsDeadSDL2 = 1 Then GameteIsDead = 1
+                ElseIf SegDistortersAreIndependant = False Then 'both death required, epistasis
+                    If SDL2ConditionalToSDL1 = True Then
+                        If GameteIsDeadSDL1 = 0 And GameteIsDeadSDL2 = 1 Then GameteIsDead = 1
+                    ElseIf SDL2ConditionalToSDL1 = False Then
+                        If GameteIsDeadSDL1 = 1 And GameteIsDeadSDL2 = 1 Then GameteIsDead = 1
+                    End If
+                End If
+    
+                If GameteIsDead = 1 Then GoTo 10
+                
+                'random traits (normally distributed)
+                For mytrait = 1 To NumberOfSimulatedRandomTraits
+                    myvalue = gauss * SD_AA_QTL1_for_Random_Backcross + Mean_QTL1_for_Random_Backcross
+                    SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, SimulatedGenotype) = Format(myvalue, "###.00")
+'                    SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+                Next 'mytrait
+                
+            Next 'SimulatedGenotype
+            
+            If RandomErrorRate > 0 And RandomErrorRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= (1 - RandomErrorRate) Then
+                                If SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A" Then
+                                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                                ElseIf SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H" Then
+                                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                                End If
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+            
+            If RandomMissingDataRate > 0 And RandomMissingDataRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= 1 - RandomMissingDataRate Then
+                                SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "-"
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+            
+            For mytrait = 1 To NumberOfSimulatedRandomTraits
+                SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, 0) = "RandomQTrait_" & mytrait
+            Next 'mytrait
+
+            SimulatedGenotypes(NumberOfSimulatedMarkers + 1, 0) = "SimulatedQTrait_1"
+            
+            
+'BACKCROSS - RARE
+        ElseIf SimulatedPopulationType = "BC1bizarre" Then 'try to mimic better meiosis behavior
+        
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+                StartRow = 1
+                For chromosome = 1 To NumberOfSimulatedChromosomes
+                    ChrSize = VirtualMap(26 + chromosome, 11)
+                    NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        MarkerPositionOld = 0
+                        Marker = StartRow
+                        If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                        SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                        If Rnd > 0.5 Then 'Recombined Gamete
+                            For Marker = 2 To NbOfMrksInChr
+                                MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                                IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld 'in cM
+                                If MapFunctionSimul = 1 Then
+                                    IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                                ElseIf MapFunctionSimul = 2 Then
+                                    temp = 2 * (IntervalWithPreviousMarker / 100)
+                                    IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                                End If
+                                'Randomize (Rnd)
+                                If Rnd > 1 - 2 * IntervalWithPreviousMarkerInRF Then 'a recombination occurred - seems to generate too many 2ble recombinations
+                                    If Genotype = "A" Then
+                                        Genotype = "H"
+                                    ElseIf Genotype = "H" Then
+                                        Genotype = "A"
+                                    End If
+                                End If
+                                SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                                MarkerPositionOld = MarkerPosition
+                            Next 'Marker
+                        Else 'parental gamete
+                            For Marker = 2 To NbOfMrksInChr
+                                SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                            Next 'Marker
+                        End If
+                    StartRow = StartRow + NbOfMrksInChr
+                Next 'chromosome
+                
+            Next 'SimulatedGenotype
+        
+        
+'BC2F1
+        ElseIf SimulatedPopulationType = "BC2F1" Then
+'   A FINIR
+
+            BC2FamilySize = 4 ' a mettre en variable
+            
+            BC2PopulationSize = SimulatedPopulationSize * BC2FamilySize
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+50
+                StartRow = 1
+                For chromosome = 1 To NumberOfSimulatedChromosomes
+                    ChrSize = VirtualMap(26 + chromosome, 11)
+                    NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        MarkerPositionOld = 0
+                        Marker = StartRow
+                        If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                        SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                        For Marker = 2 To NbOfMrksInChr
+                            MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                            IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                            If MapFunctionSimul = 1 Then
+                                IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                            ElseIf MapFunctionSimul = 2 Then
+                                temp = 2 * (IntervalWithPreviousMarker / 100)
+                                IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                            End If
+                            If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                                If Genotype = "A" Then
+                                    Genotype = "H"
+                                ElseIf Genotype = "H" Then
+                                    Genotype = "A"
+                                End If
+                            End If
+                            SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                            MarkerPositionOld = MarkerPosition
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                Next chromosome
+    '        For mytrait = 1 To NumberOfSimulatedRandomTraits
+    '            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+    '        Next mytrait
+            Next SimulatedGenotype
+    '        For mytrait = 1 To NumberOfSimulatedRandomTraits
+    '            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+    '        Next mytrait
+
+
+
+' F2 POPULATION
+        ElseIf SimulatedPopulationType = "F2" Then
+        
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+100
+                StartRow = 1
+                Gamete1_Is_Dead = 0
+                Gamete2_Is_Dead = 0
+                Zygote_Is_Dead = 0
+                For chromosome = 1 To NumberOfSimulatedChromosomes
+                    ChrSize = VirtualMap(26 + chromosome, 11)
+                    NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                    MarkerPositionOld = 0
+                    'Randomize (Rnd)
+                    If Rnd > 0.5 Then genotype1 = "A" Else genotype1 = "B"
+                    'Randomize (Rnd)
+                    If Rnd > 0.5 Then genotype2 = "B" Else genotype2 = "A"
+                    Marker = StartRow
+                    If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "A"
+                    If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "B"
+                    If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                    If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                    
+                    For Marker = 2 To NbOfMrksInChr
+                        MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                        IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                        If MapFunctionSimul = 1 Then
+                            IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                        ElseIf MapFunctionSimul = 2 Then
+                                temp = 2 * (IntervalWithPreviousMarker / 100)
+                                IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                        End If
+                        
+                        'male gamete
+                        'Randomize (Rnd)
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If genotype1 = "A" Then
+                                genotype1 = "B"
+                            ElseIf genotype1 = "B" Then
+                                genotype1 = "A"
+                            End If
+                        End If
+                        'selection on male gamete (SDL 1)
+                        If SelTypeGameticMaleSDL1 = "x" Then
+                            If chromosome = ChromosomeSDL1 Then
+                                 If Marker = PositionSDL1 Then
+                                        ''Debug.Print "SDL1 - gam sel", ViabU
+                                        myRandkill = Rnd
+                                        If ViabU < 1 Then 'selection against A gametes
+                                            If genotype1 = "A" And myRandkill > ViabU Then 'gamete A dies
+                                                Gamete1_Is_Dead = 1
+                                                ''Debug.Print "killed"
+                                            End If
+                                        ElseIf ViabU > 1 Then 'selection against B gametes
+                                            If genotype1 = "B" And myRandkill > 1 / ViabU Then 'gamete B dies
+                                                Gamete1_Is_Dead = 1
+                                                ''Debug.Print "killed"
+                                            End If
+                                        End If
+                                End If
+                            End If
+                        End If
+                        If SelTypeGameticMaleSDL2 = "x" Then
+                            'selection on male gamete (SDL 2)
+                            If chromosome = ChromosomeSDL2 Then
+                                 If Marker = PositionSDL2 Then
+                                        myRandkill = Rnd
+                                        If ViabV < 1 Then 'selection against A gametes
+                                            If genotype1 = "A" And myRandkill > ViabV Then 'gamete A dies
+                                                Gamete1_Is_Dead = 1
+                                            End If
+                                        ElseIf ViabU > 1 Then 'selection against B gametes
+                                            If genotype1 = "B" And myRandkill > 1 / ViabV Then 'gamete B dies
+                                                Gamete1_Is_Dead = 1
+                                            End If
+                                        End If
+                                End If
+                            End If
+                        End If
+                        'female gamete
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If genotype2 = "A" Then
+                                genotype2 = "B"
+                            ElseIf genotype2 = "B" Then
+                                genotype2 = "A"
+                            End If
+                        End If
+'                        Gamete2_Is_Dead = 0
+                        'selection
+                        If SelTypeGameticFemaleSDL1 = "x" Then
+                            If chromosome = ChromosomeSDL1 Then
+                                 If Marker = PositionSDL1 Then
+                                        myRandkill = Rnd
+                                        If ViabU < 1 Then 'selection against A gametes
+                                            If genotype2 = "A" And myRandkill > ViabU Then 'gamete A dies
+                                                Gamete2_Is_Dead = 1
+                                            End If
+                                        ElseIf ViabU > 1 Then 'selection against B gametes
+                                            If genotype2 = "B" And myRandkill > 1 / ViabU Then 'gamete B dies
+                                                Gamete2_Is_Dead = 1
+                                            End If
+                                        End If
+                                End If
+                            End If
+                        End If
+                        If SelTypeGameticFemaleSDL2 = "x" Then
+                            If chromosome = ChromosomeSDL2 Then
+                                 If Marker = PositionSDL2 Then
+                                        myRandkill = Rnd
+                                        If ViabV < 1 Then 'selection against A gametes
+                                            If genotype2 = "A" And myRandkill > ViabV Then 'gamete A dies
+                                                Gamete2_Is_Dead = 1
+                                            End If
+                                        ElseIf ViabU > 1 Then 'selection against B gametes
+                                            If genotype2 = "B" And myRandkill > 1 / ViabV Then 'gamete B dies
+                                                Gamete2_Is_Dead = 1
+                                            End If
+                                        End If
+                                End If
+                            End If
+                        End If
+                        'mating or selfing
+                        If genotype1 = "A" And genotype2 = "A" Then Genotype = "A"
+                        If genotype1 = "B" And genotype2 = "B" Then Genotype = "B"
+                        If genotype1 = "A" And genotype2 = "B" Then Genotype = "H"
+                        If genotype1 = "B" And genotype2 = "A" Then Genotype = "H"
+                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                        'simulate non-random trait
+                        If ChrQTL1 = chromosome Then
+                            If Marker = MarkerNumberQTL1 Then
+                                If Genotype = "H" Then
+                                    myABvalue = gauss * SD_AB_QTL1 + MeanAB_QTL1
+                                    SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myABvalue, "###.00")
+                                ElseIf Genotype = "A" Then
+                                    myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                                    SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myAAvalue, "###.00")
+                                ElseIf Genotype = "B" Then
+                                    myBBvalue = gauss * SD_BB_QTL1 + MeanBB_QTL1
+                                    SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myBBvalue, "###.00")
+                                End If
+                            End If
+                        End If
+                        MarkerPositionOld = MarkerPosition
+                        
+                    Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    
+                Next chromosome
+                
+                If Gamete1_Is_Dead = 1 Then GoTo 100
+                If Gamete2_Is_Dead = 1 Then GoTo 100
+
+                 'eliminate zygotes
+                
+                If ChromosomeSDL1 > 0 Then ' there is a segregation distorter
+                    myRandkill = Rnd
+                    If ViabU1 < 1 Then 'selection against AA genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > ViabU1 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabU1 > 1 Then 'selection against AA genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "B" And myRandkill > 1 / ViabU1 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+                    If ViabU2 < 1 Then 'selection against AB genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill > ViabU2 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabU2 > 1 Then 'selection against BB genotype versus AB genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "B" And myRandkill > 1 / ViabU2 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+                    If ViabU3 < 1 Then 'selection against AB genotype relative to AA genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill > ViabU3 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabU3 > 1 Then 'selection against AA genotype versus AB genotype
+                        If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > 1 / ViabU3 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+               End If
+                
+                If ChromosomeSDL2 > 0 Then ' there is a segregation distorter
+                    myRandkill = Rnd
+                    If ViabV1 < 1 Then 'selection against AA genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > ViabV1 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabV1 > 1 Then 'selection against AA genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "B" And myRandkill > 1 / ViabV1 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+                    If ViabV2 < 1 Then 'selection against AB genotype relative to BB genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill > ViabV2 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabV2 > 1 Then 'selection against BB genotype versus AB genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "B" And myRandkill > 1 / ViabV2 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+                    If ViabV3 < 1 Then 'selection against AB genotype relative to AA genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill > ViabV3 Then  'AA killed
+                            Zygote_Is_Dead = 1
+                        End If
+                    ElseIf ViabV3 > 1 Then 'selection against AA genotype versus AB genotype
+                        If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > 1 / ViabV3 Then 'sorry, bad luck... you're dead, gamete B!
+                            Zygote_Is_Dead = 1
+                        End If
+                    End If
+               End If
+                
+                If Zygote_Is_Dead = 1 Then GoTo 100
+                 'random traits
+                For mytrait = 1 To NumberOfSimulatedRandomTraits
+                    SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+                Next mytrait
+                
+            Next SimulatedGenotype
+            
+            SimulatedGenotypes(NumberOfSimulatedMarkers + 1, 0) = "SimulatedQTrait_1"
+            
+             'random traits names
+            For mytrait = 1 To NumberOfSimulatedRandomTraits
+                SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, 0) = "RandomQTrait_" & mytrait
+            Next mytrait
+            
+            If RandomErrorRate > 0 And RandomErrorRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= 1 - RandomErrorRate Then
+                                If SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A" Then
+                                    If Rnd >= (1 / 3) Then
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                                    Else
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B"
+                                    End If
+                                ElseIf SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B" Then
+                                    If Rnd >= (1 / 3) Then
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                                    Else
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                                    End If
+                                ElseIf SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H" Then
+                                    If Rnd >= 0.5 Then
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                                    Else
+                                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B"
+                                    End If
+                                End If
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+            If RandomMissingDataRate > 0 And RandomMissingDataRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= 1 - RandomMissingDataRate Then
+                                SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "-"
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+
+
+
+'SSD - RIL
+        ElseIf SimulatedPopulationType = "SSD" Then 'it's fake but should work. however does not allow to test for biased estimates
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+                StartRow = 1
+                For chromosome = 1 To NumberOfSimulatedChromosomes
+                    ChrSize = VirtualMap(26 + chromosome, 11)
+                    NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        MarkerPositionOld = 0
+                        Marker = StartRow
+                        If Rnd > 0.5 Then Genotype = "A" Else Genotype = "B"
+                        SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                        For Marker = 2 To NbOfMrksInChr
+                            MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                            IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                            If MapFunctionSimul = 1 Then
+                                myRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                                myR = 2 * myRF / (1 + 2 * myRF)
+                                IntervalWithPreviousMarkerInRF = myR ' fake conversion of recombination rate in SSD (the 'R')
+                                
+'                                IntervalWithPreviousMarkerInR_MartinAndHospital = 1 ' fake conversion of recombination rate in SSD (the 'R')
+                                
+                            ElseIf MapFunctionSimul = 2 Then
+                                temp = 2 * (IntervalWithPreviousMarker / 100)
+                                myRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                                myR = 2 * myRF / (1 + 2 * myRF)
+                                IntervalWithPreviousMarkerInRF = myR ' fake conversion of recombination rate in SSD (the 'R')
+                            End If
+                            'Randomize (Rnd)
+                            If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                                If Genotype = "A" Then
+                                    Genotype = "B"
+                                ElseIf Genotype = "B" Then
+                                    Genotype = "A"
+                                End If
+                            End If
+                            SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                            
+                            'simulate non-random trait
+                            If ChrQTL1 = chromosome Then
+                                If Marker = MarkerNumberQTL1 Then
+                                    If Genotype = "A" Then
+                                        myAAvalue = gauss * SD_AA_QTL1 + MeanAA_QTL1
+                                        SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myAAvalue, "###.00")
+                                    ElseIf Genotype = "B" Then
+                                        myBBvalue = gauss * SD_BB_QTL1 + MeanBB_QTL1
+                                        SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedGenotype) = Format(myBBvalue, "###.00")
+                                    End If
+                                End If
+                            End If
+                            MarkerPositionOld = MarkerPosition
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                Next chromosome
+                
+                For mytrait = 1 To NumberOfSimulatedRandomTraits
+                    SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+                Next mytrait
+                
+            Next SimulatedGenotype
+            
+            If RandomErrorRate > 0 And RandomErrorRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= 1 - RandomErrorRate Then
+                                If SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A" Then
+                                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B"
+                                ElseIf SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B" Then
+                                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                                End If
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+            If RandomMissingDataRate > 0 And RandomMissingDataRate <= 1 Then
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    StartRow = 1
+                    For chromosome = 1 To NumberOfSimulatedChromosomes
+                        NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                        For Marker = 1 To NbOfMrksInChr
+                            If Rnd >= 1 - RandomMissingDataRate Then
+                                SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "-"
+                            End If
+                        Next Marker
+                    StartRow = StartRow + NbOfMrksInChr
+                    Next chromosome
+                Next SimulatedGenotype
+            End If
+            SimulatedGenotypes(NumberOfSimulatedMarkers + 1, 0) = "SimulatedQTrait_1"
+            
+            For mytrait = 1 To NumberOfSimulatedRandomTraits
+                SimulatedGenotypes(NumberOfSimulatedMarkers + NumberOfSimulatedTraits + mytrait, 0) = "RandomQTrait_" & mytrait
+            Next 'mytrait
+'        ElseIf SimulatedPopulationType = "SSD" Then
+'
+'        'en cours a finir
+'
+'            NumberOfSelfings = 6 ' Additional selfings after the F2 = S1 so 6 gives an F8
+'
+'            For SimulatedGenotype = 1 To SimulatedPopulationSize
+'                StartRow = 1
+'                For Chromosome = 1 To NumberOfSimulatedChromosomes
+'                    ChrSize = VirtualMap(26 + Chromosome, 11)
+'                    NbOfMrksInChr = VirtualMap(26 + Chromosome, 12)
+'                    MarkerPositionOld = 0
+'                    'Randomize (Rnd)
+'                    If Rnd > 0.5 Then genotype1 = "A" Else genotype1 = "B"
+'                    'Randomize (Rnd)
+'                    If Rnd > 0.5 Then genotype2 = "B" Else genotype2 = "A"
+'                    marker = StartRow
+'                    If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(marker, SimulatedGenotype) = "A"
+'                    If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(marker, SimulatedGenotype) = "B"
+'                    If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(marker, SimulatedGenotype) = "H"
+'                    If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(marker, SimulatedGenotype) = "H"
+'
+'                    For marker = 2 To NbOfMrksInChr
+'                        MarkerPosition = VirtualMap(StartRow + marker + 4, 5)
+'                        IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+'                        If MapFunctionSimul = 1 Then
+'                            IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+'                        ElseIf MapFunctionSimul = 2 Then
+'                                temp = 2 * (IntervalWithPreviousMarker / 100)
+'                                IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+'                        End If
+'                        'gamete 1
+'                        'Randomize (Rnd)
+'                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+'                            If genotype1 = "A" Then
+'                                genotype1 = "B"
+'                            ElseIf genotype1 = "B" Then
+'                                genotype1 = "A"
+'                            End If
+'                        End If
+'                        'gamete 2
+'                        'Randomize (Rnd)
+'                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+'                            If genotype2 = "A" Then
+'                                genotype2 = "B"
+'                            ElseIf genotype2 = "B" Then
+'                                genotype2 = "A"
+'                            End If
+'                        End If
+'                        If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "AA"
+'                        If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "BB"
+'                        If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "AB"
+'                        If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "BA"
+'                        MarkerPositionOld = MarkerPosition
+'                    Next marker
+'
+'                    StartRow = StartRow + NbOfMrksInChr
+'                Next Chromosome
+'
+'                For mytrait = 1 To NumberOfSimulatedRandomTraits
+'                    SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+'                Next mytrait
+'
+'                For SelfingGeneration = 1 To NumberOfSelfings
+'                StartRow = 1
+'
+'                    For Chromosome = 1 To NumberOfSimulatedChromosomes
+'                        ChrSize = VirtualMap(26 + Chromosome, 11)
+'                        NbOfMrksInChr = VirtualMap(26 + Chromosome, 12)
+'                        MarkerPositionOld = 0
+'
+'                        If Rnd > 0.5 Then genotype1 = "A" Else genotype1 = "B"
+'
+'                        If Rnd > 0.5 Then genotype2 = "B" Else genotype2 = "A"
+'                        marker = StartRow
+'                        If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(marker, SimulatedGenotype) = "A"
+'                        If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(marker, SimulatedGenotype) = "B"
+'                        If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(marker, SimulatedGenotype) = "H"
+'                        If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(marker, SimulatedGenotype) = "H"
+'
+'                        For marker = 2 To NbOfMrksInChr
+'                            MarkerPosition = VirtualMap(StartRow + marker + 4, 5)
+'                            IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+'                            If MapFunctionSimul = 1 Then
+'                                IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+'                            ElseIf MapFunctionSimul = 2 Then
+'                                    temp = 2 * (IntervalWithPreviousMarker / 100)
+'                                    IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+'                            End If
+'                            'gamete 1
+'                            'Randomize (Rnd)
+'                            If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+'                                If genotype1 = "A" Then
+'                                    genotype1 = "B"
+'                                ElseIf genotype1 = "B" Then
+'                                    genotype1 = "A"
+'                                End If
+'                            End If
+'                            'gamete 2
+'                            'Randomize (Rnd)
+'                            If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+'                                If genotype2 = "A" Then
+'                                    genotype2 = "B"
+'                                ElseIf genotype2 = "B" Then
+'                                    genotype2 = "A"
+'                                End If
+'                            End If
+'                            If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "AA"
+'                            If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "BB"
+'                            If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "AB"
+'                            If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(StartRow + marker - 1, SimulatedGenotype) = "BA"
+'                            MarkerPositionOld = MarkerPosition
+'                        Next marker
+'
+'                        StartRow = StartRow + NbOfMrksInChr
+'                    Next Chromosome
+'
+'                Next SelfingGeneration
+'
+'            Next SimulatedGenotype
+'
+'            For mytrait = 1 To NumberOfSimulatedRandomTraits
+'                SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+'            Next mytrait
+
+
+        End If
+        
+        Application.ScreenUpdating = True: DoEvents: Application.ScreenUpdating = False
+        
+        'Writing results
+       If NumberOfSimulatedPopulations > 1 And AutoBatchMapOn = 0 Then  ' write pop in a .txt file
+       
+            FileSaveName = "SimulPop_" & SimulatedPopulation & ".txt"
+            Open FileSaveName For Output As #1
+            If SimulatedPopulationType = "BC1" Then Print #1, "data type f2 backcross"
+            If SimulatedPopulationType = "F2" Then Print #1, "data type f2 intercross"
+            Print #1, Format(SimulatedPopulationSize, "###"); Spc(1); Format(NumberOfSimulatedMarkers, "###"); Spc(1); "0"
+            For Ligne = 1 To NumberOfSimulatedMarkers
+                NomLocus = VirtualMap(Ligne + 5, 4)
+                TexteLocus = ""
+                For Colonne = 1 To SimulatedPopulationSize
+                    TexteLocus = TexteLocus & SimulatedGenotypes(Ligne, Colonne)
+                Next
+                Print #1, ; "*" & NomLocus; Spc(1); TexteLocus
+            Next
+            Close #1
+            
+        ElseIf AutoBatchMapOn = 1 Then  ' write pop in a array
+            For i = 1 To NumberOfSimulatedMarkers
+                For j = 1 To SimulatedPopulationSize
+                    Data(j - 1, i - 1) = SimulatedGenotypes(i, j) 'stay in memory
+                Next
+            Next
+
+        ElseIf NumberOfSimulatedPopulations = 1 And AutoBatchMapOn = 0 Then  ' write pop in a sheet
+            Application.StatusBar = "Writing results"
+            ClearSimulatedPopulation
+            OutPutSheetName = "SimulPop"
+            'MsgBox SimulatedPopulationSize
+'            If MaxNumberOfColumns = 256 Then
+'                If SimulatedPopulationSize <= 254 And SimulatedPopulationSize <= 508 Then
+'                    ThisWorkbook.Sheets(OutPutSheetName).Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits, SimulatedPopulationSize + 1) = SimulatedGenotypes
+'                ElseIf SimulatedPopulationSize > 254 And SimulatedPopulationSize <= 508 Then
+'                    ThisWorkbook.Sheets(OutPutSheetName).Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits, 255) = SimulatedGenotypes
+'                    ReDim SimulatedGenotypesPlus(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits, SimulatedPopulationSize - 254 + 1) As Variant
+'                    For i = 1 To NumberOfSimulatedMarkers
+'                        For j = 255 To SimulatedPopulationSize
+'                            SimulatedGenotypesPlus(i, j - 254) = SimulatedGenotypes(i, j)
+'                        Next
+'                    Next i
+'                    ThisWorkbook.Sheets("SimulPop+").Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits, SimulatedPopulationSize - 253) = SimulatedGenotypesPlus
+'                End If
+'            ElseIf MaxNumberOfColumns = 16384 Then
+                ThisWorkbook.Sheets(OutPutSheetName).Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedRandomTraits + 1, SimulatedPopulationSize + 1) = SimulatedGenotypes
+'            End If
+            Application.ScreenUpdating = False
+            ThisWorkbook.Sheets("SimulMap").Activate
+            Range("D6").Select
+            Range(Selection, Selection.End(xlDown)).Select
+            Selection.Copy
+            [N22].Select
+            ThisWorkbook.Sheets(OutPutSheetName).Activate
+            Range("B16").Select
+            Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
+                :=False, Transpose:=False
+            [c15].Resize(1, SimulatedPopulationSize) = genotypeNames
+            Range("D4") = SimulatedPopulationType
+            Range("D5") = SimulatedPopulationSize
+            Range("D6") = NumberOfSimulatedMarkers
+            Range("D7") = NumberOfSimulatedRandomTraits + 1
+            Range("D8") = nbOfSupPops
+            Range("H4") = "A"
+            Range("H5") = "B"
+            Range("H6") = "C"
+            Range("H7") = "D"
+            Range("H8") = "H"
+            Range("H9") = "-"
+            Range("K4") = "2"
+            [b10] = "Simulated " & SimulatedPopulationType & " data " & Date & " - " & Time
+            Application.CutCopyMode = False
+            [K6] = ""
+            [B15].Select
+        End If
+Next SimulatedPopulation
+
+GoTo 1000
+900
+''Debug.Print  Chromosome
+''Debug.Print  StartRow
+''Debug.Print  marker
+''Debug.Print  SimulatedGenotype
+''Debug.Print  genotype
+''Debug.Print  StartRow + marker - 1
+1000
+'Application.ScreenUpdating = True
+'Next i
+
+If AutoBatchMapOn = 0 Then
+    Sheets("SimulMap").[K6] = ""
+    Application.StatusBar = False
+End If
+
+End Sub
+
+
+
+Public Sub Simulate_Population_From_Map_Save()
+Read_Parameters
+
+'On Error GoTo 900
+Application.ScreenUpdating = False
+If AutoBatchMapOn = 0 Then If AutoBatchMapOn = 0 Then Application.StatusBar = "Starting..."
+
+ThisWorkbook.Sheets("SimulMap").Activate
+
+'Population Parameters
+NumberOfSimulatedPopulations = [N15]
+SimulatedPopulationType = [N16] 'ONLY BC1 & F2 pops at this stage
+If SimulatedPopulationType <> "BC1" And SimulatedPopulationType <> "F2" And SimulatedPopulationType <> "BC2F1" Then MsgBox "BC1 or F2 populations only, sorry.": GoTo 1000
+SimulatedPopulationSize = [N17] 'NOT more than 508 individuals at this stage
+If AutoBatchMapOn = 0 And SimulatedPopulationSize > 508 And NumberOfSimulatedPopulations = 1 Then
+    SimulatedPopulationSize = 508
+    MsgBox "I can't simulate more than 508 individuals. If you need more indivuduals, please simulate 2 or more populations; this will generate Mapmaker/EXP compatible files that you can analyze using the BatchMap command."
+End If
+
+MapFunctionSimul = [N18]
+
+ChromosomeSDL1 = [R9]
+PositionSDL1 = [R10]
+ViabU = [R15] 'differential viability fA/fB for SDL1
+SDL1AbsoluteNumber = 0
+For i = 1 To ChromosomeSDL1 - 1
+    SDL1AbsoluteNumber = SDL1AbsoluteNumber + Cells(26 + i, 12)
+Next
+SDL1AbsoluteNumber = SDL1AbsoluteNumber + PositionSDL1
+
+ChromosomeSDL2 = [R11]
+PositionSDL2 = [R12]
+ViabV = [R17] 'differential viability fA/fB for SDL2
+SDL2AbsoluteNumber = 0
+For i = 1 To ChromosomeSDL2 - 1
+    SDL2AbsoluteNumber = SDL2AbsoluteNumber + Cells(26 + i, 12)
+Next
+SDL2AbsoluteNumber = SDL2AbsoluteNumber + PositionSDL2
+
+NumberOfSimulatedTraits = [N19]
+
+If NumberOfSimulatedPopulations > 1 Then
+    If AutoBatchMapOn = 0 Then
+        OperatingSystem = ThisWorkbook.Sheets("hide").Cells(8, 1)
+        OperatingSystem = UCase(Left(OperatingSystem, 3))
+        If OperatingSystem = "WIN" Then
+            FileSaveName = Application.GetSaveAsFilename(Title:="Locate the folder to export the data (Don't rename the MyPop.txt file please): ", InitialFileName:="MyPop.txt", fileFilter:="Text Files (*.txt), *.txt")
+        ElseIf OperatingSystem = "MAC" Then
+            FileSaveName = Application.GetSaveAsFilename(Title:="Locate the folder to export the data (Don't rename the MyPop.txt file please): ", InitialFileName:="MyPop.txt", fileFilter:="TEXT")
+        End If
+        If FileSaveName = False Then GoTo 1000
+        
+        LengthOfPath = Len(FileSaveName)
+        OutputPath = Left(FileSaveName, LengthOfPath - 9)
+        'MsgBox "The output folder will be " & Chr(13) & OutputPath
+    End If
+End If
+
+'Reading map parameters
+If AutoBatchMapOn = 0 Then Application.StatusBar = "Reading map parameters"
+
+VirtualMap = ThisWorkbook.Sheets("SimulMap").UsedRange.Value
+'TotalMapSize = [N22]
+NumberOfSimulatedMarkers = [N23]
+NumberOfSimulatedChromosomes = [R22]
+
+If NumberOfSimulatedPopulations > 1 And AutoBatchMapOn = 0 Then  ' simulate genotypes and write in txt files
+    For SimulatedPopulation = 1 To NumberOfSimulatedPopulations
+        If AutoBatchMapOn = 0 Then Application.StatusBar = "Simulating population " & SimulatedPopulation & " of " & NumberOfSimulatedPopulations
+        
+        ReDim SimulatedGenotypes(NumberOfSimulatedMarkers + 1, SimulatedPopulationSize + 1) As Variant
+        
+        StartRow = 1
+        
+        If SimulatedPopulationType = "BC1" Then
+            For chromosome = 1 To NumberOfSimulatedChromosomes
+                ChrSize = VirtualMap(26 + chromosome, 11)
+                NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    MarkerPositionOld = 0
+                    'Randomize (Rnd)
+                    If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                    Marker = StartRow
+                    SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                    For Marker = 2 To NbOfMrksInChr
+                        MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                        IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                        If MapFunctionSimul = 1 Then
+                            IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                        ElseIf MapFunctionSimul = 2 Then
+                            temp = 2 * (IntervalWithPreviousMarker / 100)
+                            IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                        End If
+                        'Randomize (Rnd)
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If Genotype = "A" Then
+                                Genotype = "H"
+                            ElseIf Genotype = "H" Then
+                                Genotype = "A"
+                            End If
+                        End If
+                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                        MarkerPositionOld = MarkerPosition
+                    Next Marker
+                Next SimulatedGenotype
+                StartRow = StartRow + NbOfMrksInChr
+            Next chromosome
+            
+        ElseIf SimulatedPopulationType = "F2" Then
+            For chromosome = 1 To NumberOfSimulatedChromosomes
+                ChrSize = VirtualMap(26 + chromosome, 11)
+                NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                For SimulatedGenotype = 1 To SimulatedPopulationSize
+                    MarkerPositionOld = 0
+                    'Randomize (Rnd)
+                    If Rnd > 0.5 Then genotype1 = "A" Else genotype1 = "B"
+                    'Randomize (Rnd)
+                    If Rnd > 0.5 Then genotype2 = "B" Else genotype2 = "A"
+                    Marker = StartRow
+                    If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "A"
+                    If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "B"
+                    If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                    If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                    
+                    For Marker = 2 To NbOfMrksInChr
+                        MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                        IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                        If MapFunctionSimul = 1 Then
+                            IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                        ElseIf MapFunctionSimul = 2 Then
+                                temp = 2 * (IntervalWithPreviousMarker / 100)
+                                IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                        End If
+                        'gamete 1
+                        'Randomize (Rnd)
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If genotype1 = "A" Then
+                                genotype1 = "B"
+                            ElseIf genotype1 = "B" Then
+                                genotype1 = "A"
+                            End If
+                        End If
+                        'gamete 2
+                        'Randomize (Rnd)
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If genotype2 = "A" Then
+                                genotype2 = "B"
+                            ElseIf genotype2 = "B" Then
+                                genotype2 = "A"
+                            End If
+                        End If
+                        If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                        If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B"
+                        If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                        If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                        MarkerPositionOld = MarkerPosition
+                    Next Marker
+                Next SimulatedGenotype
+                StartRow = StartRow + NbOfMrksInChr
+            Next chromosome
+        End If
+        
+        'Writing in .txt file
+       
+        FileSaveName = "SimulPop_" & SimulatedPopulation & ".txt"
+        Open FileSaveName For Output As #1
+        If SimulatedPopulationType = "BC1" Then Print #1, "data type f2 backcross"
+        If SimulatedPopulationType = "F2" Then Print #1, "data type f2 intercross"
+        Print #1, Format(SimulatedPopulationSize, "###"); Spc(1); Format(NumberOfSimulatedMarkers, "###"); Spc(1); "0"
+        For Ligne = 1 To NumberOfSimulatedMarkers
+            NomLocus = VirtualMap(Ligne + 5, 4)
+            TexteLocus = ""
+            For Colonne = 1 To SimulatedPopulationSize
+                TexteLocus = TexteLocus & SimulatedGenotypes(Ligne, Colonne)
+            Next
+            Print #1, ; "*" & NomLocus; Spc(1); TexteLocus
+        Next
+        Close #1
+
+    Next SimulatedPopulation
+    
+ElseIf NumberOfSimulatedPopulations = 1 Or AutoBatchMapOn = 1 Then 'if only one population, write in the SimulPop sheet
+    If AutoBatchMapOn = 0 Then Application.StatusBar = "Simulating population"
+    ReDim SimulatedGenotypes(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedTraits, SimulatedPopulationSize + 1) As Variant
+    StartRow = 1
+    
+    If SimulatedPopulationType = "BC1" And ChromosomeSDL1 = 0 Then 'there is no SDL (Segregation Distortion Locus) in the genome
+        For chromosome = 1 To NumberOfSimulatedChromosomes
+            ChrSize = VirtualMap(26 + chromosome, 11)
+            NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+                MarkerPositionOld = 0
+                Marker = StartRow
+                If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                For Marker = 2 To NbOfMrksInChr
+                    MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                    IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                    If MapFunctionSimul = 1 Then
+                        IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                    ElseIf MapFunctionSimul = 2 Then
+                        temp = 2 * (IntervalWithPreviousMarker / 100)
+                        IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                    End If
+                    If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                        If Genotype = "A" Then
+                            Genotype = "H"
+                        ElseIf Genotype = "H" Then
+                            Genotype = "A"
+                        End If
+                    End If
+                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                    MarkerPositionOld = MarkerPosition
+                Next Marker
+                For mytrait = 1 To NumberOfSimulatedTraits
+                    SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+                Next mytrait
+            Next SimulatedGenotype
+            StartRow = StartRow + NbOfMrksInChr
+        Next chromosome
+        For mytrait = 1 To NumberOfSimulatedTraits
+            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+        Next mytrait
+    
+    ElseIf SimulatedPopulationType = "BC2F1" And ChromosomeSDL1 = 0 Then 'there is no SDL (Segregation Distortion Locus) in the genome
+        BC2FamilySize = 4
+        BC2PopulationSize = SimulatedPopulationSize * BC2FamilySize
+        For SimulatedGenotype = 1 To SimulatedPopulationSize
+50
+            StartRow = 1
+            For chromosome = 1 To NumberOfSimulatedChromosomes
+                ChrSize = VirtualMap(26 + chromosome, 11)
+                NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                    MarkerPositionOld = 0
+                    Marker = StartRow
+                    If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                    SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                    For Marker = 2 To NbOfMrksInChr
+                        MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                        IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                        If MapFunctionSimul = 1 Then
+                            IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                        ElseIf MapFunctionSimul = 2 Then
+                            temp = 2 * (IntervalWithPreviousMarker / 100)
+                            IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                        End If
+                        If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                            If Genotype = "A" Then
+                                Genotype = "H"
+                            ElseIf Genotype = "H" Then
+                                Genotype = "A"
+                            End If
+                        End If
+                        SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                        MarkerPositionOld = MarkerPosition
+                    Next Marker
+                StartRow = StartRow + NbOfMrksInChr
+            Next chromosome
+'        For mytrait = 1 To NumberOfSimulatedTraits
+'            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+'        Next mytrait
+        Next SimulatedGenotype
+'        For mytrait = 1 To NumberOfSimulatedTraits
+'            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+'        Next mytrait
+    
+    ElseIf SimulatedPopulationType = "BC1" And ChromosomeSDL1 > 0 Then 'there is a SDL (Segregation Distortion Locus) in the genome
+    ''Debug.Print "zzz"
+        For SimulatedGenotype = 1 To SimulatedPopulationSize
+100
+            StartRow = 1
+            For chromosome = 1 To NumberOfSimulatedChromosomes
+                ChrSize = VirtualMap(26 + chromosome, 11)
+                NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+                MarkerPositionOld = 0
+                Marker = StartRow
+                If Rnd > 0.5 Then Genotype = "A" Else Genotype = "H"
+                SimulatedGenotypes(Marker, SimulatedGenotype) = Genotype
+                For Marker = 2 To NbOfMrksInChr
+                    MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                    IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                    If MapFunctionSimul = 1 Then
+                        IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                    ElseIf MapFunctionSimul = 2 Then
+                        temp = 2 * (IntervalWithPreviousMarker / 100)
+                        IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                    End If
+                    If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                        If Genotype = "A" Then
+                            Genotype = "H"
+                        ElseIf Genotype = "H" Then
+                            Genotype = "A"
+                        End If
+                    End If
+                    SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = Genotype
+                    MarkerPositionOld = MarkerPosition
+                Next Marker
+                StartRow = StartRow + NbOfMrksInChr
+            Next chromosome
+            
+            'fA/fH = u
+            'frequency of A = fA = u /(1+u)
+            'frequency of B = fB = 1 /(1+u)
+            myRandkill = Rnd
+            If ViabU < 1 Then 'selection against AA genotypes or A gametes
+                If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > (ViabU / (1 + ViabU)) Then  'sorry, bad luck... you're dead, gamete A!
+                'kill individual/gamete and simulate again
+                    GoTo 100
+                End If
+            ElseIf ViabU > 1 Then 'selection against H genotypes or B gametes
+                If SimulatedGenotypes(SDL1AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill < (1 / (1 + ViabU)) Then  'sorry, bad luck... you're dead, gamete B!
+                'kill individual/gamete and simulate again
+                    GoTo 100
+                End If
+            End If
+            If ChromosomeSDL2 > 0 Then
+                If ViabV < 1 Then 'selection against AA genotypes or A gametes
+                    If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "A" And myRandkill > (ViabV / (1 + ViabV)) Then  'sorry, bad luck... you're dead, gamete A!
+                    'kill individual/gamete and simulate again
+                        GoTo 100
+                    End If
+                ElseIf ViabV > 1 Then 'selection against H genotypes or B gametes
+                    If SimulatedGenotypes(SDL2AbsoluteNumber, SimulatedGenotype) = "H" And myRandkill < (1 / (1 + ViabV)) Then  'sorry, bad luck... you're dead, gamete B!
+                    'kill individual/gamete and simulate again
+                        GoTo 100
+                    End If
+                End If
+            End If
+            For mytrait = 1 To NumberOfSimulatedTraits
+                SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+            Next mytrait
+        Next SimulatedGenotype
+        For mytrait = 1 To NumberOfSimulatedTraits
+            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+        Next mytrait
+    
+    ElseIf SimulatedPopulationType = "F2" Then
+        For chromosome = 1 To NumberOfSimulatedChromosomes
+            ChrSize = VirtualMap(26 + chromosome, 11)
+            NbOfMrksInChr = VirtualMap(26 + chromosome, 12)
+            For SimulatedGenotype = 1 To SimulatedPopulationSize
+                MarkerPositionOld = 0
+                'Randomize (Rnd)
+                If Rnd > 0.5 Then genotype1 = "A" Else genotype1 = "B"
+                'Randomize (Rnd)
+                If Rnd > 0.5 Then genotype2 = "B" Else genotype2 = "A"
+                Marker = StartRow
+                If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "A"
+                If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "B"
+                If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(Marker, SimulatedGenotype) = "H"
+                
+                For Marker = 2 To NbOfMrksInChr
+                    MarkerPosition = VirtualMap(StartRow + Marker + 4, 5)
+                    IntervalWithPreviousMarker = MarkerPosition - MarkerPositionOld
+                    If MapFunctionSimul = 1 Then
+                        IntervalWithPreviousMarkerInRF = 0.5 * (1 - Exp(-2 * (IntervalWithPreviousMarker / 100))) ' Haldane inverse function
+                    ElseIf MapFunctionSimul = 2 Then
+                            temp = 2 * (IntervalWithPreviousMarker / 100)
+                            IntervalWithPreviousMarkerInRF = 0.5 * ((Exp(temp) - Exp(-temp)) / (Exp(temp) + Exp(-temp))) 'Kosambi inverse funciton
+                    End If
+                    'gamete 1
+                    'Randomize (Rnd)
+                    If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                        If genotype1 = "A" Then
+                            genotype1 = "B"
+                        ElseIf genotype1 = "B" Then
+                            genotype1 = "A"
+                        End If
+                    End If
+                    'gamete 2
+                    'Randomize (Rnd)
+                    If Rnd > 1 - IntervalWithPreviousMarkerInRF Then  'a recombination occurred
+                        If genotype2 = "A" Then
+                            genotype2 = "B"
+                        ElseIf genotype2 = "B" Then
+                            genotype2 = "A"
+                        End If
+                    End If
+                    If genotype1 = "A" And genotype2 = "A" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "A"
+                    If genotype1 = "B" And genotype2 = "B" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "B"
+                    If genotype1 = "A" And genotype2 = "B" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                    If genotype1 = "B" And genotype2 = "A" Then SimulatedGenotypes(StartRow + Marker - 1, SimulatedGenotype) = "H"
+                    MarkerPositionOld = MarkerPosition
+                Next Marker
+            Next SimulatedGenotype
+            For mytrait = 1 To NumberOfSimulatedTraits
+                SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, SimulatedGenotype) = Format(100 + Rnd * 20, "###.00")
+            Next mytrait
+            StartRow = StartRow + NbOfMrksInChr
+        Next chromosome
+        For mytrait = 1 To NumberOfSimulatedTraits
+            SimulatedGenotypes(NumberOfSimulatedMarkers + mytrait, 0) = "RandomQTrait_" & mytrait
+        Next mytrait
+    End If
+    
+    If AutoBatchMapOn = 0 Then
+        Application.StatusBar = "Writing results"
+        ClearSimulatedPopulation
+        OutPutSheetName = "SimulPop"
+        'MsgBox SimulatedPopulationSize
+        If SimulatedPopulationSize <= 254 And SimulatedPopulationSize <= 508 Then
+            ThisWorkbook.Sheets(OutPutSheetName).Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedTraits, SimulatedPopulationSize + 1) = SimulatedGenotypes
+        ElseIf SimulatedPopulationSize > 254 And SimulatedPopulationSize <= 508 Then
+            ThisWorkbook.Sheets(OutPutSheetName).Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedTraits, 255) = SimulatedGenotypes
+            ReDim SimulatedGenotypesPlus(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedTraits, SimulatedPopulationSize - 254 + 1) As Variant
+            For i = 1 To NumberOfSimulatedMarkers
+                For j = 255 To SimulatedPopulationSize
+                    SimulatedGenotypesPlus(i, j - 254) = SimulatedGenotypes(i, j)
+                Next
+            Next i
+            ThisWorkbook.Sheets("SimulPop+").Range("B15").Resize(NumberOfSimulatedMarkers + 1 + NumberOfSimulatedTraits, SimulatedPopulationSize - 253) = SimulatedGenotypesPlus
+        End If
+        Application.ScreenUpdating = False
+        ThisWorkbook.Sheets("SimulMap").Activate
+        Range("D6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Selection.Copy
+        [N22].Select
+        ThisWorkbook.Sheets(OutPutSheetName).Activate
+        Range("B16").Select
+        Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
+            :=False, Transpose:=False
+            
+        Range("D4") = SimulatedPopulationType
+        Range("D5") = SimulatedPopulationSize
+        Range("D6") = NumberOfSimulatedMarkers
+        Range("D7") = NumberOfSimulatedTraits
+        Range("H4") = "A"
+        Range("H5") = "B"
+        Range("H6") = "C"
+        Range("H7") = "D"
+        Range("H8") = "H"
+        Range("H9") = "-"
+        Range("K4") = "2"
+        If SimulatedPopulationType = "BC1" Then
+            [b10] = "Simulated BC1 data " & Date & " - " & Time
+        ElseIf SimulatedPopulationType = "F2" Then
+            [b10] = "Simulated F2 data " & Date & " - " & Time
+        End If
+        Application.CutCopyMode = False
+        [B15].Select
+    ElseIf AutoBatchMapOn = 1 Then
+'    MsgBox SimulatedPopulationSize
+        For i = 1 To NumberOfSimulatedMarkers
+            For j = 1 To SimulatedPopulationSize
+                Data(j - 1, i - 1) = SimulatedGenotypes(i, j) 'stay in memory
+            Next
+        Next
+        OutPutSheetName = "Data"
+        
+        ThisWorkbook.Sheets(OutPutSheetName).Range("D4") = SimulatedPopulationType
+        ThisWorkbook.Sheets(OutPutSheetName).Range("D5") = SimulatedPopulationSize
+        ThisWorkbook.Sheets(OutPutSheetName).Range("D6") = NumberOfSimulatedMarkers
+        ThisWorkbook.Sheets(OutPutSheetName).Range("D7") = "0"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H4") = "A"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H5") = "B"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H6") = "C"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H7") = "D"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H8") = "H"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("H9") = "-"
+        ThisWorkbook.Sheets(OutPutSheetName).Range("K4") = "2"
+    End If
+    
+End If
+
+GoTo 1000
+900
+''Debug.Print  Chromosome
+''Debug.Print  StartRow
+''Debug.Print  marker
+''Debug.Print  SimulatedGenotype
+''Debug.Print  genotype
+''Debug.Print  StartRow + marker - 1
+1000
+'Application.ScreenUpdating = True
+'Next i
+If AutoBatchMapOn = 0 Then Application.StatusBar = False
+End Sub
+
+Public Sub ClearSimulatedPopulation()
+    Application.ScreenUpdating = False
+    ThisWorkbook.Sheets("SimulPop").Activate
+    Range("B15:xfd1048576").ClearContents
+    Range("D4").ClearContents
+    Range("D5").ClearContents
+    Range("D6").ClearContents
+    Range("D7").ClearContents
+    Range("H4").ClearContents
+    Range("H5").ClearContents
+    Range("H6").ClearContents
+    Range("H7").ClearContents
+    Range("H8").ClearContents
+    Range("H9").ClearContents
+    Range("K4").ClearContents
+    [b10].ClearContents
+'    ThisWorkbook.Sheets("SimulPop+").Activate
+'    Range("B15:IV60000").ClearContents
+'    Range("D4").ClearContents
+'    Range("D5").ClearContents
+'    Range("D6").ClearContents
+'    Range("D7").ClearContents
+'    Range("H4").ClearContents
+'    Range("H5").ClearContents
+'    Range("H6").ClearContents
+'    Range("H7").ClearContents
+'    Range("H8").ClearContents
+'    Range("H9").ClearContents
+'    Range("K4").ClearContents
+'    [B13].ClearContents
+    ThisWorkbook.Sheets("SimulPop").Activate
+If AutoBatchMapOn = 0 And ClearAllResultsMode = 0 Then Application.ScreenUpdating = True
+End Sub
+
+
+Public Sub Clear_BatchMap_Results()
+Application.ScreenUpdating = False
+If ThisWorkbook.Sheets("BatchMap").[h8] <> "" Then
+    If ClearAllResultsMode = 0 Then
+        msg = "Are you sure you want to clear the batch results?"
+        Style = vbYesNo + vbCritical + vbDefaultButton1
+        Title = "Clear results"
+        Ctxt = 1000
+        response = MsgBox(msg, Style, Title)
+    End If
+    
+    If response = vbYes Or ClearAllResultsMode = 1 Then
+        If MaxNumberOfColumns = 256 Then
+            ThisWorkbook.Sheets("BatchMap").Range("A5:IV60000").ClearContents
+            ThisWorkbook.Sheets("BatchMap").Range("H5:IV6").Interior.ColorIndex = xlNone
+        Else
+            ThisWorkbook.Sheets("BatchMap").Range("A5:XFD60000").ClearContents
+            ThisWorkbook.Sheets("BatchMap").Range("H5:XFD6").Interior.ColorIndex = xlNone
+        End If
+        ThisWorkbook.Sheets("BatchMap").[D4].ClearContents
+        ThisWorkbook.Sheets("BatchMap").[f4].ClearContents
+        ThisWorkbook.Sheets("BatchMap").Range("C8:G60000").Interior.ColorIndex = xlNone
+        
+    Else
+        GoTo 1000
+    End If
+End If
+'[H8].Select
+1000
+'If AutoBatchMapOn = 0 And ClearAllResultsMode = 0 Then Application.ScreenUpdating = True
+End Sub
+
+Public Sub Clear_BigBatchMap_Results()
+Application.ScreenUpdating = False
+If [c6] <> "" Then
+    If ClearAllResultsMode = 0 Then
+        msg = "Are you sure you want to clear the batch results?"
+        Style = vbYesNo + vbCritical + vbDefaultButton1
+        Title = "Clear results"
+        Ctxt = 1000
+        response = MsgBox(msg, Style, Title)
+    End If
+    If response = vbYes Or ClearAllResultsMode = 1 Then
+        If MaxNumberOfColumns = 256 Then
+            ThisWorkbook.Sheets("BigBatch").Range("A4:IV60000").ClearContents
+        Else
+            ThisWorkbook.Sheets("BigBatch").Range("A4:XFD60000").ClearContents
+        End If
+    Else
+        GoTo 1000
+    End If
+End If
+[c6].Select
+1000
+If AutoBatchMapOn = 0 And ClearAllResultsMode = 0 Then Application.ScreenUpdating = True
+End Sub
+
+Public Sub Compute_Stats_Batch()
+'On Error GoTo 900
+Application.ScreenUpdating = False
+Application.StatusBar = "Computing statistics..."
+
+ThisWorkbook.Sheets("BatchMap").Activate
+'Expected cM values from simulated map
+NumberOfMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+For i = 2 To NumberOfMarkers
+    ExpectedRelativeDistance = ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 5) - ThisWorkbook.Sheets("SimulMap").Cells(i + 4, 5)
+    ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 6) = ExpectedRelativeDistance
+    If ThisWorkbook.Sheets("SimulMap").Cells(i + 4, 2) <> ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 2) Then 'new chr is detected
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 6).ClearContents
+    End If
+Next i
+
+NumberOfProcessedFiles = [D4]
+'Average etc on batch results
+For i = 1 To NumberOfMarkers
+    mkchr = ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 2)
+    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 1) = mkchr
+    SumOfDistances = 0
+    MinDist = 1000000
+    MaxDist = 0
+    ExpectedRelativeDistance = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6)
+    If ExpectedRelativeDistance <> "" Then
+        For SimPop = 1 To NumberOfProcessedFiles
+            BatchDist = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7)
+            SumOfDistances = SumOfDistances + BatchDist
+            If BatchDist < MinDist Then MinDist = BatchDist
+            If BatchDist > MaxDist Then MaxDist = BatchDist
+        Next
+        AverageDistance = SumOfDistances / NumberOfProcessedFiles
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 3) = MinDist
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 4) = MaxDist
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 5) = AverageDistance
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7) = AverageDistance - ExpectedRelativeDistance
+    End If
+Next i
+
+For SimPop = 1 To NumberOfProcessedFiles
+    TotalMapSize = 0
+    For i = 1 To NumberOfMarkers
+        BatchDist = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7)
+        If BatchDist <> "" Then
+            TotalMapSize = TotalMapSize + BatchDist
+        End If
+    Next
+    ThisWorkbook.Sheets("BatchMap").Cells(5, SimPop + 7) = TotalMapSize
+    ThisWorkbook.Sheets("BatchMap").Cells(6, SimPop + 7) = TotalMapSize / NumberOfMarkers 'mk density
+Next
+     
+SumAverageSize = 0
+SumAverageDensity = 0
+MinSize = 1000000
+MaxSize = 0
+MinDensity = 1000000
+MaxDensity = 0
+
+For SimPop = 1 To NumberOfProcessedFiles
+    Mapsize = ThisWorkbook.Sheets("BatchMap").Cells(5, SimPop + 7)
+    mapdensity = ThisWorkbook.Sheets("BatchMap").Cells(6, SimPop + 7)
+    If Mapsize > MaxSize Then MaxSize = Mapsize
+    If Mapsize < MinSize Then MinSize = Mapsize
+    If mapdensity > MaxDensity Then MaxDensity = mapdensity
+    If mapdensity < MinDensity Then MinDensity = mapdensity
+    SumAverageSize = SumAverageSize + Mapsize
+    SumAverageDensity = SumAverageDensity + mapdensity
+Next
+AverageSize = SumAverageSize / NumberOfProcessedFiles
+AverageDensity = SumAverageDensity / NumberOfProcessedFiles
+
+[C5] = MinSize
+[D5] = MaxSize
+[c6] = MinDensity
+[D6] = MaxDensity
+[E5] = AverageSize
+[E6] = AverageDensity
+
+[G5] = "Total map size"
+[G6] = "Marker density"
+[A7] = "Chr"
+[b7] = "Marker"
+[C7] = "Min"
+[D7] = "Max"
+[E7] = "Average"
+[F7] = "Expected"
+[G7] = "Diff."
+
+Range(Cells(8, 3), Cells(7 + NumberOfMarkers, 7)).Select
+Selection.Interior.ColorIndex = 15
+Range(Cells(5, 8), Cells(6, 256)).Select
+Selection.Interior.ColorIndex = 15
+
+GoTo 1000
+900
+MsgBox "An error occurred."
+1000
+[h8].Select
+Application.StatusBar = False
+'Application.ScreenUpdating = True
+End Sub
+
+Public Sub Compute_Stats_Batch_Chi2()
+'On Error GoTo 900
+Application.ScreenUpdating = False
+Application.StatusBar = "Computing statistics..."
+
+ThisWorkbook.Sheets("BatchMap").Activate
+'Expected cM values from simulated map
+NumberOfMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+For i = 2 To NumberOfMarkers
+    ExpectedRelativeDistance = ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 5) - ThisWorkbook.Sheets("SimulMap").Cells(i + 4, 5)
+    'keep it or it won't work
+    'ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 6) = ExpectedRelativeDistance
+    If ThisWorkbook.Sheets("SimulMap").Cells(i + 4, 2) <> ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 2) Then 'new chr is detected
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 6).ClearContents
+    End If
+Next i
+
+NumberOfProcessedFiles = ThisWorkbook.Sheets("BatchMap").[D4]
+
+MyBatchresults = ThisWorkbook.Sheets("BatchMap").Range(Cells(1, 1), Cells(NumberOfMarkers + 7, NumberOfProcessedFiles + 7)).Value
+
+'Average etc on batch results
+For i = 1 To NumberOfMarkers
+    mkchr = ThisWorkbook.Sheets("SimulMap").Cells(i + 5, 2)
+    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 1) = mkchr
+    SumOfChi2 = 0
+    MinChi2 = 1000000
+    MaxChi2 = 0
+    'ExpectedRelativeDistance = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6)
+    If ExpectedRelativeDistance <> "" Then
+        For SimPop = 1 To NumberOfProcessedFiles
+            BatchChi2 = MyBatchresults(i + 7, SimPop + 7)
+            SumOfChi2 = SumOfChi2 + BatchChi2
+            If BatchChi2 < MinChi2 Then MinChi2 = BatchChi2
+            If BatchChi2 > MaxChi2 Then MaxChi2 = BatchChi2
+        Next
+        AverageChi2 = SumOfChi2 / NumberOfProcessedFiles
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 3) = MinChi2
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 4) = MaxChi2
+        ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 5) = AverageChi2
+    End If
+Next i
+
+'ranks
+SDL1on = ThisWorkbook.Sheets("SimulMap").[R9]
+SDL2on = ThisWorkbook.Sheets("SimulMap").[R11]
+If SDL1on = 1 Then SDL1number = ThisWorkbook.Sheets("SimulMap").[R10]
+If SDL2on = 1 Then SDL2number = ThisWorkbook.Sheets("SimulMap").[R12]
+'For i = 2 To NumberOfMarkers - 1
+'    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6) = 0
+'        'ExpectedRelativeDistance = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6)
+'    For SimPop = 1 To NumberOfProcessedFiles
+'        BatchChi2 = MyBatchresults(i + 7, SimPop + 7)
+'        BatchChi2Before = MyBatchresults(i + 6, SimPop + 7)
+'        BatchChi2After = MyBatchresults(i + 8, SimPop + 7)
+'        If BatchChi2 < BatchChi2Before Or BatchChi2 < BatchChi2After Then
+'            ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6) = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6) + 1
+'            If i = SDL1number Or i = SDL2number Then
+'                ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7).Select
+'                Selection.Interior.ColorIndex = 15
+'            End If
+'        End If
+''        If BatchChi2 < BatchChi2Before Then
+''            If i = SDL1number Or i = SDL2number Then
+''                ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 7) = ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 7) + 1
+''                ThisWorkbook.Sheets("BatchMap").Cells(i + 6, SimPop + 7).Select
+''                Selection.Interior.ColorIndex = 15
+''            End If
+''        End If
+''        If BatchChi2 < BatchChi2After Then
+''            If i = SDL1number Or i = SDL2number Then
+''                ThisWorkbook.Sheets("BatchMap").Cells(i + 8, 7) = ThisWorkbook.Sheets("BatchMap").Cells(i + 8, 7) + 1
+''                ThisWorkbook.Sheets("BatchMap").Cells(i + 8, SimPop + 7).Select
+''                Selection.Interior.ColorIndex = 15
+''            End If
+''        End If
+'    Next
+'    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6) = 100 * (ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6) / NumberOfProcessedFiles)
+'    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7) = 100 - ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 6)
+''    If i = SDL1number Or i = SDL2number Then
+'''        ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 7) = 100 * (ThisWorkbook.Sheets("BatchMap").Cells(i + 6, 7) / NumberOfProcessedFiles)
+'''        ThisWorkbook.Sheets("BatchMap").Cells(i + 8, 7) = 100 * (ThisWorkbook.Sheets("BatchMap").Cells(i + 8, 7) / NumberOfProcessedFiles)
+''    End If
+'Next i
+
+
+For SimPop = 1 To NumberOfProcessedFiles
+    SDL1isNotaPeak = 0
+    mychi2SDL1 = ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, SimPop + 7)
+    For i = 1 To NumberOfMarkers
+        myChi2 = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7)
+        If myChi2 > mychi2SDL1 Then
+            'it's a peak, not on SDL1
+            ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7) = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7) + 1
+                ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7).Select
+                Selection.Interior.ColorIndex = 15
+            SDL1isNotaPeak = 1
+        End If
+        Next
+        If SDL1isNotaPeak = 1 Then 'in this pop
+            'SDL1 is not a peak
+            ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 6) = ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 6) + 1
+        End If
+Next
+
+For i = 1 To NumberOfMarkers
+    ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7) = 100 * (ThisWorkbook.Sheets("BatchMap").Cells(i + 7, 7)) / NumberOfProcessedFiles
+Next
+
+ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 6) = 100 * (ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 6) / NumberOfProcessedFiles)
+ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 7) = 100 - ThisWorkbook.Sheets("BatchMap").Cells(SDL1number + 7, 6)
+
+'For SimPop = 1 To NumberOfProcessedFiles
+'    TotalMapSize = 0
+'    For i = 1 To NumberOfMarkers
+'        BatchDist = ThisWorkbook.Sheets("BatchMap").Cells(i + 7, SimPop + 7)
+'        If BatchDist <> "" Then
+'            TotalMapSize = TotalMapSize + BatchDist
+'        End If
+'    Next
+'    ThisWorkbook.Sheets("BatchMap").Cells(5, SimPop + 7) = TotalMapSize
+'    ThisWorkbook.Sheets("BatchMap").Cells(6, SimPop + 7) = TotalMapSize / NumberOfMarkers 'mk density
+'Next
+     
+SumAverageChi2 = 0
+SumAverageDensity = 0
+MinChi2 = 1000000
+MaxChi2 = 0
+MinDensity = 1000000
+MaxDensity = 0
+
+For SimPop = 1 To NumberOfProcessedFiles
+    MapChi2 = ThisWorkbook.Sheets("BatchMap").Cells(5, SimPop + 7)
+    mapdensity = ThisWorkbook.Sheets("BatchMap").Cells(6, SimPop + 7)
+    If MapChi2 > MaxChi2 Then MaxChi2 = MapChi2
+    If MapChi2 < MinChi2 Then MinChi2 = MapChi2
+    If mapdensity > MaxDensity Then MaxDensity = mapdensity
+    If mapdensity < MinDensity Then MinDensity = mapdensity
+    SumAverageChi2 = SumAverageChi2 + MapChi2
+    SumAverageDensity = SumAverageDensity + mapdensity
+Next
+AverageChi2 = SumAverageChi2 / NumberOfProcessedFiles
+AverageDensity = SumAverageDensity / NumberOfProcessedFiles
+
+[C5] = MinChi2
+[D5] = MaxChi2
+[c6] = MinDensity
+[D6] = MaxDensity
+[E5] = AverageChi2
+[E6] = AverageDensity
+
+'[G5] = "Total map size"
+'[G6] = "Marker density"
+[A7] = "Chr"
+[b7] = "Marker"
+[C7] = "Min"
+[D7] = "Max"
+[E7] = "Average"
+[F7] = "% no peak"
+[G7] = "% peak"
+
+Range(Cells(8, 3), Cells(7 + NumberOfMarkers, 7)).Select
+Selection.Interior.ColorIndex = 15
+Range(Cells(5, 8), Cells(6, 256)).Select
+Selection.Interior.ColorIndex = 15
+
+GoTo 1000
+900
+MsgBox "An error occurred."
+1000
+[h8].Select
+Application.StatusBar = False
+'Application.ScreenUpdating = True
+End Sub
+
+Public Sub ShowDistribMapSizesFromBigBatch()
+    FromBigBatch = 1
+    ShowDistributionMapSizes
+    FromBigBatch = 0
+End Sub
+
+Public Sub ShowDistributionMapSizes()
+
+Application.ScreenUpdating = False
+Application.StatusBar = "cleaning..."
+ThisWorkbook.Sheets("Distrib").Activate
+
+NumberOfIntervals = ThisWorkbook.Sheets("Distrib").[M516]
+ThisWorkbook.Sheets("Distrib").Cells.ClearContents
+ThisWorkbook.Sheets("Distrib").Cells.Interior.ColorIndex = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlDiagonalDown).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlDiagonalUp).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlEdgeLeft).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlEdgeTop).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlEdgeBottom).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlEdgeRight).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlInsideVertical).LineStyle = xlNone
+    ThisWorkbook.Sheets("Distrib").Cells.Borders(xlInsideHorizontal).LineStyle = xlNone
+    
+popType = ThisWorkbook.Sheets("SimulMap").[N16]
+popsize = ThisWorkbook.Sheets("SimulMap").[N17]
+EmapSize = ThisWorkbook.Sheets("SimulMap").[N22]
+
+If FromBigBatch = 0 Then
+    MyBatchMatrix = ThisWorkbook.Sheets("BatchMap").Range("A5:IV6").Value
+    valeurMax = ThisWorkbook.Sheets("BatchMap").[D5]
+    valeurMin = ThisWorkbook.Sheets("BatchMap").[C5]
+    NumberOfValues = ThisWorkbook.Sheets("BatchMap").[D4]
+Else
+    MyBatchMatrix = ThisWorkbook.Sheets("BigBatch").Range("C6:D10010").Value
+    
+    ThisWorkbook.Sheets("BigBatch").Activate
+    NumberOfValues = ThisWorkbook.Sheets("SimulMap").[n3]
+    '
+    Range("C6:D10005").Select
+    Range("D6").Activate
+    Selection.Sort Key1:=Range("D6"), Order1:=xlAscending, Header:=xlGuess, _
+        OrderCustom:=1, MatchCase:=False, Orientation:=xlTopToBottom
+
+    [D6].Select
+    valeurMin = [D6]
+    valeurMax = Cells(5 + NumberOfValues, 4)
+End If
+
+Application.StatusBar = "writing intervals"
+If NumberOfIntervals = 0 Or NumberOfIntervals = "" Then NumberOfIntervals = 40
+[M516] = NumberOfIntervals
+IntervalSize = Int(valeurMax - valeurMin) / NumberOfIntervals
+
+ReDim VectorFrequencies(1, NumberOfIntervals + 1) As Long
+ReDim VectorBounds(2, NumberOfIntervals + 16) As Double
+
+mininterval = valeurMin
+For Interval = 1 To NumberOfIntervals + 1
+    ThisWorkbook.Sheets("Distrib").Cells(510, Interval + 15) = mininterval
+    VectorBounds(0, Interval + 15) = mininterval
+    maxinterval = mininterval + IntervalSize
+    ThisWorkbook.Sheets("Distrib").Cells(509, Interval + 15) = maxinterval
+    VectorBounds(1, Interval + 15) = maxinterval
+    mininterval = maxinterval
+    'ThisWorkbook.Sheets("Distrib").Cells(515, interval + 15) = "0"
+    VectorFrequencies(0, Interval - 1) = 0
+Next
+
+Application.StatusBar = "counting frequencies"
+For Colonne = 8 To NumberOfValues + 7
+    If FromBigBatch = 0 Then Mapsize = MyBatchMatrix(1, Colonne) Else Mapsize = MyBatchMatrix(Colonne - 7, 2)
+    'ThisWorkbook.Sheets("Distrib").Cells(5, Colonne) = MapSize
+    For Interval = 1 To NumberOfIntervals + 1
+'        If MapSize >= ThisWorkbook.Sheets("Distrib").Cells(510, interval + 15) Then
+'            If MapSize < ThisWorkbook.Sheets("Distrib").Cells(509, interval + 15) Then
+'                VectorFrequencies(0, interval + 15) = VectorFrequencies(0, interval + 15) + 1
+'                ThisWorkbook.Sheets("Distrib").Cells(515, interval + 15) = ThisWorkbook.Sheets("Distrib").Cells(515, interval + 15) + 1
+'            End If
+'        End If
+        If Mapsize >= VectorBounds(0, Interval + 15) Then
+            If Mapsize < VectorBounds(1, Interval + 15) Then
+                VectorFrequencies(0, Interval - 1) = VectorFrequencies(0, Interval - 1) + 1
+            End If
+        End If
+    Next
+Next
+
+[P515].Resize(1, NumberOfIntervals + 1) = VectorFrequencies
+
+Application.StatusBar = "drawing..."
+For Interval = 1 To NumberOfIntervals + 1
+    effectif = VectorFrequencies(0, Interval - 1)
+    For Ligne = 1 To effectif
+        ThisWorkbook.Sheets("Distrib").Cells(500 - Ligne, Interval + 15).Interior.ColorIndex = 25
+    Next Ligne
+Next
+Application.StatusBar = "finishing..."
+
+ThisWorkbook.Sheets("Distrib").Activate
+If FromBigBatch = 0 Then
+    [M517].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.005)": [K517] = "Percent. 0.005"
+    [M518].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.995)": [K518] = "Percent. 0.995"
+    [M519].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.025)": [K519] = "Percent. 0.025"
+    [M520].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.975)": [K520] = "Percent. 0.975"
+    [M521].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.05)": [K521] = "Percent. 0.05"
+    [M522].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.95)": [K522] = "Percent. 0.95"
+    [M523].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.1666)": [K523] = "Percent. 0.167"
+    [M524].Formula = "=PERCENTILE(BatchMap!H5:IV5,0.8333)": [K524] = "Percent. 0.833"
+    [M525].Formula = "=MEDIAN(BatchMap!H5:IV5,0.8333)": [K525] = "Median"
+    [M526].Formula = "=Average(BatchMap!H5:IV5)": [K526] = "Average"
+    
+    ActiveSheet.Shapes("cbDrawAgain").Select
+    Selection.OnAction = "ShowDistributionMapSizes"
+
+    ActiveSheet.Shapes("cbBack").Select
+    Selection.OnAction = "ActiverBatchMap"
+Else
+
+    [M517].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.005)": [K517] = "Percent. 0.005"
+    [M518].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.995)": [K518] = "Percent. 0.995"
+    [M519].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.025)": [K519] = "Percent. 0.025"
+    [M520].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.975)": [K520] = "Percent. 0.975"
+    [M521].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.05)": [K521] = "Percent. 0.05"
+    [M522].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.95)": [K522] = "Percent. 0.95"
+    [M523].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.1666)": [K523] = "Percent. 0.167"
+    [M524].Formula = "=PERCENTILE(BigBatch!D6:D60000,0.8333)": [K524] = "Percent. 0.833"
+    [M525].Formula = "=MEDIAN(BigBatch!D6:D60000,0.5)": [K525] = "Median"
+    [M526].Formula = "=Average(BigBatch!D6:D60000)": [K526] = "Average"
+    
+    ActiveSheet.Shapes("cbDrawAgain").Select
+    Selection.OnAction = "ShowDistribMapSizesFromBigBatch"
+
+    ActiveSheet.Shapes("cbBack").Select
+    Selection.OnAction = "ActiverBigBatchMap"
+    
+End If
+
+[M516] = NumberOfIntervals
+[o516] = "Histogram of map sizes"
+[K516] = "Number of intervals:"
+[o518] = "Based on [P1 x P2] simulated map "
+[o519] = NumberOfValues & " " & popType & " populations of " & popsize & " individuals simulated"
+[o520] = "Expected map size: " & Format(EmapSize, "####.##") & " cM"
+
+1000
+[M516].Select
+Application.StatusBar = False
+
+End Sub
+
+
+Public Sub Draw_Simulated_Map()
+Application.StatusBar = "Starting"
+Application.ScreenUpdating = False
+Dim MyTime
+MyTime = Time
+
+'Init
+'If Probleme = 1 Then GoTo 1000
+
+' Sort_Frameed_Map Macro
+ThisWorkbook.Sheets("SimulMap").Activate
+
+If [b6] = "" Then MsgBox "No simulated map to draw": GoTo 1000
+
+'remise ö z_ro de la feuille de r_sultats ("Draw")
+Application.DisplayAlerts = False
+    ThisWorkbook.Sheets("DrawSimulMap").Delete
+    ThisWorkbook.Sheets("DrawSimulMap0").Select
+    ThisWorkbook.Sheets("DrawSimulMap0").Copy Before:=ThisWorkbook.Sheets(11)
+    ThisWorkbook.Sheets("DrawSimulMap0 (2)").Select
+    ThisWorkbook.Sheets("DrawSimulMap0 (2)").Name = "DrawSimulMap"
+    ThisWorkbook.Sheets("DrawSimulMap0").Activate
+Application.DisplayAlerts = True
+
+ResultSheetName = "DrawSimulMap"
+
+ThisWorkbook.Sheets("SimulMap").Activate
+ThisWorkbook.Sheets("SimulMap").Range("B6", Range("B1048576").End(xlUp)).Select
+NbOfRowsMap = Selection.Rows.Count
+[b6].Select
+
+TotalNumberOfMarkers = NbOfRowsMap
+
+firstRow = 6
+LastRow = firstRow + TotalNumberOfMarkers - 1
+
+Read_Map_Drawing_Parameters
+MyMinLabelsSpacing = MyLabelsSpacing
+
+If EchelleVert = "" Or EchelleVert = 0 Then EchelleVert = Application.InputBox(prompt:="Drawing scale? (Suggestion: choose from 1 to 5)", Title:="Draw Map", Type:=1)
+If EchelleVert = False Then GoTo 1000
+Application.StatusBar = "Resetting window"
+
+Application.StatusBar = "Loading data"
+'VirtualRawData = ThisWorkbook.Sheets("SimulMap").UsedRange.Value
+VirtualRawData = ThisWorkbook.Sheets("SimulMap").Range(Cells(1, 1), Cells(firstRow + TotalNumberOfMarkers, 26)).Value
+ThisWorkbook.Sheets(ResultSheetName).Activate
+
+ChrSpace = MyChromosomeSpacing * 15
+SizeSmallSecondLine = EchelleHorizont
+LargCarte = 7
+
+Xorigin = 0
+Yorigin = 100
+
+MarkerChromosomeOld = 0
+FirstChromosome = ThisWorkbook.Sheets("DrawSimulMap").[b6]
+
+Application.StatusBar = "Drawing map..."
+
+For Row = firstRow To LastRow
+    MarkerChromosome = VirtualRawData(Row, 2)
+    If MarkerChromosome <> MarkerChromosomeOld Then ' new LG
+        Application.StatusBar = "Drawing LG " & MarkerChromosome
+        
+        For Row2 = Row To TotalNumberOfMarkers + firstRow 'look for chr size
+            If VirtualRawData(Row2, 2) <> MarkerChromosome Then
+                ChromosomeSize = VirtualRawData(Row2 - 1, 5)
+                NbOfMarkersInChr = Row2 - Row
+                GoTo 100
+            End If
+        Next Row2
+100
+        If LabelPositionsOption = 2 Then
+            MarkerLabelPositionOld = 0 - Yorigin
+        ElseIf LabelPositionsOption = 1 Then 'even mark label positions
+            MyLabelsSpacing = ChromosomeSize * EchelleVert / (NbOfMarkersInChr - 1)
+            If MyMinLabelsSpacing < MyLabelsSpacing Then MyLabelsSpacing = MyMinLabelsSpacing
+            MarkerLabelPositionOld = 0 - MyLabelsSpacing / EchelleVert
+        End If
+        
+        With ActiveSheet.Shapes.AddLine(Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 16, _
+            Yorigin + 5, _
+            Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 16, _
+            Yorigin + ChromosomeSize * EchelleVert + 5)
+            .Name = "MarkChr_" & MarkerChromosome
+        End With
+        With ActiveSheet.Shapes.AddShape(msoShapeRectangle, Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 16, Yorigin - 20, 100#, MyFontSize + 10)
+            .TextFrame.Characters.Text = "LG " & MarkerChromosome
+            .Fill.Visible = False
+            .Line.Visible = False
+            .Name = "MarkChrName_" & MarkerChromosome
+            .TextFrame.Characters.Font.Name = MyFontName
+            .TextFrame.Characters.Font.Size = MyFontSize + 5
+        End With
+        MarkerChromosomeOld = MarkerChromosome
+    End If
+    
+        MarkerName = VirtualRawData(Row, 4)
+        MarkerPosition = VirtualRawData(Row, 5) ' absolute positions
+        
+        If LabelPositionsOption = 2 Then
+            If MarkerPosition - MarkerLabelPositionOld < MyLabelsSpacing / EchelleVert Then 'if < 5 cm
+                MarkerLabelPosition = MarkerLabelPositionOld + MyLabelsSpacing / EchelleVert
+                MarkerLabelPositionOld = MarkerLabelPosition
+            Else
+                MarkerLabelPosition = MarkerPosition
+                MarkerLabelPositionOld = MarkerLabelPosition
+            End If
+        ElseIf LabelPositionsOption = 1 Then
+            If MarkerPosition - MarkerLabelPositionOld > MyLabelsSpacing Then
+                MarkerLabelPositionOld = MarkerLabelPositionOld + (MarkerPosition - MarkerLabelPositionOld) / 1.5
+            End If
+            MarkerLabelPosition = MarkerLabelPositionOld + MyLabelsSpacing / EchelleVert
+            MarkerLabelPositionOld = MarkerLabelPosition
+        End If
+        MarkerNumber = Row
+        
+        With ActiveSheet.Shapes.AddLine(Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 16, Yorigin + MarkerPosition * EchelleVert + 5, Xorigin + _
+            (MarkerChromosome - FirstChromosome) * ChrSpace - 9, Yorigin + MarkerPosition * EchelleVert + 5)
+            '.Name = "Markline_" & MarkerNumber
+        End With
+        With ActiveSheet.Shapes.AddLine(Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 9, _
+            Yorigin + MarkerPosition * EchelleVert + 5, _
+            Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace - 2 + SizeSmallSecondLine, _
+            Yorigin + MarkerLabelPosition * EchelleVert + 5)
+            '.Name = "Markline_" & MarkerNumber
+        End With
+        With ActiveSheet.Shapes.AddShape(msoShapeRectangle, _
+            Xorigin + (MarkerChromosome - FirstChromosome) * ChrSpace + SizeSmallSecondLine, _
+            Yorigin + MarkerLabelPosition * EchelleVert, _
+            120#, MyFontSize * 2)
+            .TextFrame.Characters.Text = MarkerName & " (" & Format(MarkerPosition, "###0.0") & ")"
+            .Fill.Visible = False
+            .Line.Visible = False
+            '.TextFrame.Characters.Font.Name = MyfontName
+            '.TextFrame.Characters.Font.Size = MyfontSize
+            .Name = "Mark_" & MarkerNumber
+        End With
+Next Row
+
+Dim MyTime2: MyTime2 = Time: TempsEcoule = MyTime2 - MyTime
+ActiveWindow.Zoom = 100
+'Application.CommandBars("Map Zoom").Visible = True
+ThisWorkbook.Sheets("DrawSimulMap").Activate
+ResetFontIndicator
+ResetFontSize
+1000 Application.StatusBar = False
+
+End Sub
+
+
+Public Sub Read_Chr_Simulated()
+Application.ScreenUpdating = False
+ThisWorkbook.Sheets("SimulMap").Activate
+
+Range("B6", Range("B1048576").End(xlUp)).Select
+NbOfRows = Selection.Rows.Count
+
+Range("J27:N1000").ClearContents
+
+NumbOfMarkInChrTotal = 0
+TotalMapSize = 0
+NumberOfChr = 0
+For i = 6 To NbOfRows + 6
+    If Cells(i + 1, 2) <> Cells(i, 2) Then 'new chr is detected
+        NumberOfChr = NumberOfChr + 1
+        Cells(26 + Cells(i, 2), 10) = Cells(i, 2) 'chr number
+        ChrSize = Cells(i, 5)
+        Cells(26 + Cells(i, 2), 11) = ChrSize 'chr size
+        TotalMapSize = TotalMapSize + Cells(i, 5)
+        NumbOfMarkInChr = i - 5 - NumbOfMarkInChrTotal
+        Cells(26 + Cells(i, 2), 12) = NumbOfMarkInChr
+        MarkerDensity = ChrSize / NumbOfMarkInChr
+        Cells(26 + Cells(i, 2), 14) = MarkerDensity
+        NumbOfMarkInChrTotal = NumbOfMarkInChrTotal + NumbOfMarkInChr
+    End If
+Next
+[N22] = TotalMapSize
+[N23] = NumbOfMarkInChrTotal
+[N24] = TotalMapSize / NumbOfMarkInChrTotal
+[R22] = NumberOfChr
+
+[N22].Select
+Application.ScreenUpdating = True
+End Sub
+
+Public Sub Add_Simulated_To_Map_1()
+Application.ScreenUpdating = False
+    Sheets("SimulMap").Activate
+        Range("B6:F6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Copy_Values_MD
+        Range("B6").Select
+    Sheets("CompMap").Activate
+        Range("B6:f1048576").Select
+        Selection.ClearContents
+        Range("B6").Select
+'        ActiveSheet.Paste
+        Paste_Values_MD
+        Range("B6").Select
+End Sub
+
+Public Sub Add_Simulated_To_Map_2()
+Application.ScreenUpdating = False
+    Sheets("SimulMap").Activate
+        Range("B6:F6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Copy_Values_MD
+        Range("B6").Select
+    Sheets("CompMap").Activate
+        Range("H6:L1048576").Select
+        Selection.ClearContents
+        Range("H6").Select
+'        ActiveSheet.Paste
+        Paste_Values_MD
+        Range("H6").Select
+End Sub
+
+
+Public Sub BatchMap()
+'Reads a series of mapmaker/EXP compatible data files, computes the map and reports a summary
+'Assumes that the linkage group are already defined and ordered.
+'Needs to have the sequences pasted in My sequences
+'Max 254 genotypes
+
+'On Error GoTo 900
+Dim MyTime
+MyTime = Time
+BatchMode = 1
+AutoBatchMapOn = 0
+
+Application.ScreenUpdating = False
+
+msg = "Are you sure you want to continue? This will clear the previous batch results and the current sequences."
+Style = vbYesNo
+Title = "Confirmation"
+Ctxt = 1000
+response = MsgBox(msg, Style, Title)
+If response = vbNo Then GoTo 1100
+
+If ThisWorkbook.Sheets("SimulMap").[b6] = "" Then MsgBox ("Please simulate a map first"): GoTo 1100
+
+ThisWorkbook.Sheets("Menus").Range("N24:IV1023").ClearContents
+ThisWorkbook.Sheets("Menus").Range("J24:J1023").ClearContents
+
+CopySequencesFromSimulatedMap
+
+nbseq = ThisWorkbook.Sheets("Menus").Cells(23, 19)
+If nbseq = 0 Then
+    MsgBox ("Please simulate a map first")
+    ActiverSimulate
+    GoTo 1100
+End If
+
+Application.StatusBar = "reading Menu"
+MenuMatrix = ThisWorkbook.Sheets("Menus").UsedRange.Value
+
+Read_Mapping_Options
+
+If ThisWorkbook.Sheets("Data").[D4] <> "" And ThisWorkbook.Sheets("Data").[c16] <> "" Then
+    msg = "Do you want to store the data that are in the 'Data' window before importation?"
+        Style = vbYesNo + vbCritical + vbDefaultButton1
+        Title = "Import Mapmaker data"
+        Ctxt = 1000
+        response = MsgBox(msg, Style, Title)
+        If response = vbYes Then Store_Data
+End If
+
+Application.StatusBar = "Erasing current data"
+ThisWorkbook.Sheets("Data").Range("B15:IV60000").ClearContents
+'ThisWorkbook.Sheets("Data+").Range("B15:IV60000").ClearContents
+
+NumberOfFilesToProcess = ThisWorkbook.Sheets("SimulMap").[n3]
+If NumberOfFilesToProcess > MaxNumberOfColumns - 7 Then BigBatchMode = 1 Else BigBatchMode = 0
+
+NbLocus = ThisWorkbook.Sheets("SimulMap").[N23]
+Indiv = ThisWorkbook.Sheets("SimulMap").[N17]
+
+myfilename = Application.GetOpenFilename(Title:="Locate the SimulPop_1.txt file")
+If myfilename = False Then GoTo 1100
+Workbooks.OpenText Filename:=myfilename, DataType:=xlDelimited, Tab:=True, Space:=True, ConsecutiveDelimiter:=True
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Range("A5:IV60000").ClearContents
+    ThisWorkbook.Sheets("BatchMap").[D4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").[f4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").Range("C8:G60000").Interior.ColorIndex = xlNone
+    ThisWorkbook.Sheets("BatchMap").Range("H5:IV6").Interior.ColorIndex = xlNone
+    ReDim BatchResults(NbLocus + 1, NumberOfFilesToProcess) As Variant
+Else
+    ThisWorkbook.Sheets("BigBatch").Range("C6:D60000").ClearContents
+    ReDim BigBatchResults(NumberOfFilesToProcess, 2) As Variant
+End If
+
+LengthOfPath = Len(myfilename)
+OutputPath = Left(myfilename, LengthOfPath - 5)
+myfilename = ActiveWorkbook.Name
+Workbooks(myfilename).Close savechanges:=False
+
+For CurrentPopulation = 1 To NumberOfFilesToProcess
+    'Application.StatusBar = "Population " & CurrentPopulation & ":" ' Importing data file
+    Import_Mapmaker_data_For_Batch
+    Application.StatusBar = "Population " & CurrentPopulation & ": Computing map"
+    ComputeMapForAutoBatch
+    'If BigBatchMode = 0 Then ThisWorkbook.Sheets("BatchMap").Cells(7, CurrentPopulation + 7) = "Pop" & CurrentPopulation
+    If BigBatchMode = 0 Then BatchResults(0, CurrentPopulation - 1) = "Pop" & CurrentPopulation
+Next
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+    ThisWorkbook.Sheets("BatchMap").[D4] = NumberOfFilesToProcess
+    NumberOfMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+    ThisWorkbook.Sheets("BatchMap").[f4] = NumberOfMarkers
+    [h7].Resize(NbLocus + 3, NumberOfFilesToProcess) = BatchResults
+    ThisWorkbook.Sheets("SimulMap").Activate
+        Range("D6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Selection.Copy
+        [N22].Select
+    ThisWorkbook.Sheets("BatchMap").Activate
+        Range("B8").Select
+        Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
+            :=False, Transpose:=False
+            
+    Compute_Stats_Batch
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+    [c6].Resize(NumberOfFilesToProcess, 2) = BigBatchResults
+End If
+
+Dim MyTime2
+MyTime2 = Time
+TempsEcoule = MyTime2 - MyTime
+Application.StatusBar = "Batchmap total time was " & Format(TempsEcoule, "h:mm:ss")
+
+GoTo 1000
+900 MsgBox "An error occurred during the BatchMap procedure."
+GoTo 1100
+1000
+BatchMode = 0
+'Application.StatusBar = False
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+End If
+
+Beep
+[h8].Select
+1100
+End Sub
+
+
+Public Sub Auto_Batch_Map()
+'Same than BatchMap but simulates populations at the same time so saves CPU time
+'Simulates a series of populations, computes the map and reports a summary
+'Assumes that the linkage groups are already defined and ordered.
+'Needs to have the sequences pasted in My sequences
+
+'On Error GoTo 900
+Dim MyTime
+MyTime = Time
+BatchMode = 1
+AutoBatchMapOn = 1
+
+Application.ScreenUpdating = False
+
+msg = "Are you sure you want to continue? This will clear the previous batch results and the current sequences."
+Style = vbYesNo
+Title = "Confirmation"
+Ctxt = 1000
+response = MsgBox(msg, Style, Title)
+If response = vbNo Then GoTo 1100
+
+If ThisWorkbook.Sheets("SimulMap").[b6] = "" Then MsgBox ("Please simulate a map first"): GoTo 1100
+
+ThisWorkbook.Sheets("Menus").Range("N24:IV1023").ClearContents
+ThisWorkbook.Sheets("Menus").Range("J24:J1023").ClearContents
+
+CopySequencesFromSimulatedMap
+
+nbseq = ThisWorkbook.Sheets("Menus").Cells(23, 19)
+If nbseq = 0 Then
+    MsgBox ("Please simulate a map first")
+    ActiverSimulate
+    GoTo 1100
+End If
+
+Application.StatusBar = "reading Menu "
+    MenuMatrix = ThisWorkbook.Sheets("Menus").UsedRange.Value
+Application.StatusBar = "reading simulation parameters"
+    VirtualMap = ThisWorkbook.Sheets("SimulMap").UsedRange.Value
+    Read_Simulation_Parameters
+    
+If SimulatedPopulationType <> "BC1" And SimulatedPopulationType <> "F2" And SimulatedPopulationType <> "BC2F1" And SimulatedPopulationType <> "SSD" Then MsgBox "BC1 or F2 populations only, sorry.": GoTo 1000
+
+Application.StatusBar = "reading mapping options"
+    Read_Mapping_Options
+
+NumberOfFilesToProcess = ThisWorkbook.Sheets("SimulMap").[n3]
+
+BigBatchMode = 0
+If MaxNumberOfColumns = 256 Then
+    If NumberOfFilesToProcess > 249 Then BigBatchMode = 1 Else BigBatchMode = 0
+End If
+
+NbLocus = ThisWorkbook.Sheets("SimulMap").[N23]
+Indiv = ThisWorkbook.Sheets("SimulMap").[N17]
+
+ReDim Data(Indiv, NbLocus) As Variant
+Application.StatusBar = "Clearing"
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Range("A5:IV60000").ClearContents
+    ThisWorkbook.Sheets("BatchMap").[D4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").[f4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").Range("C8:G60000").Interior.ColorIndex = xlNone
+    ThisWorkbook.Sheets("BatchMap").Range("H5:XFD60000").Interior.ColorIndex = xlNone
+    ReDim BatchResults(NbLocus + 1, NumberOfFilesToProcess) As Variant
+Else
+    ThisWorkbook.Sheets("BigBatch").Range("C6:D60000").ClearContents
+    ReDim BigBatchResults(NumberOfFilesToProcess, 2) As Variant
+End If
+
+For CurrentPopulation = 1 To NumberOfFilesToProcess
+    Application.StatusBar = "Processing population " & CurrentPopulation
+    Simulate_Population_From_Map
+    ComputeMapForAutoBatch
+    'If BigBatchMode = 0 Then ThisWorkbook.Sheets("BatchMap").Cells(7, CurrentPopulation + 7) = "Pop" & CurrentPopulation
+    If BigBatchMode = 0 Then BatchResults(0, CurrentPopulation - 1) = "Pop" & CurrentPopulation
+Next
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+    [D4] = NumberOfFilesToProcess
+    NumberOfMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+    [f4] = NumberOfMarkers
+    [h7].Resize(NbLocus + 3, NumberOfFilesToProcess) = BatchResults
+    ThisWorkbook.Sheets("SimulMap").Activate
+        Range("D6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Selection.Copy
+        [N22].Select
+    ThisWorkbook.Sheets("BatchMap").Activate
+        Range("B8").Select
+        Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
+            :=False, Transpose:=False
+    
+    Compute_Stats_Batch
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+    [c6].Resize(NumberOfFilesToProcess, 2) = BigBatchResults
+End If
+
+Dim MyTime2
+MyTime2 = Time
+TempsEcoule = MyTime2 - MyTime
+Application.StatusBar = "Autobatchmap total time was " & Format(TempsEcoule, "h:mm:ss")
+GoTo 1000
+900  MsgBox "An error occurred during the AutoBatchMap procedure."
+GoTo 1100
+1000
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+End If
+AutoBatchMapOn = 0
+BigBatchMode = 1
+'Application.StatusBar = False
+Beep
+[h8].Select
+1100
+End Sub
+
+
+Public Sub Auto_Batch_Map_Chi2()
+'Same than BatchMap but simulates populations at the same time so saves CPU time
+'Simulates a series of populations, computes the map and reports a summary
+'Assumes that the linkage groups are already defined and ordered.
+'Needs to have the sequences pasted in My sequences
+
+'On Error GoTo 900
+Dim MyTime
+MyTime = Time
+BatchMode = 1
+AutoBatchMapOn = 1
+
+Application.ScreenUpdating = False
+
+msg = "Are you sure you want to continue? This will clear the previous batch results and the current sequences."
+Style = vbYesNo
+Title = "Confirmation"
+Ctxt = 1000
+response = MsgBox(msg, Style, Title)
+If response = vbNo Then GoTo 1100
+
+If ThisWorkbook.Sheets("SimulMap").[b6] = "" Then MsgBox ("Please simulate a map first"): GoTo 1100
+
+ThisWorkbook.Sheets("Menus").Range("N24:IV1023").ClearContents
+ThisWorkbook.Sheets("Menus").Range("J24:J1023").ClearContents
+
+CopySequencesFromSimulatedMap
+
+nbseq = ThisWorkbook.Sheets("Menus").Cells(23, 19)
+If nbseq = 0 Then
+    MsgBox ("Please simulate a map first")
+    ActiverSimulate
+    GoTo 1100
+End If
+
+Application.StatusBar = "reading Menu "
+    MenuMatrix = ThisWorkbook.Sheets("Menus").UsedRange.Value
+Application.StatusBar = "reading simulation parameters"
+    VirtualMap = ThisWorkbook.Sheets("SimulMap").UsedRange.Value
+    Read_Simulation_Parameters
+    
+If SimulatedPopulationType <> "BC1" And SimulatedPopulationType <> "F2" And SimulatedPopulationType <> "BC2F1" And SimulatedPopulationType <> "SSD" Then MsgBox "BC1 or F2 populations only, sorry.": GoTo 1000
+
+Application.StatusBar = "reading mapping options"
+    Read_Mapping_Options
+
+NumberOfFilesToProcess = ThisWorkbook.Sheets("SimulMap").[n3]
+
+BigBatchMode = 0
+If MaxNumberOfColumns = 256 Then
+    If NumberOfFilesToProcess > 249 Then BigBatchMode = 1 Else BigBatchMode = 0
+End If
+
+NbLocus = ThisWorkbook.Sheets("SimulMap").[N23]
+Indiv = ThisWorkbook.Sheets("SimulMap").[N17]
+
+ReDim Data(Indiv, NbLocus) As Variant
+
+Application.StatusBar = "Clearing"
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Range("A5:IV60000").ClearContents
+    ThisWorkbook.Sheets("BatchMap").[D4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").[f4].ClearContents
+    ThisWorkbook.Sheets("BatchMap").Range("C5:XFD60000").Interior.ColorIndex = xlNone
+    'ThisWorkbook.Sheets("BatchMap").Range("H5:XFD7").Interior.ColorIndex = xlNone
+    ReDim BatchResults(NbLocus + 1, NumberOfFilesToProcess) As Variant
+Else
+    ThisWorkbook.Sheets("BigBatch").Range("C6:D60000").ClearContents
+    ReDim BigBatchResults(NumberOfFilesToProcess, 2) As Variant
+End If
+
+For CurrentPopulation = 1 To NumberOfFilesToProcess
+    Application.StatusBar = "Processing population " & CurrentPopulation
+    Simulate_Population_From_Map
+    'ComputeMapForAutoBatch
+    
+    ComputeChiSquareForAutoBatch
+    
+    'If BigBatchMode = 0 Then ThisWorkbook.Sheets("BatchMap").Cells(7, CurrentPopulation + 7) = "Pop" & CurrentPopulation
+    If BigBatchMode = 0 Then BatchResults(0, CurrentPopulation - 1) = "Pop" & CurrentPopulation
+Next
+
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+    [D4] = NumberOfFilesToProcess
+    NumberOfMarkers = ThisWorkbook.Sheets("SimulMap").[N23]
+    [f4] = NumberOfMarkers
+    
+    [h7].Resize(NbLocus + 3, NumberOfFilesToProcess) = BatchResults
+    
+    ThisWorkbook.Sheets("SimulMap").Activate
+        Range("D6").Select
+        Range(Selection, Selection.End(xlDown)).Select
+        Selection.Copy
+        [N22].Select
+    ThisWorkbook.Sheets("BatchMap").Activate
+        Range("B8").Select
+        Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
+            :=False, Transpose:=False
+    
+    Compute_Stats_Batch_Chi2
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+    [c6].Resize(NumberOfFilesToProcess, 2) = BigBatchResults
+End If
+
+Dim MyTime2
+MyTime2 = Time
+TempsEcoule = MyTime2 - MyTime
+Application.StatusBar = "Autobatchmap total time was " & Format(TempsEcoule, "h:mm:ss")
+GoTo 1000
+900  MsgBox "An error occurred during the AutoBatchMap procedure."
+GoTo 1100
+1000
+If BigBatchMode = 0 Then
+    ThisWorkbook.Sheets("BatchMap").Activate
+Else
+    ThisWorkbook.Sheets("BigBatch").Activate
+End If
+AutoBatchMapOn = 0
+BigBatchMode = 1
+'Application.StatusBar = False
+Beep
+[h8].Select
+1100
+End Sub
+