simulation of reduction of variance and Bulmer effect

As it is so difficult to understand reduction of variance due to selection and Bulmer effect I just simulated it in julia (code below). This is a rather crude simulation,  but it gives the following plot according to expectattions (the plot using Plotly is actually quite poor).

# simulation of Bulmer effect

using LinearAlgebra
using Distributions
using Random
#using Plotly
using PyPlot

Random.seed!(1234);

nqtl=300
nanim=10000
nselect=Int(nanim/20)
h2=0.8
#additive genetic variance is ~1
varp=1/h2
vare=varp-1
ngen=100

print("Base Population")
# a~ N(0,1)
a=randn(nqtl);

# p~ unif
p=rand(nqtl);
sum2pq=sum(2 .* p .* (1 .-p))
a=a/sqrt(sum2pq);

# storage for QTL genotypes
geno=Array{Int64,3}(undef,nanim,nqtl,2);
Z=Array{Int64,2}(undef,nanim,nqtl);

# simulate base animals
for i in 1:nanim
    for j in 1:nqtl
        for k in 1:2
            geno[i,j,k]=rand(Binomial(1,p[j]),1)[1]
        end
    end
end

vargen=Array{Float64,1}(undef,ngen)

for igen in 1:ngen

    println("-----------------------")
    println("generation number=", igen)
    println("-----------------------")
    for i in 1:nanim
        Z[i,:]=geno[i,:,1]+geno[i,:,2]
    end
    tbv=Z*a;
    println("var tbv= ",var(tbv))
    vargen[igen]=var(tbv)

    # generate phenotype with h2
    y=tbv+randn(nanim)*sqrt(vare);
    println("realized var y= ",var(y))
    println("realized h2= ",var(tbv)/var(y))

    # returns position in ranking based on y
    ranking=sortperm(y,rev=true);
    # pick best animals
    best=ranking[1:nselect];

    tbvselect=tbv[best];
    println("var tbvselect= ",var(tbvselect))

    # now create next generation sampling 2 alleles at random from 
    # selected individuals (possibly with selfing)
    genonew=Array{Int64,3}(undef,nanim,nqtl,2);
    for i in 1:nanim;
        for j in 1:nqtl
            for k in 1:2
                ancestor=rand(1:nselect);
                orig=rand(1:2);
                genonew[i,j,k]=geno[ancestor,j,orig];
            end
        end
    end
    geno[:,:,:]=genonew #not the same as geno=genonew !!!
end
plot(vargen)

handling dates in python

I have this data set with elementary test-day records from Manech Tete Noire also known as burubeltza

64422506;01;45000130047;052;03;2014;1;26FEB2014;0;02APR2014

64422506;01;45000130047;052;03;2014;1;26FEB2014;0;05MAY2014

64422506;01;45000130047;052;03;2014;1;26FEB2014;0;05JUN2014

64422506;01;45000130047;052;03;2015;2;07DEC2014;0;30DEC2014

and I need to compute difference between date of lambing and date of milk recording. I thought it would be easy using python(3).

import copy
import datetime

First, read and parse the data

fname=sys.argv[1]
fhand=open(fname)
for i,line in enumerate(fhand):
    Dacont = line.split(sep=";")[9]
    Datago = line.split(sep=";")[7]

Then, I need to convert these strings 26FEB2014 into something meaningful using module datetime, and compute the difference:

    daysfrommisebas= datetime.datetime.strptime(Dacont.lower(),'%d%b%Y') - datetime.datetime.strptime(Datago.lower(), '%d%b%Y')
and finally, print it, not forgetting to paste the strings using "+" and removing the EOL in "line":


    print(line.rstrip('\n')+";"+str(daysfrommisebas.days))

renumf90 anf preGSf90 with pedigree and genotypes only

Right, so you have genotypes and pedigree and you would like to use preGSf90 to build G and A22, or for quality control, or for whatever reason. Then you are told to use renumf90 first, but how can you use it with no phenotypes?

Well you can just create your own, fake ones. So, if your genotyped animals are in file marker.geno.clean with this aspect

  64000670990546 12012020210211

     45214790003 12111120211101

     45214790004 22111020110102

you can simply create a random number for each animal, and a file with "animal" and overall mean, using awk:

awk '{print 1,$1,rand()}' marker.geno.clean > phenotypes.txt

and then use the basic renumf90 parameter file in https://artadia.blogspot.com/2019/03/basic-renumf90-parameter-file.html.

After that, you edit the renf90.par file to put the required OPTIONS. For instance, these are the options to get lots of information, check heritability of gene content 

OPTION msg 100

OPTION h2_gene_content

basic renumf90 parameter file

DATAFILE
phenotypes.txt
TRAITS
3
WEIGHT(S)

RESIDUAL_VARIANCE  
 0.5 
EFFECT
1 cross alpha
EFFECT
2 cross alpha  #animal
RANDOM
animal
OPTIONAL

FILE
pedigree
FILE_POS
1 2 3 0 0
SNP_FILE
marker.geno.clean
PED_DEPTH
0
(CO)VARIANCES
    0.5 

How to split a large SNP file into individual chunks

To convert this (one animal per line, UGA format as explained here and here)

         1 012202001001000020202000

         2 012202001001000021211110

         3 012202001001000020202000

         4 012211001100010020202000

         5 111211101100010110111001

         6 002202000000000021211110

         7 012211001100010021211110

         8 022220002200020020202000

         9 012202001001000020202000

        10 111211101100010110111001

into as many individual SNP files as loci,e.g.

zcat singlemarker/x000001.gz | head

        1   1

        2   1

        3   1

        4   1

        5   2

        6   1

        7   1

        8   1

        9   1

       10   2

zcat singlemarker/x000002.gz | head

         1   2

         2   2

         3   2

         4   2

         5   2

         6   1

         7   2

         8   3

         9   2

        10   2

How to do it efficiently?

1-transpose the file (in my case via Fortran program) to:

0000100001

1111101211

2222122221

2222222222

0001101201

e.g. the first line is the first marker, and so on

2-use the extraordinary GNU split to split into files, one line (=one marker) at a time:

split -l 1 -d BB.700K.gen_transposed -a 6 --numeric-suffixes=1 --filter='gzip >$FILE.gz'

This command splits one line at a time (-l 1) creates a series of files with numeric suffixes starting in 1 (--numeric-suffixes=1) of width 6 (-a 6) like x000002.gz and finnaly "piped" through gzip to obtain compressed files (--filter='gzip >$FILE.gz')

Still the files look

0000100001

1111101211

2222122221

3-use  GNU fold to insert a newline after each character:

zcat  singlemarker/x000001.gz | fold -w 1 

format bash script

from https://linuxconfig.org/bash-printf-syntax-basics-with-examples

for i in $( seq 1 10 ); do printf "%03d\t" "$i"; done
001     002     003     004     005     006     007     008     009     010

submitting parallel jobs in genotoul / genologin with slurm

In another post I found out how to submit parallel jobs in Genotoul. Genotoul has a new cluster (genologin) and a new job scheduler: slurm, instead of "genotoul" with sge.

Here is how to submit a parallel job:

alegarra@genologin2 ~/work/mice_tpb $ cat run.sh

#!/bin/bash

#SBATCH --mem=1G

#SBATCH --cpus-per-task=10 

#SBATCH --time=10:20:00 

module load mpi/openmpi-2.1.2-intel2018.0.128

module load compiler/intel-2018.0.128

export OMP_NUM_THREADS=10

preGSf90 mice_G.par

Note the "module load" that before was not needed. This is because libraries are now mostly dynamic.

Edit look at: https://hpcrcf.atlassian.net/wiki/spaces/TCP/pages/7287288/How-to+Submit+an+OpenMP+Job