Estimating Allele Frequencies using ILINK

You need copies of an "ALLELE SPEED 0" version of the two programs "unknown" and "ilink". On helix, these programs are called unknown.as0 and ilink.as0. If you are not using helix, and are not set up to compile your own "ALLELE SPEED 0" versions of unknown and ilink, please contact BIMAS for assistance.
Select your pedigree and locus parameter (often called pedin.dat and datain.dat) files for your allele frequency estimation.
Determine the locus number for the locus you want to estimate the allele frequencies of (see below for an example).
Run the lcp program. Choose:
- general pedigrees
- ilink
- specific order
- no sex difference
- locus order 1 x (x = the number of the locus to estimate) e.g. 1 2, 1 5, 1 9, etc. wherever it is in your files
- control-z to write a pedin script
Edit the pedin script you just created - remove everything but the portion that runs the lsp program (see example below).
Run the pedin script that you just edited (i.e. run lsp).
Look in the pedfile.dat file that lsp just wrote. You are looking for the largest allele number that appears in the second set of 2 columns. Make a note of that number.
Edit the datafile.dat file that lsp just wrote:
- Change the line that has a 3 and another integer in it to look like: 3 y (where y is the largest allele number you just found in the above step). (see example below to help identify which line to change)
- Change the next line below that line: modify this line to contain y equal numbers that all add up to 1 (see table below)
- Change the next to last line to contain a 2.
- Change the last line: it should become a 1 followed by y more 1s (see example below). This tells ilink to estimate the allele frequencies for every allele.
Run unknown.as0 (or whatever you have called the ALLELE SPEED 0 version of unknown).
Run ilink.as0 (or whatever you have called the ALLELE SPEED 0 version of ilink).
Look in the final.dat file that ilink just wrote to see the estimated allele frequencies.

Example of using ILINK to estimate allele frequencies

Pedigree file "ch17test1.ped":

16  1  0  0 18  0  0 1 1  0 1  0 1  0 1      0 0    0 0    0 0    0 0    0 0    248  Ped: 16  Per: 1
16  2  0  0 18  0  0 2 0  2 6  2 6  2 6      0 0    0 0    0 0    0 0    0 0    250  Ped: 16  Per: 2
16  3 17 11  0  4  4 1 0  2 2  2 2  2 2      0 0    3 4    6 10   4 9    4 5    094  Ped: 16  Per: 3
16  4 17 11  0  5  5 1 0  2 3  2 3  2 3      0 0    3 4    6 8    1 9    5 6    095  Ped: 16  Per: 4
16  5 17 11  0  6  6 1 0  1 2  1 2  1 2      0 0    0 0    0 0    0 0    0 0    096  Ped: 16  Per: 5
16  6 17 11  0  7  7 2 0  2 3  2 3  2 3      0 0    1 4    6 10   4 6    4 6    144  Ped: 16  Per: 6
16  7 17 11  0  8  8 2 0  1 3  1 3  1 3      0 0    3 4    6 10   4 9    4 5    145  Ped: 16  Per: 7
16  8 17 11  0  9  9 2 0  1 3  1 3  1 3      0 0    3 3    6 8    1 6    6 6    146  Ped: 16  Per: 8
16  9 17 11  0 10 10 2 0  1 2  1 2  1 2      0 0    3 4    6 10   4 9    4 6    147  Ped: 16  Per: 9
16 10 17 11  0  0  0 2 0  1 2  1 2  1 2      0 0    3 4    6 10   4 9    4 5    148  Ped: 16  Per: 10
16 11 18 19  3 12 12 2 0  2 6  2 6  2 6      0 0    1 3    6 6    6 9    5 6    001  Ped: 16  Per: 11
16 12 18 19  0 13 13 1 0  1 6  1 6  2 6      0 0    1 3    6 6    6 9    6 6    004  Ped: 16  Per: 12
16 13 18 19  0 14 14 1 0  1 6  1 6  1 6      0 0    1 3    6 6    6 9    6 6    005  Ped: 16  Per: 13
16 14 18 19  0 15 15 1 0  1 6  1 6  1 6      0 0    3 4    4 6    9 9    4 5    006  Ped: 16  Per: 14
16 15 18 19  0 16 16 2 0  1 6  1 6  1 6      0 0    3 4    4 6    9 9    4 5    029  Ped: 16  Per: 15
16 16 18 19  0  0  0 2 0  1 5  1 5  1 5      0 0    3 4    4 6    9 9    4 5    030  Ped: 16  Per: 16
16 17  0  0  3  0  0 1 0  1 5  1 5  1 5      0 0    3 4    8 10   1 4    4 6    002  Ped: 16  Per: 17
16 18  1  2 11  0  0 1 0  1 6  1 6  1 6      0 0    0 0    0 0    0 0    0 0    244  Ped: 16  Per: 18
16 19  0  0 11  0  0 2 0  1 6  1 6  1 6      0 0    0 0    0 0    0 0    0 0    245  Ped: 16  Per: 19

Parameter file "ch17r395.par":>

8 0 0 5  << NO. OF LOCI, RISK LOCUS, SEXLINKED (IF 1) PROGRAM
0 0.0 0.0 0   << MUT LOCUS, MUT RATE, HAPLOTYPE FREQUENCIES (IF 1)
  1  2  3  4  5  6  7  8
1   2 << AFFECTION, NO. OF ALLELES
 0.200000 0.800000  << GENE FREQUENCIES
 6 << NO. OF LIABILITY CLASSES
 0.1700 0.0010 0.0010
 0.3400 0.0010 0.0010
 0.4800 0.0050 0.0050
 0.6400 0.0050 0.0050
 0.7600 0.0100 0.0100
 0.8500 0.0100 0.0100 << PENETRANCES (RECESSIVE)
1   2 << AFFECTION, NO. OF ALLELES
 0.200000 0.800000  << GENE FREQUENCIES
 6 << NO. OF LIABILITY CLASSES
 0.1700 0.0010 0.0010
 0.3400 0.0010 0.0010
 0.4800 0.0050 0.0050
 0.6400 0.0050 0.0050
 0.7600 0.0100 0.0100
 0.8500 0.0100 0.0100 << PENETRANCES (RECESSIVE)
1   2 << AFFECTION, NO. OF ALLELES
 0.200000 0.800000  << GENE FREQUENCIES
 6 << NO. OF LIABILITY CLASSES
 0.1700 0.0010 0.0010
 0.3400 0.0010 0.0010
 0.4800 0.0050 0.0050
 0.6400 0.0050 0.0050
 0.7600 0.0100 0.0100
 0.8500 0.0100 0.0100 << PENETRANCES (RECESSIVE)
3    2 << ALLELE #S, # OF ALLELES  PVU2.D18S21
 0.464882943 0.535117057 << GENE FREQS
3    4 << ALLELE #S, # OF ALLELES  PCR.D18S37
 0.019886364 0.134943182 0.589488636 0.255681818 << GENE FREQS
3   10 << ALLELE #S, # OF ALLELES  PCR.D18S53
 0.007042254 0.069014085 0.023943662 0.28028169 0.069014085 0.307042254 0.123943662 0.084507042 0.023943662 0.0112676040000001 << GENE FREQS
3   13 << ALLELE #S, # OF ALLELES  PCR.JO@72A
 0.045454545 0.039312039 0.066339066 0.084766585 0.071253071 0.097051597 0.100737101 0.076167076 0.243243243 0.128992629 0.030712531 0.00982801 0.00614250700000012 << GENE FREQS
3    9 << ALLELE #S, # OF ALLELES  PCR.JO@57
 0.102439024 0.02804878 0.031707317 0.140243902 0.108536585 0.323170732 0.193902439 0.062195122 0.00975609900000009 << GENE FREQS
0 0  << SEX DIFFERENCE, INTERFERENCE (IF 1 OR 2)
.1 .1 .1 .1 .1 .1 .1  << RECOMBINATION VALUES
1 0.10000 0.45000 << REC VARIED, INCREMENT, FINISHING VALUE

I happen to want to estimate alleles for locus 7.

I run lcp and select general pedigrees, ilink, specific order, no sex difference, locus order 1 7, and type control-z to write a pedin script file.

The pedin script that lcp just created looks like:


:
#
#	**
#	**  LCF  Generated  Command  File
#	**
#	**  Version : V4.0-00 ( 7-Jan-1991) - UNIX/Bourne
#	**  Created : 23-Apr-97  11:34:01
#	**
#

trap "
      if [ $1x != 'nodelete'x -a $1x != 'NODELETE'x -a $2x != 'nodelete'x -a $2x != 'NODELETE'x ]
      then
        rm -f datafile.dat final.dat ipedfile.dat lsp.log lsp.stm lsp.tmp\
              outfile.dat pedfile.dat recfile.dat speedfile.dat stream.dat tempdat.dat\
              tempped.dat *.tpd *.tdt *.mpd *.mdt voutfile.dat vstream.dat
      fi;
      exit
     " 0 1 2 3 15

rm -f lsp.log
rm -f lsp.stm
rm -f lsp.tmp
rm -f final.dat
rm -f stream.dat
rm -f vstream.dat
rm -f outfile.dat
rm -f voutfile.dat
rm -f recfile.dat

if [ -f final.out ]
then
  mv final.out final.bak
fi
cp /dev/null final.out
if [ $? != '0' ]
then
  exit 1
fi

if [ -f stream.out ]
then
  mv stream.out stream.bak
fi
cp /dev/null stream.out
if [ $? != '0' ]
then
  exit 1
fi

echo
echo
echo '******************************************************************************'
echo 'Run   1 - ILINK : p1 p7'
echo '******************************************************************************'
echo
cp ch18test1.ped ch18test1.tpd
cp ch18r395.par ch18r395.tdt
lsp \
	  ilink \
	  ch18test1.tpd ch18r395.tdt \
	  2 \
	   1  7 \
	  0 \
	  0 \
	  0.1
if [ $? = '0' -o $? = '1' ]
then
  cat lsp.log >> final.out
  cat lsp.stm >> stream.out
  unknown
  if [ $? = '0' ]
  then
    ilink
    if [ $? = '0' ]
    then
      cat final.dat  >> final.out
      cat stream.dat >> stream.out
    fi
  fi
fi

I edit the pedin file by removing most of it's lines to make it look like:


cp ch18test1.ped ch18test1.tpd
cp ch18r395.par ch18r395.tdt
lsp \
	  ilink \
	  ch18test1.tpd ch18r395.tdt \
	  2 \
	   1  7 \
	  0 \
	  0 \
	  0.1

I have removed everything before the "cp ch18test1.ped..." line, and everything after the 0.1, the last parameter into lsp.

The idea is that I only want to run the lsp program, and not do anything else that the pedin script contained.

I then run the pedin script by typing 'pedin'.

The pedin script runs, and creates these files:

ch18r395.tdt
ch18test1.tpd
datafile.dat
lsp.log
lsp.stm
pedfile.dat

I look in pedfile.dat, which looks like:

   16   1   0   0  18   0   0 1 1     0   1    0   0
   16   2   0   0  18   0   0 2 0     2   6    0   0
   16   3  17  11   0   4   4 1 0     2   2    4   9
   16   4  17  11   0   5   5 1 0     2   3    1   9
   16   5  17  11   0   6   6 1 0     1   2    0   0
   16   6  17  11   0   7   7 2 0     2   3    4   6
   16   7  17  11   0   8   8 2 0     1   3    4   9
   16   8  17  11   0   9   9 2 0     1   3    1   6
   16   9  17  11   0  10  10 2 0     1   2    4   9
   16  10  17  11   0   0   0 2 0     1   2    4   9
   16  11  18  19   3  12  12 2 0     2   6    6   9
   16  12  18  19   0  13  13 1 0     1   6    6   9
   16  13  18  19   0  14  14 1 0     1   6    6   9
   16  14  18  19   0  15  15 1 0     1   6    9   9
   16  15  18  19   0  16  16 2 0     1   6    9   9
   16  16  18  19   0   0   0 2 0     1   5    9   9
   16  17   0   0   3   0   0 1 0     1   5    1   4
   16  18   1   2  11   0   0 1 0     1   6    0   0
   16  19   0   0  11   0   0 2 0     1   6    0   0

I want to see what the largest number is in the last set of two columns. In this case, the largest number is a 9.

I now want to edit datafile.dat. The datafile.dat file starts out looking like this, as it was produced by lsp:

2 0 0 3
0 0.00000000 0.00000000 0
 1 2
 
1 2
 0.20000000 0.80000000
6
 0.17000000 0.00100000 0.00100000
 0.34000000 0.00100000 0.00100000
 0.48000000 0.00500000 0.00500000
 0.64000000 0.00500000 0.00500000
 0.76000000 0.01000000 0.01000000
 0.85000000 0.01000000 0.01000000
 
3 13
 0.04545454 0.03931204 0.06633907 0.08476659 0.07125307 0.09705160 0.10073710 0.07616708 0.24324324 0.128
99263 0.03071253 0.00982801 0.00614251
 
0 0
 0.10000000
0
 1

I want to change the line that says "3 13" to say "3 9", 9 being the largest number that I identified in the above step.

I also want to change the line below that line to have 9 instances of equal numbers. I consult the below table to see what 9 numbers that are equal add up to 1.

I also want to change the next-to-last line to contain a 2 instead of a 0.

And I want to change the last line to be a 1 with nine 1's after it.

The resulting file looks like:

2 0 0 3
0 0.00000000 0.00000000 0
 1 2
 
1 2
 0.20000000 0.80000000
6
 0.17000000 0.00100000 0.00100000
 0.34000000 0.00100000 0.00100000
 0.48000000 0.00500000 0.00500000
 0.64000000 0.00500000 0.00500000
 0.76000000 0.01000000 0.01000000
 0.85000000 0.01000000 0.01000000
 
3 9 
0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.12
0 0
 0.10000000
2
 1 1 1 1 1 1 1 1 1 1

I now run unknown.as0. I get this as output:

Program UNKNOWN version 5.20
The following maximum values are in effect:
      20 loci
     325 single locus genotypes
      25 alleles at a single locus
    2000 individuals in one pedigree
       3 marriage(s) for one male
       3 quantitative factor(s) at a single locus
      60 liability classes
      25 binary codes at a single locus
       8 maximum number of loops
Opening DATAFILE.DAT
YOU ARE USING LINKAGE (V5.20) WITH  2-POINT
YOU ARE USING FASTLINK (V3.0P)
 AUTOSOMAL DATA
Opening PEDFILE.DAT
Ped. 16

I then run ilink.as0. It runs for a little while and produces this output:

YOU ARE USING LINKAGE (V5.1 (1-Feb-1991)) WITH  2-POINT
YOU ARE USING FASTLINK (V3.0P (29-Sep-1995)) AUTOSOMAL DATA
 
Program ILINK version  5.10 (1-Feb-1991)
 
 
FASTLINK (slow) version 3.0P (29-Sep-1995)
 
The program constants are set to the following maxima:
     8 loci in mapping problem
    25 alleles at a single locus (maxall)
 50000 maximum of censoring array (maxcensor)
  1000 individuals in all pedigrees combined
   150 pedigrees (maxped)
     3 quantitative factor(s) at a single locus
    15 liability classes
    25 binary codes at a single locus
    3.00 base scaling factor for likelihood (scale)
    2.00 scale multiplier for each locus (scalemult)
 0.00000 frequency for elimination of heterozygotes (minfreq)
 
 
See README.allele
FASTLINK will renumber alleles for you automagically
ALLELE DIAGNOSIS: pedigree 16 does not use all alleles at locus index 2
 
Number of alleles at locus 1 is 2
Number of alleles at locus 2 is 9
The number of haplotypes used is 18
Maxcensor can be reduced to       -32767
ITERATION     1 T =      0.100 NFE =    10 F =  7.97493e+01
ITERATION     2 T =      0.002 NFE =    21 F =  7.31730e+01
ITERATION     3 T =      0.000 NFE =    32 F =  7.31607e+01
ITERATION     4 T =      0.000 NFE =    43 F =  7.31607e+01
ITERATION     5 T =      0.024 NFE =    53 F =  7.05223e+01
ITERATION     6 T =      0.009 NFE =    63 F =  7.02004e+01
ITERATION     7 T =      1.000 NFE =    73 F =  6.81615e+01
ITERATION     8 T =      0.975 NFE =    83 F =  6.71899e+01
ITERATION     9 T =      1.000 NFE =    93 F =  6.68298e+01
ITERATION    10 T =      1.000 NFE =   103 F =  6.67557e+01
ITERATION    11 T =      1.000 NFE =   113 F =  6.67420e+01
ITERATION    12 T =      2.000 NFE =   123 F =  6.67355e+01
ITERATION    13 T =      1.000 NFE =   133 F =  6.67296e+01
ITERATION    14 T =      1.000 NFE =   143 F =  6.67289e+01

I then look in the file final.dat and see this:

CHROMOSOME ORDER OF LOCI : 
  1  2
****************** FINAL VALUES **********************
PROVIDED FOR LOCUS   2 (CHROMOSOME ORDER)
******************************************************
GENE FREQUENCIES :
 0.167370 0.000691 0.000691 0.167007 0.000691 0.175027 0.000691 0.000691 0.487144
******************************************************
THETAS:
 0.500
******************************************************
-2 LN(LIKE) =  6.67289e+01
LOD SCORE =-6.35973e-09
NUMBER OF ITERATIONS =    14
NUMBER OF FUNCTION EVALUATIONS =   143
PTG = -1.62539e-04
******************************************************
******************************************************

So my estimated allele frequencies are: 0.167370 0.000691 0.000691 0.167007 0.000691 0.175027 0.000691 0.000691 0.487144

I'm done!

A table of (almost) equal numbers that add up to 1

 2: 0.50 0.50

 3: 0.33 0.33 0.34

 4: 0.25 0.25 0.25 0.25

 5: 0.20 0.20 0.20 0.20 0.20

 6: 0.16 0.16 0.17 0.17 0.17 0.17 

 7: 0.14 0.14 0.14 0.14 0.14 0.15 0.15

 8: 0.12 0.12 0.12 0.12 0.13 0.13 0.13 0.13

 9: 0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.12

10: 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

11: 0.09 0.09 0.09 0.09 0.09 0.09 0.09 0.09 0.09 0.09 0.10

12: 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.09 0.09 0.09 0.09

13: 0.07 0.07 0.07 0.07 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08

14: 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.08 0.08

15: 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07

16: 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.07

17: 0.05 0.05 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06

18: 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06

19: 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.06 0.06 0.06 0.06

20: 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05