||||
|\/|
|/\|
||||

David L. Duffy, MBBS PhD.
Queensland Institute of Medical Research,
300 Herston Road,
Herston, Queensland 4029, Australia.
Email: davidD@qimr.edu.au

CONTENTS

• Introduction
• Methods
• Usage
• Datasets
• Tips and tricks
• Documentation of routines
• Limitations
• References
• Licence
• Program history
• Appendix: Embedded Scheme

INTRODUCTION

Program Sib-pair performs a number of simple analyses of family data that tend to be "nonparametric" or "robust" in nature. It is modelled to some extent on the Genetic Analysis System [Young, 1995] in terms of the command language and types of analysis. Included are routines for:

• Imputation of genotypes.
• Estimation of allele frequencies in codominant genetic systems.
• Simple and complex segregation analysis of a binary trait.
• Estimation of familial correlations and sibship variances for a quantitative trait. Variance components analysis of quantitative and binary traits using a variety of likelihoods.
• Haseman-Elston sib-pair regression of a quantitative (or binary) trait using full and half-sib data, and variance components linkage analysis for normally distributed quantitative traits.
• Carrying out multiple versions of the transmission-disequilibrium test.
• Testing allelic association with a binary or quantitative trait -- Monte Carlo simulation of null distribution of simple tests, or "measured genotype" familial analysis including combined association and linkage analysis.
• Single locus Affected Pedigree Method identity-by-state and identity-by-descent linkage analysis. This includes Wards [1993] extensions to include unaffected pedigree members.
• Writing out of locus and pedigree files in the formats used by the programs APM, Arlequin, ASPEX, CRIMAP, FISHER, GAS, GDA, Genehunter, LINKAGE, LOKI, MENDEL, MERLIN, PAP and SAGE.

An example of use

Sib-pair is command line oriented, and writes output to the standard output (the screen if you are using it interactively). Output can be saved to a file by diverting it to a file using the "out" command (in which case you will no longer see it coming to the screen). Alternatively, you may use another program that can copy from the standard output, such as tee, emacs (the powerful editor usually found on Unix systems), syn (a powerful Windows editor) or the Tcl/Tk based GUI program isp which can found on the same site as Sib-pair itself.

Bearing this in mind, here is a sample interactive session. We start at the command line of your operating system shell:

> sib-pair

|||| SIB-PAIR: A program for simple genetic analysis
|\/| Version : Version 1.00.beta (12-Sep-2008)
|/\| Author : David L Duffy (c) 1995-2008
|||| Job run : Tue Sep 16 12:25:34 2008 (moonboom)

Type "help" for help, "quit" to quit, "ctrl-C" to interrupt.


>>

A double arrow command prompt appears. We read in a previously prepared script:

>> include ex.in

The contents of ex.in are a series of Sib-pair commands:

# declare four loci
set loc a affection
set loc b quantitative
set loc m1 marker 0.0 cM
set loc m2 marker 5.1 cM
# read the pedigree data
read ped inline
ex1 1a x x m n x 1 3 1 2
ex1 1b x x f n x 1 2 3 4
ex1 2a 1a 1b m n 3.5 1 2 1 3
ex1 2b x x f n 1.1 2 2 2 3
ex1 3a 2a 2b m y 4.3 1 2 1 2
ex1 3b 2a 2b m n 2.0 2 2 2 3
ex1 3c 2a 2b f n 0.8 2 2 3 3
ex1 4a 3c 3d f y x 1 2 2 3
ex1 3d x x m n x x x x x
ex1 4b 3b 3e m y 4.7 1 2 3 4
ex1 4c 3b 3f m n 1.6 2 2 1 3
;;;;
# The four semicolons ends the in-line data
run

The "run" command actually starts the initial processing of the pedigree.

NOTE:  Imputation level 1.  Imputing untyped parental
        genotypes where unequivocal.

Pedigree file        = inline.ped                                                                                                                                      
Number of loci       =     4

Locus      Type Position
---------- ---- --------
a           a     6
b           q     7
m1          m     8--  9
m2          m    10-- 11

Number of marker loci=    2
Bonferroni corr. 5%  =    0.025321
Bonferroni corr. 1%  =    0.005013
Bonferroni corr. 0.1%=    0.000500

NOTE:  Creating dummy record for ex1-3e.
NOTE:  Creating dummy record for ex1-3f.

NOTE: Father and mother of person ex1-4a appear to be swapped around.  Reordering.
NOTE:  Person ex1-3e appears as a mother and sex was unspecified.
       Setting sex to female.
NOTE:  Person ex1-3f appears as a mother and sex was unspecified.
       Setting sex to female.
Pedigree: ex1        No. members:   13 No. founders:    6 No. sibships:    5

Total number of pedigrees  =     1
Number with only 1 member  =     0
Total number of sibships   =     5
Total number of subjects   =    13
Total subjects genotyped   =    10 ( 76.9%)
Total number of genotypes  =     20
Largest pedigree (members) =    13 (Pedigree ex1)
Deepest pedigree (genrtns) =     4 (Pedigree ex1)

Mean size of pedigrees     =    13.0
Mean pedigree depth        =     4.0
Mean size where >1 members =    13.0
Mean depth where >1 geners =     4.0
Number of imputed genotypes=     0

Number of pedigree errors  =     0
Number of deleted records  =     0

We obtain a few routine messages and summary statistics. The small table of Bonferroni corrections is a reminder that Sib-pair does not usually correct P-values for multiple tests; it is up to the user to decide what the appropriate significance thresholds are.

Blank records are created for named but missing parents. This is necessary, as a formal pedigree should have neither or both parents present for each person.

The "ls" commands shows trait loci (with an asterisk appended), and marker loci.

The "drop" commands drops loci from the scope of subsequent commands, while the "undrop" command returns access to all loci.

------------------------------------------------
Allele frequencies for locus "m1        "
------------------------------------------------
   Allele  Frequency   Count  Histogram
       1     0.3000        6  ******
       2     0.6500       13  *************
       3     0.0500        1  *

Number of alleles    =    3
Heterozygosity (Hu)  =    0.5105
Poly. Inf. Content   =    0.4064
Number persons typed =   10 ( 76.9%)

Marker         NAll  Allele(s)    Freq   Het    Ntyped
-------------- ---- -----------   ------ ------ ------
m1                3    1 ..   3    -     0.5105    10
m2                4    1 ..   4    -     0.7211    10

The "describe" command gives summary information about loci. For marker loci, it tabulates simple allele counts and proportions in the dataset. Given the keyword "snp", it gives a summary for all markers, one line per locus. For traits, "describe" gives familial correlations or recurrence risks. The "tabulate" gives simpler tables:

>> tab a

a          x:   2      y:   3      n:   8

>> tab a m1

------------------------------------------------
Cross-tabulation of "a         " ... "m1        "
------------------------------------------------

                      m1        
a                 1/2          1/3          2/2   
           ---------------------------------------
     n         2 (.286)     1 (.143)     4 (.571)
     y         3 (1.00)     0 (.000)     0 (.000)

    No. complete observations =    10
    LR contingency chi-square =   5.5
           Degrees of freedom =   2
                      P-value =0.0643

There are a few other simple descriptive commands. The "count" command gives information about families and individuals:

Count where "b > 1":
 
Pedigree   Con=T   Num  ASPs Trios    4+
---------- ----- ----- ----- ----- -----
ex1            6    13     1     0     0
Total          6    13     1     0     0

The "select" command selects pedigrees containing a specified number (or one or more) individuals meeting the criterion. The "print" command is individual oriented:

Print where "b > 1":

ped=ex1 id=2b fa=x mo=x  sex=f b=1.1000 m1=2/2 m2=2/3
ped=ex1 id=4c fa=3b mo=3f sex=m b=1.6000 m1=2/2 m2=1/3
ped=ex1 id=4b fa=3b mo=3e sex=m b=4.7000 m1=1/2 m2=3/4
ped=ex1 id=3a fa=2a mo=2b sex=m b=4.3000 m1=1/2 m2=1/2
ped=ex1 id=3b fa=2a mo=2b sex=m b=2.0000 m1=2/2 m2=2/3
ped=ex1 id=2a fa=1a mo=1b sex=m b=3.5000 m1=1/2 m2=1/3

Number of matched persons   =    6 out of    13 ( 46.2%)
Number of matched pedigrees =    1 out of     1 (100.0%)

If we had instead issued:

we would obtain:

Print where "b > 1":

ped=ex1 id=2b fa=x mo=x  sex=f m1=2/2 m2=2/3
ped=ex1 id=4c fa=3b mo=3f sex=m m1=2/2 m2=1/3
...

As of 2006-Mar-01, we can write expressions containing genotypes:

Print where "m2 = = 3/4":

ped=ex1 id=1b fa=x mo=x  sex=f a=n b=x m1=1/2 m2=3/4
ped=ex1 id=4b fa=3b mo=3e sex=m a=y b=4.7000 m1=1/2 m2=3/4

The "associate" command gives results from family based tests of allelic association for a trait versus all active marker loci:

--------------------------------------------------
Allelic association testing for trait "a         "
--------------------------------------------------

Marker     Typed  Allels Chi-square Asy P  Emp P  Iters
---------- ------ ------ ---------- ------ ------ ------
m1             10      3        1.9 0.3930 0.1575    127 AssX2-HWE .  
m1              1      3        1.1 0.5978 0.5978     23 RC-TDT    .  
m2             10      4        0.9 0.8192 0.7692     26 AssX2-HWE .  
m2              1      4        2.1 0.4471 0.4471    160 RC-TDT    .  

The "sibpair" command gives results from regression based tests of linkage for a trait versus all active marker loci. The P-values for these tests can be "empirically" estimated by gene-dropping marker alleles under the null hypothesis of no linkage:

---------------------------------------------------------
Sham S+D H-E for trait "b         " v. all markers
---------------------------------------------------------

Marker     FSibs  HSibs  t-value    Asy P  Emp P  Iters
---------- ------ ------ ---------- ------ ------ ------
m1              3      1        2.6 0.1180 0.0288    695 H-E +  
m2              3      1        1.5 0.1908 0.2128     94 H-E .  

Finally, we log transform the quantitative trait values and write out a MERLIN type pedigree file, so we can further examine our data using a multipoint program:

If we like, Merlin can be run from within Sib-pair:

The "$" command shells out to run another program. When we exit from Merlin, we will be returned to Sib-pair:

This job took  1.8 minutes

There are a number of operations useful for manipulating pedigree data prior to analysis. These allow you to prune ("prune") pedigrees down to selected individuals and only those relatives needed to connect the index people, to reduce families to unrelated cases and controls ("case"), to break up large pedigrees into component nuclear families ("nuclear") or into unrelated cliques ("subped"), if the pedigree file does not specify all the connecting relatives (between different branches say). One can also select ("select" and "unselect") particular groups of pedigrees, and specify different values of a variable in each group:

Sib-pair also can be used as a calculator, and has a few genetics utilities, notably the "sml" and "grr" commands which gives expected recurrence risks, ibd's and genotype frequencies for a specified diallelic model:

>> sml 0.01 0.5 0.2 0.1

------------------------------------------------
Single Major Locus Recurrence Risk Calculation
------------------------------------------------

 Frequency(A): 0.010000; Pen(AA): 0.500; Pen(AB): 0.200; Pen(BB): 0.100
 Trait Prev  : 0.102020; Pop AR:   2.0%; Var(Add): 0.000206; Var(Dom): 0.000004

 Measure      MZ Twin       Sib-Sib        Par-Off      Second    
----------   ----------    ----------     ---------    ----------
 Rec Risk         0.104         0.103         0.103         0.103
 Rel Risk         1.023         1.011         1.011         1.006
 Odds Rat         1.025         1.012         1.012         1.006
 PRR              1.020         1.010         1.010         1.005
 Tet Corr         0.007         0.003         0.003         0.002
 ibd|A-A          1.000         0.502         0.500         0.251
 ibd|A-U          1.000         0.500         0.500         0.250

 Freq of A if Affected: 0.019898 (0.000,0.039,0.961) 
 Freq of A if Unaffctd: 0.008875 (0.000,0.018,0.982)

 Mating       Proportion    Risk to offspring
----------   -----------   ------------------ 
UnA x UnA         0.806         0.102
Aff x UnA         0.183         0.103
Aff x Aff         0.010         0.104

Finally, most commands can be interrupted by typing "ctrl-C" (that is holding the control and C keys down simultaneously). This returns you to the Sib-pair prompt, after finishing the current command as quickly as possible. To interrupt the program completely, you need to strike "ctrl-C" six times.

METHODS

Analytic Methods

Imputation of unobserved genotypes. This is performed using the algorithm described by Lange & Goradia [1987]. Firstly, (0) A phenoset (all possible genotypes) for one locus is generated for each individual in the pedigree. Then, iterate by nuclear families, repeating the next two steps until no further updates: (1) Parental genotypes inconsistent with their offspring are removed; (2) child genotypes inconsistent with their parents are removed. Finally, (3) If zero genotypes remain, report an inconsistency; if one genotype remains, this becomes the imputed genotype; if the joint spouse genotype is unambiguous, but the specific genotype each spouse carries is ambiguous, if requested, randomly assign a genotype to each parent. These latter genotypes might be used only for calculation of statistics for offspring of the pair, but not for the parents themselves. A further extension is to sequentially (founders then nonfounders) impute the remaining missing genotypes as the most likely member of the phenoset. This or a faster randomised algorithm is always run (unless the imputation flag is set to "-1", see below) to give starting values for the MCMC methods, but these simulated genotypes will not be saved unless the imputation level is set to "3" or "full".

Sex imputation. Likelihood of observed homozygosity at multiple sex-linked loci is calculated under hypothesis of male and female sex, assuming 0.1% of male and female genotypes are miscalled as heterozygotes.

Allele frequencies. These are for codominant systems only. Either a straight allele count is used, or the contribution of each pedigree is weighted by the number of founders it contains. Alternatively, the imputed and observed genotypes in the founders can be counted, or the MLE of the founder allele frequencies calculated by MCMC (see below), or the BLUE of the founder allele frequencies via the approach of McPeek et al [2004].

Hardy-Weinberg proportions. These are tested by a Pearson chi-square, with the P-value estimated via "gene-dropping" (Monte-Carlo,MC) simulation. These are based on the genotypes in the founders of that pedigree, where typed, or on the gene frequencies in the total sample where the founder is untyped, and the structure of the pedigree. For diallelic markers, the exact test of Hardy-Weinberg equilibrium (assuming all individuals is unrelated) is also calculated.

Binary trait correlation. The ANOVA type estimate of the within-pedigree intraclass correlation is calculated, along with Tarone's score test [1979] for extra-binomial variation (testing whether r>0).

Quantitative and categorical trait correlation. The pairwise (all-possible-pairs) type estimates of the intraclass and interclass correlations are calculated for MZ twins, parent-offspring, full-sib, half-sib, grandparent, avuncular, cousin and marital pairs. In the case of categorical traits, this correlation is the (unweighted) Cohen kappa. The polychoric correlation (a faster two-stage estimator) is available for ordinal traits. Standard errors are estimated using a random subset delete-d jackknife method [Shao and Tu 1995]. The draw size is max(1, min(10, Nobs/10)), and iter pseudosamples are evaluated. The standard error of the test statistic is:

JSE=((n-d)/dm * Sum(r_i-mean(r_i, i=1, m))²)^1/2

where d is the number of observations dropped for each pseudosample, r_i is the test statistic value for the ith pseudosample, m is the number of pseudosamples, and n is the number of observations [Shao and Tu 1995].

Segregation ratios. The default analysis assumes the pedigrees are unascertained, and gives the naive estimators. The ascertainment corrected analysis is for nuclear families and follows Davie [1979]: p=(r-j)/(t-j), where p is the risk, r is the total number of affected children, j, the number of sibships containing exactly one proband, t, the number of children. A fairly efficient approximate sampling variance is also given.

Haplotype reconstruction. This routine is not currently available in the Fortran 95 version of Sib-pair. This is performed on a nuclear family by family basis, though incorporating grandparental information where available. Initial reconstruction is performed using a simulated annealing algorithm that maximizes a sharing based criterion based on length of runs of the same alleles on a putative chromosome among sib pairs, parent-offspring pairs, and grandparent-child pairs. The order of loci in the pedigree file is treated as the linkage order, and map distance information is not used. Missing parental haplotypes are filled in using the childrens' haplotypes in a simple fashion. Mendelian inconsistencies are flagged in the printout.

The second algorithm is more ambitious and attempts to construct recombination minimized haplotypes, again on a nuclear family by family basis, and dealing more intelligently with missing data. A simulated annealing algorithm using multiple restarts is used. Recombination events and Mendelian errors are flagged in the output.

Admixture analysis. This refers to testing for a mixture of specified distributions in the empirical distribution of a quantitative trait, usually a mixture of normals, though Sib-pair also offers mixtures of exponential and Poisson distributions. Information from the relationship between family members is not utilised. The usual EM approach is used.

Test for normality. The Filliben correlation [Filliben 1975] is calculated as a test for normality. This is the correlation between the observed data and the rankits expected under the assumption of normality. The P-value for this statistic is approximate, and is produced using an approach modelled on that used by Royston [1993] for the related Shapiro-Francia W':

log(1-r_f) ~ N(m,s)
m:=1.0402 (log(log(n))-log(n)) - 1.99196
s:=0.788392/log(n)) + 0.31293

where n is the sample size. This performs reasonably well versus the empirical percentiles:

	Coverage
n	5	10	20	50	100	500	1000	5000
Nominal P=0.05	0.021	0.044	0.055	0.057	0.059	0.052	0.045	0.033
Nominal P=0.01	0.00	0.006	0.009	0.014	0.019	0.007	0.007	0.007

A common use of the Filliben correlation is as the criterion for selecting an optimal transformation of the data.

The J_0.02 statistic is a variance-corrected order-statistic based skewness measure:

[(P_0.02+P_0.98)/2-P_0.50]/[P_0.75-P_0.25]

[David & Johnson 1956]. A test of normality can be constructed, using the standard error of J_0.02. This is based on interpolation of results from Monte-Carlo simulations (SE(J_0.02):=1.36/sqrt(N) cf Resek [1974]).

Sibship variance test. This is the linear model suggested by Fain [1977] for the detection of the phenotypic effects of quantitative trait loci. Briefly, if parental trait values are at the extreme of the population distribution, then they will be carrying multiple increasing or decreasing alleles at the QTLs. As a result, the trait variance among their offspring is decreased compared to sibships whose parents have trait values close to the population mean. This U-shaped relationship between midparent value and sibship variance can be detected by fitting linear or quadratic curves.

Variance components analysis. This is the usual mixed effects analysis of quantitative traits assuming multivariate normality. The log-likelihood for a quantitative trait:

LL = -0.5 [log(det(S)) + (y-u)' inv(S) (y-u)]

is maximized, where S is the variance-covariance matrix for the trait values for each phenotyped pedigree member, and y and u are the trait values and their expected values:

u= B'X,

reducing to the grand mean in the absence of covariates. The main diagonal elements of S takes the value:

V_A+V_Q+V_D+V_E,

and the off-diagonal elements:

R_i,jV_A+ibd_i,jV_Q+K_i,jV_D,

where,

R_i,j is the coefficient of relationship for the i-jth pair of relatives
K is the coefficient of fraternity
ibd is the average ibd sharing at the marker location being tested for linkage to a QTL.

The maximization is performed using AS319 (variable metric minimizer with numerically estimated gradients) or BOBYQA (bound optimization by quadratic approximation). The phenometric models (ADE, AE, E) are fitted to the intact entire pedigree, but the models including V_Q can be fitted to the entire pedigree, or to sibships only.

Standard errors for the fixed effects coefficients are the usual GLS estimates,

(X' S^-1 X)^-1 .

For binary traits, one can fit the above likelihood, which often gives quite adequate parameter estimates but dubious hypothesis testing, or move to the probit-normal mixed model, using the Multifactorial Threshold Model type formulation.

The log-likelihood to be maximized is then:

LL = log[cdf_MVN(u,S)]

where u is now the threshold value, and the probability is integrated from u to +Inf when y=1, and to -Inf when y=0. Integration of the multivariate normal is carried out using either the code provided by Alan Genz [Genz 1992] or the Mendell-Elston approximation code in TOMS717 [Bunch et al 1993]. The former uses quasi-Monte Carlo integration, and is limited to less than 500 dimensions (individuals per pedigree). Because this is a Monte Carlo algorithm, the maximizer can have problems settling on an exact solution. The Mendell-Elston approach is deterministic and faster, but less accurate (biased).

Categorical trait association analysis. This is the Pearson goodness-of-fit based test for equality of allelic gene frequencies at a marker locus in individuals expressing different trait values (RxC table). Both the nominal (ignoring relatedness of the sample) and empirical P-values for the test are output. The empiric P-value is estimated via gene-dropping simulation. These are based on the gene frequencies in the total sample, and the structure of the pedigree, and are conditional on the observed occurrence of the binary trait in the families. The gene-dropping can also be further conditioned on identity by descent at another marker (which may in fact be the marker being tested for association). This therefore gives a "model-free" test of association conditional on linkage.

This can be compared to results from the quasi-likelihood score methods (WQLS and MQLS) of Bourgain et al [2003], and Thornton and McPeek [2007]. These use the kinship matrix to adjust for the resulting score for familial correlations in marker genotype (and trait phenotype by using this as the predictor). These methods are implemented as per those papers for analysis of binary traits:

T = (p_a-p_e) (SD_diff)^-1

p = Σ(w_ig_i)/Σ(w_i)

g_i is half the allele count for the ith individual.

For the WQLS,

wa =	A-1a
we =	A-11
Vardiff =	0.5p(1-p) (aA-1a) (aA-11)-2 - (1A-11)-1

where A is the numerator relationship matrix, a the phenotype vector, and 1 a unit vector.

It is relatively straightforward to extend the WQLS approach to multiple categories and to quantitative traits. In the former case, the score contains terms for the multiple indicators for different levels of the trait variable.

Formerly, a sibship permutation based test was provided, calculating the same chi-square statistic for members of sibships that contain at least one affected and one unaffected typed individual, and generating a (within-sibship) permutation P-value. This is now replaced by a more powerful score (FBAT) test combining the within-sibship association and transmission-disequilibrium tests after Knapp [1999] and Laird et al [2000]. In this approach, the complete or partial genotype of untyped parents is reconstructed from the genotypes of the affected and unaffected children. The transmission of alleles from these reconstructed parents is conditioned on the genotypes of the children used in the reconstruction. The appropriate conditional distribution under the null hypothesis is approximated via Monte Carlo simulation and rejection sampling. For binary traits, the sibship permutation test in still available in the form of a conditional logistic regression stratifying on sibship. Finally, it is also possible to construct permutation based tests by using the permute or get commands.

In the case of multiple generation pedigrees, Sib-pair scans the pedigree upwards nuclear family by nuclear family, saving imputed genotypes, so that these contribute as children to subsequent nuclear families. I don't think this will be biased, but have yet to check. This behaviour can be turned off.

Population genetic F statistics (F_IS, F_IT and F_ST) are also calculated for each marker assuming individuals with different indicator trait values come from separate related demes. These are estimated following the approach of Pons and Chaouche [1995], as described by Excoffier [2001]. Results for each locus are combined to give multilocus summaries for the sample as:

F*IS = (Σ HS - Σ HO)/Σ HS

F*IT = (Σ HT - Σ HO)/Σ HT

F*ST = (Σ HT - Σ HS)/Σ HT

Estimation of empirical (marker identity by state based) kinship uses an EM algorithm described by Choi et al [2009].

Multidimensional scaling (MDS) analysis of interindividual genetic distances is performed on the mean IBS sharing matrix for all pairs of individuals. The algorithm is Classical MDS, and so uses EISPACK routines [Smith et al 1976] to extract the first D eigenvectors from the doubly-centred distance matrix (1-mean(IBS)).

Quantitative trait association analysis. This fits an additive (allelic means) model predicting an individual's trait value from his/her genotype at a marker locus. The residual sum-of-squares is compared to those obtained via a gene-dropping simulation of the pedigrees, giving an empiric P-value.

Two-locus and N-locus linkage disequilibrium estimation. One algorithm finds informative founder matings, or informative matings where all the grandparents are untyped, and imputes the two-locus haplotypes transmitted to the offspring. The loci are assumed to be tightly linked, so that four parental haplotypes are counted, as opposed to the more usual two haplotypes from the child. This is the default for markers where there are more than 6 alleles per marker. Both D and D' measures are calculated.

For two locus analyses, the second algorithm combines phased and unphased genotype data in an EM fitted log-linear model. This can handle X-linked loci. If the markers are diallelic, then the model is solved directly as a cubic equation (note that for speed, the resulting chi-square is the Pearson chi-square, based on the r²). For more than two markers, and for testing haplotype association to a trait, phase information is discarded (ie all individuals are treated as unrelated), and only autosomal markers can be used.

Homozygosity analysis. This tests for an increase in observed homozygosity at a marker locus in individuals expressing a binary trait, comparing this to the predicted homozygosity based on the allele frequencies in the total sample. This may occur in the presence of allelic association with a recessive trait locus, and/or deletional loss of heterozygosity (the parents would be untyped for such an individual not to have been flagged as a Mendelian inconsistency of course). A one-tailed empiric P-value is estimated via "gene-dropping" simulation, based on the gene frequencies in the total sample, and the structure of the pedigree. The multipoint homozygosity analysis uses the mean maximum marker homozygosity run length in the set of cases. Again, gene dropping is used to produce the distribution of this statistic under the null given the marker map, allele frequencies and observed pedigrees.

The runs-of-homozygosity estimate of the inbreeding coefficient can also be calculated. The sum of the lengths (bp) of all runs exceeding a minimum threshold (default 1.5 Mbp) is divided by the total genome map length.

Transmission-disequilibrium test. The original formulation of the TDT is for a diallelic marker [Spielman et al 1993]. The TDT statistic calculated by the program is the Pearson goodness-of-fit based test of symmetry in the square table of transmitted versus nontransmitted alleles to each affected child [Haberman, 1979]. Empiric P-values are produced by randomization of the table. This global allelic form does not correct for the (usually small) correlation in parental genotypes induced by linkage disequilibrium (absent of course under the null hypothesis). Another allelic test provided is the marginal allelic test suggested by Spielman and Ewens [1996], which is probably slightly more powerful than the global symmetry test [Kaplan et al 1997].

The genotypic TDT P-value is estimated via gene-dropping based on the genotypes of typed ancestors of probands (where both parents of the proband and all antecedents must be typed), and the structure of the pedigree. The test statistic compares the observed number of each genotype transmitted with the number expected based on the parental genotypes. Pairs of cells whose total count is less than a given cutoff may be excluded from the analysis to increase power. The P-value for the TDT testing each allele versus all others in turn ("allele-by-allele") is the exact two-sided binomial probability (via the beta distribution). When the plevel is zero, only the best of the allele-by-allele test results are printed, and the P-value is Bonferroni corrected for the number of alleles at that marker.

Note that the default option is to use probands for the TDT only one parent is typed. For the diallelic marker case with unequal allele frequencies, using one parent families does lead to biased results. The unified transmission test available via the "assoc" command does not suffer from this problem (see above).

The Schaid and Sommer [1993] genotypic risk ratio tests for familial association under the assumption of Hardy-Weinberg equilibrium, or conditional on parental genotypes, is also offered. This uses log-linear modelling (implemented as iteratively reweighted least squares) of a biallelic locus with both parents genotyped. The attributable risk is also produced.

Haplotype Relative Risk analysis. This is the original familial association statistic comparing transmitted and nontransmitted allele frequencies unmatched on family (Falk and Rubinstein 1987; Knapp et al 1993). As usual, an empiric P-value is estimated via "gene-dropping" simulation, based on the gene frequencies in the total sample, and the structure of the pedigree.

Sib-pair analysis. Identity by descent estimation is based on the sib pair and parental genotypes when available. In the case of untyped parents, the full-sib sharing is the sum of sharing for each possible set of parental genotypes weighted by their likelihood based on all children in the sibship. Half-sib sharing is estimated based only on known genotypes, whether observed or unequivocally imputed. The effective degrees of freedom for the t-test of the slope of the regression line is given as the number of individuals in the sample (counted once only) who are in a sib pair where both members are typed at trait and marker, minus two. For a sample made up of nuclear families (no halfsibs), this will be equivalent to the Σ(sibship size-1)-2 value used by SIBPAL 2.6, and originally suggested by Hodge [1984]. For binary traits, the same ordinary-least-squares analysis is performed -- the t-statistics from these results are only really applicable to large samples, and tend to be too liberal. The quantity regressed is not the usual Haseman-Elston squared trait difference, but a function of the squared trait sums and differences following Sham and Purcell [2001]. This approach is supposed to approach the power of the variance components approach according to those authors, and gives appropriate Type 1 error rates.

For the Fulker & Cardon methods, the expected ibd through the interval between the two markers is estimated using the equation given in Olson [1995]. Haseman-Elston regressions are performed at a series of points across the interval using the ibd sharing of the two flanking markers, and the given size of the interval. The Haldane mapping function is used.

Affected sib pair analysis. This is the original IBS based approach described by Lange [1986a], extended to half-sibs as per Bishop and Williamson [1990]. No correction for sibship size is made -- that is all possible pairs are treated as independent. The usual two d.f. chi-square is calculated, with expected counts being calculated based on the observed gene frequencies in the total sample. The IBD based mean test is also calculated.

Affected Pedigree Method. This uses the measures of genetic similarity described by Weeks and Lange [1988; see also, Lange, 1986a, 1986b], Whittemore and Halpern [1994], and Ward [1993, 1995]. The expected mean and variance for each pedigree is estimated via gene-dropping simulation. These are based on the observed gene frequencies in the total sample, and the structure of the pedigree. Both ibs and ibd based family scores can be estimated. The original APM statistics, the APM statistic of Whittemore and Halpern and the T(AB) and GPM statistics of Ward are calculated.

Multilocus IBS sharing statistics. These are used to confirm the pedigree structure using marker data. One approach calculates the overall probability of sharing two alleles IBS for full and half sibs, summing over all loci, and ignoring any linkage between markers. The second approach uses gene dropping based on the given marker map. Two lists are generated, one by individual, the other by relative pair.

Martingale residuals. The elegant approach of Commenges [1994] to genetic analysis of age-at-onset is to analyse the residuals obtained from a nonparametric or semiparametric survival analysis. Sib-pair implements the former, calculating the martingale residuals using the Nelson-Aalen estimator for the integrated hazard [eg Andersen et al 1993], which are then transformed following Therneau et al [1990] to give a more symmetrical distribution.

Generalized linear models and survival regression. These are the usual IRLS algorithms (using AS 164, [Stirling 1981]). The exponential, Weibull and Extreme Value Distribution regressions are implemented as Poisson regressions (with log time as offset) as per Aitken and Clayton [1980]. Conditional logistic regression calls the original logccs subroutine of Breslow and coworkers [Smith et al 1981].

Bivariate survival analysis. This is a ranks-based approach, using Kendall's tau as the measure of association. Both the original estimator [Oakes 1982] and the "normalized" estimator suggested by Oakes [2008] to get around high censoring are used. The within pair correlations in age at onset are calulated for twin and sib pairs. A bootstrap is used to estimate the standard error for tau, and the zygosity labels for the pairs are permuted to give a test of the MZ:DZ difference in correlation. This currently assume pairs are independent.

Jonckheere-Terpstra trend test. This is currently for a marker versus a quantitative trait, where the genotypes are assumed to be naturally ordered eg diallelic loci. An empiric P-value can either be calculated by permutation of the quantitative phenotype, or by gene dropping.

Stratified proportional odds test of Whitehead and Whitehead. This is a test [Whitehead and Whitehead 1991] in the spirit of the Cochran-Armitage-Berry approach, and comes in fixed and random effects flavours, the latter following Der Simonian and Laird (1987). It can be applied to binary and ordinal data.

Other standard statistical tests. A number of classical tests for independent data (eg unrelated cases) are implemented, such as contingency chi-square test (with Monte Carlo "exact" P-values, see below), ordinal by ordinal trend test for contingency tables [Yates 1948], and nonparametric one way analysis of variance via the Kruskal-Wallis test. Sib-pair can calculate the Pearson correlation coefficient for between-trait association, the Kaplan-Meier estimator for the survivor function, and Nelson-Aalen estimator for the hazard function, and measures of agreement for contingency tables.

Deterministic Estimation of Pedigree Likelihoods

Rather than the original Elston-Stewart algorithm [Elston and Stewart 1971] augmented by pedigree traversal, Sib-pair uses the iterative peeling approach described by Wang et al [1996] (following Janss et al [1992]). This has the advantage of "automatically" dealing with loops, although the resulting likelihood in that case is only approximate: this approximation can be improved by several methods, but these are not currently implemented. The algorithm follows the detailed description in Chapter 2.1 of Schelling [2004], though at present only allows evaluation at one or two codominant loci.

Briefly, the iterative approach peels up and down simultaneously by calculating anterior and posterior values for each individual, where the anterior and posterior values represent the scaled likelihood contributions for ancestors and siblings, and mates and descendants respectively. The values for any individual rely on those for the other relatives, so these are reciprocally updated over multiple iterations until they converge. Usually a maximum of 10-20 iterations suffices, and if an ideal ordering of evaluations is used, a single iteration in unlooped pedigrees is sufficient ie it is equivalent to the usual pedigree traversal algorithms eg Lange and Boehnke [1983]; Wang et al [1996], Fernandez and Fernando [2001].

Monte Carlo Algorithms

Gene-dropping. The "unconditional" algorithm producing null distributions is as follows. Repeat the following 1-3 steps a large number of times. (1) Founder genotypes are assigned using the allele frequencies in the observed sample, assuming panmixia and Hardy- Weinberg equilibrium (HWE). Iterate, until all genotypes are assigned: (2) If both parental genotypes are nonmissing, randomly assign the index a genotype based on Mendelian autosomal inheritance (ie if parental genotypes are {1/2} and {3/4}, a child's genotype is randomly selected from {{1/3}, {1/4}, {2/3}, {2/4}}, with each genotype having a probability of selection of 0.25). Once complete, (3) calculate the test statistic based on the family's simulated genotype. Following completion of the outer loop, (4) summarize the distribution of the resulting test statistic. This procedure is used to generate null distributions for the association Pearson chi-square.

For the share command, this also allows for recombination between markers based on the given linkage map

The null distributions for the genotypic marginal TDT is generated using a "conditional on founders" algorithm, that takes observed founder/ancestor genotypes as given. Only typed nonfounders genotypes where both parents were typed are simulated.

ibd estimation. The modification to calculate ibd distributions gives each founder (two) unique alleles in his/her "typing genotype". A simple gene-drop gives the null distributions for the ibs and ibd statistics of the APM method.

I based the "conditional on observed genotypes" (gene-drop with rejection sampling) algorithm for calculating ibd on that described by Blangero et al [1995]. As before, each founder is assigned a typing genotype made up of two unique alleles. Offspring are only assigned ibd-typing genotypes that are consistent with the observed genotype at the observed locus of interest. For example, say the observed genotypes in the parents are 100/102 and 100/102, and the typing genotypes associated with these are set to {1/2} and {3/4} respectively. If the child is 100/102, assignment of typing genotypes {1/3} and {2/4} will be rejected. This is equivalent to a child's genotype being randomly selected from {{1/4}, {2/3}}, with each genotype having a probability of selection of 0.5. The resulting ibd statistics based on the typing genotype will approximate ibd for the marker locus of interest.

Multipoint ibd estimation. This is a simplification of the Cardon-Fulker approach as extended by Almasy and Blangero [1998]. The ibd estimates from sets of markers in complete linkage are combined by calculating the variance-weighted mean ibd for each pair of relatives over the set. The estimates of the ibd variances for each marker arise as a byproduct of the MCMC algorithm.

Gene-dropping conditional on ibd. The previous algorithms can be combined by gene dropping a marker conditional on "ibd indicator" genotypes that represent segregation at either that same locus, or a neighbouring more informative marker. This allows us to simulate the null distribution assuming linkage but no association.

Missing genotype simulation by Monte-Carlo Markov Chain (MCMC). The calculation of ibd in the presence of missing genotypes is performed via a Metropolis algorithm. This algorithm is a multiallelic extension of that described by Lange and Matthysse [1989]. One iteration of the generation of a legal constellation of imputed and observed genotypes is produced by:

(1) Perform (a) (b) (c) or (d):

(a) Simulate ibd, then "mutating" up to four imputed founder alleles. These propagate through the pedigree using the current pattern of ibd transmission as indicated by the ibd-typing alleles, and are rejected and resimulated if an (unordered) inconsistency with an observed genotype occurs.

(b) Simulate ibd, then switch the parent of origin for an individual heterozygote. Propagate this change up through the pedigree to the originating founder(s), but not below the chosen "pivot" individual.

(c) A "conditional on observed genotypes" dropping of ibd-typing alleles, with the refinement that this ignores imputed genotypes. Inconsistencies thus generated for imputed genotypes are resolved by changing the imputed genotypes. This procedure will be slow in the presence of a large number of untyped nonfounders.

For (a)-(c), additional local proposals (as below) are compounded to these, and the resulting proposed constellation is accepted or rejected via the Metropolis criterion.

(d) Resimulate all untyped x untyped founder mating joint genotypes conditional on their offspring and other spouses, then other pedigree members singly, again conditional on surrounding genotypes. This is a simple Gibbs sampler, and is more efficient than the above when there are many missing genotypes in larger pedigrees.

MCMC burn-in. In releases prior to version 0.96.0, there was no burn-in for this Metropolis algorithm, as preliminary empirical tests had found the results from this program agreed well with "exact" results from programs such as GENEHUNTER. Subsequently, I have found some pedigrees where using the starting genotypes from the Lange-Goradia approach does lead to biased ibd estimates for certain pairs of relatives. Therefore, the program now performs a number of burn-in iterations (default 100) prior to those used to estimate ibd. The required number of such iterations depends on the number of missing genotypes in the pedigree.

Metropolis generalized linear mixed model and finite polygenic model sampler. This is either a "standard" or "slice" Metropolis sampler, where the simulated variables include diallelic QTL genotypes, Gaussian breeding values, a single QTL allele frequency (shared by all QTLs in the FPM), up to three genotypic means (shared by all QTLs in the FPM), polygenic and environmental variances (including pedigree ("VC") and maternally-derived sibship ("VS") variances).

The trait model can be gaussian, binomial (with identity, probit or logit link), poisson (including log link), weibull or MFT.

Proposals for diallelic QTLs genotypes are straightforward to generate, and do not usually give rise to noncommunication between sets of legal genotype proposals. Proposals for continuous variables are generated from random normal deviates, and a tuning parameter can be set that alters the variance of these proposal distributions.

The likelihood contribution from the ith individual to the Metropolis criterion for these models is (see for example, Guo and Thompson [1993]):

LL = F*Σ(log(P(G_j))) + F*log(f(a|V_A)) + (1-F)*log(f(a|a_FA, a_MO)) + log(c|V_A) + I*log(f(y|G₁,...G_j, a, c, V_E))

where,
P(x) denotes the probability of x,
f(x) denotes the density of x,
y is the trait value,
a is the breeding value,
c is the pedigree-specific intercept,
G_j is the genotype at the jth QTL,
V_A is the additive polygenic variance,
V_C is the familial environmental variance,
V_E is the error variance,
F=1 when a founder, 0 when a nonfounder
I=1 when phenotype observed, 0 when unobserved.

The conditional density for the breeding values of offspring includes the correction for inbreeding (the segregation variance being 1-0.5*(F_FA+F_MO). The random effects are modelled as zero-mean gaussian.

In the case of the logistic-normal GLMM, the heritability is estimated as:

h² = V_A/(V_A+pi²/3).

The realizations of the parameters are summarized as means, and approximate standard errors produced by batching (default B=iter^1/2 [Jones et al 2005]). The interbatch lag-1 serial correlation is calculated as a diagnostic for the appropriate number of values to simulate [Ripley 1987].

The implementation of the generalized linear mixed models is quite straightforward in the chosen Metropolis paradigm (it would be more work to produce a Gibbs sampler, I believe), but for the "standard" sampler, it is important to check that the proposal acceptance rates are in the optimum range (usually stated as 0.2-0.6, Ripley 1987). This is less critical for the slice sampler, where the tabulated acceptance rates are actually the ratio of accepted proposals to the number of function evaluations (and so are just a measure of algorithm efficiency). Models fitting V_C and V_S or V_A are two-level GLMMs and so I have fitted a number of test datasets from the literature. There are surprising differences between results from standard software for some of these datasets, so although Sib-pair sometimes does not give identical results to that from non-simulation-based maximum likelihood methods, this may reflect approximations used by other programs.

Increasing the number of random effects chains is realized by duplicating the families the appropriate number of times and correcting the likelihood and standard errors. One is essentially averaging over multiple estimates of the random effects for each individual, as global parameters such the fixed effects regression coefficients and overall variances are the same over the replicate chains at that iteration. This seems to reduce bias in the estimation of the random effects, but with the side effect of increasing the between-batch correlation, and so slowing estimation. The tabulated results below generally used 4 chains run for 10000 iterations after a 1000 iteration burnin.

Binomial GLMM analysis of seed germination dataset of Crowder et al (1978) using different approaches. PQL1 is the penalised quasilikelihood approach implemented as glmmPQL() in the MASS package [Venables and Ripley 2002], while PQL2, AGQ are results from lmer() in the lme4 package of Bates and Sarkar [2005] using penalized quasilikelihood, adaptive Gaussian Quadrature respectively. The BUGS results are from the BUGS Examples manual.

	Parameter Estimate (SE)
Method	Sib-pair	AGQ	PQL1	PQL2	BUGS
SD of Plate Effect	0.28 (0.17)	0.24 (0.09)	0.23	0.24	0.29 (0.15)
Intercept	-0.51 (0.13)	-0.54 (0.17)	-0.54 (0.17)	-0.54 (0.16)	-0.51
Seed	0.06 (0.32)	0.10 (0.28)	0.09 (0.27)	-0.09 (0.28)
Root	1.31 (0.26)	1.33 (0.24)	1.32 (0.23)	1.32 (0.24)
Seed x Root	-0.79 (0.44)	-0.81 (0.38)	-0.81 (0.38)	-0.81 (0.38)

Binomial GLMM analysis of "bacteria" dataset from the R MASS package [Venables and Ripley 2002] using 5 different approaches. PQL1 is the penalised quasilikelihood approach implemented as glmmPQL() by Ripley in the MASS package, while PQL2, AGQ and Laplace are results from lmer() in the lme4 package of Bates and Sarkar [2005] using penalized quasilikelihood, adaptive Gaussian Quadrature and the Laplace approximation respectively.

	Parameter Estimate (SE)
Method	Sib-pair	AGQ	PQL1	PQL2	Laplace
RE Variance	2.05 (1.23)	1.70 (1.05)	1.98	3.27	1.66
Intercept	3.70 (0.76)	2.86 (0.48)	2.74 (0.38)	2.75 (0.48)	2.81 (0.48)
Low dose	-1.46 (0.74)	-1.36 (0.82)	-1.25 (0.64)	-1.25 (0.82)	-1.35 (0.82)
High dose	0.60 (0.74)	0.58 (0.85)	0.49 (0.67)	0.49 (0.85)	0.58 (0.85)
Week>2	-1.66 (0.51)	-1.63 (0.46)	-1.61 (0.36)	-1.61 (0.46)	-1.57 (0.46)

Binomial GLMM analysis of contraception usage data from the 1988 Bangladesh Fertility Survey [Steele et al 1996].

	Parameter Estimate (SE)
Method	Sib-pair	PQL
RE Variance	0.25 (0.06)	0.22
Intercept	-1.67 (0.16)	-1.66 (0.15)
Age	-0.03 (0.01)	-0.03 (0.01)
Urban	0.72 (0.07)	0.72 (0.12)
1 child	1.10 (0.12)	1.09 (0.16)
2 children	1.36 (0.15)	1.35 (0.17)
3+ children	1.32 (0.17)	1.32 (0.18)

Poisson GLMM analysis of epileptic seizure count data of Thall and Vail [1990] using different approaches. PQL1 is the penalised quasilikelihood approach implemented as glmmPQL() in the MASS package [Ripley and Venables 2002], while PQL2, AGQ are results from lmer() in the lme4 package of Bates and Sarkar (2005).

	Parameter Estimate (SE)
Method	Sib-pair	AGQ	PQL1	PQL2
RE variance	0.268 (0.101)	0.252	0.197	0.101
Intercept	1.834 (0.112)	1.833 (0.074)	1.870 (0.106)	1.870 (0.074)
Progabide	-0.346 (0.152)	-0.334 (0.105)	-0.310 (0.149)	-0.309 (0.105)
log basal rate	0.861 (0.119)	0.883 (0.091)	0.882 (0.129)	0.882 (0.091)
Base:therapy	0.394 (0.157)	0.339 (0.143)	0.342 (0.203)	0.342 (0.143)
log age	0.513 (0.330)	0.481 (0.244)	0.534 (0.346)	0.533 (0.244)
Period 4	-0.159 (0.054)	-0.160 (0.055)	-0.160 (0.077)	-0.160 (0.143)

Poisson GLMM analysis of European male melanoma death rate dataset of Langford et al (1998) using different approaches. PQL1 is the penalised quasilikelihood approach implemented as glmmPQL() by Ripley and Venables [2002] in the MASS package, while PQL2, AGQ are results from lmer() in the lme4 package of Bates and Sarkar (2005). The STATA result used the xtpois command, and comes from the review article at http://www.mlwin.com/softrev/revstata.html.

	Parameter Estimate (SE)
Method	Sib-pair	AGQ	PQL1	PQL2	STATA
Region variance	0.188 (0.037)	0.170 (-)	0.161	0.125	0.102 (-)
Intercept	-0.151 (0.058)	-0.139 (0.043)	-0.129 (0.049)	-0.129 (0.043)	-0.138 (0.017)
UVB insolation	-0.035 (0.011)	-0.034 (0.009)	-0.038 (0.010)	-0.038 (0.009)	-0.056 (0.004)

The final results are sometimes sensitive to the choice of starting values for the random effects (the fixed effects are started automatically from the marginal model parameter estimates using "reg"), and to the proposal step size. Because of the correlation between random and fixed effects in GLMM's other than Gaussian (since the intercept affects variance), differences in the estimated random effect size do alter the fixed effects regression coefficients.

The output (with plevel set to 1) allows plotting of parameter estimates to assess convergence of the chain.

Multiple genotype imputation. This is one approach to dealing with association in families where trait data is available for some individuals without marker genotype data. The genotype sampler is used to simulate missing marker genotypes conditional only on observed marker genotypes in the pedigree and the population allele frequencies for that marker (which can be fixed to a specified value). It is not, unfortunately, simulated conditional on values at a given trait, so it assumes data is missing-completely-at-random (MCAR). Multiple cycles of imputation followed by association analysis of the augmented data are carried out, then averaged results are produced. The estimate of the mean effect (with standard error) for an allele is produced following Rubin [1987].

Randomized TDT. The randomization test for the global allelic TDT permutes the transmission table by randomly selecting a single proband-parent pair and reversing the transmitted and nontransitted alleles. One "shuffle" of the table involves N such permutations, where N is the number of such informative parent-proband pairs in the observed pedigrees (this reduces the correlation between successive tables in the random walk to close to zero).

Sequential empiric P-values. The Monte-Carlo P-values provided for the various MC-based tests are produced using the sequential approach described by Besag and Clifford [1991]. In this refinement, we only generate as many pseudosamples as is necessary to give a P-value numerator of size mincount; the denominator is the number of pseudosamples. The practical effect of this procedure is that if the true P-value is large, then relatively few pseudosamples are generated to give a less precise estimate of this uninteresting value. Besag and Clifford suggest a value for mincount of 10-20. It is necessary to set a maximum denominator to avoid excessive computation for "highly significant" results. For cases where none of the simulated test statistics exceeds the observed statistic, I provide an estimate of the P-value based on extreme value theory following Davis & Resnick [1984]. If the observed statistic is sufficiently large, it is in the tail of the distribution function which will tend to one of the classical extreme value distributions, the Pareto if the tail index is finite. The tail index a can be estimated following Hill [1975] as:

a = m^-1 Σ (log(X_(i))-log(X_(m+1)),

where X_(i) is the m'th order statistic. The estimate of the tail is then:

P = (m/n) (x/X_(m+1))^1/a.

This estimator seems to be conservative, but not as conservative as the previous estimate 1/(n+1).

Other empiric P-values. An exception to this is the algorithm used for empirical P-values for the APM. Here, a P-value for each family is simulated at the same time as the mean and variance. These P-values were previously combined using the procedure due to Fisher [Hedges and Olkin 1985], that is, twice the sum of the natural logarithms of the P-values was treated as a chi-square variate with 2*N degrees of freedom, where N is the number of contributing families. This does not seem to be particularly powerful, so each P-value is now inverse-normal transformed to a Z-score, and these combined in an unweighted fashion [Hedges and Olkin 1985].

Shortest paths. A shortest path through a pedigree between two individuals is generated using Djikstra's algorithm. The presence of multiple paths is not flagged, and the path is not purely genetic relationships.

Pedigree loops. This is a depth-first search with backtracking. Currently prints only one loop per pedigree.

USAGE

The program reads commands from standard input, and writes results to standard output. Therefore, the program can be run interactively, or if a series of commands is to be found in a file, in batch mode. If the input file was test.in, entering "sib-pair <test.in >test.out" would perform the commands in test.in, and write results to test.out.

A command is a single line of keywords, locus names and/or variable values. If a "\" character is the last word of a line, the next line is interpreted as a continuation of the previous command. Sib-pair is case-sensitive, so that the keyword "READ" is not equivalent to "read". Commands are either global, which can be entered at any time; descriptive (set impute, set locus, read pedigree), which must precede the run statement; the run statement, that causes the dataset to be read and processed; or analytic, which act only after the run statement.

One command, set plevel, controls verbosity of output. Some useful descriptive tables are only printed if plevel is at least 1.

Sib-pair's parser can evaluate simple algebraic and logical expressions for each record in a datafile, but does not directly allow complex programming. It does offer simple macro facilities including a loop construct. The eval command accesses a Scheme interpreter that can carry out more sophisticated computing.

There are command line options that may save a few keystrokes. Running "sib-pair -i test.in" starts Sib-pair and "includes" the contents of test.in. The command "sib-pair -l test.in" starts Sib-pair and runs the "locus" command on the contents of test.in. The command "sib-pair -h" (or --help) starts Sib-pair and runs the help command. If a file with suffix ".in" ".bin", or ".bin.gz" is the first argument, then it will automatically be included, while the ".bin", or ".bin.gz" is taken to indicate the file is a Sib-pair binary format pedigree file, and Sib-pair will attempt to read this into memory using the "read binary" command. This behaviour can be used to associate Sib-pair with that file type in a graphical directory browser, such as the Microsoft Windows Explorer.

If the program has been compiled with g95, the --g95 flag prints useful information about environmental variables that can be tweaked to alter performance, and run time error messages.

Global commands

1. !|#. The rest of the line is a comment, and is echoed to standard output.
2. echo. Prints the rest of the line to standard output.
3. $ <operating system command>. The rest of the line (up to position 80) is a command, and is passed to the shell for execution.
4. dir [<OS specific args>]. Prints list of files in current directory.
5. [set] pwd [<new directory>]. Prints name of current directory, or changes directory to that specified.
6. file rename <name> [<new_name>]. Renames a file.
7. file delete|query|cat|head <file1> [...<fileN>]. Deletes, tests for existence, prints contents or first 10 lines of a file or files.
8. file print [/<search string>/] [(<fmt>)] [+] [NR] [<col1> ...<colN>] <filename>. Prints or searches contents of a text file. Can select subset of fields to print from each record, and can use a Fortran format statement to format the output. A + allows one to print lines below a line matching the search string. Search recognizes wild cards "*" and ".".
9. file transpose <file1> [...<fileN>]. Transposes rows and columns forming contents of a file or files. If a row is short, then the corresponding column will have missing value tokens making up the total.
10. file inverse <file1> [...<fileN>] [<ridge_constant>] Inverts a numerical matrix read in sparse form (each line of file give row and column indices, followed by the element value.
11. file vcf <filename> ([<start_position> [<end_position>]]) | ([<locus_name> [... <locus_nameN>]]). Prints locus names and positions within a VCF files. Can select subset of loci to print by physical position or locus name.
12. file vcf order|liftover <VCF_infile> <VCF_outfile>. Subsets and orders records of a VCF file to match the order of the active locus map. Writes a new VCF file. Sorting a VCF file thus requires three steps: read the loci into Sib-pair, use the order command, then file vcf order. The liftover modifier leads to the VCF locus genomic position being altered to that on the Sib-pair map.
13. file tbi <filename> [(<position> | <locus_name>)] [ann] . Prints summary of tabix index file (list of indexed sequences and range of map positions) or prints the record from the indexed VCF file that matches those coordinates (directly specified or current map position for a given locus). If the annotation modifier is added, the values of the INFO variables for that record are printed one line per variable-value pair.
14. file fasta <filename> ([<start_position> [<span>]]) | ([<index>. Prints summary for each sequence in the FASTA file, a specified subsequence, or produces an index file of the (samtools) faidx format.
15. clear. Restarts the program, closing all workfiles and zeroing all arrays.
16. help [All|Examples|Globals|Operators|Data| Analysis|<search_string>]. Prints a brief description of the commands -- either all, a subset, or all matching the search string.
17. info. Information about program settings and the current dataset. For the latter, gives counts of active and inactive pedigrees, individuals, and loci and a table of numbers observed for every trait and marker.
18. list|ls|which [<loc1> [[to] <locN>]] [$(a|q|m|h|x|A|D)[m|r]]. List of loci in current analysis. The "$<letter>" keyword denotes the list of of one class of locus, which can be reordered according to genetic map position (m) or in reverse (r) to the dataset ordering. The which command gives the position/rank of each locus argument in the total list of loci. If the plevel is set below -1, an unadorned list of locus names is printed (suitable for directly being read by other programs).
19. list|ls|which where <search_string>...|(position <pos1> [--] [<pos2>] ...) List of loci in current analysis subsetted by annotation matching a given search string or for a range of map positions. Positions take the form <chr>:<coordinate> or chrXX <coordinate>. If a chromosome is not specified for a given position, then the chromosome for the previous position is used.
20. show chromosomes. List of active marker loci listed by chromosome.
21. show loci. List of active loci along with table of numbers available (as per info).
22. show pedigrees. Same as gener, with print level 0.
23. show ids [<search-string1>]. List of indidivual IDs tabulated versus pedigree(s), or list of IDs meeting a search criterion.
24. show map (position <pos1> [--] [<pos2>]...)| [<loc1>...]|(where (chromosome <chr>...)|<search_string>...). Shows the current marker map, or a subset of the map selected by map position or range of map positions, or by locus name or annotation string. Positions are given <chr>:<pos>, where position is in the current map units.
25. show macros. Shows the current macro definitions.
26. show missing. Same as show loci but with number of missing values per locus.
27. show sex. Table of counts of each sex
28. show spectra. Table of nucleotide allele sprectra of active loci.
29. time. Print time elapsed since start of the program.
30. set timer [on|off]. Show time taken by each command.
31. include <command_file>. Read in Sib-pair commands from a file. If the command file is not specified, a directory browser will start up.
32. output [<output_file>]. Divert Sib-pair output from standard output (the screen) to a file. Issuing the output command a second time closes the output file. If the command file is not specified on the first call, a directory browser will start up.
33. macro <macro_name> [ (= <macro value>) | (<- all|fre <marker> | bur | che | epo | imp | ite | las | lik | ls | min | ple | pri| pva | pwd| sex | twi)]. Create a macro variable or function. If an equals sign is present, the rest of the line is taken as the macro variable body. If the <- (setting) sign is present, the macro variable takes the current value of the named program setting (or list of alleles or markers).

    Otherwise, multiple lines of the body of a macro function are expected to follow on subsequent lines, terminating with an empty line or ";;;;". The function body contains Sib-pair commands and parameters for substitution, represented as %1, %2 etc. The macro is called either as a function or a variable.

    To call a macro as a function, the entire macro name must be the first word on the command line, followed by any required parameters (which can be any legal string). The function body is evaluated, and the resulting Sib-pair command(s) is then acted upon. For example:

    # Draw a pedigree (need to have dot and gv):
    # First set up macro
    >> macro draw
    draw> select pedigree %1
    draw> write dot %%.dot %2
    draw> $ dot -Tps -o %%.ps %%.dot
    draw> $ gv --scale=-4 %%.ps
    draw> $ rm %%.dot %%.ps
    draw> unselect
    draw>
    # Call macro
    >> draw ped1 trait
    

    The %% keyword is replaced by the macro call ID, a unique 5 character string. The %0 keyword is replaced by all the macro parameters, while %+2, %+3... is replaced by that parameter and all subsequent parameters.

    To call the macro as a variable, the entire macro name is prepended by the "%" character. It can then appear anywhere in a command. The macro variable is replaced by the macro body, but there is no evaluation of macro function parameters or %%. The variable name can be delimited by brackets, so that it can be embedded in another string. The resulting command is then parsed:

    >> macro a=1
    >> macro b=+
    >> macro b1=2
    >> %a%b%a
    => 2.
    >> 1%b1
    => 12.
    >> 1%(b)1
    => 2.
    >> macro a=D14S52 D14S43
    >> tab %a
    ------------------------------------------------
    Cross-tabulation of "D14S52" ... "D14S43"
    ------------------------------------------------
    [...]
    

34. eval [<Scheme expression>]. Accesses the Scheme read-eval-print-loop. If called without arguments, presents the Scheme prompt "%%", otherwise evaluates the rest of the line as if it were a Scheme expression. Sib-pair macro variables and functions are stored within the Scheme environment, and so can read and written to. The Sib-pair specific extensions to Scheme are "(nloci)", "(ls)", "(loc)", "(loctyp)", "(locord)", "(locnotes)", "(locnotes-set!)", "(read-line)", "(run)" and "(help)".

    (ls ["m" | "x" | "h" | "a" | "q" | "d[mxhqa]"])	List locus names
    (nloci)	returns total number of loci
    (loc <index>)	returns locus at that position in the locus list
    (locord <loc>)	returns position of a locus in the locus list
    (locnotes <loc>)	returns notes for a locus
    (locnotes-set! <loc> <string>)	rewrites notes for a locus
    (loctyp <loc>)	evaluates type of a locus ("adhmqx")
    (read-line <port>)	Read next line from port as a string
    (run "<Sib-pair command>")	Run a Sib-pair command
    (string-split <string> [<sep>]	Convert from string to list of words split on white-space (or a specified character separator)
    (substring? <sub> <str>)	returns start of substring in string.
    (system "<operating system command>")	Run an OS command
    (help)	Information about this Scheme implementation

    See the appendix for more details.

35. last [<line_number>]. With no argument, displays the command history, otherwise submits that line of the history for reevaluation. A negative line number counts backwards from the current line. The command history is saved to a file "sib-pair.log".
36. quit|bye. Halts the program.
37. set prompt [on|off]. Displays a prompt, and activates/resets the command line history.
38. set gui [on|off]. Activates or stops the GUI. The GUI is either Java based (using the japi library to interface with the Java AWT library), or GTK2.0 based (using the pilib library to interface), depending on the compile time choice.
39. set ndecimal_points [<nwid>] <ndec>. The total width (number of characters) of a quantitative variable written to a new pedigree file defaults to 9 (and is fixed to 8 for some files, notably MENDEL and FISHER) but can be set as high as 20 for GAS and LINKAGE format files. The number of decimal places can be set to ndec.
40. set epoch [iso|jul|<epoch>]. Set or show the epoch used for julian dates. Defaults to "iso" epoch of 1970-01-01.
41. set out|plevel <level>|verbose|on|off. Increasing the print level causes more information to be printed by almost all procedures. Print level 1 prints out the identities and genotypes of parents imputed where the genotype was missing, raw counts of genotypes for the hwe procedure, expected ibs probabilities for the asp procedure etc. Print level 2 (or verbose) writes out the statistics for each simulated dataset for the MC based procedures, the intrapair variance and ibd sharing for each pair in the sib pair analysis, etc. Print level -1 omits outputting a list of pedigrees.
42. set printstyle [rectangular|pairwise]. Sets the format used by the print command. The default type is rectangular, so that the selected records are printed on one line with values at each variable in regular columns. If the pairs style has been set, then each value is printed as a variable_name=value pair.
43. set categorical_trait levels|labels. Sets the output representation of values of a categorical trait to either the level (1..N) or the label (as seen when first read in, or from the annotation.
44. set table_separator <separator_character>. Sets the character used to separate columns in summary tables of test results. Defaults to " " (a space). This affects the output from the summary, hwe, frequency, assoc commands.
45. set weight founders|imputed. Weights contribution of each pedigree to the allele frequencies by the number of typed founders, or alternatively gives the count of the founder alleles, observed and imputed.
46. set analysis [imputed | observed]. Includes imputed genotypes in generalized linear models fitted using the regress command. Genotypes are automatically reimputed each time regress is called. The idea of this is to allow multiple imputation association analysis. Imputation is performed using Sib-pair's single locus genotype MCMC sampler. The regress command with imputed genotypes respects the set frequencies command, so that the population allele frequencies can be specified in advance.
47. set burn-in <number of iterations>. Controls the number of Markov Chain Monte-Carlo iterations used by the apm algorithm discarded before estimation commences. Default is 100 iterations. Setting bur to zero means no burn-in is performed (the old default).
48. set iteration<number of iterations>. Controls the number of iterations used by the various Monte Carlo algorithms. Default is 200 iterations. Setting ite to zero means the Monte Carlo procedures are not performed.
49. set starting_trials <number of iterations>. Controls the number of Monte-Carlo trials used to generate a starting genotype configuration for the MCMC algorithms. Defaults to 5000. Sib-pair uses either this gene dropping approach or a sequential imputation approach using Lange-Goradia elimination to produce a legal set for ungenotyped individuals in a pedigree.
50. set emit <number of iterations>. Controls the number of (Monte-Carlo) Expectation Maximization iterations used by the mcf algorithm.
51. set batch <number of batches>. Controls the number of batches used by the fpm algorithm for the estimation of parameter standard errors.
52. set chain <number of MCMC chains>. Controls the number of chains used by the fpm algorithm.
53. set tune <MCMC tuning parameter>. Controls the multiplier for the MCMC proposal step size used by the fpm algorithm. The base step size is usually the fixed effects model standard error for that parameter, and tune defaults to 0.3.
54. set jackknife <number of jackknife draws>|off. Controls the number of deleted observations used by the describe command's algorithm for the estimation of standard errors for familial correlations for a quantitative trait. If off modifier chosen, reverts to default behaviour.
55. set mft (mendell-elston|genz <number of function evaluations>. [<absolute error> [<relative error>]]). Controls the multivariate normal CDF algorithm (TOMS717 Mendell-Elston or Genz). The number of function evaluations and absolute or relative tolerances used for integration by the Genz mft algorithm can also be specified.
56. set mincount <minimum numerator of P-value>. Controls the number used for Monte Carlo simulation of a P-value. Default is 20 pseudosamples with a test statistic more extreme than that for the observed statistic. Set mincount equal to iter if this is not desired.
57. set order_statistics <number of order statistic>. Controls the number of Monte Carlo simulated statistics used to extrapolate P-values for even more extreme observed statistic values using the approach of Davis & Resnick [1984]. Defaults to highest 10.
58. set chi-square [gibbs|pearson]. Set whether the Pearson or Gibbs chi-square is used for categorical trait association. Defaults to Gibbs.
59. set seeds <seed1> <seed2>< seed3>. Initializes random number generator seeds to given values, rather than via system time.
60. set optimizer varmet|bobyqa. Choose BOBYQA or VARMET optimizer for variance components model fitter.
61. set fbatimpute on|off. Enables (default) or disables imputation of missing child alleles based on their own offspring.
62. set vcf [<Alternate_allele_count_variable> <Total_allele_count_variable>]. Defines names of the VCF INFO variables from which reference counts of alternate and total alleles will be read by the "ass vcf" command. Defaults to "AC" and "AN".
63. set tdt bot[h parents]|one [parent]|first. Limits TDT statistic to cases where either both parents or at least one parent is typed, or one proband per family where both parents typed.
64. set hre zero|children. Assume zero recombinants between markers for dis command where parents genotyped, thus counting four imputed parental haplotypes. Alternatively, only utilize two haplotypes from children.
65. set roh [<minimum_length>]. Sets the minimum length of a run of homozygosity to be used to contribute to the F_roh estimator of inbreeding. Defaults to 1.5 Mbp (or 1.5 cM).
66. set model [allelic|genotypic]. Select allelic or genotypic encoding for marker loci in regression models. Default is allelic encoding.
67. set map function kosambi|haldane. Set the mapping function used by multipoint analytic and locus file output routines.
68. set map units [cM|M|Mbp|kbp|bp]. Set the map units for reading scripts and writing to the screen.
69. set sml <Frequency of A allele> <Penetrance of AA genotype> <Penetrance of AB genotype> <Penetrance of BB genotype>. Set default single major locus model written to locus files (usually p=0.05, pen=(0.5, 0.5, 0.05)).
70. set liability <binary trait> <liability class indicator> <number of classes>. Declare a quantitative trait to indicate liability class for the named binary trait. Used to generate Linkage format locus and pedigree files.

Utilities

71. pchisq <chi-square> <degrees of freedom> [(ncp <ncp>) | <df2>]. Calculate P-value for central or noncentral (if the ncp keyword is present) Chi-square distribution (or F-distribution).
72. qchisq <P-value> <degrees of freedom> Calculate quantile corresponding to the given P-value for the central Chi-square distribution.
73. chisq <nrows> <ncols>. Calculate contingency chi-square and permutation P for flat table entered via keyboard.
74. polychoric <nrows> <ncols>. Calculate polychoric correlation and association P for flat table entered via keyboard.
75. proportion <numerator> <denominator> <confidence interval width>. Calculate accurate confidence interval following Wilson (as described by Agresti and Coull) for a proportion.
76. tetrachoric <prevalence> <recurrence risk ratio> Calculate tetrachoric correlation equivalent to given recurrence risk ratio and trait prevalence.
77. sml <Frequency of A allele> <Penetrance of AA genotype> <Penetrance of AB genotype> <Penetrance of BB genotype>. Calculates recurrence risks and segregation ratios under a specified diallelic generalized single major locus model.
78. sml <Frequency of A allele> <Mean for AA genotype> <Mean for AB genotype> <Mean for BB genotype> <standard deviation for AA genotype>. [<AB SD> [<BB SD>]]. Calculates mean, variance components and parent-offspring regression results under a specified diallelic generalized single major locus model.
79. grr <trait prevalence> (<Frequency of A allele> <genotypic risk ratio> [add|dom|rec]) | (<Frequency of A allele in cases> <Frequency of A allele in controls> case-control [population-controls). Calculates recurrence risks and segregation ratios under a diallelic generalized single major locus model specified via trait prevalence, ratio of penetrances and pattern of inheritance (codominant multiplicative, dominant or recessive). If the case-control modifier is present, instead it expects the prevalence, risk allele frequency in cases, and risk allele frequency in controls who can either be unaffected, or unknown status (population).
80. ito <Frequency of A allele> [<Penetrance of AA genotype> <Penetrance of AB genotype> <Penetrance of BB genotype>]. Calculates conditional genotype frequencies for a relative of a proband of given genotype (if only allele frequency provided) or affection status under the specified diallelic generalized single major locus model (if penetrances were also specified). ITO refers to the matrices used to perform the calculation for pairs of relatives [Li & Sacks 1954].

Algebraic operators and functions

81. "<allele1>/<allele2>". Double quotes mark the contained text for special evaluation by the parser. A constant genotype is written as two numbers (1-999) or letters (a-zA-Z) separated by a slash and surrounded by quotes. Other quoted items are passed intact to be read, either as a reserved command or as a single Fortran real, so "1+3" is evaluated as 1000, and "1 1" as 11.
82. (<value>|<locus>) *|/|mod|+|-|^ (<value>|<locus>). Arithmetic operations combining numerical constants and/or trait values. The result of an operation involving constants is a single constant, but an operation involving a trait value results in nobs results (where nobs is the number of individuals in the pedigree file.
83. (<value>|<locus>) <|>|=<|=>|ge|le| eq|==|ne|^=|and|or (<variable>|<locus>). Logical operations comparing numerical constants and/or trait values. when operating on genotypes, the equality and inequality operators require both pairs of alleles to meet the criterion, but the comparison operators test true if either pair of alleles meets the criterion. That is "2/2">"1/3" evaluates to True, but "1/2"=="2/2" evaluates to False.
84. if <logical expression> then <action> [else if <logical expression> then <action>]... [else <action>]. Conditional evaluation of expressions. The if statements can be nested.
85. log|log10|sqrt|exp|sin|cos|tan|asin|acos|atan|abs|int|round (<variable>|<locus>). Functions acting on numerical constants and/or trait values.
86. rand|rnorm. Produce a random value from U(0..1) or N(0,1).
87. istyp|untyp <marker>. Test if genotyped at given marker. Necessary since if imputation is higher than -1, all untyped individuals have a genotype containing negative allele numbers (used to start MCMC algorithm).
88. ishom|ishet <marker>. Test if homozygous (or heterozygous) at given marker.
89. alla|allb <marker>. Return the first or second allele for each individual at the given marker.
90. marcom. Show the maximum of the number of markers an individual and any of his relatives are both genotyped at.
91. numtyp|anytyp|alltyp|protyp. Show number of markers an individual is genotyped at, or indicate whether genotyped at any one or all marker loci, or the proportion typed out of all marker loci .
92. male|female|isfou|isnon. Test sex and founder status of individual.
93. num|nfound. Number of members and number of founders of the pedigree containing an individual.
94. famnum|index. Position of pedigree and of individual in the active dataset.
95. chosen. Indicates individuals who were affected by or contributed to the last operation eg hash, merge, table.
96. <expr> : <expr> [:...]. The colon separates the members of a sequence of algebraic expressions. Its main function is to allow multiple expressions in the branches of an if-then-else statement, but it will minimize overheads in repetitive calculations. If multiple bracket-enclosed expressions are encountered, then a : is automatically inserted:

    >> (a=rnorm) (b=a^2)
    

    Command iteration

97. ... { (<val1> <val2> ...<valN>)|{<val1> : <val2>} | ($[mxqa]) } .... Repeat the associated command, placing each value of the iterator list into the position the list currently occupies. There may be multiple iterators, and iterator lists may be nested. This macro extension allows much of the functionality of proper computer languages:

    Iteration over a range of numerical values of a function	sml {0.1 0.3 0.5} 1 0 0
    Iteration over a range of loci of a function	tab a1 m1 m2
    Iteration over a range of functions	{tdt ass} a1
    Combinatorial generation of strings eg locus names	set loc a{1 : 3} aff
    Compound statements	if (male) then m{1 2}=x

    Data Declaration commands

98. set datadirectory <pathname>. Sets directory to be searched for pedigree files.
99. set workdirectory <pathname>. Sets directory to which temporary files are written.
100. set impute off|on|full|sequential|lange-goradia|nil|<level>. Toggles imputation routine.
     • Imputation level 0 (off) does not impute genotypes. It does generate legal genotypes via gene-dropping for all untyped individuals, that are used to start the MCMC genotype sampler. These genotypes are "hidden" to other routines. The gene-dropping approach is slowed by, and can fail (stochastically) in large pedigrees. In this case, the imputation level can be set to -1 or nil (completely off!).
     • Imputation level 1 (on) imputes single individual's genotypes if unambiguous using results of the Lange-Goradia imputation algorithm. All other missing genotypes are silently imputed via gene-dropping.
     • Imputation level 2 (full) carries out the same imputation as level 1, but the imputed values for all missing genotypes are printable and saveable using write. If imputation has already been carried out, changing the imputation level to 2 affects just the printing of "hidden" imputed values of unobserved genotypes.
     • Imputation level 3 (lange-goradia or sequential) imputes values for all missing genotypes via sequential application of the Lange-Goradia genotype elimination algorithm. That is, the missing genotype is imputed to be the most likely genotype conditional on the typed members of the pedigree, and those genotypes imputed prior to that iteration. The imputed genotypes are not visible to non-MCMC commands (eg write) unless the imputation level is changed to 2.
101. set errordrop off|on|<level>. Toggles automatic deletion of genotypes that give rise to a Mendelian inconsistency, either an entire nuclear family (level 1), or an entire pedigree (level 2, the default).
102. set checking off|on. Toggles the first level testing for Mendelian inconsistencies within nuclear families.
103. set locus <locus name> <locus type> [<map position> [<description...>]]. Declares position (by order within list), name and type of locus within pedigree file. Locus type may be either:
     
     

     marker	an autosomal (fully) codominant marker
     xmarker	X-linked codominant marker
     haploid	Y or mitochondrial codominant marker
     quantitative	quantitative (or interval or ordinal) trait
     affection	binary trait
     categorical	categorical trait

     
     
     It is best to avoid a locus name containing reserved characters (eg "+-*/()^"), if algebraic manipulation of that variable will be required (otherwise quotation of the name is required). Names identical to commands also cause trouble unless protected by brackets.

     The fifth column (optionally) contains the genetic map position. All subsequent words (up to a total of 40 characters) are stored as an annotation. The annotation is appended to the long form of output of some commands (eg show loci or list), and is searchable by some commands (currently keep|drop where). If you wish to annotate, but do not have a map position, a "." for the fifth column is unobtrusive and accepted by Sib-pair.

     The annotation for a categorical variable can also include labels for the levels of the trait. These are of the form "level=label". If the categorical variable is string-valued in the input dataset, then Sib-pair will convert to levels and write an appropriate annotation.

104. declare loci <number>(m|x|q|a) [<number>(m|x|q|a)... ]. Declare a batch of loci, automatically generating names: "trait1", "trait2"..."traitN", "mar1", "mar2"..."marN" as appropriate.
105. rename (<locus name> [to] <new name>)|[bim | vcf] map_file_name [INFO_variable]. Change name of previously declared loci. Either the new name is specified on the command line, or it is obtained from a map file where the appropriate replacement name is associated with a locus at the same map position as the target locus. If an INFO variable in a VCF file gives a locus name, this can be specified as an additional field, and used instead of the VCF file ID variable.
106. loci <command file>. Read in Sib-pair locus and pedigree file declarations from a file. If the command file is not specified, a directory browser will start up.
107. read locus linkage <locus file name>. Read locus names, types and map positions from a Linkage-format locus (.dat) file. Does not recognise factor coding of genotypes, but does create a new quantitative trait for liability class
108. read locus merlin <locus file name> [xli] [snp]. Read locus names, types from a Merlin-format locus (.dat) file. Markers are treated as sex-linked when the xli modifier is used. Markers are stored as SNPs when the snp modifier is used, with reduced memory consumption.
109. read locus plink <locus file name> [append] [human]. Read locus names, types and map positions from a PLINK locus (.map) file. The append modifier stops a dummy binary trait being declared as the first locus, and appends to rather than overwrites any existing locus information. The human modifier causes chromosomes 23...26 to be interpreted as X, Y, XY, MIT.
110. read locus file <locus file name> <locus_name> <chromosome> <map_position> <reference_allele> <alternate_allele> [human] . Read locus names, types and map positions from a rectangular locus (.map) file where the first line contains names for each column. The column names corresponding to the locus name, chromosome, map position (bp), and referance and alternate alleles are specified as arguments. These can be numeric if names are not provided.
111. read locus vcf <VCF file name> [<start> [<end>]] [human]. Read locus names, types from a VCF. Reading can be retricted to a subset of loci within an interval (specified in base pairs, prefixed optionally by a chromosome identifier). The human modifier recognizes chromosomes 23-26 as sex and mitochrondrial chromosomes.
112. set append [skip|version]. Controls how repeat declarations of a locus are treated. Default behaviour is to modify the name of the repeated locus ("name_v2", "name_v3"...), but this can be changed so that such declarations can be skipped.
113. read pedigree <pedigree file name>|inline [skip <lines_at_beginning>] [nopedid] [nosex]. Reads a GAS type pedigree file either from an external file, or inline following the command. The inline data is terminated by a line containing ";;;;". The skip keyword leads to the skipping past that number of lines at the beginning of the file. The nopedid and nosex modifiers allow the pedigree file to be missing a pedigree ID column (in which case all records are taken to belong to a single kindred), or a sex column (where sex defaults to missing for all records, and is usually imputed in a mating-consistent fashion)
114. read linkage <pedigree file name>|inline. Reads a pre-MAKEPED LINKAGE type pedigree file.
115. read ppd <pedigree file name>|inline. Reads a post-MAKEPED LINKAGE type pedigree file.
116. read cases <data file name>|inline [noid [sex|nosex] [sep <field_separator>] [skip <lines_at_beginning>]. Reads a data file containing unrelated individuals. The individual ID is the first column of data, which is followed by all the phenotypes. If the noid keyword is present, then an ID field is not expected, and row number will be used as ID. If the sex keyword is present then the second column of data is expected to be the sex. The skip keyword leads to the skipping past that number of lines at the beginning of the file.
117. read csv <data file name> [noheader] [<field_separator>]. Reads a CSV format file containing data for unrelated individuals. The default is to expect a header row giving the column names. The type of each column ("a", "q", "m", "c") is inferred based on their contents, and used to declare the variables to Sib-pair, and to read the dataset.
118. read vcf <VCF file name> [...VCF file N] [ped_id] [<lines_length>]. Reads only the IDs from a VCF file header and sets up the corresponding pedigree in memory. If the ped_id modifier is present, the ID is split to give a pedigree and individual ID (not applicable if multiple VCF file names are give). Neither the marker names nor genotype data are read from the VCF file.
119. read annotation <VCF file name> [<INFO_var1>[:<field>]...<INFO_varN>[:<field>] [sep] [dump]. Reads INFO variable values from a VCF file. If a variable name is not specified, lists the descriptions of all the variables from the header. Otherwise, finds the matching locus in the current dataset (by name then by position if the former fails), and prepends the value(s) to the Sib-pair locus annotation (locnotes), or dumps these to output. If an INFO variable is an array (with values separated by "|"), the field name or index can be specified.
120. show annotation <VCF file name> [<INFO_var1>[:<field>]...<INFO_varN>[:<field>] Prints cross-tabulations of INFO variable values from a VCF file. If a variable name is not specified, lists the descriptions.
121. read bin <data file name>. Reads a Sib-pair "binary" pedigree file. The "run" statement is bypassed, so this command reads in locus and pedigree data immediately. No checking is done for pedigree and Mendelian errors and reordering of pedigree members is not performed. These can be large files, but on systems where gzip is available, can be automatically compressed (see write bin) and decompressed. The default format for the file arises from a Fortran unformatted write of the locus and pedigree arrays, and so will be compiler and platform specific.
122. read hapmap <data file name>. Reads a HapMap style genotype data file. All individuals are assumed unrelated.
123. read plink <file prefix> [compress]. Reads a PLINK .bed and ancillary files: <prefix>.bim, <prefix>.fam and <prefix>.bed. The "run" statement is bypassed, so this command reads in locus and pedigree data immediately. No checking is done for pedigree and Mendelian errors or reordering of pedigree members -- this is assumed to have been performed by the program that prepared the files. If the compress modifier is present, then the genotype data is stored in a 4 bits per genotype format internally.
124. set sex marker <locus name>. Declares a sex-informative marker, such as Amelogenin.
125. set sex <report_threshold> [<hom_to_het_rate>]. Sets reporting threshold and assumed X-linked genotype heterozygote to homozygote miscall rate for diagnosis of sex from X-marker information.
126. set twin <quantitative trait> [merlin]. Declares a variable to identify twin zygosity, or specifically just monozygotic (MZ) twins. All individuals within a pedigree with the same value of this variable are taken to be part of the same twin sibship (twin pair or higher order multiple). Different values indicate different pairs within the same family. If the merlin coding has been specified, odd numbers indicate MZ twins, and even numbers DZ twins. For the default coding, nonzero values indicate MZ twins, and DZ twins can be indicated using a zero value for the indicator). This information is used to write out MZ twin indicators to the pedigree files used by MENDEL, MERLIN and SOLAR, and to test for marker inconsistencies.
127. set twin error <threshold> <minimum markers>.

     The error modifier shows and allows setting of the quick genotype discordance test (test dup). The first argument is the allowed mistyping (ie discordance) rate for genotypes compared between putative genetically identical pairs of samples. This defaults to 0.005, which is acceptable for SNP or sequence markers. The second argument determines how many markers must be present before the comparison is made, defaulting to 100 markers.

128. set skiplines <slines>. Skip first slines lines in pedigree file) when reading in.
129. order <loc1>...[<locB> to <locC>]... [$(m|x|q|a)[r|m]]...<locN>. Set order of loci. Addition of r to a class eg $mr, reverses the order of all members of that class, while the m modifier causes the order to be the genetic map order. You may have to revise the genetic map order (by set map or set dist to get sensible export files for some programs such as Linkage (Sib-pair assumes a map position lower than the preceding position implies the markers are unlinked).
130. set map <pos1>...<posN>. Set map positions for the marker loci. This will overwrite any original map positions.
131. set distances <dis1_2> <dis2_3>...<posN-1_N>. Set interlocus map distances map positions for the marker loci. Distances are in centiMorgans. This will overwrite any original map positions.
132. set chromosome <chr1> [...<chrN>]. Assign each marker locus to a specific chromosome. If there are more markers than specified chromosomes, the last specified chromosome is reused.
133. read map <map file name> [[k]bp]. Read in map positions for loci from a file, matching via names of previously declared markers. Should recognize most formats of map file automatically eg GVF, Merlin, Mendel, Solar. Tests number of columns and whether column contents are numeric or alphabetical, skipping first row as possible header, or looks for identifying strings in the header. Default units are cM or Mbp. The bp (kbp) modifier tells Sib-pair to divide the positions by 10⁶ (10³); thus the map distances become Mbp, and remain readable.
134. read chain <chain file name>. Read in the specified UCSC chain file and update the map positions of currently active loci using that information.
135. set frequencies <marker> [ [<allele_name1>] <allele_freq1>... [<allele_nameN>] <allele_freqN>)]. Sets the "population" allele frequencies for a marker to be used by MCMC procedures that simulate missing genotypes for that marker. This currently only affects the set analysis and gpe commands. Only one marker at a time can have specified frequencies. The obs option can fix the frequency to that observed for the currently active dataset. To free the frequencies (allow calculation from the dataset), call the command without specifying any frequencies after the marker name.
136. run. Reads in pedigree file and creates working pedigree file. Imputes genotypes if requested.

Data manipulation commands

137. keep|drop [<loc1>...[<locB> to <locC>]... [$(m|x|q|a)]...<locN>]. Retain or exclude loci for subsequent analysis. Consecutive loci can be summarized as a range, as can all members of a particular class of locus type (marker, quantitative, affection) via a class ($type) token. Note that dropped variables can still be used in algebraic and logical expressions.
138. undrop [<loc1>...[<locB> to <locC>]... [$(m|x|q|a)] ...<locN>] Return previously dropped loci to analysis. Default is to undrop all dropped loci. Loci can be selected on name (including wild cards), class. This is not the reverse of the delete command.
139. keep|drop|undrop where (monomorphic
     | diallelic
     | snp
     | allele_number [<c_op>] <numal>
     | spectrum <allele1-allele2-...alleleN>
     | max <frequency> [<c_op>] <allele_frequency>
     | number_typed [<c_op>] <ntyp>
     | missing [<c_op>] <ntyp>
     | chromosome <chr1> [...<chrN>] [trait]
     | distance <smallest_gap>
     | r2 <smallest_r2>
     | every number_skipped
     | position <pos1> [--] [<end_position>;] ...
     | near <loc> [...<locN>] [<N>]
     | hwe_p [<c_op>] [<critical_P>]
     | homoz [<c_op>] [<threshold>]
     | test_p [<c_op>] [<critical_P>]
     | in [<file>]
     | covers <trait> [<n_uncovered_cats>] | <search_string>). Retain or exclude loci for analysis. Note that dropped variables can still be used in algebraic and logical expressions. The where condition can be used to match the set of loci meeting that condition. Available conditions are: test that a marker is monomorphic or diallelic or a specified number of alleles, that the commonest marker allele frequency exceeds or falls below a threshold, the number (or proportion) of individuals typed or missing exceeds or falls below a threshold, the marker is closer than a set amount to the last included marker (map distance), marker is in too high linkage disequilibrium (r2) with the last included marker, every Nth marker in list, on a given chromosome or chromosomes (may be a trait eg expression level), the N nearest loci to a target locus, within an interval or intervals on the genetic map, the HWE test or most recently applied test P-value is smaller (or larger) than a critical level (defaults to 0.05/Nmarkers), the marker homozygosity, the number of categories of a specified trait that are completely ungenotyped (defaults to zero), locus name is listed in a text file, or the marker annotation matches the search string.
140. select [containing|exactly <nprobands>] [where] <a logical expression>. Select pedigrees containing one or more individuals with a trait value meeting the criterion.
141. select pedigree|id [[not] in] (<ped1>...<pedN>)|file <file>. Select pedigrees included or excluded from a list of pedigree or individual names, which can be in a file. The names can contain wildcard characters: "." (match any character in the target at that position in the search string) and "*" (match any characters zero or more times in the target at that position in the search string).
142. unselect [<Nth>]. Returns all pedigrees excluded by a select command back to the analysis. If an integer argument is given, this gives how many select statement subsettings to roll back. The argument can be negative (reversing the effect of a previous targetted unselect).
143. pack loci|pedigrees. Permanently delete all loci currently excluded by a drop command, or all pedigrees currently excluded by a select command from the work file.
144. edit <pedigree> <person>|all <trait> to <value> [<new value>]. Allows editing of trait values or genotypes. The all keyword performs the action on all members of that pedigree: since wildcards can now be used, an equivalent is "edit <ped> *".
145. copy <pedigree1> <person1> <pedigree2> <person2> [insert | merge] Copies all active trait values and genotypes from first record to second record. If the merge modifier is present, copying only occurs where the existing value of the target is missing at the variable.
146. update|merge [<file_type>] <file_name> :

     merge	[genotypes] [illumina|locus_first] [qua <c_op> <thr>] [skip lines_to_skip] |
     	[probabilities (key <mergekey>)|pedid|id [mach | (impute2 file <locus_file>)]]
     	[plink|bed <root_prefix>] [join|compress] [id] [pos]
     	[vcf <VCF_file_name>] [ped_id] | [quality <QCstat> [<c_op> <thr>]].
     	[dose <PLINK_dosage_file_name> [<thresh>]] |
     	fimpute <locus_file_name> <genotype_file_name> [compare]
     	[mac <datfile> <pedfile> [<thresh>] |
     [<locus1> ...<locusN>] <phenotype_file_name> [lines_to_skip].

     Updates phenotype/genotype data in the current dataset using values read from a file. The default file format requires the first line of the file to give the names of the variables that are included in the subsequent lines. Usually, the first column is the pedigree_ID and is named ped[igree]; the second column is the individual_ID, and is named id. Alternatively, the pedigree ID can be omitted, so the merge requires the individual IDs to be all unique. The remaining columns should have names that match locus names in the current dataset (data for nonmatching names are skipped).

     ped id locus_name1 locus_name2 locus_name3 ...
     1    1  A/A         12.434      y ...
     1    2  A/B           x         n ...
     ...
     

     Where the pedigree and individual IDs for a record in the update file match that of an active individual in the current dataset, the corresponding phenotype and genotype values for that individual are updated using the values read from the file. If merge was called, only missing values in the current dataset are updated, but all values are overwritten if the call was to update. By specifying locus names on the command line, you can further control the loci that are updated from the new file.

     The merge genotypes option reads files containing one record per individual locus genotype (ie fields are individual_ID, locus name, genotype):

     1    locus1  A/A
     1    locus2  1/2
     ...
     

     or if the locus_first modifier is present:

     locus1  1 A A 0.85 ...
     locus2  1 1 2 0.99 ...
     ...
     

     or if the illumina modifier is present:

     [Header]
     ....
     [Data]
     SNP Name  Sample ID GC Score  Allele1 - Forward Allele2 - Forward
     locus1  1 0.85 A A  
     locus2  1 0.91 2 9
     ...
     

     If a quality score is present, this can be used to filter out poor quality genotypes. This is automatically enabled for the Illumina report format, with a threshold of 0.6. This can be altered using the quality_score option. For other formats, the quality score is assumed to be in column 5.

     The merge plink option reads the supplied PLINK .bed, .bim, .fam files and merges the new SNP genotypes to the existing pedigree data. The loci merged are those active in the current SIb-pair dataset and present in the .bim data under the same locus name, unless the pos modifier is present, when matching is by map position. If the join modifier is present, it also automatically reads the accompanying .bim file and declares all these loci before merging -- this is silently called by the "read plink" command. The id keyword sets the merge key to be only the individual ID. The compress option stores genotypes as 4 bits per genotype, but is less flexible than the "snp" storage mode.

     The merge probabilities option reads SNP genotypic probability files from Beagle, Impute2, and Mach, converting to the most likely genotype. One can specify the name of a mergekey variable using the key modifier (giving the column number that corresponds to each individual) or the pedid or id keywords to specify that the probabilities file column headers contain ids of the form <pedigree_ID>-<individual_ID> or <individual_ID> respectively. If headers are absent, as in the case of Impute2, these can be read from a locus file, as specified by the file modifier.

     The merge dose option reads SNP allelic dosages from a dosage file prepared for PLINK. A "hard" threshold is used to convert dosages to genotypes, which can be set by appending a <threshold> value to the command. The default threshold value is 0.5, giving cutpoints of 0.5 and 1.5 (the dosage scores take values 0-2).

     The merge mach option reads SNP allelic dosages from a pedigree prepared using MaCH or minimac. The names of the pedigree file and datfile specifying the locus names need to be specified. A "hard" threshold is used to convert dosages to genotypes, which can be set by appending a <threshold> value to the command. The default threshold value is 0.5, giving cutpoints of 0.5 and 1.5 (the MaCH dosage scores take values 0-2).

     The merge vcf option reads SNP genotypes from a VCF file. The ped_id modifier tells Sib-pair to match on pedigree and individual ID, where the VCF IDs are expected to take the form <pedigreeID>_<individualID>. The name of a genotype quality score present in the VCF file can be given, and optionally a comparison that will be used to filter the genotypes to be tested. The mean quality score for each locus is saved, and is accessible to the summary command.

     The merge fimpute option merges in genotypes from an FImpute genotype file (both FImpute locus and genotype files need to be specified). The compare option compares genotypes in the specified genotype file to those in the current dataset matching on individual and locus, including enumerating the updateable genotypes from FImpute.

147. delete <pedigree> <person>|all Sets all data to missing for a specified individual. The all keyword performs the action on all members of that pedigree.
148. delete [<locus1>...<locusN>] [whe|where] <a logical expression>. Sets specified data to missing for all individuals meeting particular criteria.
149. get <relationship> <statistic> <trait> [<newtrait>]

     get	all	mean	<trait>	[<newtrait>]
     	offspring|children	minimum
     	sons	maximum
     	daughters	sum
     	parents	count
     	father	sample
     	mother
     	spouse
     	husband
     	wife
     	siblings
     	brothers
     	sisters
     	mztwins

     
     Summarizes trait values of all relatives of the specified class of all individuals, saving the result to a trait locus if requested. The sample option samples with replacement unless the all option is used, when the sampling is without replacement. Statistics for siblings include ego ie give sibship mean etc; sampling for siblings excludes ego.

150. recode (<marker>|$(m|x)) [frequencies|letter|number|nucleotides]. Recodes alleles at that marker or set of markers to 1..N, where the ordering defaults to the allele size (or collation order for letter alleles). If the freq modifier is present, the numbering is by ascending allele frequency. If the let modifier is present, the numbering is "1..4" for "ACGT", and the reverse for num. The nuc option recodes for genotypes using "A/B" allele coding to the appropriate nucleotides provided in the locus annotation, where it should be in the format "[<A1>/<A2>]".
151. recode (<marker>|$(m|x)) [reference <filter_trait>]. Recodes missing genotypes at that marker or set of markers to the wild type homozygote. An indicator trait can be used filter this so this restricted to a subset.
152. recode <marker> <all1|value1>...<allN|valueN> to <new allele|new value> [...<newN>]. Allows pooling and/or recoding of marker alleles prior to subsequent analysis. If there are fewer new values than old values, the last new value is recycled.
153. combine <marker1> [...<markerN>] [<threshold>]. Pool rare alleles for a marker into one new allele. "Rare" defaults to a frequency of 5%, but can be changed via the last parameter on the command line.
154. swap <marker> [to <marker2> ] [...<markerN>]. Swap diallelic marker alleles names eg "[12]" -> "[21]".
155. swap <trait> [lod <swap_lod_threshhold>]. Tests for and repairs likely allele swaps with particular levels of a categorical trait. Calculates likelihood for table of association between marker and the target trait, and then for each table where the allele labels are swapped for a single level of the trait. If the ratio of likelihoods for each ordering expressed as a lod exceeds the threshhold - default 2 - then the allele labels will be swapped for individuals where the value of the target trait matches the most deviant level.
156. flip <marker> [to <marker2> ] [...<markerN>]. Recode SNP marker alleles to the complementary strand coding ie "[ACGT]" -> "[TGCA]".
157. flip [map|fasta <file>] Compares observed alleles for a SNP to reference alleles, either from the locus annotation (where they will be written [A/B] with A the consensus or reference allele) or those in an external VCF, GFF file, and flips or swaps the alleles to match. It also alters the locus annotation reference alleles to match the external file. If allele frequencies are provided in the external file (eg VCF INFO AF variable), this is used to test for strand flips for ambiguous SNP ie A/T and G/C. The flip fasta option changes only the locus annotation. The FASTA file must be indexed, ie have a .fai file as produced by the samtools faidx program or Sib-pair, and flip fasta changes only the locus annotation. Results are sent to locstat, so one can keep/drop/undrop on the value of the comparison for each SNP.
158. date <quantitative_trait> [julian|gregorian|year]. Convert a numeric variable from Julian to Gregorian, Gregorian to Julian, or Gregorian to "decimal" year. The "chronological" Julian date is the number of days since the epoch, usually 1970-01-01 or -4712-01-01. Gregorian dates are represented as 8 (or 9) digit integers of the form of (-)YYYYMMDD. The decimal years are YYYY.x, where the decimal part is the day of year number (from 1...366) divided by the length of that year (365 or 366).
159. date [(<yyyymmdd> julian)| (<juldate> gregorian)]. Convert a single date from Julian to Gregorian or Gregorian to Julian.
160. test age <age_variable> [<threshold>]. Test that parent and offspring ages are consistent. The threshold controls the minimum age of parents at birth of offspring, and defaults to zero, since the units are unknown.
161. test dob <DOB_variable> [gregorian] [<threshold>]. Test that parent and offspring DOBs are consistent. The threshold controls the minimum age of parents at birth of offspring, and defaults to 4380 days, if gregorian has been used to declare the DOB variable to represent Gregorian dates.
162. test age <age_variable> [<threshold>]. Test that parent and offspring ages are consistent. The threshold controls the minimum age of parents at birth of offspring, and defaults to zero, since the units are unknown.
163. test sex. Uses X-chromosome and/or Amelogenin to test the putative sexes of genotyped individuals. If Y-chromosome markers are declared, then presence of data for an individual is taken as further evidence of male sex; absence is ambiguous, as may simply be missing.
164. test haploid [mitochondrial]. Uses Y-chromosome or mitochondrial markers to test the putative relationships of genotyped individuals.
165. test map [(merge [<thresh>]) | <map_file>]. Test for multiple loci dclared to have the same map position. If the merge modifier has been specified, Sib-pair will merge the genotype data for the two loci if the allelic spectrum is compatible ie same alleles, or compatible with a strand mixup. If too many incompatibilities between genotypes at the same individual are detected (default threshold 0.005), then the merge will not proceed. If a file name is specified, the current map is compared to that contained in that file, which is read as if by "read map".
166. test ids <file>. Compare IDs in current dataset to those in a specified file, which may be either a simple list of pedigree and individual IDs, or a Sib-pair binary or VCF file.
167. test vcf <VCF_file> [ped_id] [qua <QCstat> [<c_op> <thr>]]. Compare genotypes for the matching IDs and loci in current dataset to those in the specified VCF file. Flips strand as necessary. The name of a genotype quality score present in the VCF file can be given, and optionally a comparison that will be used to filter the genotypes to be tested. The mean quality score for each locus is saved, and is accessible to the summary command.
168. test strand [<marker 1> [...<marker N>]]. Test if mix of genotypes for a locus are compatible with a strand mixup.
169. test flips map|fasta <file>|<source-indicator> [lod]. Compare reference alleles given in locus annotations to those in an external VCF, GFF or FASTA file. The FASTA file must be indexed, ie have a .fai file such as that produced by the samtools faidx program or Sib-pair. Results are sent to locstat, so one can keep/drop/undrop on the value of the comparison for each SNP. If the name of a stratifying variable is given, then differences in allele frequency between strata are tested to see if they are compatible with an allele swap in one stratum, and the identity of the most deviant level and associate lod printed.
170. test duplicate [ids [merge]] [<indicator>]. Test if genetically identical individuals are present in the dataset. Write an indicator value for each individual tested to a phenotypic variable if specified. If the ids modifier is specified, then testing is limited to persons with the same individual ID in different pedigrees. If found, and the merge modifier has been specified, then Sib-pair will reciprocally fill in missing data in one record using data from the other record. A duplication indicator can be written to a phenotype variable if this is specified.
171. test locus [<marker 1> [...<marker N>]]. Carry out the usual Sib-pair Mendelian error screen on the specified locus and only on active pedigrees. Useful if checking was turned off when the pedigree was first read in.
172. test lange-goradia [<marker 1> [...<marker N>]]. Carry out the usual Sib-pair Lange-Goradia elimination algorithm test for Mendelian errors on the specified locus and only on active pedigrees. Useful if checking was turned off when the pedigree was first read in.
173. test <pedigree> <individual [<ped2> <id2>]. Default is to test for the most highly related individual in the entire dataset. If a second ID is given, prints the genotypes at loci where they are discordant, along with the proportion discordance - good for detection of somatic mutations or genotyping errors

Analysis commands

174. transform <xtrait> <divisor> <subtractand> <power> <lower threshold> <higher threshold>. This transform the quantitative trait xtrait as:
     boxcox({xtrait-subtractand}/divisor)
     where boxcox() is (a slightly altered) Box-Cox transformation, so that:
     
     • if power=0, the transformation is log({x-s}/d);
     • if power=1, it is {x-s}/d;
     • and otherwise [({x-s}/d)^p-1]/p.
     
     The resulting transformed value can then be truncated above or below using a specified low or high threshold.
175. normality_test [<trait> [...<traitN>]]. Calculate Filliben correlation (correlation between value and expected quantile under assumption of normality) for each quantitative trait. The test statistics are accessible to subsequent summary, keep, and drop commands.
176. standardize <trait> [familywise]. Replace each trait value with its Z-score, ie (x-xbar)/sd, where xbar is the total sample mean, and sd the total sample standard deviation. This will be performed using the individual's family mean and standard deviation, if the fam keyword is included.
177. adjust <ytrait> on <xtrait> [to <adjustment value of xtrait|m|f>]. Performs linear regression of quantitative trait ytrait on quantitative or binary trait xtrait (or sex, if sex is set to on), calculates residuals, and adds adjustment value or, if not specified, the mean value of xtrait. The residuals then replace the original values of ytrait. A multiple regressive adjustment of Y on X₁ and X₂ requires sequential adjustment of Y on X₂,X₁ on X₂, and then Y on the adjusted X₁. Has been superceded by the more powerful residuals command.
178. residuals <ytrait> on <loc1>...[to]...<locN> [complete_obs]. Replace quantitative trait with the residuals from the multivariate regression on the list of predictors (which may include the average allele length of a marker locus). The com option means only individuals with no missing values for any of the listed traits will be updated. Otherwise, missing values are replaced with the sample mean for that phenotype when calculating the predicted value.
179. predict <ytrait> on <loc1>...[to]...<locN> [complete_obs] Replace quantitative trait with the predicted value from the multivariate regression on the list of predictors (which may include the average allele length of a marker locus). The com option means only individuals with no missing values for any of the listed traits will be updated. Otherwise, missing values are replaced with the sample mean for that phenotype. when calculating the predicted value
180. impute <ytrait>. Replace missing quantitative trait values with the predicted values from a regression on the spouse, sibling and offspring observed values. Designed mainly for imputing missing age or date of birth.
181. impute <ytrait> on <loc1>...[to]...<locN> [complete_obs] Replace missing quantitative trait values with the predicted value from the multivariate regression on the list of predictors (which may include the average allele length of a marker locus). The com option means only individuals with no missing values for any of the listed predictor traits will be updated. Otherwise, missing values are replaced with the sample mean for that phenotype when calculating the predicted value.
182. kaplan-meier <age-at-onset> <censor> [residuals]. Prints the product-limit estimator for the survivor function for the quantitative trait age-at-onset, where censor is the binary outcome trait, which is affected when age-at-onset represents the age at which the individual first expressed the trait. The age-at-onset is replaced by a nonparametric residual when requested. If affected, this is:
     sgn(1-H(t)).(-2(1-H(t)+log(H(t))))^1/2
     If unaffected:
     -(-2H(t))^1/2
     where H(t) is the Nelson-Aalen estimate of the integrated hazard function at that age t.
183. lifetable (<start>|0) <end> <censor> [<cohort stratum width> [<period stratum width>]] [time] [covariate <covariate>]. Prints the life table for the quantitative trait end, where censor is the binary outcome trait, which is affected when end represents the age at which the individual first expressed the trait. Start represents the time of entry into the study of the individual. If the start trait is "0", the start of observation is taken to be 0 for all individuals. The <cohort stratum width> is the bin width used to divide up the entry times into the study. The <period stratum width> is the bin width used to divide up person time of followup (and defaults to years). The <time> modifier causes the start and end values to be treated as a continous measure, rather than a Julian date (the default). The covariate keyword declares the following named variable to be a categorical covariate which will be used to stratify the dataset.
184. surv <age-at-onset> <censor> [<covariate1> [..<covariateN>]]. Carries out the logrank nonparametric test for equality of survival curves across the different strata formed by crossclassifying one or more covariates. The variable censor is the binary outcome trait, which is affected when age-at-onset represents the time at which the individual first expressed the trait. If no covariate is specified, tests all active markers in turn, carrying out a gene-dropping test of significance.
185. cumulative_position <trait> <index>. Generates a sequential index for each individual in a pedigree that meets a logical condition, actually if the trait is non-missing. In other words, numbers the probands.
186. rank <trait> <rank>. Write the ranks of a quantitative trait to the quantitative variable rank.
187. blom <trait> <blom_score>. Write the approximate inverse normal scores for a quantitative trait to the quantitative variable blom_score.
188. quantile_normalization <trait1>...<traitN>. Carry out quantile normalization over a set of quantitative traits, replacing each observed value with the mean value at the corresponding order statistics of all the chosen traits. Missing values are dealt with by reusing the nearest value rather than linear interpolation.
189. simulate <trait> [<h2> [<linked extant marker>]]. The data for the named trait is replaced by simulated data for a trait under the control of a QTL (or polygenes) with a total heritability of h2 (defaulting to 50%). If a second marker name is given, the controlling QTL is simulated as being completely linked to the second marker.
190. simulate <marker> [<linked extant marker>] [<number of equifrequent marker alleles> | <allele 1 frequency>...<allele N frequency>]. The data for the named autosomal marker is replaced by simulated data. If a second marker name is given, the new marker is simulated as being completely linked to the second marker. Either a set of allele frequencies, or the number of (equifrequent) alleles, can be given for the simulation. If the sum of the given allele frequencies is less than 1, an extra allele will be added automatically.
191. simulate qtl <trait> <marker> [<h2>]. The data for the named autosomal marker is replaced by simulated genotype data that is generated conditional on the trait values and the genetic model for the trait. The model for the QTL is taken from the results of the set sml command, and the residual heritability is specified as the last argument to the command. Binary traits are modelled under the multifactorial threshold model, and quantitative traits are not currently supported.
192. simulate pedigrees [<nped> [<ngen> [<min_number_of offspring> [<max_number_of_offspring> ] [<pedigreeID_prefix>]]]]. Generate a set of nped (default 100) random pedigrees, each of ngen (default 2) generations. The component nuclear families each contain between min_number_of offspring and max_number_of_offpsring offspring (defaulting to a range 0-2). The generated pedigrees are each descended from a single founder couple (with marry-ins). Multiple calls to simulate pedigree can be made in an incremental fashion. A string to be prefixed to the pedigree ID can be given as argument 5.
193. permute <trait> The phenotype values for the named trait are permuted within pedigrees.
194. nuclear [maxsibs] [grandparents]. Split pedigrees into component nuclear families, duplicating individuals as necessary. If maxsibs is set, then sibships containing more than maxsibs members are truncated. The gra option includes the grandparents as well.
195. subpedigrees. Split nominal pedigrees into component true pedigrees. Sib-pair normally can analyse a group of individuals with the same pedigree ID, even if they are not all related. This command splits such groups into uniquely named formal pedigrees.
196. join <pedname1> [<pedname2>...<pednameN>] Join or rejoin pedigrees into a single pedigree, appropriately dealing with shared individual IDs and their associated data. Note that Sib-pair allows noncontiguous pedigree blocks to use the same pedigree name, so multiple pedigrees of the same name will be collected up.
197. prune [<binary trait> [|<quantitative trait> over|under <threshold>]]. Reduce pedigree to contain probands and minimum number of connecting relatives.
198. cases <locus>. Reduce pedigree to unrelated individuals with non-missing values at the trait i.e. the informative founders, and any informative nonfounders who are not directly related to any individuals already selected.
199. unique_id [sequential]. Generate unique consecutive (within family) numerical IDs for all individuals (as well as new numeric pedigree IDs). The sequential gives IDs from 1...total_records, instead of 10001, 10002...20001...
200. hash [ <file_of_IDs> | file <file_of_IDs> [<col1> [<col2>]] | <pedigree_ID> <individual_ID> | id <individual_ID> | show [locus ]| delete | size <percent_table_size>] | vcf <VCF_file> . Sets up, show, deletes or utilizes a hashed index for searching out IDs. Does not allow use of wild cards cf print. If the single argument is a file name, each line in the file is expected to contain a pedigree and individual ID as the first two words, or just an individual ID as the single field. These will be searched for in the current dataset. If the file keyword was used, then the file name can be followed by one or column names specifying which contains the individual ID, or pedigree and individual IDs. Increasing the plevel gives lists of unmatched and matched IDs. If two arguments are given, these are taken to be a pedigree and individual ID to be searched for, while id followed by a string is taken to be an individual ID to be searched for. The show, delete, and size are for tuning, and will not be needed for ordinary use. The vcf option tabulates the overlap in IDs between the present dataset and a VCF file
201. print [where] <a logical expression>. Print trait values for individuals, with a combination of trait values meeting the criterion.
202. print ped <Ped1>...<PedN> [id <Id1>...<IdN>] Print trait values for individuals, with specified combination of pedigree and individuals IDs. The pedigree and ID names can contain wildcard characters: "." (match any character at that position in the search string) and "*" (match zero or more characters).
203. write [<pedigree file name>]. Writes a GAS type pedigree file from the current dataset. Default is to screen.
204. head|tail [map|loci|kinships] [<nrec>| (<skip> <nrec>)]. Writes the first or last nrec records (default 10) of the current dataset, loci, marker map or kinship matrix to the screen. If two arguments are present, then the first represents the number of records to skip over before writing nrec records.
205. more [<nrec>]. Pages through the current dataset, nrec records (default 20) per page. Allows paging backwards and forwards by full or half pages.
206. write pap. Writes the required pedigree files trip.dat and phen.dat (note that you may have to sort trip.dat).
207. write bin <pedigree file name> [compress]. Writes a Sib-pair "binary" pedigree file. The file actually contains both locus descriptions and pedigree/genotype/phenotype data, and so saves the state of the program at that time. An image of Scheme's memory and the twinship and sex marker indicators are also saved. These can be large files, but on systems where gzip is available, will be compressed if the compress modifier is present. Such files are automatically decompressed by the read bin command (if the filename has a .gz suffix. The default format for the file arises from a Fortran unformatted write of the locus and pedigree arrays, and so will be compiler and platform specific.
208. write <format> <pedigree file name> <modifier>

     write	pedigree|gas	<pedigree file name>	[header]
     	arl		<population>
     	asp|tcl
     	beagle		[fou|tri]
     	crimap|tcl
     	csv		[<delimiter> [<missing_value_token>]
     	dot		[[<trait>] [(<marker>|<trait>) [<colour-y> [<colour-n> [<colour-x> [<marker-background-colour>]]]]]]
     	fimpute		[chip <chip_variable>] [<filter_trait>]
     	fisher
     	gda		[all]
     	haploview
     	linkage|ppd|gh		[dummy] [numbered_alleles]
     	mendel
     	merlin
     	phe
     	plink		[<trait>]
     	sage
     	sas
     	sib-pair		[header]
     	roadtrips
     	snp
     	snap
     	solar		[phe] [nopedigreeID]
     	structure		[fou]
     	vcf		[ped_id] [<trait>]

     Use of the keywords pedigree or gas writes a GAS type pedigree file from the current dataset. Quantitative values are written as F9.x or F8.4. If the modifier header is present, a line containing the names of all the variables is prepended. The keyword dot writes a script that will draw all the active pedigrees using the dot program as a marriage node graph (a binary trait locus and one marker locus can be represented by colour and the written genotype in the node respectively); the keyword gda writes a GDA Nexus datafile containing all current marker genotypes for founders. If the keyword all is added, nonfounders will be included as well, but the "gdatype" format will not differentiate between relatives. Similarly, arlequin writes a data file for the program Arlequin containing genotypes from all active individuals; the data may be separated into populations, based on the categories of a specified population membership indicator trait. The keywords linkage and ppd write a pedigree file from the current dataset suitable for use by the LINKAGE (and FASTLINK) programs, the former type requires preprocessing by the Makeped program (note that if a quantitative trait value is zero -- that is nonmissing -- it is recoded to 0.0001); aspex (or tcl) writes a linkage style pedigree file but with the marker locus names as the first line, as the ASPEX programs prefer; gh writes a linkage style pedigree file with a dummy affection trait as the first trait and all the quantitative traits last, with "-" for missing quantitative trait values. The dummy option added to linkage or gh writes a dummy affection locus as the the first trait (everybody affected). The numbered_allele option skips recoding alleles to numbered alleles. The haploview option is a linkage style file with recoding of letter alleles from "ACGT" to "1234". The sage keyword writes a pedigree file from the current dataset suitable for use by the program FSP included in the SAGE package; mendel writes a pedigree file from the current dataset suitable for use by the programs MENDEL or SIMWALK; merlin writes a pedigree file suitable for Merlin -- actually a LINKAGE "pre" format file with zygosity included as the first trait, if the "set twin" command has been previously issued; fisher writes a pedigree file from the current dataset suitable for use by the program FISHER; cri writes the ".gen" style file required by CRI-MAP; phe writes the "pheno.dat" style file required by Mapmaker-Sibs; both csv and solar give a comma delimited file, with header naming columns, from which the pedigree ID column can be dropped via the nop option, and the SOLAR phenotype (or genotype, depending on a prior keep statement) file written by the phe option. The SOLAR pedigree file has two additional columns: MZ twin indicator (requiring a previous "set twin") and a household (actually pedigree) indicator. The sas command writes a datastep with the pedigree data as inline "cards". The structure, roa, and beagle commands write genotype data files for Structure, ROADTRIPS, and Beagle respectively

     (and can be restricted to writing just founder data using the founders option). The plink option writes PLINK format .bed, .bim and .fam files. If a (trailing) trait name is given, this is written as the trait in the .fam file. The snp modifier writes a PLINK .tped file (row-major pedigree file), while snap writes a similar file for WOMBAT. A VCF file containing the genotypes of the currently active markers is written when the vcf modifier is present. The ped_id modifier tells Sib-pair to write a combination pedigree and individual ID, where the VCF IDs take the form <pedigreeID>_<individualID>. If a trait name is given, then eonly individuals where this trait is non-missing will have data written to the output file. The fimpute writes genotype data files for the FImpute imputation program of Sargolzaei et al. The chip modifier gives the trait identifying which SNP array (or genotyping method) was used for that individual. A trait name at the end of the command restricts the outputted data to those where the trait is non-missing.
209. write map mendel|merlin|loki <map file name>. Writes out the map file required by MENDEL 4.0, MERLIN or LOKI.
210. write locus pap. Writes the required locus files header.dat and popln.dat.

211. write locus	aspex|tcl	<locus file name>
     	beagle
     	fisher
     	fimpute		[chip< chip>]
     	gas
     	haploview
     	linkage|gh		[dummy] [xlinked]
     	loki
     	mendel		[trait]
     	merlin
     	relpair
     	sage
     	sib-pair		[<pedigree file name>]
     	structure		<pedigree file name>
     	superlink		[dummy] [xlinked]

     Use of the keyword gas writes a GAS type locus file from the current dataset; haploview writes an "info" file, giving the marker position as the first word of the annotation if it is numerical, or the map position multiplied by 10⁶; linkage writes a locus file from the current dataset suitable for use by the LINKAGE (and FASTLINK) programs; gh writes the same as linkage save that map distances are in cM. The dummy option is used when the first trait is a dummy trait generated by write linkage <file> dummy, while the xlinked option declares the markers to be all X-linked. The keyword loki writes a control file for LOKI's prep program; sage writes a locus file from the current dataset suitable for use by the program FSP included in the SAGE package; mendel writes a locus file from the current dataset suitable for use by the programs MENDEL or SIMWALK (with binary traits defaulting to a factor, but given as a diallelic locus if the trait modifier is present); fisher writes a locus file from the current dataset suitable for use by the program FISHER; merlin for MERLIN; tcl or aspex writes the tcl command file required by ASPEX programs such as SIB_PHASE; sib-pair writes a Sib-pair style script; superlink writes a LINKAGE style locus file modified for the SUPERGH option of SUPERLINK; relpair writes the RELPAIR format (modified from that for MENDEL): it infers chromosome number from the map position (as multiples of 1000) or from the locus name (if it takes the form of "DxSxxx"); beagle writes a Beagle marker list (marker name, position in bp, allele names). fimpute writes a FImpute locus file. The chip keyword preceded the trait indicating the SNP array or genotyping method for each marker.
212. write var [mendel] <var file name>. Writes out the var file (list of quantitative traits) required by MENDEL.
213. write grm <GRM file name prefix> [<filter_trait>]. Writes out the current dataset kinship matrix in the GCTA GRM binary file format (actually 3 files, one with IDs, N used to estimate each kinship, and the actual kinships). An indicator trait can be used filter this so this restricted to a subset.
214. read grm <GRM file name prefix>. Read in kinships from a GCTA binary format GRM matrix binary file format that match the IDs of the current Sib-pair dataset. GCTA actually provides 3 files, one with IDs (pedigree, individual), N used to estimate each kinship (which Sib-pair discards), and the actual kinships. These replace the current dataset "big" kinship matrix.
215. generations [<quantitative trait> [reverse]]. List founders/marry-ins and sibships by generation number for all pedigrees, (over)writing the generation number to a quantitative trait if requested. If the reverse modifier is present, the generation number counts up from the bottom of the pedigree, rather than from the top.
216. loops [<binary trait>]. Prints marital or inbreeding loops in the active pedigrees. If a binary trait is specified, the members of the loop are flagged as "y" at that trait, with nonmembers set to missing at the trait.
217. gpe <codominant marker> [mcmc] [<allele dose estimate>]. Gives iterative peeling or MCMC (Monte-Carlo Markov Chain) estimates of the genotype probability estimates for the given marker for each individual: a vector of probabilities corresponding to the possible genotypes at that marker. For an observed genotype, this is 100% for the observed value and 0% for all other possible genotype values. For an unobserved genotype, it gives the probability distribution of possible genotypes conditional on the sample allele frequencies (assuming Hardy-Weinberg Equilibrium) and the observed genotypes in the individual's relatives. If the name of an extant quantitative trait is appended to the command, the expected gene dose for the first (ie lowest in the collation order) allele will be written to this variable. The gpe command respects the set frequencies command, so that the population allele frequencies can be specified in advance.
218. peel <codominant marker>. Calculate the pedigree likelihood for the given marker.
219. haplotypes <marker1> <marker2> <newmarker> [<threshold>]. This infers phased genotypes when the two markers are in complete (or near-complete) LD. The threshold sets the maximum number of the rare haplotype that is acceptable when LD is not complete.
220. triads This routine lists haplotypes inferred from fully typed parent-offspring triads, along with counts of obligate recombinants.
221. relatives <ped> <id>. This routine lists relatives of an index individual: parents, sibs, spouses, offpsring and descendants. If the pedigree is small (<12 members) or the plevel is set to 1 or higher, then a list of shortest paths from each pedigree member to the index individual is also printed.
222. ancestors <binary trait> |(<quantitative trait> >|>=|<|<=|==|^= < threshold>). This prints the IDs of the ancestor (and ancestral mating) shared by the greatest possible number of probands in a family. The mean intrafamilial inbreeding coefficient for the probands is also output.
223. typed [<binary trait>|<quantitative trait>]. Prints number of individuals phenotyped at each locus. If a stratifying variable is specified, these counts are versus each level of the trait.
224. frequencies|describe [[<codominant marker>| <binary trait>| <quantitative trait>]...[to]...<trait>] | [snp|polychoric] . Print allele frequencies for marker loci, segregation ratios for a binary trait, relative-pair pairwise agreement for a categorical trait (Cohen kappa), polychoric correlation for an ordinal trait, or means, variances, familial correlations and a sibship variance test for a quantitative trait. Default is to describe all loci. The snp option prints minor allele frequencies and number typed for all diallelic marker loci.
225. mcfrequency <codominant marker> | $m. Print MCEM (Monte-Carlo Expectation-Maximization) estimates of the founder allele frequencies for marker loci. A fixed number of EM iterations are carried out, usually 20. This can be set higher if desired using the set emit command.
226. count [where] <a logical expression>. Count individuals, and sibships and pedigrees containing such individuals, with a combination of trait values meeting the criterion.
227. print [where] <a logical expression>. Print phenotype data for individuals with a combination of trait values meeting the criterion.
228. print [where] <a logical expression>. Print phenotype data for individuals with a combination of trait values meeting the criterion.
229. hist <quantitative trait> [<number_of_bins>]. Produce Alternative interface to mix with one distribution. The number of bins in the lineprinter histogram can be set.
230. plot <quantitative trait> <quantitative trait> [<categorical trait>] [<file>]. Produce an Encapsulated Postscript scatterplot for two quantitative traits. If a third binary or quantitative trait is specified, then this controls the plot symbol used for each point (10 different symbols are available for the trait values 1 to 10). Graphic file name defaults to "sib-pair.eps".
231. tabulate [ped <trait>) | ((([ordered] [showmiss] [polychoric] [sampleweights <trait>]) | allelic) <trait 1>...[<trait N>])]. Print contingency table for one, two or N traits, along with contingency chi-squares, Kruskal-Wallis test or odds ratio if appropriate. For RxC contingency tables where the second variable is a diallelic marker locus, allele frequencies and exact P-values testing Hardy-Weinberg Equilibrium are printed for each level of the first trait. If the showmiss keyword is present, then an additional level indicating that the trait value is missing is included in the table, while the polychoric option will calculate the polychoric correlation for a two-way table. One can specify for a oneway table to be printed for each pedigree using the pedigree option, which is supplemented by the Tarone score test for extrabinomial variation if the trait is binary. The default behaviour is to print a one-way table for each active variable. The allelic option requires the last locus to be a marker, in which case it tabulates allele frequency BLUEs versus the other strata variables.
The ordered option produces a table ordered by cell count rather than variable value collation order.
232. llm [(] <trait 1> [[+] <trait 2>...[)^<N>]] [[+] <trait 1> * <trait 2>...][-1] [allelic]. Carry out log-linear modelling of a multidimensional contingency table under the specified model. The allelic option causes an allelic model to be fitted for genotypes of the first listed marker in the model. Formula exponentiation expands to all interactions up to the level of the exponents for the selected terms. The "lrt" command can be used to compare sequentially fitted models.
233. kruskal-wallis <quantitative trait> <trait>. Print table of means for the quantitative trait for each level of factor, along with the Kruskal-Wallis chi-square.
234. means|correlations [<trait 1> ...[<trait N>]]. Calculates means, standard deviation and correlation matrix for a list of traits.
235. pca [<quantitative trait 1> <quantitative trait 2> [...<quantitative trait N>]]. Carries out principal components analysis for the listed traits. If a single integer is given, this allows entry of a covariance matrix for that number of variables via the keyboard.
236. regress <ytrait> = <x1>...[to]... <xN> [offset <offset>] [poisson|(exponential|weibull|evd [<censoring_trait>)] [shape <shape>] [sim] [rep <nreplicates>]. Performs linear or logistic or poisson or weibull or EVD regression of trait ytrait on set of loci x1...xN. If an x variable is a marker genotype, that independent variable is the mean allele size in the genotype, with the exception of the first marker locus encountered in the list, which is fully allelic effect coded. The offset option reads an offset for the linear predictor from the specified trait. Addition of a binary trait name to the end of the keyword list when the regression is weibull or exponential declares this as the censoring indicator. The shape keyword declares a starting value for the solution of the Weibull distribution shape parameter. The sim keyword gives a gene-dropped P-value for the first marker locus in the list. The rep keyword specifies a number of replicates for multiple imputation of the test marker locus genotypes, and is usually used when set analysis imputed has already been issued.
237. clreg <ytrait> = <x1>...[to]... <xN> [ped|stratum <stratum_indicator>]. Performs conditional logistic regression of trait ytrait on set of loci x1...xN. If an x variable is a marker genotype, that independent variable is the mean allele size in the genotype, with the exception of the first marker locus encountered in the list, which is fully allelic effect coded. The default stratification is on sibship, but the addition of the ped keyword gives an analysis stratified on pedigree.
238. mixture <quantitative trait> [<Number of distributions> [normal|pooled_normal|exponential|poisson]]. Estimate mixing proportions, means and standard deviations for a 1..5 component mixture model describing the specified quantititative trait. The default is a mixture of Normal (Gaussian) distributions with different means and variances, but a common variance can alternatively be specified. Other distributions available are the exponential and Poisson. A line-printer type histogram is produced.
239. kinship [inbreeding [mc] | pairwise | dominance | roadtrips | <binary trait> [|<quantitative trait> [>|>=|<|<=|==|^= <threshold>]]. Write the numerator relationship matrix (matrix of coefficients of relationship) for each pedigree in a lower triangular form or as a list of pairs (in the latter case, the coefficient of fraternity is also printed, along with an indicator as to whether a pair are full sibs. When the dominance keyword is present, a pairwise list is printed of bilineally-related pairs (defined as K>0) that are not full sibs. Alternatively, if requested, print a list of individuals with a non-zero inbreeding coefficient, using a Monte Carlo estimator if mc modifier is present. If a binary trait is specified, the NRM is only for the affecteds if plevel=1; for plevel=0, only a summary for each pedigree is printed: number of affecteds, number of "sporadic" cases ie cases unrelated to any other affected family members (eg marry-ins with no affected descendants), mean coefficients of relationship for affected relative pairs and of inbreeding for cases, and permutation test of significance of observed mean coefficient. For use with the ROADTRIPS programs, a pairs format where the kinship, rather than coefficient of relationship, and inbreeding coefficients are printed for self-pairs (rather than 1+F) can be printed.
240. ibd <codominant marker> i [ ... <codominant marker N>][pairwise]. Write the estimated mean identity-by-descent sharing at a marker or set of closely linked markers for all relative pairs in each pedigree as a lower triangular matrix or a list of pairs.
241. ibs [ibd | moment | <binary trait> | <quantitative trait> >|>=|<|<=|==|^= <threshold>]. Prints the estimated mean identity-by-state sharing at all active markers for all pairs in the dataset as a list of pairs, or if the ibd modifier is present, then the MLEs for the kinship coefficients. The moment modifier gives the shrinkage estimate of the kinship matrix of Endelman and Jannink [2012]. If a binary trait is specified, sharing is only calculated for pairs where both are affected. Ungenotyped individuals are skipped completely, but genotyped individuals sharing no markers with any other individual will be evaluated (so the mean IBS will be printed as missing).
242. set kinship A|C|ridge_constant <constant>. A kinship matrix for the entire dataset (stored internally as matrix kinmat) can be specified by this command, or read in from an external file. While A and C set it to the pedigree-based NRM or a block diagonal family environment correlation matrix, the ridge_constant option can be used to modify an existing matrix.
243. hwe [founders] [<mar1> ..[to].. <marN>] [$(m | x)]. Prints chi-square statistic for Hardy-Weinberg equilibrium for all marker loci. Analysis may be restricted to founders, and if the marker is diallelic, an exact test is carried out. If nonfounders are included, then a gene-dropping simulated P-value is produced. The mean IBS sharing for all typed matings is also calculated, and compared to its expected value. This latter test may allow detection of homogamy or assortment.
244. cksib. Lists all sib pairs, and the mean of IBS at all marker loci where both members of the pair are typed at the marker, comparing this to that expected if related as specified by the pedigree structure. The output is to the standard output.
245. share [pairs]. Lists all relative pairs, and the mean of IBS at all marker loci where both members of the pair are typed at the marker, comparing this to that expected if related as specified by the pedigree structure, allele frequencies and linkage map. The output is to the standard output. The default lists individuals whose Z score measuring deviation from expected exceeds 1.65 with any other relative. The pairs option prints the statistic for each deviant pair, or all pairs if output is set to verbose.
246. mds <dim1>[... <dimN>]. Performs classical (metric) multidimensional scaling of interindividual genetic distances. These are based on identity by state (IBS) sharing at all active marker loci across all pairs of active individuals in the dataset. Individuals have to be typed at at least 50% of the active markers to be included in the analysis. The number of dimensions output is controlled by the number of quantitative traits named in the command: each will contain the estimated coordinates for every individual on one of the dimensions.
247. mztwin <monozygosity_indicator> |(<zygosity_score> even|odd|(>|>=|<|<=|==|^= <threshold>)) [clean|delete|find]. Using a binary or quantitative trait which indicates which sib pairs are monozygotic twin pairs, list markers at which the twins carry discordant genotypes. Gives proportion discordance for each marker. This is useful for estimating genotyping error rates. The clean option deletes genotypes for pairs where there is an inconsistency, and fills in missing genotypes where that for the cotwin is available. The delete option drops the member of the pair with the fewest nonmissing phenotypes, and averages (across the pair) quantitative phenotypes where both are observed. The find option uses the currently active markers to identify likely MZ twin pairs as pairs of relatives (in the same pedigree) sharing 99.5% or greater marker concordance. A value indicating membership of an MZ twin pair will be written to the currently active twin indicator trait or to the named variable. The even and odd zygosity score test require the score to be positive, and are designed for use with MERLIN's coding system.
248. davie <binary trait> [<proband indicator>]. Print segregation ratios and standard errors for a binary trait, adjusted for the ascertainment scheme. Probands are indicated by being affected at the proband indicator "locus". If no proband indicator variable is given, complete ascertainment is assumed (equivalent to "davie trait trait").
249. twin <trait> [<zygosity_score> even | odd | ( >|>=|<|<=|==|^= <threshold>)]. Calculate, and test for equality of, MZ and DZ twin correlations or concordances. If a zygosity indicator is not specified, the indicator declared by the "set twin" command is used to divide siblings into MZ twins and non-MZ full siblings (DZ twins and "singletons"). If a zygosity indicator is specified, then this is used to define three groups: MZ twins (indicator expression TRUE or greater than zero), DZ twins (indicator FALSE or equal to zero), and nontwin siblings (indicator value is missing). The Z test for equality of the Pearson correlation coefficients is actually the likelihood ratio test [Brandner 1933]
.
250. kendall <quantitative trait> <censoring trait> [<zygosity_score> even | odd | ( >|>=|<|<=|==|^= <threshold>)]. Bivariate survival analysis for twins. Kendall's tau is used as the measure of association, standard errors are produced by bootstrapping, and a permutation P for the MZ:DZ difference.
251. varcomp|mft <quantitative trait> [ae] [ce] [ace] [ade] [covariate <covariate trait 1> [ ...+ <covariate trait N>]]. Performs MVN or MFT variance components analysis for a quantitative trait using all phenotyped individuals. Currently fits ADE, ACE, AE and CE models. The ae option fits AE and E only, and so forth. Multiple covariates can be included in the fixed effects part of the model via adding the cov keyword and a list of covariates at the end of the command line. Only the first marker locus in the list of covariates is fully allelic effect coded, with subsequent markers included as the mean of their allele values (i.e. 1="1/1", 1.5="1/2", 2="2/2" for a diallelic marker).
252. lrt. Constructs likelihood ratio test comparing the last two variance components (var), log-linear models (llm), generalized linear (reg) or generalized linear mixed (fpm) models fitted. The compared likelihood statistics from fpm are actually the mean model loglikelihoods. This allows one to carry out tests of linkage after conditioning out the effects of genotype at a candidate polymorphism, and traditional "measured genotype" association analysis in pedigrees, including non-normal data (currently binomial and poisson distibuted traits).
253. blup (<quantitative trait> <h2>) | (<codominant marker> [<binary trait> | (<quantitative trait> [(>|>=|<|<=|==|^= <threshold>]))). Calculates BLUPs (best linear unbiased predictions of the breeding value) for a quantitative trait assuming the given heritability, or if a codominant marker is specified, the BLUEs for the marker allele frequencies, either for the entire sample, or for one subgroup.
254. fpm <quantitative trait> [>|>=|<|<=|==|^= <threshold>]|<binary trait>
     [nqtl <number simulated QTLs>]
     [link logit|probit|mft|ln]
     [likelihood_family gaussian|binomial|poisson]
     [p|fre] [a|va] [d|vd] [g|vg|h2] [c|vc] [s|vs]
     [fix a|c|d|e|g|mu|s|var]
     [aval|cval|dval|eval|gval|mu|pval|sval| var|AA|AB|BB|SD <value>]
     [cov <x1> [+ <x2>...]]
     [save <target locus].
     Uses a Markov Chain Monte Carlo (Metropolis) algorithm to perform Generalized Linear Mixed Model analysis (when nqtl set to zero), complex segregation analysis (nqtl set to one), "genetic" mixed model analysis (nqtl set to one, g estimated), or a finite polygenic model analysis for a quantitative or binary trait using all phenotyped individuals. The a, d, g, c and s options respectively include additive QTL effects, dominance QTL effects, additive Gaussian polygenic random effects, pedigree environmental effects and maternally derived sibship effects in the model. The fix option allows that variable to be held fixed. The aval (dval etc) keyword allows the starting value of that parameter to be set. Either a Gaussian, Binomial, Weibull or Poisson likelihood can be used. In addition to the canonical links for these three types, one can also fit the multifactorial threshold model to binary data. Because the identity link function for a binomial likelihood is of genetic interest, notably in the case of the single major locus model, special checks to reject models where predictions lie outside 0-1 are made, so that this model can be successfully fitted. "sav" modifier keyword for "fpm" can save BLUPs or genotypes to a variable.
255. dis [<marker locus 1> [<marker locus 2> [.. <marker locus N> [<trait>]]]. Estimates frequencies of haplotypes based on individuals with two informative parents or no informative parents. For pairs of multiallelic markers (6 or more alleles), only the former group is used. In the informative parent-offspring triads, the loci are assumed to be tightly linked, so that four parental haplotypes may be inferred, but this maybe turned off. If no markers are specified, the sequence of pairs of markers in the list of loci is produced (ie marker1 with marker2, marker2 with marker3...). If one marker is specified, then the pairing of all other markers with this index marker is analysed. For pairs of markers with few alleles, information from both phased and unphased genotypes is combined using an EM-fitted log-linear model, which can also deal with X-linked data. If three or more markers are specified, or two or more markers along with a trait, the genotypes are all treated as if unphased (and X-linked data is not handled). The same log-linear modelling framework is used. When a trait is specified, the haplotype frequencies within each trait category (two if a binary trait, as many levels as observed for a quantitative trait) are calculated, along with a chi-square testing equality of marker haplotype frequencies across categories. If exactly two numbers are given, then these are assumed to be the number of alleles at two loci, and genotype counts may be entered at the command line.
256. neff [<Map window> [<critical P value>]]. Calculates the effective number of association tests using the set of active (autosomal) markers, allowing for the observed intermarker linkage disequilibrium. Implements the approach of Moskvina & Schmidt [2008]. The optional parameters are the width of the sliding map window containing markers to be tested versus the current iterated marker (currently 1 cM/Mbp), and the genome-wide critical P-value (defaults to 0.05).
257. assoc (<binary trait> |<quantitative trait> [(>|>=|<|<=|==|^= <threshold>]) | categorical) [founders] [vcf <VCF_file_name>] [covariate <covariate>] [genotypic] [ibd <ibd_marker>] [snp|freq|maf|risk]. For a binary or dichotomised quantitative trait, prints chi-square statistics for association for all marker loci versus the trait, either affected versus unaffected if the trait is binary, or above or below the threshold if the trait is quantitative. If the categorical modifier is present, then each unique value of a quantitative trait is taken to represent a category of a multinomial trait. For dichotomous traits, a second table in the output shows the results from the reconstruction-combined TDT within informative sibships. Also prints F statistics assuming trait is indicator of membership of different subpopulations. For a quantitative trait, prints the model and residual sums-of-squares and allelic means with naive standard errors from an additive allelic ANOVA model. Monte-Carlo empiric P-values are produced for either analysis via gene-dropping, which can either be unconditional, or completely linked to that of another marker. Analysis may be restricted to founders. Genotypic rather than allelic analyses can also be specified, using the genotype flag, and if the trait prevalence has been specified (by the set prev command), the attributable risks and sibling recurrence risk ratios are also calculated (plevel > 0). Covariates can be added to some analyses. The snp modifier limits the output to an allelic odds ratio and 95% confidence intervals (along with the P-value) for a binary trait, and allelic regression coefficient and SE for quantitative traits; the freq, maf and risk modifiers limit output to case and control frequencies (plus test statistic and P-values), where the frequencies are that of the base allele, minor (in control) allele, and risk allele respectively. The vcf option, by contrast, performs a naive chi-square test comparing case and control allele frequencies to those recorded in a VCF file as INFO variables defaulting to "AC" and "AN". These can be altered by the "set vcf" command
258. mito|yhaplo <quantitative trait> [<marker 1> [...<marker N>]]. Performs ANOVA analysis of the association between a quantitative trait and Y chromosome or mitochondrial haplotypes, defaulting to each eligible marker in turn if markers making up a haplotype are not explicitly listed. Haplotypes are imputed where unequivocal, and Monte-Carlo empiric P-values are produced by gene-dropping.
259. rare (<trait>|<marker>) [<MAF_threshold>]. Pools rare alleles in all currently active markers below a given threshold allele frequency. Default maximum MAF for inclusion is 0.01. If a trait is specified, lists carrier individuals by trait locus status, and provide table of number of alleles carried versus trait locus status, along with a contingency chi-square test. If a marker is specified, the genotypes for this are replaced by dummy genotypes representing the number of rare alleles carried by the individual ("1/1"=0, "1/2"=1, "2/2"=2 or more).
260. skat <trait> [madsen_browning|beta_1_25_weights) Performs SKAT combining all currently active markers. Analysis can be unweighted, or using two different allele frequency based weights.
261. wqls <trait> [kin <kinship_file_name>] [ridge_constant <ridge_constant>]. Calculates the WQLS score test of allelic association for binary, categorical (more than 2 categories), or quantitative traits.
262. mqls <binary trait> [|<quantitative trait> >|>=|<|<=|==|^= <threshold>] [prevalence <prevalence>] [kin <kinship_file_name>] [ridge_constant <ridge_constant>]. [hwe|robust]. Calculates the MQLS, WQLS and modified-chi-square quasi-likelihood score tests of allelic association where the trait is binary, or above or below the given threshold if the trait is quantitative, or categorical (corrected chi-square only). The MQLS test requires an estimate of the trait prevalence, which can be specified by appending a value to the command following the prevalence keyword, or preset using the set prev command.
Rather than using the given pedigree, a kinship matrix can be read from a file, containing one element per line, and fields corresponding to a header line:

ped1 id1 ped2 id2...kin

A ridge constant can be specified to be added to the diagonal of this matrix, and either "naive" (hwe) or robust (robust) forms of the tests used.

263. homoz [<binary trait> [|<quantitative trait> >|>=|<|<=|==|^= <threshold>]]. Prints the asymptotic Z statistic and a one-sided MC P-value for whether homozygosity at each marker locus is increased in probands, either affected if the trait is binary, or above or below the given threshold if the trait is quantitative.
264. multihomoz [(<binary trait> [|<quantitative trait> >|>=|<|<=|==|^= <threshold>]) | save <quantitative trait>]. Prints the asymptotic Z statistic and a one-sided MC P-value for whether the maximum length of runs of homozygosity at marker loci along the specified map is increased in probands, either affected if the trait is binary, or above or below the given threshold if the trait is quantitative. The simulated P-value assumed linkage equilibrium so the markers must be appropriately thinned. If trait is not specified, or the save modifer used, the observed and expected multilocus homozygosity and the runs-of-homozygosity inbreeding coefficient estimator (F_roh) is written for each individual. The latter can be save to a quantitative trait.
265. fstats [<population_indicator> [founders]. Prints the F statistics comparing identity-by-state sharing of alleles within loci, within populations (demes), and between populations. Population membership of an individual is indicated by the value of a specified binary or quantitative trait. There can be more than two populations.
266. tdt <binary trait>|(<quantitative trait> >|>=|<|<=|==|^= <threshold>) [cutoff <cutoff>] [mat|pat]. Prints transmission-disequilibrium statistics for all marker loci versus the trait, where an index person is either affected with a binary trait, or whose value for a quantitative trait exceeds the given threshold. Since binary traits are coded internally as 2=y and 1=n, an analysis using unaffecteds as proband can be performed as tdt <binary trait> under 2. Similarly, in unascertained families, tdt <binary trait> over 0 tests for segregation distortion. Calculation of the TDT statistic can be restricted to pairs of cells whose total is greater than cutoff, eg 5, and to the maternal or paternal contributions, if parent-of-origin effects are suspected.
267. hrr <binary trait>|(<quantitative trait> >|>=|<|<=|==|^= <threshold>). Performs the Haplotype Relative Risk test, comparing case offspring allele frequencies to those in their parents. A gene-dropping P-value is generated, so that it is applicable to general pedigrees.
268. schaid <binary trait> [<marker> [<allele>]]. Performs the Schaid and Sommer [1993] genotypic risk ratio test for familial association under the assumption of Hardy-Weinberg equilibrium, as well as the "Conditional on Parental Genotypes" version that is equivalent to the genotypic TDT. Only one allele (defaulting to the commonest) is tested versus all others, and two penetrance ratios (GRR₂=f₂/f₀ and GRR₁=f₁/f₀) are estimated, along with the LR chi-square test that GRR₂=GRR₁=1.
269. sdt <binary trait> [ped|stratum <stratum_indicator>]. Prints sibship-matched case-control conditional logistic regression results for all marker loci versus the trait, where sibships contain affected and unaffected individuals at the binary trait. If the ped or stratum modifier is present, then the stratification is by pedigree of level of the stratifying variable.
270. stratified_association <trait> [<marker>] <stratum_indicator> Performs fixed and random effects score tests on a specified locus or all marker loci versus the trait, stratifying on the provided stratum indicator.
271. interaction_association <trait> [<marker>] <stratum_indicator> Performs the fixed and random effects score tests on a specified locus or all marker loci versus the trait, stratifying on the provided stratum indicator, but returns the heterogeneity test result as the summary statistic.
272. trend <quantitative trait> <marker> [permute] Prints Jonckheere-Terpstra trend test result for the specified marker loci versus the trait. Most suitable for SNP markers, where genotype ordering is clear. Monte-Carlo empiric P-values are produced either via gene-dropping or permutation.
273. asp <binary trait>|(<quantitative trait> >|>=|<|<=|==|^= <threshold>). Prints IBS-based affected (full and half) sib-pair statistics for all marker loci versus the trait. It also prints the mean IBD sharing for full sibs, along with the exact (binomial) two-tailed P-value for the "mean" test. All possible sib-pairs are used, and are treated as independent.
274. pen <locus1> <locus2>). Performs Penrose [1935, 1937] sib-pair linkage statistics for two loci. Quantitative traits are treated as if categorical. The output consists of a two-by-two table of the sibling concordances at the two loci.
275. apm <binary trait>|(<quantitative trait> >|>=|<|<=|==|^= <threshold>) [ibd|ibs]. Prints APM statistics for all marker loci versus the trait.
276. sibpair|he1|he2|vis <quantitative trait>| <binary trait> [<Weight variable>] [sim] [mean<trait mean>] [var | sd <trait variance or SD>] [cor<trait sibling correlation>]. Performs Haseman-Elston regressions (Sham & Purcell [2000] as the default, but Visscher-Hopper [Visscher & Hopper 2001], traditional and "new" Haseman-Elston also available) for all marker loci versus the trait using full and half-sib relative pairs. The contribution of each pair can be weighted by the mean of their values at a quantitative trait. Empirical P-values can be simulated, if requested. For the S+P regression, the "true" population trait mean, variance and sibling correlation can be specified, to facilitate analysis of selected samples.
277. twopair <quantitative trait>| <binary trait> <marker locus 1> <marker locus 2> <theta12>. Performs Fulker & Cardon's Haseman-Elston interval regression for first and second marker loci versus the trait using full- and half-sib relative pairs. The recombination distance between the markers theta12 must be given. Haseman-Elston regression is performed using ibd estimated at ten points in the interval.
278. qtlpair <quantitative trait> [full] [cqe] [covariate <covariate trait 1>]. Performs variance components linkage analysis for all marker loci versus the trait using full-sib relative pairs, or if the full option is active, all genotyped individuals. Both the polygenic background and the QTL are modelled as additive genetic. Covariates, which can include codominant marker loci, are added using as cov keyword-trait pairs. Only the first marker locus is fully allelic effect coded, with subsequent markers included as the mean of their allele values (i.e. 1="1/1", 1.5="1/2", 2="2/2" for a diallelic marker).
279. linkage [<marker locus 1> [<marker locus 2>]]. Performs Elston and Keats sib pair linkage analysis for codominant markers. Default is adjacent pairs of markers (ie marker 1 with marker 2, marker 2 with marker 3...). If one marker is named, then gives estimate of recombination distance to all other markers.
280. lod <marker locus 1> <marker locus 2> [<recombination fraction>]. Performs two-point lod score linkage analysis for codominant markers. The algorithm is fairly slow in the presence of many missing genotypes.
281. summary [<N highest results> | plot [qq] [<EPS file>] | table | dump [<output file>]| get <variable_name>]. List the N (defaulting to 5) most significant P-values from the last linkage or association analysis. If the plot keyword is issued, a Postscript plot of the -log₁₀ P-values versus map position of the active marker loci tested for association or linkage in the last analysis, unless the qq modifier keyword is present. In that case, a quantile-quantile plot of the -log₁₀ P-values is produced. This is saved to a file (default name "sib-pair.eps"). For the table modifier keyword, a table of binned P-values or test statistics is produced, usual bin sizes are powers of ten. The dump modifier prints all the P-values to a file. The get <varname> modifier extracts numerical values from locus annotations taking the form varname=value, or from the i'th column.
282. summary combine [<marker locus 1> [... <marker locus N>]]. Combine P-values from test statistics for the specified loci using the Fisher method, and new Cauchy method of Liu and Xie [2018].

The following script performs a number of analyses on a dataset containing four loci.

Test.in

set work c:\tmp\
set weight founders
set out verbose
set impute on
set locus quant quantitative 
set locus trait affection 
set locus marker1 marker 
set locus marker2 nam
read pedigree test.ped
run 
freq
mix quant 2
assoc trait
tdt trait 
ass quant
sibpair quant
recode marker1 126 128 to 999 
freq marker1
tdt trait
apm trait
sibpair quant
! create a new trait
set locus new_quant qua
if (quant le 0) then new_quant= -sqrt(-2*quant) else new_quant=log(quant)
! adjust the binned allele sizes of marker
if (marker1 ne 999) then marker1=marker1+1
drop trait quant
write pedigree testout.ped

DATASETS

pedigree-id person-id father-id mother-id sex-of-person locus-value-1...locus-value- N

A pedigree ID may be up to 20 alphanumeric characters, and an individual's personal ID up to 14 characters. Missing values are denoted x (or .), and represented internally as a trait value of -9999. Locus values for a binary trait are y (expresses trait), n (does not express trait). Sex takes the values m (male) and f (female), and may be missing. Alleles at a marker locus are integers between 1 and 999 or single letters. Slashes dividing alleles of a genotype are optional. A pedigree file may contain a comment at any time, prefaced by ! or #, and may contain a locus header of the form (though this has no function and is included to allow compatibility with the GAS pedigree format):

pedigree locus <locus-name-1>...<locus-name-N>.

If only one parent of an individual is specified in the pedigree file, a dummy record and ID number for the other parent is generated by the program.

Here is the data set analysed by the script test.in:

Test.ped

! test pedigrees including one halfsib in pedigree 1000
!                                          Marker 1  Marker 2
! The seven mating types:  1000-1 x 1000-2 Type VII  Type III
!                          1000-1 x 1000-3 Type VI   Type II
!                          1001-1 x 1001-2 Type IV   Type V
!                          
1000 1   x   x   m   10  y   126/132   1/1
1000 2   x   x   f   10  n   128/130   1/2
1000 3   x   x   f   25  n   128/132   2/2
1000 4   1   2   f   20  y   126/128   1/1
1000 5   1   2   m   30  y   130/132   1/1
1000 6   1   2   m   40  n   128/132   1/2
1000 7   1   2   f   50  n   126/130   1/2
1000 8   1   3   f   60  n   126/128   1/2
1000 9   1   3   m   40  y   132/132   1/2
1001 1   x   x   m   20  y   124/124   1/2
1001 2   x   x   f   30  n   126/128   1/2
1001 3   1   2   f   40  n   124/128   1/1
1001 4   1   2   m   30  n   124/126   1/2
1001 5   1   2   m   40  n   124/126   2/2
1001 6   1   2   m   40  n   124/128   1/2
1002 1   x   x   m   10  y     x/x     x/x
1002 2   x   x   f   40  n     x/x     x/x
1002 3   1   2   m   30  n   126/126   1/2
1002 4   1   2   m   60  n   126/126   1/1
1003 1   x   x   m   20  n   126/126   1/2
1003 2   x   x   f   25  n     x/x     x/x
1003 3   1   2   m   40  n   126/126   1/2
1003 4   1   2   m   15  y   126/126   1/2
! end-of-pedigrees

The pedigree file written by Sib-pair contains the original records plus any additional dummy records for missing parents of nonfounders. It is sorted by founder/nonfounder status, generation number, and the collation order of the parental IDs and individual ID.

TIPS AND TRICKS

How do I?

• Get data into Sib-pair: Sib-pair can now read a few file formats other than its native formats (notably PLINK, MERLIN, pre- and post-Makeped Linkage format). Sib-pair cannot handle multiple datasets simultaneously, though you can use the "update" command to merge new genotype data into your existing pedigree for example. I have also written a number of external utilities to do that kind of work -- see the website for programs such as mergeped. Often I do preprocessing in R.
• Get more output: Set the print level higher.

  #
  # Get full TDT tables but summary results for association analysis
  #
  set ple 1
  tdt trait
  set ple 0
  ass trait
  

• Analyse only selected traits: Use the drop then undrop commands to specify loci or ranges of loci to be included or exclude from particular analyses:

  # Drop out year of birth and don't show any allele frequencies
  drop yob $m
  describe
  undrop
  

• Test for duplicate individuals Several different commands are useful. The test duplicates option is the most general, but one can limit the individuals to be compared by either test ped1 id1 ped2 id2 or mztwin find.
• Write out only the markers: Use keep $m, then write pedigree. Note that this will not return previously dropped markers to the analysis (to this, you would need to first undrop.
• Delete marker data for particular individuals: Use keep $m, then delete <pedigree> <person>, followed by und.

  ...
  set loc errprob qua
  keep $m
  delete where errprob>0.9
  und
  

• Use the parser: This hopefully reduces the number of steps in data manipulation required for a complete analysis.

  #
  # Create a new variable that is a function of
  # three existing quantitative variables
  #
  set loc b1 aff
  set loc q1 qua
  set loc q2 qua
  set loc q3 qua
  read ped ex.ped
  run
  set loc new_var qua
  new_var=log(q1+q2)/q3^2
  if (male) then new_var=new_var+10
  #
  # Select a subset of pedigrees where two or more probands
  # meeting multiple criteria
  #
  select containing 2 where new_var>35 and q1 le q2 and isnon
  write newped.ped
  

• Use wildcard selection: This allows selection on pedigree or person names in a flexible fashion:
  

  #
  # Select first six pedigrees in file
  #
  >> select famnum<=6
  >> wri
  

  ! Pedigree  Person   Father   Mother  x
  !
  a                 a        x        x x
  and             and        x        x x
  are             are        x        x x
  as               as        x        x x
  at               at        x        x x
  be               be        x        x x
  

  
  print ped * id a.
  
  id=as-as sex=x
  id=at-at sex=x
  
  >> print ped a*e
  id=are-are sex=x
  
  >>print ped *s*
  id=as-as sex=x
  

• Use variable names in formulae when the variable name shares the first three letters with a command eg "trait" and "transform": surround the variable name with brackets, or precede the command with let.
• Analyse multiple pedigree files: Jobs can be stacked, providing each begins with the clear command. The loci will have to be declared each time however.
• Delete a marker that is giving too many mendelian inconsistencies: Use the drop command on that marker before the run command.
• Ignore error messages from a marker that is giving too many mendelian inconsistencies: Drop the marker out before error checking as before, then use the undrop command before the hopefully robust type of analysis chosen. Alternatively set checking off and set impute -1 will turn off checking for all markers.
• Test for segregation distortion: Do a TDT with everyone affected, for example,

  set loc dummy aff
  dummy=y
  tdt dummy

• Log transform a quantitative trait: Use tra trait 1 0 0 or trait=log(trait) (slower).
• Get a histogram for a quantitative trait: Use the hist command (this is a synonym for mixture trait 1). Setting the print level higher for mixture also prints out posterior probabilities of membership of the different distributions -- useful for choosing thresholds.
• Print genotype frequencies: Set plevel to 1 so that the full table is printed when the hwe command is used.
• Remove unrelated individuals or nuclear families that have become disconnected from the main pedigree, although they have the same pedigree ID: Use the subped command, followed by select num > 20, or some other suitable number that will keep only the main pedigree.
• Do a multipoint ASP linkage analysis: Write the appropriate format pedigree and locus file, and call another program like SIB_PHASE:

  # Write a locus file
  write locus aspex batch.tcl
  # Write out the pedigrees as nuclear families (if multigenerational)
  nuclear
  write aspex batch.ped
  # The resulting output will be included in the Sib-pair output
  $ sib_phase -f batch.tcl batch.ped
  $ rm batch.tcl batch.ped
  

• Which TDT result should I trust? The genotypic TDT is currently a more experimental test. In the case of nuclear families with typed parents, it reduces to a simulation based CPG GRR test. In larger pedigrees with multiple generations of affecteds, matings must have both parents genotyped before they are included in the simulation, but this is done over the complete pedigree.
• What does the "founder" option for the allele frequencies give (updated)? Since a simple count of alleles in typed founders can miss alleles segregating in the pedigree (and thus inherited from untyped founders), I have provided a method that enumerates all alleles in each pedigree, but weights the contribution of the pedigree by the number of founders it carries (that is, the number of representatives from the population whose allele frequencies one is trying to estimate). So, in a simple example,

  1 1   x x               1 1   x x                
  |     |                 |     |
  +--+--+                 +--+--+            
   |                       |
  +-----+-----+                 |
  1 2   1 2   1 2               1 3
  

  the contributions from each family would be weighed equally, as each contains two founders. For allele 1 for example,

  Family 1      Family 2
  2 * 5/8  +  2 * 3/4       
  --------------------  =   11/16
  4 (total founders)
  

  	Allele 1	Allele 2	Allele 3
  The naive estimate is:	8/12 (.67)	3/12 (.25)	1/12 (.08)
  The weighted estimate:	11/16 (.69)	3/16 (.19)	2/16 (.13)
  The imputation estimate:	6/8 (.75)	1/8 (.12)	1/8 (.13)
  The MLEs (MENDEL USERM13):	.6254	.2182	.1564
  The MCEM MLEs (Sib-pair mcf):	.6250	.2201	.1549
  The BLUEs (Sib-pair blu):	.6818	.2045	.1136

  In this example, the frequencies of the 2 and 3 alleles are better estimated by the weighted method than the naive method. In the imputation estimation approach, the untyped parents were imputed as 1/2 and 1/3 respectively. The MCEM estimate gives the MLEs within the limits of stochastic error. In general, providing there are enough pedigrees, the naive estimate is as good as any.

• What do I do if I have covariates or liability classes? (updated) Several quantitative trait analyses allow the addition of covariates to the analysis. Eventually the various binary trait analyses Sib-pair does will allow for covariates. At the moment, the best thing you can do is create multiple phenotypes eg male diabetes and female diabetes, with the phenotype set to missing appropriately. Then one can do the analysis within each stratum, and pool test statistics in various ways used in meta-analysis. In the case of quantitative traits, analysis of residuals formed by adjusting for covariates will carry you a long way. Using the "set liab" allows you to write liability classes to programs such as MLINK.
• Why can't I have variables called d1 or e1 etc? Following Fortran conventions, d1 and e1 are read to mean 0d1 and 0e1 and evaluate to a real number: 0. You can have a1, f1, 1d, dd1 etc as names.
• Does the fpm command actually work? It does seem to in examples, but multiple runs should be performed, and those with high coefficient of variation of log likelihood looked at with a jaundiced eye. I have applied it now to a number of nongenetic generalized linear mixed model problem datasets, where it seems to give answers that roughly agree with those from other programs ;).

DOCUMENTATION OF ROUTINES

Regression (multiple) is performed by AS (Applied Statistics) 164 [Stirling 1981], which uses modified Givens rotations to perform weighted least squares regression including linear constraints. It is also used to give generalised linear models (poisson and binomial regression) via IRLS. The random number generator is the well known AS 183 [Wichmann and Hill 1982]. The approximate randomization routine is styled after general templates described by Noreen [1989]. Mixtures of distributions are fitted using AS 203. Various standard distributions are evaluated using AS 3, 66, [111, retired], 241, 275. Matrix inversion for likelihood evaluation of the variance components model is performed using the LINPACK [Dongarra et al 1979] routines DGEDI and DGEFA (which replace AS6 and AS7 as of April 2008). Extraction of eigenvectors and eigenvalues is done using the appropriate EISPACK routine (RS). Likelihood maximization is performed using either AS319 ("varmet") [Koval 1997] or BOBYQA [Powell 2009], both of which seems to do a good job of it. Conditional logistic regression is performed using AS162 [Smith et al 1981]. Evaluation of multivariate normal cumulative probabilities uses routines from MVNDSTPACK [Genz 1992], or alternatively from TOMS717 [Bunch et al 1993]

The Scheme interpreter is a port of the Minischeme and Tinyscheme interpreters. Minischeme is the work of Atsushi Moriwaki and Akira Kida, following Matsuda and Saigo (Programming of LISP, archive No 5 1987, 6-42). Tinyscheme is based on Minischeme and is more fully featured. It was written by Dimitrios Souflis (dsouflis@acm.org). The Sib-pair Scheme (like Minischeme) does not support vectors, bignums and complex numbers. It does offers double precision floating point and 64 bit integer arithmetic, and provides most of the standard arithmetic and string handling functions.

LIMITATIONS

The Fortran 77 old versions of Sib-pair had limits on pedigree size, numbers of loci and allele numbers per loci. The new version uses allocatable arrays, and so is limited only by memory. Note that rewrites in 2012 allow work data to be stored as files, so one can analyse datasets that do not fit into available memory.

Locus names, and pedigree and individual identifiers are restricted to a maximum length (20, 20, 14 characters respectively). Locus annotations are also restricted to 40 characters each. These can be increased by recompiling after changing the constants in module idstring_widths or locstring_widths (this may make old binary images unreadable).

The Monte-Carlo based routines are computationally intensive. Generally speaking 200-300 iterations of such a routine are sufficient to give a good estimate of a mean or variance (as in the apm routine), but 1000 iterations or more are advised for an accurate P-value. Using "set iter 0" will provide only the parametric estimators, e.g. for screening purposes. The MCMC algorithms need even longer runs to be ensure that they have reached stationarity.

There are only a few multipoint procedures: diseq, qtl, multihomoz and haplotype. The twopoint command has not been ported (yet) to the Fortran 95 Sib-pair.

COMPILATION

The program is (fairly) standard Fortran 95, and has run successfully on PCs (under DOS, NT or Linux), SUN Sparcstations, DEC Alpha, Macintoshs running OSX, and HP9000 workstations. The code compiles using g95, ifort, gfortran, and sun f95, and currently I distribute Windows and Linux binaries. The code has sometimes encountered problems compiling using the Linux ifort with optimization (-O) on.

I have done some comparisons of different Fortran compilers.

Job 0:
read locus plink hapmap3_r1_b36_fwd.CEU.qc.poly.recode.mapfile
read pedigree hapmap3_r1_b36_fwd.CEU.qc.poly.recode.ped
set imp -1; set che off
run

Job 1:
read locus plink hapmap3_r1_b36_fwd.CEU.qc.poly.recode.mapfile
read pedigree hapmap3_r1_b36_fwd.CEU.qc.poly.recode.ped
run

Job 2:
read locus plink hapmap3_r1_b36_fwd.CEU.qc.poly.recode.mapfile.gz
read pedigree hapmap3_r1_b36_fwd.CEU.qc.poly.recode.ped.gz
run

Program	Job 0	Job 1	Job 2
i7
Sib-pair gfortran -O2	6.4m	6m24s	6m38s/6m45s
Sib-pair g95 -O2	11m39s	-	11m57s
Sib-pair sun f95 -O2	3m15s	-	3m36s/3.7m
gunzip	19.9s	-	-
plink	2m15s	-	-
merlin	-	3m50s	3m46s

Running the standard testsuite (testsuite.in) gives:

64bit
g95	8.52 s
gfortran 4.4.3	6.64 s
Pathscale 4.0.10 gcc version	6.64 s
32bit
g95	11.20 s
gfortran 4.6.0	7.38 s
gfortran 4.4.3	8.56 s
open64	9.32 s
sunf95	9.62 s

ACKNOWLEDGEMENTS

This program was developed while the author was an Australian National Health and Medical Research Council Neil Hamilton Fairley Postdoctoral Fellow and later a Research Fellow. This included a period working in the Genetic Epidemiology Division of the Johns Hopkins University School of Public Health and in the Epidemiology Unit at the Queensland Institute of Medical Research.

REFERENCES

• Aitkin MA, Clayton D (1980): The fitting of exponential, Weibull and extreme value distributions to complex censored survival data using GLIM. Appl Statist 29: 156-163.
• Almasy L, Blangero J (1998): Multipoint quantitative trait linkage analysis in general pedigrees. Am J Hum Genet 62: 1198-1211.
• Andersen PK, Borgan O, Gill RD, Keiding N (1993): Statistical models based on counting processes. New York: Springer Verlag.
• Besag J, Clifford P (1991): Sequential Monte Carlo p-values. Biometrika 78: 301-304.
• Bishop DT, Williamson JA (1990): The power of identity-by-state methods for linkage analysis. Am J Human Genet 46: 254-265.
• Blangero J, Samollow PB, Rocha MB, Hixson JE, Rogers J (1995): The IGF1 locus is a major determinant of serum osteocalcin levels in Mexican Americans. Fourth Annual Meeting of the International Genetic Epidemiology Society, Snowbird, Utah, June 20- 22, 1995.
• Bourgain C, Hoffjan S, Nicolae R, Newman D, Steiner L, Walker K, Reynolds R, Ober C, McPeek MS (2003): Novel case- control test in a founder population identifies P-selectin as an atopy-susceptibility locus. Am J Hum Genet 73: 612-626.
• Brandner FA (1933): A Test of the Significance of the Difference of the Correlation Coefficients in Normal Bivariate Samples. Biometrika 25: 102-109.
• Bunch DS, Gay DM, Welsch RE (1993): Algorithm 717: Subroutines for maximum likelihood and quasi-likelihood estimation of parameters in nonlinear regression models ACM Transactions on Mathematical Software (TOMS) 19: 109-130.
• Canning C, Thompson EA, Skolnick E (1978): Probability functions on complex pedigrees. Adv Appl Prob 10: 26-61.
• Choi Y, Wijsman EM, Weir BS (2009): Case-Control Association Testing in the Presence of Unknown Relationships. Genetic Epidemiology 33: 668-678.
• David F, Johnson NL (1956): Some tests of significance with ordered variables. J R Statist Soc B 18: 1-20.
• Davie AM (1979): The 'singles' method for segregation analysis under incomplete ascertainment. Ann Hum Genet 42: 507-10.
• Davis R, Resnick S (1984): Tail Estimates Motivated by Extreme Value Theory. Ann Stat 12: 1467-87.
• Davis S, Schroeder M, Goldin LR, Weeks DE (1996): Nonparametric simulation-based statistics for detecting linkage in general pedigrees. Am J Hum Genet 58: 867-80.
• Dongarra JJ, Moler CB, Bunch JR, Stewart GW (1979): LINPACK Users' Guide. Philadelphia: SIAM.
• Elston RC, Stewart J (1971): A general model for the genetic analysis of pedigree data. Hum. Hered. 21: 523-542.
• Endelman JB, Jannink J-L (2012): Shrinkage Estimation of the Realized Relationship Matrix. G3 2: 1405-1413.
• Excoffier L (2001): Analysis of population subdivision. In: Balding DJ et al. Handbook of Statistical Genetics. London: Wiley and Sons. 271-307.
• Fain PR (1977): Characteristics of simple sibship variance tests for the detection of major loci and application to height, weight and spatial performance. Am J Hum Genet 42: 109-20.
• Falk CT, Rubinstein P (1987). Haplotype relative risks: an easy reliable way to construct a proper control sample for risk calculations. Ann Hum Genet 51: 227-233.
• Fernandez SA, Fernando R (2002): Technical Note: Determining Peeling Order Using Sparse Matrix Algorithms. J Dairy Sci 85: 1623-1629.
• Filliben J (1975): The probability plot correlation coefficient test for normality. Technometrics 17: 111-117.
• Gauderman WJ (2003). Candidate gene association analysis for a quantitative trait, using parent-offspring trios. Genet Epidemiol 25: 327-338.
• Genz A (1992). Numerical computation of multivariate normal probabilities. J Comp Graph Stat 1: 141-149.
• Guo SW, Thompson EA (1994): Monte Carlo estimation of mixed models for large complex pedigrees. Biometrics 50: 417-432.
• Haberman SJ (1979): Analysis of quantitative data. Volume 2. New developments. New York: Academic Press.
• Haseman JK, Elston RC (1972): The investigation of linkage between a quantitative trait and a marker locus. Behav Genet 2: 3-19.
• Hastings WK [1970] Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57:97-109.
• Hedges LV, Olkin I (1985): Statistical methods for metaanalysis. San Diego: Academic Press.
• Hill HM (1975): A simple general approach to inference about the tail of a distribution. Ann Statist 3: 1163-1174.
• Janss LLG, Van der Werf JHJ, Van Arendonk JAM (1992): Detection of a major gene using segregation analysis in data from several generations, Proceedings of the 43rd Annual Meeting of the European Association of Animal Production, Madrid (Spain), Vol. 1 of Session 5a Free Communications, p. 144.
• Jones GL, Haran M, Caffo BS, Neath R (2005): Fixed-width output analysis for Markov Chain Monte Carlo [Preprint]. http://www.stat.umn.edu/~galin/mcse_rev.pdf
• Kaplan NL, Martin ER, Weir BS (1997): Power studies for the transmission/disequilibrium tests with multiple alleles. Am J Human Genet 60: 691-702.
• Keats BJ, Elston RC (1986): Determination of the order of loci on the short arm of chromosome 11 using two and three locus linkage analyses of pedigree and sib pair data. Genet Epidemiol Suppl 1:147-52.
• Knapp M, Seuchter SA, Baur MP (1993). The haplotype-relative-risk (HRR) method for analysis of association in nuclear families. Am J Hum Genet 52: 1085-1093.
• Knapp M, Wassmer G, Baur MP (1995): The relative efficiency of the Hardy-Weinberg Equilibrium-likelihood and the Conditional on Parental Genotype-likelihood methods for candidate-gene association studies. Am J Hum Genet 57: 1476-1485.
• Knapp M (1999): The transmission/disequilibrium test and parental-genotype reconstruction: the reconstruction-combined transmission/disequilibrium test. Am J Hum Genet 64: 861-870.
• Koval JJ (1997): Algorithm AS 319: Variable Metric Function Minimization. Applied Statistics 46: 515-521.
• Kruglyak L, Daly MJ, Reeve-Daly MP, and Lander ES (1996): Parametric and nonparametric linkage analysis: a unified multipoint approach. Am J Hum Genet 58: 1347-1363.
• Kuonen D (1999): Saddlepoint approximations for distributions of quadratic forms in normal variables. Biometrika 86: 929¿935
• Laird N, Horvath S, and Xu X (2000): Implementing a unified approach to family based tests of association. Genetic Epi 19(Suppl 1): S36-S42.
• Lange K (1986a): The affected sib-pair method using identity by state relations. Am J Hum Genet 39: 148-150.
• Lange K (1986b): A test statistic for the affected-sib-set method. Ann Hum Genet 50: 283-290.
• Lange K (1997): Mathematical and statistical methods for genetic analysis. New York: Springer-Verlag.
• Lange K, Boehnke M (1983): Extensions to pedigree analysis. V. Optimal calculation of Mendelian likelihoods. Hum Hered 33: 291-301.
• Lange K, Goradia T (1987): An algorithm for automatic genotype elimination. Am J Hum Genet 40: 250-256.
• Lange K, Matthysse S (1989): Simulation of pedigree genotypes by random walks. Am J Hum Genet 45: 959-970.
• Li CC, Sacks L (1954): The derivation of joint distribution and correlation between relatives by the use of stochastic matrices. Biometrics 10:347-360.
• Madsen BE, and Browning SR (2009): A groupwise association test for rare mutations using a weighted sum statistic. PLoS Genetics. 5:e1000384.
• McPeek MS, Wu X, Ober C (2004): Best Linear Unbiased Allele-Frequency Estimation in Complex Pedigrees Biometrics 60: 359-367.
• Meyer MC (2001): An alternative unimodal density estimator with a consistent estimate of the mode. Statistica Sinica 11: 1159-1174.
• Moskvina V, Schmidt KM (2008). On multiple-testing correction in genome-wide association studies. Genet Epidemiol 32: 567-573.
• Noreen EW (1989): Computer-intensive methods for testing hypotheses: an introduction. New York: Wiley.
• Oakes D (1982): A concordance test for independence in the presence of censoring. Biometrics 38:451-5.
• Oakes D (2008): On consistency of Kendall's tau under censoring. Biometrika 95: 997­1001.
• Olson J (1995): Robust multipoint linkage analysis. An extension of the Haseman-Elston approach. Genet Epidemiol 12: 177-194.
• Olson J, Rao S, Jacobs K, Elston RC (1998): Linkage of chromosome 1 markers to alcoholism-related phenotypes by sib-pair linkage analysis of principal components. Genetic Analysis Workshop 11. September 8-10, Arcachon, France.
• Penrose LN (1935): The detection of autosomal linkage in data that consists of pairs of brothers or sisters of unspecified parentage. Ann Eugen 8: 133-138.
• Penrose LN (1937): Genetic linkage in graded human characters. Ann Eugen 8: 233-237.
• Plummer M, Best N, Cowles K, Vines K [2006]: Coda: Output analysis and diagnostics for MCMC. R package version 0.10-7.
• Pons O, Chaouche K (1995): Estimation, variance and optimal sampling of gene diversity. II. Diploid locus. Theor Appl Genetics 90: 122-130.
• Powell MJD (2009): The BOBYQA algorithm for bound constrained optimization without derivatives. Department of Applied Mathematics and Theoretical Physics Technical Report DAMTP 2009/NA06
• Ripley BD (1987): Stochastic Simulation. New York: Wiley.
• Resek RW (1974): Alternative tests of skewness: Efficiency comparisons under realistic alternative hypothesis. Proc Bus Econ Statist Section Am Statist Assoc 1974: 546-551.
• Royston P (1993): A pocket-calculator algorithm for the Shapiro-Francia test for non-normality: An application to medicine. Statist Med 12: 181-184.
• Rubin DB (1987): Multiple imputation for nonresponse in surveys. New York: John Wiley.
• Shao J, Tu D (1995): The jackknife and bootstrap. New York: Springer
• Schaid DJ, Sommer SS (1993): Genotype relative risks: methods for design and analysis of candidate-gene association studies. Am J Hum Genet 53: 1114-1126.
• Schelling M (2004): Deterministic calculation and stochastic simulation in multi-point linkage analysis. Doctor of Technical Science Dissertation, Swiss Institute of Technology, Zurich.
• Sham PC, Purcell S (2001): Equivalence between Haseman-Elston and Variance-Components Linkage Analyses for Sib Pairs. Am J Hum Genet 68:1527-1532.
• Sheehan N [1990]. Genetic reconstruction on pedigrees. PhD thesis, University of Washington, Seattle.
• Smith B, Boyle J, Dongarra J, Garbow B, Y Ikebe, V Klema, C Moler [1976]. Matrix Eigensystem Routines, EISPACK Guide, Lecture Notes in Computer Science, Volume 6, Springer Verlag, 1976.
• Smith PG, Pike MC, Hill AP, Breslow NE and Day NE (1981): Algorithm AS 162: Multivariate Conditional Logistic Analysis of Stratum-Matched Case-Control Studies. Appl Statist 30: 190-197.
• Spielman RS, McGinnis RE, Ewens WJ (1993): Transmission test for linkage disequilibrium: the insulin gene region and insulin-dependent diabetes mellitus. Am J Hum Genet 52: 506-516.
• Spielman RS, Ewens WJ (1996): The TDT and other family-based tests for linkage disequilibrium and association [editorial]. Am J Hum Genet 59: 983-989.
• Steele F, Diamond I, Amin S (1996): Immunization uptake in rural Bangladesh: a multilevel analysis. Journal of the Royal Statistical Society, Series A 159: 289-299.
• Stirling WD (1981): Algorithm AS 164: Least Squares Subject to Linear Constraints. Applied Statistics 30: 204-212.
• Tarone RE (1979). Testing the goodness of fit of the binomial distribution. Biometrika 66: 585-590.
• Therneau TM, Grambsch PM, Fleming TR (1990): Martingale-based residuals for survival models. Biometrika 77: 147-160.
• Thompson EA (1991): Probabilities on complex pedigrees: the Gibbs sampler approach. In: Computer science and statistics: 23d Symposium on the interface, pp 371-378.
• Thornton T, McPeek MS (2007): Case-Control Association Testing with Related Individuals A More Powerful Quasi-Likelihood Score Test. American Journal of Human Genetics 81: 321-337.
• Thornton T, McPeek MS (2010): ROADTRIPS: Case-Control Association Testing with Partially or Completely Unknown Population and Pedigree Structure. American Journal of Human Genetics 86: 172-184.
• Visscher PM, Hopper JL (2001): Power of regression and maximum likelihood methods to map QTL from sib-pair and DZ twin data. Ann Human Genet 65: 583-601.
• Wang T, Fernando RL, Stricker C and Elston R C (1996): An approximation to the likelihood for a pedigree with loops, Theor Appl Genet 93: 1299-1309.
• Weeks DE, Lange K (1988): The affected-pedigree-member method of linkage analysis. Am J Hum Genet 42: 315-326.
• Ward PJ (1993): Some developments on the affected-pedigree-member method of linkage analysis. Am J Hum Genet 52: 1200-1215.
• Ward PJ, Bonaiti-Pellie C (1995): Measuring gene-disease association using a general pair method. Genet Epidemiol 12: 681-686.
• Wei GCG, Tanner MA [1990]. A Monte Carlo implementation of the EM algorithm and the Poor Mans's Data Augmentation algorithms. J Am Statist Assoc 85: 699-704.
• Whitehead A, Whitehead J (1991: A general parametric approach to the meta-analysis of randomized clinical trials Statist Med 10: 1665-1677.
• Whittemore AS, Halpern J (1994a): Probability of identity by descent: computation and applications. Biometrics 50: 113-117.
• Whittemore AS, Halpern J (1994b): A class of tests for linkage using affected pedigree members. Biometrics 50: 118-127.
• Wichmann BA, Hill ID (1982): Algorithm AS 183: An Efficient and Portable Pseudo-Random Number Generator. Applied Statistics 31: 188-190.
• Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X (2011): Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test. American Journal of Human Genetics 89: 82-93.
• Yates F (1948): The analysis of contingency tables with groupings based on quantitative characters. Biometrika 38: 176-181.
• Yazdi MH, Visscher PM, Ducrocq V, Thompson R (2002): Heritability, reliability of genetic evaluations and response to selection in proportional hazard models J Dairy Sci 85: 1563-1577.
• Young A (1995): Genetic Analysis System, version 2.0 [Computer program]. Oxford: Oxford University.

LICENCE

I am happy for others to download, use, redistribute, and adapt my own code as they please, with appropriate acknowledgment of the original authorship. Users should note that I make no warranty as to fitness of the software for any purpose ;).

A number of subroutines included in Sib-pair come from other sources (see above), and are identified as such in the comments in the source code. Some of the statistical routines are from the journal Applied Statistics, and were, as I understand it, made available subject to the restriction that no fee is charged for redistribution.

Contributed code is of course welcome.

PROGRAM HISTORY

28-Sep-2020 (1.00b)

09-Aug-2020 (1.00b)

16-Jul-2020 (1.00b)

03-Jul-2020 (1.00b)

30-Jun-2020 (1.00b)

24-Apr-2020 (1.00b)

22-Apr-2020 (1.00b)

10-Apr-2020 (1.00b)

31-Mar-2020 (1.00b)

23-Mar-2020 (1.00b)

19-Mar-2020 (1.00b)

14-Mar-2020 (1.00b)

13-Mar-2020 (1.00b)

10-Mar-2020 (1.00b)

19-Feb-2020 (1.00b)

18-Feb-2020 (1.00b)

11-Nov-2019 (1.00b)

05-Nov-2019 (1.00b)

25-Oct-2019 (1.00b)

26-Sep-2019 (1.00b)

20-Sep-2019 (1.00b)

12-Sep-2019 (1.00b)

02-Sep-2019 (1.00b)

02-Aug-2019 (1.00b)

31-Jul-2019 (1.00b)

24-Jul-2019 (1.00b)

09-Jul-2019 (1.00b)

09-May-2019 (1.00b)

30-Apr-2019 (1.00b)

27-Apr-2019 (1.00b)

26-Apr-2019 (1.00b)

03-Apr-2019 (1.00b)

01-Apr-2019 (1.00b)

23-Mar-2019 (1.00b)

01-Mar-2019 (1.00b)

24-Feb-2019 (1.00b)

22-Feb-2019 (1.00b)

05-Feb-2019 (1.00b)

01-Feb-2019 (1.00b)

12-Dec-2018 (1.00b)

13-Nov-2018 (1.00b)

Now can "show map where <string> -- finally extended search for this and for list to allow multiple search strings. Incorporated routines to read and index BGZF files. Can (finally) use tildes in file paths for all commands (some would not expand). Can declare j2000 as epoch using "set epoch". Bumped up pedigree and individual ID string length to 25. New version can read old binary pedigree files, but not vice versa.

29-Aug-2018 (1.00b)

07-Jul-2018 (1.00b)

20-Jun-2018 (1.00b)

16-Jun-2018 (1.00b)

22-May-2018 (1.00b)

21-May-2018 (1.00b)

08-May-2018 (1.00b)

24-Apr-2018 (1.00b)

19-Apr-2018 (1.00b)

Sib-pair now can read and write string values for categorical traits. The "set cat" command determines whether levels or labels of a categorical trait are printed.

The anova macro makes a nice ANOVA table for the fixed effects in a mixed model analysis (var, mft). No fancy corrections I'm afraid.

08-Apr-2018 (1.00b)

23-Jan-2018 (1.00b)

12-Jan-2018 (1.00b)

21-Nov-2017 (1.00b)

20-Oct-2017 (1.00b)

05-Apr-2017 (1.00b)

03-Feb-2017 (1.00b)

16-Jan-2017 (1.00b)

05-Dec-2016 (1.00b)

30-Oct-2016 (1.00b)

24-Oct-2016 (1.00b)

17-Oct-2016 (1.00b)

13-Oct-2016 (1.00b)

07-Oct-2016 (1.00b)

22-Jul-2016 (1.00b)

13-Jul-2016 (1.00b)

15-Apr-2016 (1.00b)

28-Jan-2016 (1.00b)

21-Dec-2015 (1.00b)

10-Dec-2015 (1.00b)

08-Dec-2015 (1.00b)

06-Dec-2015 (1.00b)

27-Nov-2015 (1.00b)

25-Nov-2015 (1.00b)

18-Nov-2015 (1.00b)

15-Nov-2015 (1.00b)

Found a long-term bug in "residuals" and "predict" command, where allelic coding was not used for first marker - would instead use mean of alleles, so usually no effect.

12-Nov-2015 (1.00b)

01-Nov-2015 (1.00b)

27-Oct-2015 (1.00b)

19-Oct-2015 (1.00b)

08-Oct-2015 (1.00b)

02-Oct-2015 (1.00b)

21-Sep-2015 (1.00b)

Part of my current pipeline for imputation is:

macro chr=1
read bin genos.bin.gz
keep where chromosome %chr
order $mm
# drop indels
drop where spectrum "ID"
# test strand problems using allelic spectrum and HWE
test strand
# merge data where multiple assays same SNP
# allows for strand swaps
test map merge
# drop the superfluous duplicates
drop whe test >= 2
# drop triallelic markers
drop whe all > 2
pack
# match strand to HRC reference data (record in annotation)
flip map /home/davidD/Genetics/Maps/HRC.r1.GRCh37.autosomes.mac5.sites.tab.gz
# filter for VCF output
set loc use aff
if (protyp > 0.9) then use = y
write vcf chr%chr.vcf use
$ bgzip chr%chr.vcf

15-Sep-2015 (1.00b)

11-Sep-2015 (1.00b)

24-Aug-2015 (1.00b)

22-Aug-2015 (1.00b)

21-Aug-2015 (1.00b)

19-Aug-2015 (1.00b)

31-Jul-2015 (1.00b)

27-Jul-2015 (1.00b)

10-Jul-2015 (1.00b)

29-Jun-2015 (1.00b)

11-Jun-2015 (1.00b)

12-May-2015 (1.00b)

17-Apr-2015 (1.00b)

17-Mar-2015 (1.00b)

12-Mar-2015 (1.00b)

11-Mar-2015 (1.00b)

09-Mar-2015 (1.00b)

04-Mar-2015 (1.00b)

26-Feb-2015 (1.00b)

19-Feb-2015 (1.00b)

17-Feb-2015 (1.00b)

04-Feb-2015 (1.00b)

06-Jan-2015 (1.00b)

#  
# ad hoc gene-based etc test window <trait> <m1>...<mN>
# 
macro window
 keep %0 
 ass %1
 neff
 eval (let ((neff (stat-result "stat")) \
            (loci (ls "m"))) \
          (pchisq (* -2 (apply + (map log (map locstat loci)))) (* 2 neff)))

01-Jan-2015 (1.00b)

26-Aug-2014 (1.00b)

23-Jul-2014 (1.00b)

02-Jul-2014 (1.00b)

22-Jun-2014 (1.00b)

27-May-2014 (1.00b)

07-May-2014 (1.00b)

07-Apr-2014 (1.00b)

19-Mar-2014 (1.00b)

07-Feb-2014 (1.00b)

10-Jan-2014 (1.00b)

12-Dec-2013 (1.00b)

21-Aug-2013 (1.00b)

08-Jul-2013 (1.00b)

08-Jun-2013 (1.00b)

05-Jun-2013 (1.00b)

03-Jun-2013 (1.00b)

28-Feb-2013 (1.00b)

04-Jan-2013 (1.00b)

14-Nov-2012 (1.00b)

08-Nov-2012 (1.00b)

25-Oct-2012 (1.00b)

09-Oct-2012 (1.00b)

05-Oct-2012 (1.00b)

04-Oct-2012 (1.00b)

02-Oct-2012 (1.00b)

14-Sep-2012 (1.00b)

12-Sep-2012 (1.00b)

17-Aug-2012 (1.00b)

13-Aug-2012 (1.00b)

04-Aug-2012 (1.00b)

03-Aug-2012 (1.00b)

12-Jul-2012 (1.00b)

05-Jul-2012 (1.00b)

04-Jul-2012 (1.00b)

22-Jun-2012 (1.00b)

13-Jun-2012 (1.00b)

24-May-2012 (1.00b)

18-May-2012 (1.00b)

01-May-2012 (1.00b)

27-Apr-2012 (1.00b)

12-Mar-2012 (1.00b)

02-Mar-2012 (1.00b)

01-Feb-2012 (1.00b)

05-Jan-2012 (1.00b)

19-Dec-2011 (1.00b)

13-Dec-2011 (1.00b)

11-Dec-2011 (1.00b)

01-Dec-2011 (1.00b)

18-Nov-2011 (1.00b)

15-Nov-2011 (1.00b)

08-Nov-2011 (1.00b)

13-Aug-2011 (1.00b)

21-Jul-2011 (1.00b)

13-Jul-2011 (1.00b)

08-Jul-2011 (1.00b)

30-Jun-2011 (1.00b)

16-Jun-2011 (1.00b)

08-Jun-2011 (1.00b)

24-May-2011 (1.00b)

13-May-2011 (1.00b)

18-Mar-2011 (1.00b)

11-Mar-2011 (1.00b)

24-Feb-2011 (1.00b)

09-Feb-2011 (1.00b)

03-Feb-2011 (1.00b)

13-Jan-2011 (1.00b)

17-Dec-2010 (1.00b)

07-Dec-2010 (1.00b)

05-Dec-2010 (1.00b)

29-Nov-2010 (1.00b)

18-Nov-2010 (1.00b)

03-Sep-2010 (1.00b)

29-Aug-2010 (1.00b)

27-Aug-2010 (1.00b)

16-Aug-2010 (1.00b)

16-Aug-2010 (1.00b)

15-Aug-2010 (1.00b)

13-Aug-2010 (1.00b)

10-Aug-2010 (1.00b)

30-Jul-2010 (1.00b)

28-Jul-2010 (1.00b)

26-Jul-2010 (1.00b)

14-Jul-2010 (1.00b)

09-Jul-2010 (1.00b)

28-Jun-2010 (1.00b)

24-Jun-2010 (1.00b)

15-Jun-2010 (1.00b)

25-May-2010 (1.00b)

21-May-2010 (1.00b)

20-May-2010 (1.00b)

19-May-2010 (1.00b)

14-May-2010 (1.00b)

11-May-2010 (1.00b)

07-May-2010 (1.00b)

eval (define (s-repeat n cmdstring) \
(let loop ((i 1)) \
(if (<= i n) \
(begin \
(run cmdstring) (loop (+ i 1))) \
(format "# end of ~d iterations~\%" n))))
macro repeat
eval (s-repeat %1 "%+2")
;;;;

05-May-2010 (1.00b)

29-Mar-2010 (1.00b)

22-Mar-2010 (1.00b)

17-Mar-2010 (1.00b)

12-Mar-2010 (1.00b)

09-Mar-2010 (1.00b)

08-Mar-2010 (1.00b)

26-Feb-2010 (1.00b)

25-Feb-2010 (1.00b)

24-Feb-2010 (1.00b)

23-Feb-2010 (1.00b)

12-Feb-2010 (1.00b)

06-Jan-2010 (1.00b)

5-Jan-2010 (1.00b)

24-Dec-2009 (1.00b)

23-Dec-2009 (1.00b)

21-Dec-2009 (1.00b)

11-Dec-2009 (1.00b)

03-Dec-2009 (1.00b)

30-Nov-2009 (1.00b)

25-Nov-2009 (1.00b)

19-Nov-2009 (1.00b)

10-Nov-2009 (1.00b)

05-Nov-2009 (1.00b)

03-Nov-2009 (1.00b)

23-Oct-2009 (1.00b)

21-Oct-2009 (1.00b)

19-Oct-2009 (1.00b)

29-Sep-2009 (1.00b)

23-Sep-2009 (1.00b)

22-Sep-2009 (1.00b)

18-Sep-2009 (1.00b)

08-Sep-2009 (1.00b)

03-Sep-2009 (1.00b)

02-Sep-2009 (1.00b)

01-Sep-2009 (1.00b)

28-Aug-2009 (1.00b)

27-Aug-2009 (1.00b)

20-Aug-2009 (1.00b)

18-Aug-2009 (1.00b)

23-Jul-2009 (1.00b)

20-Jul-2009 (1.00b)

10-Jul-2009 (1.00b)

03-Jul-2009 (1.00b)

02-Jul-2009 (1.00b)

26-Jun-2009 (1.00b)

12-Jun-2009 (1.00b)

10-Jun-2009 (1.00b)

09-Jun-2009 (1.00b)

05-Jun-2009 (1.00b)

31-May-2009 (1.00b)

25-May-2009 (1.00b)

22-May-2009 (1.00b)

19-May-2009 (1.00b)

28-Apr-2009 (1.00b)

23-Apr-2009 (1.00b)

06-Apr-2009 (1.00b)

03-Apr-2009 (1.00b)

31-Mar-2009 (1.00b)

24-Mar-2009 (1.00b)

23-Mar-2009 (1.00b)

18-Mar-2009 (1.00b)

27-Feb-2009 (1.00b)

23-Feb-2009 (1.00b)

19-Feb-2009 (1.00b)

16-Feb-2009 (1.00b)

09-Feb-2009 (1.00b)

06-Feb-2009 (1.00b)

05-Feb-2009 (1.00b)

3-Feb-2009 (1.00b)

18-Jan-2009 (1.00b)

12-Jan-2009 (1.00b)

07-Jan-2009 (1.00b)

05-Jan-2009 (1.00b)

1-Jan-2009 (1.00b)

23-Dec-2008 (1.00b)

19-Dec-2008 (1.00b)

17-Dec-2008 (1.00b)

15-Dec-2008 (1.00b)

11-Dec-2008 (1.00b)

09-Dec-2008 (1.00b)

05-Dec-2008 (1.00b)

02-Dec-2008 (1.00b)

23-Nov-2008 (1.00b)

21-Nov-2008 (1.00b)

20-Nov-2008 (1.00b)

06-Nov-2008 (1.00b)

04-Nov-2008 (1.00b)

03-Nov-2008 (1.00b)

31-Oct-2008 (1.00b)

30-Oct-2008 (1.00b)

24-Oct-2008 (1.00b)

23-Oct-2008 (1.00b)

22-Oct-2008 (1.00b)

20-Oct-2008 (1.00b)

17-Oct-2008 (1.00b)

14-Oct-2008 (1.00b)

26-Sep-2008 (1.00b)

24-Sep-2008 (1.00b)

23-Sep-2008 (1.00b)

22-Sep-2008 (1.00b)

19-Sep-2008 (1.00b)

12-Sep-2008 (1.00b)

11-Sep-2008 (1.00b)

04-Sep-2008 (1.00b)

01-Sep-2008 (1.00b)

29-Aug-2008 (1.00b)

28-Aug-2008 (1.00b)

21-Aug-2008 (1.00b)

20-Aug-2008 (1.00b)

15-Aug-2008 (1.00b)

12-Aug-2008 (1.00b)

11-Aug-2008 (1.00b)

08-Aug-2008 (1.00b)

07-Aug-2008 (1.00b)

06-Aug-2008 (1.00b)

05-Aug-2008 (1.00b)

29-Jul-2008 (1.00b)

24-Jul-2008 (1.00b)

23-Jul-2008 (1.00b)

18-Jul-2008 (1.00b)

15-Jul-2008 (1.00b)

09-Jul-2008 (1.00b)

07-Jul-2008 (1.00b)

02-Jul-2008 (1.00b)

30-Jun-2008 (1.00b)

24-Jun-2008 (1.00b)

20-Jun-2008 (1.00b)

19-Jun-2008 (1.00b)

16-Jun-2008 (1.00b)

13-Jun-2008 (1.00b)

11-Jun-2008 (1.00b)

06-Jun-2008 (1.00b)

05-Jun-2008 (1.00b)

04-Jun-2008 (1.00b)

02-Jun-2008 (1.00b)

31-May-2008 (1.00b)

30-May-2008 (1.00b)

28-May-2008 (1.00b)

27-May-2008 (1.00b)

22-May-2008 (1.00b)

16-May-2008 (1.00b)

15-May-2008 (1.00b)

14-May-2008 (1.00b)

08-May-2008 (1.00b)

01-May-2008 (1.00b)

29-Apr-2008 (1.00b)

24-Apr-2008 (1.00b)

22-Apr-2008 (1.00b)

21-Apr-2008 (1.00b)

15-Apr-2008 (1.00b)

09-Apr-2008 (1.00b)

08-Apr-2008 (1.00b)

07-Apr-2008 (1.00b)

04-Apr-2008 (1.00b)

02-Apr-2008 (1.00b)

01-Apr-2008 (1.00b)

31-Mar-2008 (1.00b)

28-Mar-2008 (1.00b)

25-Mar-2008 (1.00b)

18-Mar-2008 No. 2 (1.00b)

18-Mar-2008 (1.00b)

17-Mar-2008 (1.00b)

16-Mar-2008 (1.00b)

12-Mar-2008 (1.00b)

05-Mar-2008 (1.00b)

03-Mar-2008 (1.00b)

28-Feb-2008 (1.00b)

21-Feb-2008 (1.00b)

20-Feb-2008 (1.00b)

19-Feb-2008 (1.00b)

15-Feb-2008 (1.00b)

14-Feb-2008 (1.00b)

12-Feb-2008 (1.00b)

06-Feb-2008 (1.00b)

05-Feb-2008 (1.00b)

31-Jan-2008 (1.00b)

31-Jan-2008 (1.00b)

30-Jan-2008 (1.00b)

29-Jan-2008 (1.00b)

25-Jan-2008 (1.00b)

17-Jan-2008 (1.00b)

14-Jan-2008 (1.00b)

07-Jan-2008 (1.00b)

02-Jan-2008 (1.00b)

10-Dec-2007 (1.00b)

24-Nov-2007 (1.00b)

20-Nov-2007 (1.00b)

19-Nov-2007 (1.00b)

16-Nov-2007 (1.00b)

09-Nov-2007 (1.00b)

08-Nov-2007 (1.00b)

05-Nov-2007 (1.00b)

02-Nov-2007 (1.00b)

25-Oct-2007 (1.00b)

19-Sep-2007 (1.00b)

17-Sep-2007 (1.00b)

24-Aug-2007 (1.00b)

20-Aug-2007 (1.00b)

17-Aug-2007 (1.00b)

16-Aug-2007 (1.00b)

10-Aug-2007 (1.00b)

06-Aug-2007 (1.00b)

03-Aug-2007 (1.00b)

02-Aug-2007 (1.00b)

31-Jul-2007 (1.00b)

30-Jul-2007 (1.00b)

25-Jul-2007 (1.00b)

12-Jul-2007 (1.00b)

25-Jun-2007 (1.00b)

20-Jun-2007 (1.00b)

19-Jun-2007 (1.00b)

15-Jun-2007 (1.00b)

14-Jun-2007 (1.00b)

13-Jun-2007 (1.00b)

01-Jun-2007 (1.00b)

21-May-2007 (1.00b)

18-May-2007 (1.00b)

17-May-2007 (1.00b)

11-May-2007 (1.00b)

08-May-2007 (1.00b)

02-May-2007 (1.00b)

27-Apr-2007 (1.00b)

21-Apr-2007 (1.00b)

20-Apr-2007 (1.00b)

19-Apr-2007 (1.00b)

18-Apr-2007 (1.00b)

17-Apr-2007 (1.00b)

14-Apr-2007 (1.00b)

13-Apr-2007 (1.00b)

05-Apr-2007 (1.00b)

03-Apr-2007 (1.00b)

28-Mar-2007 (1.00b)

26-Mar-2007 (1.00b)

22-Mar-2007 (1.00b)

15-Mar-2007 (1.00b)

09-Mar-2007 (1.00b)

06-Mar-2007 (1.00b)

01-Mar-2007 (1.00a)

23-Feb-2007 (1.00a)

20-Feb-2007 (1.00a)

16-Feb-2007 (1.00a)

13-Feb-2007 (1.00a)

12-Feb-2007 (1.00a)

02-Feb-2007 (1.00a)

31-Jan-2007 (1.00a)

The Windows (and even linux, if desired) version now can call a graphical file browser, if the file name is not specified on the command line (eg "inc"). This is the Java AWT file dialogue called via the crossplatform/language JAPI library. This is slower on older computers due to a delay in starting up the Java runtime. On linux, the java processes ("japi kernel") seem to persist sleeping. Obviously, the Java runtime must be present on your computer, along with the older AWT library.

12-Jan-2007 (1.00a)

08-Jan-2007 (1.00a)

02-Jan-2007 (1.00a)

21-Dec-2006 (1.00a)

11-Dec-2006 (1.00a)

08-Dec-2006 (1.00a)

04-Dec-2006 (1.00a)

29-Nov-2006 (1.00a)

28-Nov-2006 (1.00a)

24-Nov-2006 (nsp95)

23-Nov-2006 (nsp95)

21-Nov-2006 (nsp95)

16-Nov-2006 (nsp95)

15-Nov-2006 (nsp95)

06-Nov-2006 (nsp95)

19-Oct-2006 (nsp95)

18-Oct-2006 (nsp95)

12-Oct-2006 (nsp95)

10-Oct-2006 (nsp95)

03-Oct-2006 (nsp95)

29-Sep-2006 (nsp95)

31-Aug-2006 (nsp95)

• Testing of sex assignment using sex-linked markers including Amelogenin.
• Summary of Mendelian errors by pedigree and by locus (similar to that generated by the taberrors shell script).
• Nicer summary of monozygotic twin assignment test results.
• Penalised nonparametric ML estimation of MCMC posterior modes.
• Another quantitative trait TDT, following the simplified regression model formulation of Gauderman et al 2003.
• Recoding nucleotide letter codes to/from numerical code.
• Information about dataset memory utilisation.

07-Jul-2006 (1.00a17)

06-Jul-2006 (1.00a17)

27-Jun-2006 (1.00a17)

22-Jun-2006 (1.00a17)

21-Jun-2006 (1.00a17)

20-Jun-2006 (1.00a17)

15-Jun-2006 (1.00a17)

13-Jun-2006 (1.00a16)

07-Jun-2006 (1.00a16)

01-Jun-2006 (1.00a16)

31-May-2006 (1.00a16)

26-May-2006 (1.00a16)

25-May-2006 (1.00a16)

19-May-2006 (1.00a15)

17-May-2006 (1.00a15)

15-May-2006 (1.00a14)

08-May-2006 (1.00a13)

11-Apr-2006 (1.00a12)

17-Mar-2006 (1.00a11)

16-Mar-2006 (1.00a10)

15-Mar-2006 (1.00a10)

13-Mar-2006 (1.00a9)

28-Feb-2006 (1.00a8)

27-Feb-2006 (1.00a7)

20-Feb-2006 (1.00a7)

13-Feb-2006 (1.00a6)

6-Feb-2006 (1.00a5)

3-Feb-2006 (1.00a4)

30-Jan-2006 (1.00a3)

30-Jan-2006 (1.00a2)

20-Jan-2006 (1.00a1)

24-Dec-2005 (0.99.9)

23-Dec-2005 (0.99.9)

22-Dec-2005 (0.99.9)

16-Dec-2005 (0.99.9)

09-Dec-2005 (0.99.9)

14-Nov-2005 (0.99.9)

13-Oct-2005 (0.99.9)

5-Oct-2005 (0.99.9)

29-Sep-2005 (0.99.9)

21-Sep-2005 (0.99.9)

12-Sep-2005 (0.99.9)

31-Aug-2005 (0.99.9)

27-Jul-2005 (0.99.9)

18-Jul-2005 (0.99.9)

28-Feb-2005 (0.99.9)

25-Feb-2005 (0.99.9)

06-Jan-2005 (0.99.9)

20-Dec-2004 (0.99.9)

08-Sep-2004 (0.99.9)

21-Jun-2004 (0.99.9)

3-Jun-2004 (0.99.9)

24-May-2004 (0.99.9)

20-May-2004 (0.99.9)

19-May-2004 (0.99.9)

7-May-2004 (0.99.9)

23-Apr-2004 (0.99.9)

18-Feb-2004 (0.99.9)

12-Feb-2004 (0.99.9)

21-Jan-2004 (0.99.9)

8-Jan-2004 (0.99.9)

18-Dec-2003 (0.99.9)

24-Nov-2003 (0.99.9)

11-Nov-2003 (0.99.9)

24-Oct-2003 (0.99.9)

13-Oct-2003 (0.99.9)

9-Sep-2003 (0.99.9)

21-Aug-2003 (0.99.9)

18-Aug-2003 (0.99.9)

07-Aug-2003 (0.99.9)

30-Jul-2003 (0.99.9)

29-Jul-2003 (0.99.9)

18-Jul-2003 (0.99.9)

17-Jul-2003 (0.99.9)

11-Jul-2003 (0.99.9)

26-Jun-2003 (0.99.9)

26-May-2003 (0.99.9)

24-April-2003 (0.99.9)

17-April-2003 (0.99.9)

10-April-2003 (0.99.9)

09-April-2003 (0.99.9)

28-March-2003 (0.99.9)

20-March-2003 (0.99.9)

28-February-2003 (0.99.9)

12-February-2003 (0.99.9)

31-January-2003 (0.99.9)

23-January-2003 (0.99.9)

22-January-2003 (0.99.9)

17-January-2003 (0.99.9)

04-October-2002 (0.99.9)

18-October-2002 (0.99.9)

17-October-2002 (0.99.9)

4-October-2002 (0.99.9)

25-September-2002 (0.99.9)

18-September-2002 (0.99.9)

29-July-2002 (0.99.9)

01-July-2002 (0.99.9)

26-June-2002 (0.99.9)

20-June-2002 (0.99.9)

19-June-2002 (0.99.9)

12-June-2002 (0.99.9)

28-May-2002 (0.99.9)

17-May-2002 (0.99.9)

17-Apr-2002 (0.99.9)

05-Apr-2002 (0.99.9)

13-Mar-2002 (0.99.9)

12-Mar-2002 (0.99.9)

06-Mar-2002 (0.99.9)

19-Feb-2002 (0.99.9)

18-Feb-2002 (0.99.9)

15-Feb-2002 (0.99.9)

31-Jan-2002 (0.99.9)

23-Jan-2002 (0.99.9)

18-Jan-2002 (0.99.9)

11-Jan-2002 (0.99.9)

3-Jan-2002 (0.99.9)

29-Nov-2001 (0.99.9)

27-Nov-2001 (0.99.9)

26-Nov-2001 (0.99.9)

16-Nov-2001 (0.99.9)

15-Nov-2001 (0.99.9)

Number of pedigree errors  =     1
Number of deleted records  =     4

02-Nov-2001 (0.99.9)

11-Sep-2001 (0.99.9)

04-Sep-2001 (0.99.9)

21-Aug-2001 (0.99.9)

17-Aug-2001 (0.99.9)

14-Aug-2001 (0.99.9)

18-Jul-2001 (0.99.9)

17-Jul-2001 (0.99.9)

The legal equivalent is:
if (b eq 1 and c eq 1) then d=1 else if (b eq 1 and c ne 1) then d=2 else d=3
or
if b ne 1 then d=3 else if c eq 1 then d=1 else d=2.

14-Jul-2001 (0.99.9)

10-Jul-2001 (0.99.9)

6-Jul-2001 (0.99.9)

5-Jul-2001 (0.99.9)

28-Jun-2001 (0.99.9)

20-Jun-2001 (0.99.9)

20-Jun-2001 (0.99.9)

19-Jun-2001 (0.99.9)

18-Jun-2001 (0.99.9)

14-Jun-2001 (0.99.9)

8-Jun-2001 (0.99.9)

6-Jun-2001 (0.99.9)

24-May-2001 (0.99.9)

26-Apr-2001 (0.99.9)

12-Mar-2001 (0.99.9)

03-Oct-2000 (0.99.9)

26-Sep-2000 (0.99.9)

20-Sep-2000 (0.99.9)

26-Jul-2000 (0.99.9)

14-Jul-2000 (0.99.9)

23-Jun-2000 (0.99.9)

15-Jun-2000 (0.99.9)

13-Jun-2000 (0.99.9)

08-Jun-2000 (0.99.9)

28-Feb-2000 (0.99.9)

25-Feb-2000 (0.99.9)

15-Feb-2000 (0.99.9)

14-Feb-2000 (0.99.9)

09-Feb-2000 (0.99.9)

20-Jan-2000 (0.99.9)

08-Dec-1999 (0.99.8)

26-Nov-1999 (0.99.7)

12-Nov-1999 (0.99.7)

29-Sep-1999 (0.99.7)

16-Sep-1999 (0.99.7)

03-Sep-1999 (0.99.7)

26-Aug-1999 (0.99.6)

23-Aug-1999 (0.99.5)

12-Aug-1999 (0.99.4)

05-Jul-1999 (0.99.3)

27-Jan-1999 (0.98.9)

25-Jan-1999 (0.98.8)

21-Jan-1999 (0.98.7)

28-Aug-1998 (0.98.3)

14-Aug-1998 (0.98.2)

21-Jul-1998 (0.97.9)

2-Jul-1998 (0.97.7)

26-Jun-1998 (0.97.6)

16-Jun-1998 (0.97.5)

25-May-1998 (0.97.3)

01-Apr-1998 (0.97.2)

13-Mar-1998 (0.97.0)

03-Mar-1998

02-Mar-1998 (0.96.5)

The Monte-Carlo Markov Chain algorithm for simulating missing genotypes has been further tinkered with, following problems encountered in multigeneration pedigrees where multiple consecutive generations are untyped. I have reintroduced a "switch of origins" operation for heterozygote children of untyped x untyped matings, as the existing "mutation" operation will shuffle such parental genotypes only when the child has no offspring (a fact that had eluded me until now). This seems to fix the trouble. Reanalysis of the example datasets included with Sib-pair found little difference in ibd based statistics compared to the old algorithms. Using the Lange-Goradia algorithm for producing starting genotypes for MCMC can still fail, as described in the original 1987 paper, when loops are present. If this occurs, Sib-pair now switches to its alternative gene-dropping based algorithm.

11-Feb-1998 (0.96.0)

A "set burn-in" option has been added to allow specification of the number of iterations of the MCMC algorithms to be run prior to the iterations used to calculate statistics. Previously, this was set to zero. In the empirical tests I had done, the Lange-Goradia algorithm generated starting genotypes seemed to give unbiased results, and so I had removed the burn-in. More recently, I have found example families (with missing genotypes) where estimated ibd is biased upwards. This bias is reduced by a suitable number of burn-in iterations.

28-Jan-1998 (0.95.6)

23-Jan-1998 (0.95.5)

13-Jan-1998

12-Jan-1998 (0.95.4)

11-Dec-1997

10-Dec-1997 (0.95.3)

02-Dec-1997

26-Nov-1997

25-Nov-1997 (0.95.2)

19-Nov-1997 (0.95.1)

Elapsed time now measured by time(), rather than secnds() for DOS version. A refinement only of interest for overnight runs.

Program now warns if there are fewer fields for a pedigree record than expected -- previously these were silently padded out as missing. It still silently ignores extra fields.

Memory utilisation decreased by equivalencing several work arrays.

13-Nov-1997 (0.95.0)

06-Nov-1997

05-Nov-1997

24-Oct-1997

25-Sep-1997

21-Aug-1997

20-Aug-1997

11-Aug-1997

06-Aug-1997

03-Aug-1997

16-Jul-1997

15-Jul-1997

28-Jun-1997

26-May-1997

06-Mar-1997

03-Mar-1997 (0.94)

05-Feb-1997

Feb-1997

18-Dec-1996 (0.93)

18-Oct-1996 (0.92)

Sep-1996

Jul-1996

Mar-1996

Aug-1995

May-1995

Apr-1995

Appendix: Embedded Scheme Commands

The Scheme implemented within Sib-pair supports only a subset of the language. The only atomic data types are booleans, (long) integer and double precision real numbers, and strings. Characters (and character procedures) are handled as strings of length one.

<	Numerically less than
<=	Numerically less than or equal
=	Numerically equal
>	Numerically greater than
>=	Numerically greater than or equal to
-	Integer subtraction
/	Integer division
*	Integer multiplication
+	Integer addition
abs	Absolute value
and	Logical and
append	Append to end of list
apply	Apply function to list
apropos	List matching functions
assoc	Find matching pair in an "alist" (list of paired values) using equal?
assq	Find matching pair in an "alist" using eq?
assv	Find matching pair in an "alist" using eqv?
atom?	Is a simple number or string?
begin	Start a block
boolean?	Is a boolean?
car	First element of pair or list
cdr	Remainder of pair or list
caar	First element of first element of pair or list
cadr	Und so weiter
cdar
cddr
caaar
cdddr
caadr
cadar
caddr
cdaar
cdadr
cddar
call/cc
call-with-current-continuation
case	Conditional branching
char<?	Is char lexicographically less?
char>?	Is char lexicographically greater?
char<=?	Is char lexicographically less than or equal?
char>=?	Is char lexicographically greater than or equal?
char-alphabetic?	Is char alphabetic?
char-ci<?	Is case-insensitive char lexicographically less?
char-ci>?	Is case-insensitive char lexicographically greater?
char-ci<=?	Is case-insensitive char lexicographically less than or equal?
char-ci>=?	Is case-insensitive char lexicographically greater than or equal?
char-numeric?	Is char a digit?
char-whitespace?	Is char a tab, space, CR, NL, HR?
char->integer	Return ASCII code
char-downcase	Return lower case
char-upcase	Return upper case
closure?	Is a closure
cond	Conditional branching
cons	Create a Lisp pair of elements
cons-stream
current-second	Return seconds since epoch
define	Define a variable or function
delay
display	Print a variable or function
do	Looping construct
else
eq?
equal?
eqv?
error	Print an error message
eval	Evaluate a SEXPR (Scheme expression)
even?	Test if even
exact->inexact	Cast to real
exit	Exit Sib-pair completely
expt	Integer exponentiation
force
for-each	Like map, but used for side-effects only
format	Print formatted data
gc	Force a garbage collection
gcd	Greatest common divisor
gensym	Create a new unique symbol
get-closure-code
help	Minimal information about Scheme
if	Conditional branching
inexact->exact	Cast to integer
integer?	Is an integer
integer->char	Return character for given ASCII code
lambda	Function
lcm	Least common multiple
length	Gives length of a list
let	Declare a local variable
let*
letrec
list	Create a list
list?	Test if a list?
list-ref	Gives Nth element of a list
list-tail	Drops the first k elements of a list
ls	List Sib-pair variables
macro
make-list	Create a list of given length
make-string	Create a string of given length
map	Apply a function to each element of a list in turn
max	Maximum of arguments
member	Test if present in a list
memq	Test if present in a list using equal?
memv	Test if present in a list using eqv?
min	Minimum of arguments
modulo	Give modulo
negative?	Test if negative integer
newline	Print a newline
new-segment
not	Logical negation
null?	Test if a null
number?	Test if an integer number
number->string	Convert from number to string
odd?	Test if odd
or	Logical or
pair?	Test if a Lisp pair
peek-char	Read one character without advancing
positive?	Test if a positive integer
print-width
procedure?	Is a function?
quasiquote	Allow some of quoted expression to evaluate
quit	Leave the Scheme interpreter
quote	Quoted expression left unevaluated
quotient	Integer division
random	A random integer from U(1, N)
read
read-char	Read one character from standard input or file
read-line	Read one line of input from standard input or file
remainder	Remainder
reverse	Reverse ordering of a list
run	Run a Sib-pair command
set!	Set the contents of an existing variable
set-car!	Set contents of first element of pair or list
set-cdr!	Set contents of rest of pair or list
sqrt	Integer square root
string?	Is argument a string?
string=?	Are strings equal?
substring?	Is string 1 a substring of string 2?
string<?	Is string lexicographically less?
string>?	Is string lexicographically greater?
string<=?	Is string lexicographically less than or equal?
string>=?	Is string lexicographically greater than or equal?
string-append	Concatenate strings
string-length	Length of argument string
string-ref	Get kth character from string
string-set!	Set the value of a substring
string-split	Split a string into a list of words
string->number	Convert from a string to a number
string->symbol	Convert from a string to a symbol
substring	Get the value of a substring
symbol?	Is a symbol?
symbol->string	Convert from a symbol to a string
system	Passes commands to the operating system
unquote
unquote-splicing
write
zero?	Test if equal to zero

There are a number of Sib-pair specific builtin commands for example allowing access to locus data:

apropos <str>	lists commands containing that string.
help	Information about this Scheme implementation.
isatty? <fil>	test if interactive session.
file-exists? <fil>	test file.
file-list <nam>	list contents of a directory.
file-delete <fil>	delete a file.
open-input-file <fil>	open a port.
close-input-port <port>	close a port.
read-line [<port>]	reads in next line from stdin or open file.
format <port> "{~[:num:][@][ASD~%]}" <args>...	formatted output.
string-split <str> [<sep>]	splits string on white space or optional char.
substring? <sub> <str>	returns start of substring in string.
regexp <sstring> <str>	test if search regular expession matches string.
system <cmd>	passes command to shell.
getenv <nam>	returns value of environment variable.
date	returns current date and time.
time	returns current time as number of seconds since epoch.
seq <sta> <fin> [<step>]	generate sequence.
filter <test> <list>	filter contents of list
environment-bound? <name>	test if symbol bound to a variable

version	prints Sib-pair version.
run <cmd> ...	runs a Sib-pair command.
pass-command <cmd> ...	stores Sib-pair commands to the buffer for evaluation once you return to the usual Sib-pair prompt.

Sib-pair locus dataset accessors:

ls [<typ>]	creates a list of locus names (of given type "adhmqx".
nloci [<typ>]<	returns total number of loci.
loc <index>	returns locus at that position in the locus list.
loc-set! <idx> <name>	set name of locus at that position in the locus list.
lochash-update!	update hash of locus names
locord <loc>	returns position of a locus in the locus list.
locnotes [<loc>..<loc>]	returns notes for a locus.
locnotes-set! <loc> <str>	rewrites notes for a locus.
loctyp [<loc>..<loc>]	evaluates type of a locus ("adhmqx").
loctyp-set! <loc> <typ><	set type of a locus ("adhmqx").
map-position [<loc>..<loc>]	returns map position for locus.
map-position-set! <loc> <pos>	sets map position for locus.
chromosome [<loc>..<loc>]	returns locus chromosome.
chromosome-set! <loc> <chr>	sets locus chromosome.
locstat [<loc>..<loc>]	returns last P-value for a locus.
locstat-init! <title>	initializes all locus P-values.
locstat-set! <loc> <value>	sets P-value for a locus.
locrank <loc>	returns rank of locus test statistic.
stat-result ["pval|lik|npars|lrt|df|stat|var"]	returns result of last model.

Sib-pair phenotype dataset accessors:

nobs	number of individual records.
npeds	number of pedigrees.
nactpeds	number of active pedigrees
active-status	activity status of pedigrees
set-active-status! <idx> <lev>	set activity status of pedigree
active-pedigrees [<idx>...]	list of active pedigree names
pedigrees [<idx>...]	list of pedigree names.
pedigree-size [<idx>...]	size of ith pedigree.
pedigree-members [<idx>...]	list of indices of pedigree members
individual-pedigree [<idx>...])	give pedigree ID for index
individual-name [<idx>...]	give ID for index.
set-pedigree-name! "<ped>" <newname>	Set pedigree ID string .
set-individual-name! <idx>|"<id>" |("<ped>" "<id>") <newname>	Set ID string for individual.
individual-index [<id>|<idx>...]	give index for ID.
insert-record! <idx>|(<ped> <id>) ['after]	insert data row before index.
father [<idx>...]	father indices.
mother [<idx>...]	mother indices.
imztwin [<idx>...]	MZ twin pointer.
sex [<idx>...]	sex value for individual(s).
set-sex! [<idx>...] <sex>	set sex value for individual.
data <loc> [<idx>]	phenotype for individual(s).
set-data! <idx>|"<id>" |("<ped>" "<id>") <loc> <val>	Set phenotype for individual.
data-counts <loc> [<loc>]	tabulation of phenotype levels.
allele-freqs <loc>	tabulation of alleles.

There are a number of builtin statistical functions (20100323 onwards), calling the appropriate Fortran routines.

pnorm <Z>	Gaussian upper tail probability
qnorm <p>	Gaussian quantile for given upper tail probability
pchisq <X2> <df> [<ncp>]	Central and non central chi-square upper tail probability
qchisq <p> <df>	Chi-square quantile for given upper tail probability
fp <F> <df1> <df2>	F upper tail probability
pgamma <X> <df>	Incomplete gamma integral
lgamma <X>	Log gamma
dbeta <X> <shape1> <shape2>	Density of Beta distribution.
bivnor <Z1> <Z2> <r>	Bivariate Gaussian upper tail probability
pmvnorm <p> <dir> <cor> ['genz]	Multivariate Gaussian tail probability
pchisqsum <X> <alpha>	Upper tail probability for distributions of quadratic forms in normal variables

And a few statistical data manipulation routines, that follow the XLispStat equivalents:

sort <list>	Sort a list of numbers
order <list>	Indices of order of numbers in list
rank <list>	Rank of numbers in list
quantile <list> <p>	quantile(s) of list
sample-seq <N> <size> ['replace]	sample a sequence 1..N
sample <list> <size> ['replace]	sample a list
differences <list>	Lag-1 differences in list
which <list>	Indices of non-null elements of list
unique <list>	Indices of unique elements of list
duplicated <list>	Indices of repeat elements of list
intersect <list> <list>	Common elements of lists
setdiff <list> <list>	Elements unique to first list
union <list> <list>	Elements of both lists
list-select <list> <idx-list>	Sublist using indices
stats <list>	N, missing, mean, variance, min, max of list.

Sib-pair Scheme also now (2010-05-9) contains builtin functions to call the JAPI (java application programming interface) library (http://www.japi.de), which allows building GUIs etc. JAPI allows Fortran code to interface the Java AWT (Abstract Windowing Toolkit). The JAPI website details the commands.

j_start	Start JAPI listener
j_quit	Stop JAPI listener
j_frame	Create a frame
j_panel	Create a panel
j_borderpanel	Create a borderpanel
j_dialog	Create a dialogue
j_button	Create a button
j_radiobutton	Create a radiobutton
j_radiogroup	Create a radiogroup
j_checkbox	Create a checkbox
j_list	Create a list
j_add	Add an object to a container
j_setcolor	Set the foreground colour
j_setcolorbg	Set the background colour
j_setnamedcolorbg	Set the background colour
j_getselect	Return index of selected item
j_select	Select an item
j_deselect	Deselect an item
j_fileselect	A file dialog
j_filedialog	A file dialog
j_enable	Activate an object
j_disable	Deactivate an object
j_additem	Add an item to a list
j_seperator	Draw a separator
j_textfield	Create a textfield object
j_textarea	Create a textarea
j_setborderpos	Place an object using borderlayout
j_setrows	Set number of rows eg textarea
j_setcolumns	Set number of columns eg textarea
j_getrows	Get number of rows eg textarea
j_getcolumns	Get number of columns eg textarea
j_getlength	Get length in pixels of object
j_getselstart	Start of selected text
j_getselend	End of selected text
j_selecttext	Select text
j_gettext	Get text from label or list item
j_getseltext	Get selected text
j_getitem	Get a list item
j_label	Create a label
j_getcurpos	Get current cursor position
j_setcurpos	Set current cursor position
j_setfont	Set font
j_settext	Set text in object
j_inserttext	Insert text at position in eg textarea
j_replacetext	Replace text at position in eg textarea
j_delete	Delete text
j_dispose	Free resource
j_menubar	Create a menubar
j_menu	Create a menu
j_menuitem	Create an item for a menu
j_pack	Pack layout of frame or panel using layout manager
j_show	Make object visible
j_hide	Hide object
j_keylistener	Listen for key stroke if object active
j_getkeycode	Returns pressed key code
j_getkeychar	Returns pressed key ASCII code
j_mouselistener	Listen for mouse activation
j_getmousebutton	Return pressed mouse button
j_nextaction	Return next action of any object
j_getwidth	Get width of object (pixels)
j_getheight	Get height of object (pixels)
j_getpos	Get position of object (X.Y)
j_setpos	Set position of object (X,Y)
j_setsize	Get size of object (pixels)
j_setalign	Get layout alignment for object
j_setborderlayout	Layout manager
j_setgridlayout	Layout manager
j_setflowlayout	Layout manager

Sib-pair Scheme (2010-11-29) also contains builtin functions to call the EGGX/ProCALL graphical library ( http://www.ir.isas.jaxa.jp/~cyamauch/eggx_procall/), which allows plotting and building simple GUIs etc under X Windows.

ggetdisplayinfo	Get X display info. Returns a list containing ndepth (8,16,24 bit), nrwidth (screen width), nrheight (screen height).
gopen	Open a window. Takes two optional arguments nxwidth (width of window in pixels) and nywidth (height of window in pixels). Returns the handle (integer window index) for the opened window.
gclose	Close the specified window. Takes one argument, the window handle.
gcloseall	Close all active windows. Takes no arguments.
newcoordinate	Set up new new coordinates system
newwindow	Change the coordinates system for a window. Takes five arguments: handle, x0, y0 (bottom left), x1, y1 (top right).
layer	Set drawing and display layers (0-7 per window). Takes three arguments: handle, index of display layer, index of drawing layer
copylayer	Copy one layer to another. Takes three arguments: handle, index of source layer, index of destination layer.
gsetbgcolor	Set background colour. Takes two arguments: handle, name of colour to set to (string).
gclr	Clear the specified window. Takes one argument: handle.
tclr	Clear the current terminal window
newpencolor	Change pen colour (0-14). Takes two arguments: handle, number (0-14) or colour name (black, white, red, green, blue cyan, magenta, yellow, dimgray, gray, darkred, darkgreen, darkblue, darkcyan, darkmagenta, darkyellow).
newcolor	Specify a new colour (X name). Takes two arguments: handle, name of colour to set to (string).
newrgbcolor	Specify a new colour (RGB values). Takes four arguments: handle, R, G, B.
newlinewidth	Specify line width. Takes two arguments: handle, line width.
newlinestyle	Specify line style. Takes two arguments: handle, line style (0=solid, 1=dotted).
pset	Draw a point. Takes three arguments: handle, xcoordinate, ycoordinate.
drawline	Draw a line. Takes five arguments: handle, x0, y0, x1, y1.
moveto	Move pen to point. Takes three arguments: handle, xcoordinate, ycoordinate.
lineto	Draw line from current pen location to point Takes three arguments: handle, xcoordinate, ycoordinate.
drawrect	Draw a rectangle. Takes five arguments: handle, x, y, width, height.
fillrect	Fill a rectangle. Takes five arguments: handle, x, y, width, height.
drawcirc	Draw a circle (ellipse). Takes five arguments: handle, x, y, xradius, yradius.
fillcirc	Fill a circle Takes five arguments: handle, x, y, xradius, yradius.
drawsym	Draw a symbol Takes three to five arguments: handle, x, y, optional type (1-10), optional size.
drawstr	Print a string on the window at specified position Takes four to five arguments: handle, x, y, string, optional size.
drawnum	Print a number on the window at specified position Takes four to six arguments: handle, x, y, string, optional size, optional number of digits.
gsetnonblock	Set event handling to be nonblocking
ggetch	Get a character from the keyboard if focus on window
ggetevent	Get a keyboard or mouse event for window
ggetxpress	Get a keyboard or mouse event for window