==========================
= SPLINK (version 1.08) =
==========================
Author
======
David Clayton
(MRC Biostatistics Unit, Cambridge)
Description
===========
SPLINK is a program for linkage analysis using affected sib pairs. It
uses the method of maximum likelihood to estimate the probability that
two affected sibs share 0, 1 or 2 autosomal marker haplotypes
identically by descent (IBD). For marker haplotypes on the
X-chromosome, IBD sharing of maternal haplotypes only is
considered. Parental marker data need not be present (although, of
course, they help) and neither IBD status or phase of marker haplotype
need be known. When full parental data are not available, the
uncertainty in assigning parental haplotypes may be reduced by
including data on sibs who are either unaffected or of unknown
affection status. When all available data still permit ambiguity for
parental haplotypes, the population haplotype frequencies become
relevant, and these are estimated internally, again using the method
of maximum likelihood.
Affected sib trios (and larger groupings) are dealt with approximately
by treating all pairwise comparisons as independent (see technical note below).
Since this device does not correspond to a genuine likelihood, in this
more general case the method should be referred to as maximum pseudo-
likelihood. It has been become conventional to weight the individual
sibpair comparisons in larger sets of affected sibs by 2/A where A is
the number of affected sibs in the family. However the grounds for this
are not clear and there may be good reasons to prefer the unweighted
approach.
Denoting the IBD sharing probabilities by z0, z1 and z2, the null
hypothesis for autosomal markers corresponds to
z0 = 0.25, z1 = 0.5, z2 = 0.25
and likelihood ratio chi-squared tests are carried out against either
unrestricted alternatives, or against the "possible triangle"
alternatives
z1 > 0.5, z0 > 0.5*z1
The former is a 2 degree of freedom test while the latter is a
compromise between 1 and 2 degree of freedom tests.
The analogous "score" tests based upon gradients of the log (pseudo-)
likelihood surface at the null values of z0, z1 and z2 may also be
calculated. Additionally, 1 degree of freedom tests imposing one of
the restrictions
z1 = 0.5
or
z1 = 2*z0
may also be selected. Three variants of the score test are possible,
based upon different estimates of the variance-covariance matrix of the
score vector. These estimates are:
(a) a "naive" estimate based on the assumption that the pseudo-likelihood
behaves as a true likelihood (that is, the variance is estimated from
the curvature of the log pseudo-likelihood surface),
(b) a "theoretical" estimate, and
(c) an "empirical" estimate, derived from the empirical variance-covariance
matrix of the contributions from each nuclear family to the score vector.
There is no difference between the estimates (a) and (b) when there are
are only sib pairs or when sib pairs from larger collections of affected
sibs are unweighted. These score tests are experimental and remain to be
fully evaluated, but in principle offer the best solution to the problem
of multiple sib pair comparisons in the same nuclear family.
For all test statistics, asymptotic p-values are calculated from the
relevant chi-squared distribution(s).
It should be noted that the p-values calculated from the methods above
rely on large sample asymptotic theory and may not be valid for smaller
sample sizes. For example, we would need at least 20 fully informative
sib pairs (that is, pairs in which ibd sharing status can be assigned
unequivocally) in order to observe expected frequencies no smaller
than 5 for all ibd sharing categories. In classical contingency table
analyses, this is usually regarded as the minimum for reliable use of chi-
squared tests. The same considerations apply to traditional thresholds for
lod scores, which are based on the same large sample theory. For the data
which can be used by SPLINK, however, it can be difficult to see, in effect,
how large the sample is, since some observations may contribute little
information. Because of this, the program calculates the "effective sample
size", equivalent to the number of fully informative sib pairs which would
carry the same information. Tentatively I would urge caution if this value
falls below 20. For the technically inclined, the effective sample size
is calculated by equating the determinant of the information matrix,
evaluated at the null hypothesis (z0, z1, z2) = (0.25, 0.5, 0.25), with
that for a study consisting entirely of fully informative sib pairs.
Another attempt to address the problem of small samples is an optional
"bootstrap" testing procedure. This is a Monte Carlo method which involves
creating many datasets at random by resampling the observed families
with replacement. Thus, if there are n families, each bootstrap sample also
contains n families, sampled from the original families with replacement.
Note that the bootstrap datasets will contain repeated families. For
each bootstrap dataset, we calculate whether it would support
linkage by calculating an estimate of (z0, z1, z2) subject to the possible
triangle restriction. A p-value can be obtained by doing this a large number
of times (for example 1000 or 10000 times) and counting the proportion of
occasions on which the estimate falls on the null value (0.25, 0.5, 0.25).
These calculations are quite rapid. This procedure is based on Bayesian ideas
and has yet to be properly validated.
For markers on the X-chromosome, sharing of maternal haplotypes only is
considered so that z2=0 and, under the null hypothesis
z0 = z1 = 0.5
As in the autosomal case, likelihood ratio tests and score tests may be
carried out, the latter using any of three variance estimators outlined
above. In this case, the "possible triangle" restriction degenerates to
the one-sided alternatives z1<0.5.
Optionally an output file may be created which contains the IBD
assignment scores for each family. A program is available, written in
the S programming language, which allows these to be plotted on a
triangular coordinate system. Any user of this must have access to
the S-plus statistical package.
A second output option allows the estimated locus (haplotype)
freqencies to be written to a file for later input to other programs. This
has the correct format for inclusion in the DATAFILE for the LINKAGE
programs.
A final option is to write a brief one-line summary to an output file.
This line lists heterozygosity, PIC score, LOD score, and estimates of
z0, z1 and z2. Fields are delimited by tabs.
Data input
==========
The data input file should contain, for each person, the following
blank-delimited fields:
family family or pedigree code (alphanumeric)
id person's identifier within family (alphanumeric)
father id of father (who must have the same family code)
mother ditto for mother
sex sex (2=Female, 1=Male)
affected disease status (2=affected, 1=unaffected, 0=unknown)
marker 1 coded a/b, where a and b are the two alleles. Alleles
must be coded as consecutive integers, with 0
representing unknown. Thus 0/0 represents completely
missing data but, for a biallelic marker, 2/0
represents either 2/1 or 2/2. For markers on the X
chromosome, males should have marker phenotypes coded
a/0 or a/a, so that males and females have equal length
records.
marker 2 ditto
... etc.
Although these fields must appear in the specified order, persons need
not appear on the file in any particular order. Note that parents must
be included in the data file even if no data concerning them are
available; such entries are necessary to correctly identify sibships.
Persons who appear on the data file only as parents do not need to have
valid entries in the "mother" and "father" fields. A single period (.)
is recommended for the coding of these, but any identifyier which does
not occur in the family will have the same effect. The disease status
of parents is not used by the program and may be coded as 0.
Data input is via the standard data input stream and may be fed into
SPLINK either via a filter program, or by using the > operator on the
command line, for example:
splink >input.dat
It is envisaged that the input data will be extracted from a larger
database, and it should not prove too difficult to achieve the above
format. An alternative is to use Linkage PEDFILEs as input since these
files have the same basic structure even though they contain a few extra
fields. A filter program "ped2spl" is supplied which converts Linkage
PEDFILEs into a form suitable for input to splink (see below).
Output
======
Output is to the standard output stream, but may be saved to
a file using the < command line operator:
splink >input.dat <output.lst
The optional output of family ibd scores is controlled by the -o flag
(see below).
Flags
=====
A number of flags control program operation. In the description
that follows, the # character represents the optional value to be
assigned to the flag. The value represented by # must follow the flag
directly. If this value does not commence with the - character, there may
be blank space between a flag and its value; alternatively the flag and its value
may be run together . Logical flags are set by simply including them on the command
line. If a flag, eg -mf, is set by default, it can be unset by either writing
-nomf or -mf- .
-asp Only include families with, at least, one affected sib
pair. Default is ON. Switch this option OFF if you want
such families to contribute to estimation of the
haplotype frequencies.
-bs# Do bootstrap significance test, using # bootstrap samples.
By default, bootstrap testing is not performed.
-c Only include "complete" families in the analysis (ie,
families with both parents typed).
-d# Write locus data to file #. If the filename is preceded by the
+ sign (with no intervening space), the output is appended to
the file. If # is absent, the file datafile.dat is used
-f# Specify maximum number of (nuclear) families.
(Default -f500)
-h Help (list command line options)
-i# Specify # as an identifier string to be used to label
output. If not specified, the string ??? is used
-l# Specify number of marker loci. If missing, it is assumed
that this number is given as the first item on the input
file.
-mf Allow multiple nuclear families from one pedigree
(although the relationship between these families will be
ignored). If not set, only the first nuclear family
encountered in each pedigree is used. Default is for this
flag to be set, but it may be unset by -nomf or -mf- .
-n# Specify maximum number of persons on data file.
(Default -n2000)
-o# Specifies that the IBD assignment scores for each
family will be written to file #. By default this
option is turned off.
-O# Specifies amount of output. #=0 is minimum output (results
only), #=1 identifies subjects/pedigrees omitted for one reason
or another, and #=2 lists all pedigrees. Default is 2.
-pf# Write pedigrees used in the analysis to file #.
-pt Specifies that estimation of z0, z1, z2 will be
restricted to values compatible with linkage. For an
autosomal marker this is the "possible triangle" and for
an X marker, z1 is constrained to be no smaller than z0.
Default is for the flag to be set, but it can be unset
by -nopt or -pt- .
-rs# Specify an integer seed, #, for the pseudo-random number
generator. This is used in bootstrap testing (see -bs).
-s# Selects score test. Permitted values of # are,
for autosomal markers:
2 2 df test
pt 1~2 df "possible triangle" test
1d 1 df test with the constraint z1=0.5
1r 1 df test with the constraint z1=2*z0
for X markers (only maternal haplotype scored):
m 1 df test of null z0=z1=0.5
m1 1 df 1-sided test for z1<z0 (z2=0)
More than one -s option may be selected, but the first
three are only legal for autosomal markers and the last
two for X markers. By default, if the -pt flag is set,
the "pt" test is selected for an autosomal marker and
the "m1" test is selected for an X marker. All other
tests are unselected. If the -pt flag is unset, the default
score tests are the "2" and "m" options
-S# Write brief one-line summary to file #. If the filename is
preceded by a + sign (with no intervening space), the output
is appended to the file. If # is absent, the file splink.lst
is used
-x# Specify maximum allowable ambiguity for parental
haplotyping. If there are more than # possible
parental haplotype assignments, the family is excluded.
Note that, for speed, possible parental haplotypes are
stored in dynamically allocated memory and the -x
option will help if you run out of memory.
(Default -x1000)
-X The marker loci are on the X-chromosome.
-wt Weight multiple sibpair pair comparisons within the same
nuclear family by 2/A, where A is the number of affected
sibs. Default is for such comparisons to be weighted.
-v# Specifies variance estimator to be used for score tests.
Values of # are n, t, or e for naive, theoretical and
empirical estimates respectively. More than one -v option
can be specified.
(Default -vn)
Example:
=======
splink >infile.dat -l2 -s1d -x500 -oibdscore.dat
Differences from previous version
=================================
The main features introduced in version 1.0 are the completely revised input
data format, and the upgraded score test facilities.
1.01
The weighting of multiple sib pair comparisons was made optional and
a "theoretical" variance estimate of the score vector was included.
1.02
Internal changes only, aimed at fixing problems with integer arithmetic
encountered on short word length machines like the IBM PC.
1.03
a) -X option implemented
b) This has necessitated checking of the sex data for internal consistency.
earlier versions of the program ignored the sex field and allowed mothers
and fathers to be interchanged with impunity.
c) A bug in the "possible triangle" maximization, affecting some problems
in which the restricted maximum of the likelihood corresponds with the null
hypothesis, has been fixed.
d) Some dubious code controlling the evaluation of all possible parental
haplotype configurations has been rewritten. The function in question is
Haplotype::next(). This has had little effect on the results of the test
examples run by me.
1.04
a) Bootstrap testing introduced.
b) Equivalent sample size for fully informative pairs introduced.
c) Better checking of input data format when the number of loci is read
from the input data file.
1.05
a) Optional output via -S and -d flags
b) Labelling of output via -i flag
1.06
No useful additions, but some corrections.
1.07
Introduction of -c flag to restrict analysis to affected sib-pairs with both
parents typed. Earlier versions suggested using -x1 to achieve this result.
Peter Holmans has pointed out that this leads to a bias towards 0-sharing,
since, in the absence of parental data, 0-sharing can more frequently be
inferred from offspring data alone.
1.08
a) Introduction of -pf flag to facilitate comparison of results with other
programs.
b) Improvement of documentation to clarify the "likelihood" used (see
technical note below).
c) Improvement to command line syntax
Known (or Suspected) Bugs
=========================
(1) The "theoretical" variance estimate of the score test doesn't always
perform well when using the -wt option. Sometimes the estimated variance
matrix is not positive semidefinite and it is impossible to compute a
test statistic. This may well be related to problems with the "likelihood"
(see below) and this variance estimate should be treated with caution.
Compiling:
=========
Most of SPLINK is written in C++ and must be compiled using a suitable C++
compiler. The files splink.C (or splink.cpp) and splfun.C (or splfun.cpp)
are C++ source files and cline.c, gamma.c, and random.c are plain C source
files. Finally there are "header" files splink.h, cline.h, and myrand.h
which contain class definitions etc..
In unix, compilation would normally be by:
CC splink.C splfun.C cline.c gamma.c random.c -lm -o splink
=========================================
= Technical Note: The Pseudo-likelihood =
=========================================
1. Fully informative families
=============================
To understand the pseudo--likelihood approach taken by SPLINK, first
consider the simple case of an affected sib pair with both parents
where there are 4 different alleles in the family:
(a,b) (c,d)
| |
+-----+-----+
|
+-----+-----+
| |
(?,?) (?,?)
The conventional likelihood approachto this problem, as advocated by Risch,
is to score the sibling pair for 0-, 1- or 2-IBD and to assign a
likelihood contribution of either z0, z1 or z2 accordingly. However it
is important to recognise that these likelihood contributions are NOT
probabilities of the observed data, although they are the correct
likelihood contributions in this simple case. To deal with more
complicated problems we need to be able to calculate the likelihood
properly, by calculating the probability of the data under the assumed
model for linkage.
Conditional upon parental genotypes, there are 16 possible offspring
genotypes. Of these, 4 or 0-IBD, 8 are 1-IBD, and 4 are 2-IBD. In the
absence of association, all sibling genotypes within these IBD groups
are equally probable. Thus the (conditional) probability of a pair of
sibling genotypes that have 0-IBD is z0/4, for pairs with 1- and 2-IBD
the corresponding probabilities are z1/8 and z2/4 respectively. The
uncondiitonal probability of the observation is obtained by multiplying
these conditional probabilities by the probability of the parental
genotype (for which we assume Hardy-Weinberg equilibrium and random
mating). Thus the full likelihood contribution is
Prob(Parental genotype) * { z0/4 OR z1/8 OR z2/4 }
OR
Prob(Parental genotype) * (1/16) * { 4*z0 OR 2*z1 OR 4*z2 }
Now, of course, this turns outs to have exactly the same effect as
scoring either z0, z1 or z2 since all we have done is to multiply by
terms that do not depend on the z's. However doing it more carefully is
important when one considers data which have ambiguous IBD status. Note that
the middle term is the probability of any one inheritance vector under the
hypothesis of no linkage.
2. Homozygous parents
=====================
The simplest case in which IBD status is ambiguous is the case where
one parent is homozygous at the marker locus. For example, what is the
appropriate likelihood contribution for the following family?
(a,a) (c,d)
| |
+-----+-----+
|
+-----+-----+
| |
(a,c) (a,d)
Since the two `a' alleles are indistinguishable, the data are consistent
with any one of four possible inheritance patterns. Two of these
are 1-IBD and the other two are 0-IBD, Thus the likelihood contribution
from this family is
Prob(Parental genotype) * (1/16) * {2*4*z0 + 2*2*z1}
Note that in this case and in the fully informative case considered
above, the likelihood factorizes --- the first term depends only on
the marker allele frequencies and the second term depends only on the
IBD sharing probabilities.
3. Missing parental genotypes
=============================
When parental genotypes are missing, the likelihood contribution must
be calculated by adding terms such as those shown above over all
possible parental genotypes. The number of possibilities can be reduced
by considering other sibs. Here we get into slightly murky water --- if
an additional sibs is unaffected, then the probability of his or her
genotype depends (rather weakly) on the z's. In SPLINK this is
ignored, and this equivalent to treating the extra sib as having
unknown disease status. This can rather loosely be justified for
diseases with late onset on the grounds that an unaffected sibling
might go on to develop the disease. However this doesn't hold up to
close scrutiny; if the sibling is very young, then it makes sense to
treat him as having unknown disease status, while an older sib is
almost certainly truly unaffected. In practice the real justification
for ignoring the contributions of unaffected siblings to the likelihood
for the z's is that, if the penetrance is low, the contributions are
extremely small and can be ignored.
Even when parents are typed, the parental genotype will usually be
unknown when we are combining several close markers to form a marker
haplotype, since the phase of the markers will be unknown. The
probability of each of the possible parental genotype assignments now
depends on haplotype frequencies, and calculation of the likelihood
contribution now involves summation over all possible phases.
The above discussion explains the way by which SPLINK calculates
likelihoods for families which include an affected sib pair. With the
proviso discussed in the last paragraph, these are true likelihoods and
there should be no concern about assuming the usual asymptotic
distribution results. In particular, the variance-covariance matrix of
the score vector can be consistently estimated using the information
matrix. By dividing the determinant of the information matrix for the
z's by the theoretical contribution of one fully informative family
(that is a family of the type we considered at the start of this
discussion), we can calculate an ``effective'' sample size. Note that,
in the absence of some parental genotypes, the factorization of the
likelihood into separate parts for the marker haplotype frequencies and
for the IBD sharing probabilities is lost and the two sets of
parameters are not totally orthogonal. However, in practice the
interdependence of the inferences seems to be so small as to be
ignorable. For calculating score tests, SPLINK assumes approximate
parameter orthogonality in order to avoid the necessity for inverting
large matrices.
4. Affected sib trios etc.
==========================
Things get really difficult when we consider trios and larger groupings
of affected siblings. The existing literature on this problem really
only deals with the fully informative family situation with which we
began the discussion. That is, consider the following situation.
(a,b) (c,d)
| |
+-----+-----+
|
+----------+----------+
| | |
(?,?) (?,?) (?,?)
There are now 3 possible sib pair comparisons. The general
recommendation seems to be to pretend that they are independent and
then, possibly, downweight them to compensate for the interdependence
of the comparisons. A weight of 2/3 is usually advocated by analogy with
comparisons of a quantitative trait. This is equivalent to scoring
the likelihood contribution as
{ Product( z0 OR z1 OR z2 ) }^w
where the product is taken over all three possible sibpair comparisons,
and w is the weighting.
The case for and against weighting is interesting in this situation. It
is quite easy to show that the score contributions from each sib pair
comparisons are quasi-independent under the null hypothesis, in the
sense that they have zero product moment correlation. Thus a score test
based on an unweighted likelihood has the correct asymptotic properties.
The same argument shows that, with weights > 1, the information matrix
based on this "likelihood" overestimates the variance of the score
vector. Thus tests will be conservative.
Things really get difficult when we consider trios in the case where
there is ambiguity of the IBD status of the comparisons. To construct a
likelihood we need to add the probabilities for all possible
configurations. In general these cannot be written down as a function of
(z0, z1, z2) alone and some sort of "pseudolikelihood" argument must be used.
By analogy with the case of fully informative families,
SPLINK takes as the probability of each configuration a downweighted
product of terms from each sib pair comparison:
Prob(Parental genotype) * (1/64) * { Product( 4*z0 OR 2*z1 OR 4*z2 ) }^w
As before, the middle term is the probability of any one inheritance vector
under the hypothesis of no linkage. The last term modies the probabilities
in a manner which mimics independence of all sib-pair comparisons.
On reflection this is, perhaps, not such a good idea. It certainly does
not yield the true probability of the data given the model. Indeed it
does not even yield an expression with the properties of a probability.
To demonstrate this, return to the fully informative family with an
affected sib trio. Conditional upon parental genotypes there are 64
possible configurations of the inheritance vectors for siblings, and these
fall into just four IBD patterns. These are shown below, together with the
number of ways they can occur and the "probability" assigned to them by the
(unweighted) rule described above. The "probability" shown omits the
probability of the parental genotype, and represents a conditional
"probability" of genotype of offspring conditional upon parental genotype.
IBD pattern # Configurations "Probability"
=========== ================ =============
0, 0, 2 12 (z0^2 * z2)
0, 1, 1 24 (z0 * z1^2)/4
1, 1, 2 24 (z1^2 * z2)/4
2, 2, 2 4 (z2^3)
The total "probability" when summed over all possible configurations of the
inheritance vectors is 1 when there is no linkage (z0=z2=0.25, z1=0.5),
but not otherwise. Downweighting does not help. Also, the marginal
probabilities for any one sib pair sharing 0, 1 or 2 haplotypes IBD are not
z0, z1 and z2!
The above discussion shows clearly that the function maximized by SPLINK
cannot be properly described as a likelihood when trios and larger
groupings of affected siblings are present. Indeed, when there is
uncertainty about IBD status of the comparisons within such groupings,
the function may behave rather stangely. An extreme example is provided
by including an affected sibling with completely unknown genotype; such
individuals should not affect the analysis but do so with the present
computations.
In these situations, the function maximized is a pseudo-likelihood and
asymptotic likelihood properties cannot safely be assumed. I am aware of
this problem and it provides the main motivation for the empirical estimates
of the variance of the score statistics and for experiments with
bootstrap significance testing. I am also working on theoretical estimation
of the variance of the score statistic, but am not 100% confident of the
accuracy of the current method (see "known bugs").
================
= Availability =
================
These programs are available via anonymous ftp from
ftp.mrc-bsu.cam.ac.uk
on the directory /pub/genetics.
References
==========
Holmans, P. (1993) Asymptotic properties of affected sib-pair linkage
analysis. Am J Hum Genet 52:362-374.
Holmans, P. and Clayton, D. (1995). Efficiency of typing unaffected
relatives in an affected sib-pair linkage study with single locus and
multiple tightly-linked markers. Am J Hum Genet 57:1221-1232.