As well as the standard summary statistics described above, PLINK offers some alternative measures such as estimated inbreeding coefficients for each individual and genome-wide identity-by-state and identity-by-descent estimates for all pairs of individuals. The latter can be used to detect sample contaminations, swaps and duplications as well as pedigree errors and unknown familial relationships (e.g. sibling pairs in a case/control population-based sample). PLINK also has functions to detect specific segments shared between distantly-related individuals.

All these analyses require a large number of SNPs!

Pairwise IBD estimation

The pairwise clustering based on IBS, as outlined in the previous section is useful for detecting pairs of individuals who look more different from each other than you'd expect in a random, homogeneous sample. In this section, we consider using the same genotype data to provide a complementary analysis: using estimates of pairwise IBD to find pairs of individuals who look too similar to eachother, i.e. more than we would expect by chance in a random sample.

In a homogeneous sample, it is possible to calculate genome-wide IBD given IBS information, as long as a large number of SNPs are available (probably 1000 independent SNPs at a bare minimum; ideally 100K or more).

plink --file mydata --genome

which creates the file

     plink.genome

which has the following fields:

     FID1      Family ID for first individual
     IID1      Individual ID for first individual
     FID2      Family ID for second individual
     IID2      Individual ID for second individual
     RT        Relationship type given PED file
     EZ        Expected IBD sharing given PED file
     Z0        P(IBD=0)
     Z1        P(IBD=1)
     Z2        P(IBD=2)
     PI_HAT    P(IBD=2)+0.5*P(IBD=1) ( proportion IBD )
     PHE       Pairwise phenotypic code (1,0,-1 = AA, AU and UU pairs)
     DST       IBS distance (IBS2 + 0.5*IBS1) / ( N SNP pairs )
     PPC       IBS binomial test
     RATIO     Of HETHET : IBS 0 SNPs (expected value is 2)

This file will have as many rows as there are unique pairs of individuals in the sample -- for large samples with thousands of individuals, this file can be very large (and take considerable time to generate).

HINT Instead of --genome, using the command --Z-genome will perform the same analysos but create a compressed file, plink.genome.gz. The --read-genome command can directly read compressed files, as of v1.07. This file can be decompressed by the standard gunzip utility, or processed with Unix commands such as zgrep, zcat, etc.

To calculate IBD only for members of the same family (as designated by the PED file), add the command

     --rel-check

(i.e. this greatly speeds up analysis by not considering all possible pairs of individuals, if the goal is to validate known relationships with genetic data).

To create a more verbose version of this file, add the extra command

     --genome-full

which will appended the following extra fields to the normal genome file create a file with the following fields

     IBS0      Number of IBS 0 nonmissing loci
     IBS1      Number of IBS 1 nonmissing loci
     IBS2      Number of IBS 2 nonmissing loci
     HOMHOM    Number of IBS 0 SNP pairs used in PPC test
     HETHET    Number of IBS 2 het/het SNP pairs in PPC test

HINT To produce a smaller version of this file use the command --genome-minimal instead; however, this is only useful if the purpose is to subsequently merge the data using --read-genome-minimal (i.e. when running --cluster or --segment. A disadvantage is that multiple plink.genome.min files cannot be concatenated in the same manner for normal plink.genome files; this will be remedied in future releases of PLINK (i.e. to allow parallel computation of the genome file. Note: as of 1.07, you are advised to use --Z-genome instead of this option -- see above.

HINT In 1.05 onwards, the genome files are indexed by the header row, which must be present. When using --read-genome, the only fields extracted are the four ID fields and DST and PPC when using the --cluster or --mds-plot options. You can therefore extract just these columns, if you do not need the other fields,e.g.

   gawk ' { print $1,$2,$3,$4,$12,$13 } ' plink.genome > new.genome

As mentioned above, the IBD estimation part of this analysis relies on the sample being reasonably homogeneous -- otherwise, the estimates will be biased (i.e. individuals within the same strata will show too much apparent IBD). It is therefore important to run the other population stratification measures provided by plink and other packages before estimating pairwise IBD. In addition, see the notes on the IBS test in the previous section where it is introduced as a constrain on clustering.

HINT To reduce the file size, use the --minX option to only output to plink.genome pairs where PI_HAT is greater than X. That is,

plink --file mydata --genome --min 0.05

will only display the pairs of individuals showing reasonably high levels of IBD sharing (i.e. typically it will be these pairs that are of interest, rather than the vast majority of pairs that show no excess sharing).

Hint Calculating the average pi-hat for each individual and looking for outliers is also useful (in particular, sample contamination will lead to too many heterozygote calls, which leads to fewer IBS 0 calls, which leads to over-estimated IBD with all other people in the sample). Be sure to set --min 0 and --max 1 in this case to obtain pairs for all individuals.

Advanced hint If you have access to a cluster, use the --genome-lists option to facilitate parallelization, as described in the IBS clustering section.

Inbreeding coefficients

Given a large number of SNPs, in a homogeneous sample, it is possible to calculate inbreeding coefficients (i.e. based on the observed versus expected number of homozygous genotypes).

plink --file mydata --het

which will create the output file:

     plink.het

which contains the fields, one row per person in the file:

     FID       Family ID
     IID       Individual ID
     O(HOM)    Observed number of homozygotes
     E(HOM)    Expected number of homozygotes
     N(NM)     Number of non-missing genotypes
     F         F inbreeding coefficient estimate

This analysis will automatically skip haploid markers (male X and Y chromosome markers).

Note With whole genome data, it is probably best to apply this analysis to a subset that are pruned to be in approximate linkage equilibrium, say on the order of 50,000 autosomal SNPs. Use the --indep-pairwise and --indep commands to achieve this, described here.

Note The estimate of F can sometimes be negative. Often this will just reflect random sampling error, but a result that is strongly negative (i.e. an individual has fewer homozygotes than one would expect by chance at the genome-wide level) can reflect other factors, e.g. sample contamination events perhaps.

Runs of homozygosity

A simple screen for runs of homozygous genotypes within any one individual is provided by the commands --homozyg-snp and --homozyg-kb which define the run in terms of the required number of homozygous SNPs spanning a certain kb distance, e.g.

The algorithm is as follows: Take a window of X SNPs and slide this across the genome. At each window position determine whether this window looks 'homozygous' enough (yes/no) (i.e. allowing for some number of hets or missing calls). Then, for each SNP, calculate the proportion of 'homozygous' windows that overlap that position. Call segments based on this metric, e.g. based on a threshold for the average.

The exact window size and thresholds, relative to the SNP density and expected size of homozygous segments, etc, is obviously important: sensible default values are supplied for the context of dense SNP maps, scanning for large segments. In general, this approach will ensure that otherwise long runs of homozygosity are not broken by the occassional heterozygote. (For more accurate detection of smaller segments, one might consider approaches that also take population parameters such as allele frequency and recombination rate into account, in a HMM approach for example: but for now, PLINK only supports this basic detection of long, homozygous segments).

To run this option with default values, use the command

plink --bfile mydata --homozyg

which generates a file

     plink.hom

The plink.hom file has the following format, one row per identified homozygous region:

     FID      Family ID
     IID      Individual ID
     CHR      Chromosome
     SNP1     SNP at start of region
     SNP2     SNP at end of region
     POS1     Physical position (bp) of SNP1
     POS2     Physical position (bp) of SNP2
     KB       Length of region (kb)
     NSNP     Number of SNPs in run
     DENSITY  Average SNP density (1 SNP per kb)
     PHOM     Proportion of sites homozygous
     PHET     Proportion of sites heterozygous

The options to change the behavior of this function (along with the default values as parameters) are given below.

To change the definition of the sliding 'window', use the options

    --homozyg-window-kb 5000
    --homozyg-window-snp 50

To change the number of heterozygotes allowed in a window

    --homozyg-window-het 1

To change the number of missing calls allowed in window, e.g.

    --homozyg-window-missing 5

To change the proportion of overlapping windows that must be called homozygous to define any given SNP as 'in a homozygous segment', use

     --homozyg-window-threshold 0.05

(i.e. this number is relatively low, so that SNPs at the edge of a true segment will be called, as long as the windows are sufficiently large, such that the probability of a window being homozygous by chance is sufficiently small).

The above options define the 'window' that slides across the genome; the options below relate to the final segments that are called as homozygous or not:

    --homozyg-snp 100
    --homozyg-kb 1000

You can also specify the required minimum density (in kb, i.e. 1 SNP per 50 kb)

    --homozyg-density 50

Finally, if two SNPs within a segments are too far apart (measured in kb), that segment can be split in two:

    --homozyg-gap 1000

HINT As is, this analysis should be performed on sets of SNPs that have been pruned for strong local LD (if the goal is to find long segments that are more likely to represent homozygosity by descent (i.e. autozygosity) rather than simply by chance).

To obtain pools of overlapping and potentially matching segments, we can use --homozyg-group in addition to the above, which generates the file

     plink.hom.overlap

which contains the fields

     FID     Family ID
     IID     Individual ID
     PHE     Phenotype of individual
     CHR     Chromosome
     SNP1    SNP at start of segment
     SNP2    SNP at end of segment
     BP1     Physical position of start of segment
     BP2     Physical position of end of segment
     KB      Physical size of segment
     NS      Number of segments in the pool that match this one
     GRP     Allelic-match grouping of each segment

For example, the command

plink --file test --homozyg --homozyg-group

might result in the file plink.hom.overlap containing:

 FID  IID   PHE  CHR SNP1 SNP2        BP1      BP2    KB   NS    GRP
   1    1     2    1 snp1 snp7   1000000   7000000  6000   1     1
   6    1     1    1 snp1 snp5   1000000   5000000  4000   1     1*
   2    1     1    1 snp2 snp7   2000000   7000000  5000   0     2*
 CON    3   1:2    1 snp2 snp5   2000000   5000000  3000

This implies one pool (i.e. each pool is separated by a CON (consensus row) and a space. CON is the consensus region; the ratio is the case:control segment ratio; under IID we have the number of individuals.

When there is more than one pool, they are ordered by the number of segments in the pool, then physical position. To output only pools of a particular size, use the --pool-size N option (e.g. --pool-size 10 to only display pools with at least 10 segments).

A pool contains overlapping segments, which may or may not also allelically match. For allelic matching, segments are compared pairwise, and a match is declared if at least 0.95 of jointly non-missing, jointly homozygous sites are identical. This threshold can be changed with the option

     --homozyg-match 0.99

The number of other segments in the pool that allelically match each segment is shown in the NS field. The GRP field shows how PLINK attempts to group allelically-matching segments within the pool of overlapping segments. It works as follows:

For each segment, find the number of other segments that match (NS).
Find segment with largest NS, denote as group 1, and put a * to indicate this is the index for this group.
Denote all other segments that match this index as being in GRP 1 (i.e. but without the *)
Continue to next ungrouped segment (2*, etc)

By default, we compare all segments pairwise when asking if they match; if the --consensus-match flag is given, then for a pool of overlapping segments, matches are defined only on the basis of the consensus region (i.e. the overlapping region shared by all segments). This is probably not very sensible in many cases, as the consensus region can often be small (i.e. a single SNP).

The NS field can suggest any intransitivity in matching: e.g. if B matches A and C but A does not match C, then if B has already been grouped with A, C would not be added to this group as being an allelic match. In this case C would have NS > 0 but belong to a GRP of its own.

Internally, all pools are formed but then pruned if, for instance, a smaller pool is included in a larger pool completely. That means that in certain circumstances you will see a segment in more than one pool. For example, imagine a grid with three people A, B and C along the columns, each row representing physical position, and the presence of a letter representing a homozygous run:

     . . .
     A . .
     A B .
     A B C
     A B C
     A . C
     A . .
     . . .

In this case, {A,B} and {A,C} and {B,C} pools will not be displayed, as they appear in the super-pool {A,B,C}. However, if we instead had:

     . . .
     A . .
     A B .
     A B .
     A . .
     A . C
     A . C
     A . .
     . . .

Then you will see {A,B} and {A,C} (i.e. with A shown twice) as we have two distinct consensus regions here.

Finally, if the --homozyg-verbose option is added, the plink.hom.overlap file will then display the actual segments underneath each pool. Here each individual is listed across the page, so the 3 columns refers to the 3 segments in the pool. For example:

plink --file test --homozyg-snp 2 --homozyg-group --homozyg-verbose

now generates plink.hom.overlap as follows (with annotation added in italics):

     FID  IID      PHE  CHR SNP1 SNP2    BP1    BP2       KB   NS    GRP
       1    1        2    1 snp1 snp7      1      7    0.006   1     1
       6    1        1    1 snp1 snp5      1      5    0.004   1     1*
       2    1        1    1 snp2 snp7      2      7    0.005   0     2*
     CON    3      1:2    1 snp2 snp5      2      5    0.003

     SNP    1    6    2          <-- Family ID
            1    1    1          <-- Individual ID
            1    1    2          <-- GRP code

    snp1 [A/A] [A/A]  C/A        <-- now SNPs are listed down the page
                                 
    snp2 [A/A] [A/A] [C/C]       <-- start of consensus region
    snp3 [A/A] [A/A] [C/C]
    snp4 [A/A] [A/A] [C/C]
    snp5 [A/A] [A/A] [C/C]       <-- end of consensus region
                                 
    snp6 [A/A]  A/C  [C/C]
    snp7 [A/A]  A/C  [C/C]

A bracket indicates that that genotype is part of the homozygous segment: the consensus region is the intersection. The entire union of SNPs is displayed and the consensus region is indicated by spaces before and after. i.e. the consensus region is that where all genotypes are in [brackets].

Obviously, this file can get quite large (+wide) with real data and it is not very machine-readable.

Segmental sharing: detection of extended haplotypes shared IBD

WARNING This analysis is still in the beta development stage and is considerably more involved than many others provided by this package: currently, you should only perform these analyses if you consider yourself something of an analytic expert and are confident you will be able to interpret the output! Over time, we expect that the documentation and features supporting this analysis will improve.

There are five important steps to this analysis:

Obtain a homogeneous sample
Remove very closely related individuals
Prune SNP set
Detect segments
Associate with disease

Check for a homogenous sample

This analysis requires that all individuals belong to a single, homogeneous population. To ensure this assumption is reasonable: as described here, first run

plink --bfile mydata1 --genome

to generate a plink.genome file. This will be used subsequently in a number of steps.

Then, using the available tools, such as listed here and described more fully in the section on stratification, obtain a relatively homogeneous dataset. Some relevant options are listed here:

     --cluster    (cluster individuals)
     --matrix     (generate .mdist file, used to generate MDS plots)
     --ppc        (threshold for PPC test, not to cluster individuals)
     --mds-plot   (generate a multidimensional scaling plot)
     --ibs-test   (as case/control less similar on average?)
     --neighbour  (option to find individual outliers)

Also, remove individuals who appear to have higher levels of inbreeding than expected (see above). If you have a set of individuals you want to exclude from analysis based on these steps, for example, listed in the file outliers.txt (FID, IID) then use:

./plink --bfile mydata1 --remove outliers.txt --make-bed --out mydata2

Remove very closely related individuals

The focus of this analysis is to look for extended haplotypes shared between distantly related individuals: having very closely related individuals (siblings, first cousins, etc) will likely swamp the results of the analysis. Scan the plink.genome file for any individuals with high PIHAT values (e.g. greater than 0.05). Optionally, remove one member of the pair if you find close relatives. (Alternatively, to keep them in but just exclude this pair from the segmental analysis, see below).

Prune the set of SNPs

The segmental sharing analysis requires approximately independent SNPs (i.e. linkage equilibrium). Two options to prune are documented here.

A reasonable strategy might be as follows:

plink --bfile mydata2 --mind 1 --geno 0.01 --maf 0.05 --make-bed --out mydata3

followed by

plink --bfile mydata3 --indep-pairwise 100 25 0.2

followed by

plink --bfile mydata3 --extract plink.prune.in --make-bed --out mydata4

Detecting shared segments (extended, shared haplotypes)

With a newly pruned fileset, ideally containing only independent, high quality SNPs in individuals who are not very closely related but are from the same population, run the command

plink --bfile mydata4 --read-genome plink.genome --segment

PLINK expects the 3rd column the MAP/BIM file to contain genetic distances in Morgan units. A reasonable approximation is to scale from physical position (i.e. column 4) at 1cM=1Mb. If the genetic distances are in cM instead of Morgans, add the --cm flag.

To set threshold on who to include/exclude based on genome wide IBD use

     --min 0.01
     --max 0.10

For example, this would exclude pairs who share >10% of their genomes. Alternatively, to include all pairs, irrespective of whether we estimate any genome-wide sharing or not, add the option

     --all-pairs

instead. This will use all pairs, allowing for a small prior probability of sharing for pairs that otherwise are at the boundary of IBD sharing (i.e. sharing 0% IBD). Naturally, for a large sample, it may become prohibitive to consider all possible pairs.

The --segment option generates a file

     plink.segment

which has the fields:

     FID1       Family ID of first individual
     IID1       Individual ID of first individual
     FID2       Family ID of second individual
     IID2       Individual ID of second individual
     PHE        Phenotype concordance: -1,0,1
     CHR        Chromosome code
     BP1        Start physical position of segment (bp)
     BP2        End physical position of segment (bp)
     SNP1       Start SNP of segment
     SNP2       End SNP of segment
     NSNP       Number of SNPs in this segment
     KB         Physical length of segment (kb)

Here one row represents one segment. The PHE field is coded -1,0,1 for control/control, case/control, or case/case pairs respectively.

The option

     --segment-length 2000

means to only select segments that are at least 2000 kb in length, for example. The option

     --segment-snp 100

means only to select segments that contain at least 100 SNPs, for example.

For ease of interpretation, and to increase the probably that the segments are truly shared IBD and thus tags shared rare variation between two individuals, it makes sense to restrict ones focus to very extended segments (e.g. over 1Mb in size, for example).

Another option is the --segment-group option, which generates output similar to --homozyg-group, described above; similarly, --segment-verbose prints out the actual genotypes for the individuals that overlap. However, these can be large files that are not necessarily easy to interpret.

Association with disease

Along with the --segment option, as above, if you also add:

     --mperm N

then, for case/control data, this performs a test of whether segments stack up more in case/case pairs versus non-case/case pairs at any position, performing N permutations. Equivalently, you can use an already-created segment file:

./plink --bfile mydata4 --read-segment plink.segment --mperm 10000

This will generate two files:

     plink.segment.summary

which contains one row corresponding to one SNP; there are five fields:

     CHR    Chromosome code
     SNP    SNP identifier
     CONU   Number of control/control segments over this SNP
     DISC   Case/control segments spanning this position 
     CONA   Case/case segment count

The file

     plink.segment.summary.mperm

contains empirical significance values for each position, asking whether there is a higher rate of case/case sharing than expected. It is important to note that the test statistic is still under developemtn: in this current release, it should merely be interpreted as a rough guide to the data. Naturally, the thresholds for declaring significance will be much lower than for genome-wide association analysis; precise guidelines will be put in place presently.

This document last modified Wednesday, 25-Jan-2017 11:39:27 EST