PLINK: Whole genome data analysis toolset plink...
Latest PLINK release is v1.07 (10-Oct-2009)

Whole genome association analysis toolset

Introduction | Basics | Download | Reference | Formats | Data management | Summary stats | Filters | Stratification | IBS/IBD | Association | Family-based | Permutation | LD calcualtions | Haplotypes | Conditional tests | Proxy association | Imputation | Dosage data | Meta-analysis | Result annotation | Clumping | Gene Report | Epistasis | Rare CNVs | Common CNPs | R-plugins | SNP annotation | Simulation | Profiles | ID helper | Resources | Flow chart | Misc. | FAQ | gPLINK

1. Introduction

2. Basic information

3. Download and general notes

4. Command reference table

5. Basic usage/data formats 6. Data management

7. Summary stats 8. Inclusion thresholds 9. Population stratification 10. IBS/IBD estimation 11. Association 12. Family-based association 13. Permutation procedures 14. LD calculations 15. Multimarker tests 16. Conditional haplotype tests 17. Proxy association 18. Imputation (beta) 19. Dosage data 20. Meta-analysis 21. Annotation 22. LD-based results clumping 23. Gene-based report 24. Epistasis 25. Rare CNVs 26. Common CNPs 27. R-plugins 28. Annotation web-lookup 29. Simulation tools 30. Profile scoring 31. ID helper 32. Resources 33. Flow-chart 34. Miscellaneous 35. FAQ & Hints

36. gPLINK
 

Data management tools

PLINK provides a simple interface for recoding, reordering, merging, flipping DNA-strand and extracting subsets of data.

Recode and reorder a sample

A basic, but often useful feature, is to output a dataset:
  1. with the PED file markers reordered for physical position,
  2. with excluded SNPs (negative values in the MAP file) excluded from the new PED file
  3. possibly excluding other SNPs based on filters such as genotyping rate
  4. possibly recoding the SNPs to a 1/2 coding
  5. possibly recoding the SNPs between letters and numbers (A,C,G,T / 1,2,3,4)
  6. possibly transposing the genotype file (SNPs as rows)
  7. possibly recoding the SNP to an additive and dominant pair of components
  8. possibly listing the data with each specific genotype as a distinct row
  9. possibly listing the data one genotype per row
  10. possibly listing only minor alleles
The basic option to generate a new dataset is the --recode option:
plink --file data --recode

which will output the allele labels as they appear in the original; also, the missing genotype code is preserved if this is different from 0. Also, if --output-missing-genotype is specified (which can be as well as --missing-genotype) then this value will be used instead (i.e. so that input and output files can have different missing codes; this also applies to the phenotype with --output-missing-phenotype and --missing-phenotype).

The --make-bed option does the same as --recode but creates binary files; these can also be filtered, etc, as described below.

In contrast,

plink --file data --recode12

will recode the alleles as 1 and 2 (and the missing genotype will always be 0).

Both these commands will create two new files
     plink.ped
     plink.map
(where, as usual, "plink" would be replaced by any specified --out {filename} ).

Unless manually specified, for all these options, the usual filters for missingness and allele frequency will be set so as not to exclude any SNPs or individuals. By explicitly including an option, e.g. --maf 0.05 on the command line, this behaviour is overriden (see this page).

By default, any --recode option, and also --make-bed will preserve all genotypes exactly as they are. To set to missing Mendel errors or heterozygous haploid calls, use the options --set-me-missing and --set-hh-missing respectively. For the former, you will also need to specify --me 1 1 (i.e. to invole an evalation of Mendel errors, which does not occur by default, by not excluding any individuals or SNPs based on the results, i.e. if you only want to zero-out certain genotypes).

To recode SNP alleles from A,C,G,T to 1,2,3,4 or vice versa, use --allele1234 (to go from letters to numbers) and --alleleACGT (to go from numbers to letters). These flags should be used in conjunction with a data generation command (e.g. --make-bed), or any other analysis or summary statistic option. Alleles other than A,C,G,T or 1,2,3,4 will be left unchanged.

It is sometimes useful to have a PED file that is tab-delimited, except that between alleles of the same genotype a space instead of a tab is used. A file formatted in this way can load into Excel, for example, as a tab-delimited file, but with one genotype per column instead of one allele per column. Use the option --tab as well as --recode or --recode12 to achieve this effect.

To make a new file in which non-founders without both parents also in the same fileset are recoded as founders (i.e. pat and mat codes set both to 0), add the --make-founders flag.
Transposed genotype files
When using either --recode or --recode12, you can obtain a transposed text genotype file by adding the --transpose option. This generates two files:
     plink.tped
     plink.fam
The first contains the genotype data, with SNPs as rows and individuals as columns, for example: if the original file was
     1 1 0 0 1  1  1 1  G G
     1 2 0 0 2  1  0 0  A G
     1 3 0 0 1  1  1 1  A G
     1 4 0 0 2  1  2 1  A A
then this would generate
     1 snp1 0 10001  1 1  0 0  1 1  2 1
     1 snp2 0 20001  G G  G A  G A  A A
The first four columns are from the MAP file (chromosome, SNP ID, genetic position, physical position), followed by the genotype data. The plink.fam gives the ID, sex and phenotype information for each individual. The order of individuals in this file is the same as the order across the columns of the TPED file. The FAM file is just the first six columns of the PED file (or literally the same FAM file if the input where a binary fileset).
Additive and dominance components
The following format is often useful if one wants to use a standard, non-genetic statistical package to analyse the data, as here genotypes are coded as a single allele dosage number. To create a file with SNP genotypes recoded in terms of additive and dominant components, use the option:
plink --file data --recodeAD

which, assuming C is the minor allele, will recode genotypes as follows:
     SNP       SNP_A ,  SNP_HET
     ---       -----    -----
     A A   ->    0   ,   0
     A C   ->    1   ,   1
     C C   ->    2   ,   0
     0 0   ->   NA   ,  NA
In otherwords, the default for the additive recoding is to count the number of minor alleles per person. The --recodeAD option produces both an additive and dominance coding: use --recodeA instead to skip the SNP_HET coding.

The --recodeAD option saves the data to a single file
     plink.raw
which has a header row indicating the SNP names (with _A and _HET appended to the SNP names to represent additive and dominant components, respectively).

For example, consider the following PED file, which has two SNPs:
     1 1 0 0 1  1  1 1  G G
     1 2 0 0 2  1  0 0  A G
     1 3 0 0 1  1  1 1  A G
     1 4 0 0 2  1  2 1  A A
Using the --recodeAD option generates the file plink-recode.raw:
     FID IID PAT MAT SEX PHENOTYPE snp1_2 snp1_HET snp2_G snp2_HET
     1 1 0 0 1 1  0  0   2 0
     1 2 0 0 2 1  NA NA  1 1
     1 3 0 0 1 1  0  0   1 1
     1 4 0 0 2 1  1  1   0 0
The column labels reflect the snp name (e.g. snp1) with the name of the minor allele appended (i.e. snp1_2 in the first instance, as 2 is the minor allele) for the additive component. The dominant component ( a dummy variable reflecting heterozygote state) is coded with the _HET suffix.

This file can be easily loaded into R: for example:
     d <- read.table("plink.raw",header=T)
For example, for the first SNP, the individuals are coded 1/1, 0/0, 1/1 and 2/1. The additive count of the number of common (1) alleles is therefore: 2, NA, 2 and 1, which is reflected in the field snp1_2. The field snp1_HET is coded 1 for the fourth individual who is heterozygous -- this field can be used to model dominance effect of the allele.

The behavior of the --recodeA and --recodeAD commands can be changed with the --recode-allele command. This allows for the 0, 1, 2 count to reflect the number of a pre-specified allele type per SNP, rather than the number of the minor allele. This command takes as a single argument the name of a file that lists SNP name and allele to report, e.g. if the file recode.txt contained
     snp1   1
     snp2   A
then
plink --file data --recodeAD --recode-allele recode.txt

would now report in the LOG file
     Reading allele coding list from [ recode.txt ] 
     Read allele codes for 2 SNPs
and the plink.raw file would read
     FID IID PAT MAT SEX PHENOTYPE snp1_1 snp1_HET snp2_A snp2_HET
     1 1 0 0 1 1   2  0   0 0
     1 2 0 0 2 1   NA NA  1 1
     1 3 0 0 1 1   2  0   1 1
     1 4 0 0 2 1   1  1   2 0
If the SNP is monomorphic, by default the allele code out will be 0 and all individuals will have a count of 0 (or NA). If an allele is specified in --recode-allele that is not seen in the data, similarly all individuals will receive a 0 count (i.e. rather than an error being given).

NOTE For alleles that have exactly 0.50 minor allele frequency, as for the second SNP in the example above, then which allele is labelled as minor will depend on which was first encountered in the PED file.

Listing by minor allele count
The command
     --recode-rlist
will generate a files
 
     plink.rlist
     plink.fam
     plink.map
where the plink.rlist file format is
     SNP
     GENOTYPE (BOTH ALLELES)
     FID/IID PAIRS ...
For example, consider a particular SNP, rs2379981 has a minor allele (G) seen twice (in two heterozygotes) and two individuals with a missing genotpe; all other individuals are homozygous for the major allele. In this case, we would see two rows in the pink.rlist file:
     rs2379981 HET G A CH18612 NA18612  JA18998 NA18998
     rs2379981 NIL 0 0 JA18999 NA18999  JA19003 NA19003
indicating, for example, that individual FID/IID CH18612/NA18612 has a rare heterozygote.

This command could be used in conjunction with the --reference command and --freq to list all instances of rare non-reference alleles, e.g. from resequencing study data.
Listing by long-format (LGEN)
To output a file in the LGEN format, use the command
     --recode-lgen
which generates files
      plink.lgen
      plink.fam
      plink.map
that can be read with the --lfile command. The
     --with-reference
with generate a fourth file
     plink.ref
that can be read back in with the --reference command when using --lfile.
Listing by genotype
Another format that might sometimes be useful is the --list option which genetes a file
     plink.list
that is ordered one genotype per row, listing all family and individual IDs of people with that genotype. For example, if we have a file with two SNPs rs1001 and rs2002 (both on chromosome 1):
     A 1 0 0 1  2  A A  1 1
     B 2 0 0 1  2  A C  0 0
     C 3 0 0 1  1  A C  1 2
     D 4 0 0 1  1  C C  1 2
then then option
plink --file mydata --list

will generate the file plink.list
     1 rs1001 AA A 1
     1 rs1001 AC B 2 C 3
     1 rs1001 CC D 4
     1 rs1001 00
     1 rs2002 22
     1 rs2002 21 C 3 D 4
     1 rs2002 11 A 1
     1 rs2002 00 B 2
which has columns
     Chromosome
     SNP identifier
     Genotype
     Family ID, Individual ID for 1st person
     Family ID, Individual ID for 2nd person
     ...
     Family ID, Individual ID for final person
Obviously, different rows will have a different number of columns. Here, we see that individual A 1 has the A/A genotype for rs1001, etc. This option is often useful in conjunction with --snp, if you want an easy breakdown of which individuals have which genotypes.

Write SNP list files

To output just the list of SNPs that remain after all filtering, etc, use the --write-snplist command, e.g. to get a list of all high frequency, high genotyping-rate SNPs:
plink --bfile mydata --maf 0.05 --geno 0.05 --write-snplist

which generates a file
     plink.snplist
This file is simply a list of included SNP names, i.e. the same SNPs that a --recode or --make-bed statement would have produced in the corresponding MAP or BIM files.

Update SNP information

To automatically update either the genetic or physical positions for some or all SNPs in a dataset, use the --update-map command, which takes a single parameter of a filename, e.g.
plink --bfile mydata --update-map build36.txt --make-bed --out mydata2

where, for example, the file build36.txt contains new physical positions for SNPs, based on dbSNP126/build 36, in the simple format of SNP/position per line, e.g.
     rs100001  1000202
     rs100002  6252678
     rs100003  7635353
     ...
To change genetic position (3rd column in map file) add the flag --update-cm as well as --update-map. There is no way to change chromosome codes using this command. Normally, one would want to save the new file with the changed positions, as in the example above, although one could combine other commands instead (e.g. association testing, etc) although the updated positions would then be lost (i.e. the changes are not automatically saved).

The file with new SNP information does not need to feature all of the SNPs in the current dataset: SNPs not in this file will be left unchanged. If a SNP is listed more than once in the file, an error will be reported.

NOTE When updating the map positions, it is possible that the implied ordering of SNPs in the dataset might change. If this is the case, a message will be written to the LOG file. Although the positions are updated, the order is not changed internally: as SNPs might be out of order, it is important to correct this by saving and reloading the file. For example, the if the original contains
     ...
     rs10001   500000
     rs10002   520000
     rs10003   540000
     rs10004   560000
     ...
but we update rs10002 to position 580000, the data will be
     ...
     rs10001   500000
     rs10002   580000
     rs10003   540000
     rs10004   560000
     ...
Only after saving and reloading (e.g. --make-bed / --bfile ) will the file be in the correct order
     ...
     rs10001   500000
     rs10003   540000
     rs10004   560000
     rs10002   580000
     ...
This will only be an issue for commands which rely on relative SNP positions (e.g. --hap-window, --homozyg, etc). If the LOG file does not show a message that the order of SNPs has changed after using --update-map, one need not worry.

The name and chromosome code of a SNP can also be changed, by adding the modifiers --update-name or --update-chr, e.g.
./plink --bfile mydata --update-map rsID.lst --update-name --make-bed --out mydata2

or
./plink --bfile mydata --update-map chr-codes.txt --update-chr --make-bed --out mydata2

In both case, the format of the input file should be two columns per line, e.g.
    SNP_A-1919191   rs123456
    SNP_A-64646464  rs222222
    ...
or, for chromosome codes (use numeric values and codes X, Y, etc)
   rs123456     1
   rs987654     18
   rs678678     X
   ..
You cannot update more than one attribute at a time for SNPs.

Update allele information

To recode alleles, for example from A,B allele coding to A,C,G,T coding, use the command --update-alleles, for example
./plink --bfile mydata --update-alleles mylist.txt --make-bed --out newfile

where the file mylist.txt contains five columns per row listing,
     SNP identifier
     Old allele code for one allele
     Old allele code for other allele
     New allele code for first allele
     New allele code for other allele
For example,
    rs10001  A B   G T
    rs10002  A B   A C
    ...
will change allele A to G and allele B to T for rs10001, etc.

Force a specific reference allele

It is possible to manually specify which allele is the A1 allele and which is A2. By default, the minor allele is assigned to be A1. All odds ratios, etc, are calculated with respect to the A1 allele (i.e. an odds ratio greater than 1 implies that the A1 allele increases risk).

To set a particular allele as A1, which might not be the minor allele, use the command --reference-allele, which can be used with any other analysis or data generation command, e.g.
./plink --bfile mydata --reference-allele mylist.txt --assoc

where the file mylist.txt contains a list of SNP IDs and the allele to be set as A1, e.g.
     rs10001 A
     rs10002 T
     rs10003 T
     ...
This command can make comparing results across studies easier, so that odds ratios reported can be made to be in the same direction as the other study, for example.

Update individual information

Rather than try to manually edit PED or FAM files (which is not advised), use these functions to change ID codes, sex and parental information for individuals in a fileset. The command
plink --bfile mydata --update-ids recoded.txt --make-bed --out mydata2

changes ID codes for individuals specified in recoded.txt, which should be in the format of four columnds per row: old FID, old IID, new FID, new IID, e.g.
     FA 1001        F0001   I0001
     FA 1002.dup    F0002   I0002
     ...
will, for example find the person FA/1001 and change their FID/IID values to F0001/I0001. Not all people need be listed in the file (they will not be changed; the order of the file need not match the original dataset.

Two simular commands (but that cannot be run at the same time as --update-ids) are
--update-sex myfile1.txt

that expects 3 columns per row:
     FID
     IID
     SEX    Coded 1/2/0 for M/F/missing
and
--update-parents myfile2.txt

that expects 4 columns per row:
     FID
     IID
     PAT   New paternal IID code
     MAT   New maternal IID code
PLINK does not check see whether the new parents actually exist in the current file.

With all of these commands, you need to issue a data output command (--make-bed, --recode, etc) for the changes to be preserved.

Write covariate files

If a covariate file is specified along with any of the above --recode options or with --make-bed, then that covariate file will also be written, as plink.cov by default. This option is useful if the covariate file has a different number of individuals, or is ordered differently, to produce a set of covariate values that line up more easily with the newly-created genotype and phenotype files.
plink --file data --covar myfile.txt --recode

creates also plink.cov. If you want just to create a revised version of the covariate file, but without creating a new set of genotype files, then use the --write-covar option. This can be used in conjunction with filters, etc, to output, for example, only covariates for high-genotyping (99%) cases, as in this example:
plink --file data --write-covar myfile.txt --filter-cases --mind 0.01

will output just the relevant lines of myfile.txt to plink.cov, sorted to match the order of data.ped.

To also include phenotype information in the plink.cov file add the flag --with-phenotype. This can be useful, for example, when used in conjunction with --recodeA to generate the files needed to replicate an analysis in R (e.g. extracting the appropriate genotype data, and applying filters, etc).

To recode a categorical variable to a set of binary dummy variables, add the command
     --dummy-coding
for example
./plink --bfile mydate --covar cdata.raw --write-covar --dummy-coding

If the original covariate had two fields, a categorical variable with 8 levels (coded 0 to 7, although it could have any numeric coding, e.g. 100, 150, 200, 250, etc), and a second variable that was continuous, e.g.
     A8504 1   5  0.606218
     A8008 1   1  0.442154
     A8542 1   7  0.388042   
     A8022 1   2  0.286125
     A8024 1   3  0.903004
     A8026 1   4  0.790778  
     A8524 1   -9 0.713952
     A8556 1   0  0.814292
     A8562 1   1  0.803336
     ...
then the command above will create mynewfile.cov, with added header row, with the fields:
     FID       Family ID
     IID       Individual ID
     COV1_2    Dummy variable for first covariate, coded  1/0 for 2/other
     COV1_3    Dummy variable for first covariate, coded  1/0 for 3/other
     COV1_4    etc
     COV1_5 
     COV1_6 
     COV1_7 
     COV1_0 
     COV2     Unchanged continuous covariate
Thus mynewfile.cov is as follows (spaces added for clarity):
     FID IID   COV1_2 COV1_3 COV1_4 COV1_5 COV1_6 COV1_7 COV1_0 COV2 
     A8504 1   0 0 0 1 0 0 0          0.606218 
     A8008 1   0 0 0 0 0 0 0          0.442154 
     A8542 1   0 0 0 0 0 1 0          0.388042 
     A8022 1   1 0 0 0 0 0 0          0.286125 
     A8024 1   0 1 0 0 0 0 0          0.903004 
     A8026 1   0 0 1 0 0 0 0          0.790778 
     A8524 1   -9 -9 -9 -9 -9 -9 -9   0.713952 
     A8556 1   0 0 0 0 0 0 1          0.814292 
     A8562 1   0 0 0 0 0 0 0          0.803336 
That is, for a variable with K categories, K-1 new dummy variables are created. This new file can be used with --linear and --logistic, and a coefficient for each level would now be estimated for the first covariate (otherwise PLINK would have incorrectly treated the first covariate as an ordinal/ratio measure). For covariate Y, each new dummy variable for level X is named Y_X, e.g. COV1_2, etc.

Note that one level is automatically excluded (1 in this case, i.e. there is no COV1_1), which implicitly makes 1 the reference category in subsequent analysis. If PLINK detects more than 50 levels, it assumes the variable is not categorical (i.e. like COV2) and so leaves it unchanged. The command can operate on multiple covariates in a single file at the same time. Note that missing values are correctly handled (i.e. left as missing).

NOTE Note that, unlike cluster files (see below) PLINK cannot handle any string information in covariate files.

Write cluster files

Similar to --write-covar, the --write-cluster will output the single selected cluster from the file specified by --within. Unlike covariate files, this allows string labels to be used.
plink --bfile mydata --within clst.dat --write-cluster --out mynewfile

which writes a file
     mynewfile.clst
Use --mwithin to select which of multiple clusters is selected. The --dummy-coding can not currently be used with --write-cluster however.

Flip DNA strand for SNPs

This command will read the list of SNPs in the file list.txt and flip the strand for these SNPs, then save a new PED or BED fileset (i.e. by using either the --recode or --make-bed commands):
plink --file data --flip list.txt --recode

The list.txt should just be a simple list of SNP IDs, one SNP per line.

Flipping strand means changing alleles
   A -> T
   C -> G
   G -> C
   T -> A
so, for example, a A/C SNP will become a T/G; alternatively, a A/T SNP will become a T/A SNP (i.e. in this case, the labels remain the same, but whether the minor allele is A or T will still depend on strand).

To flip strand for just a subset of the sample (e.g. if two samples have already been merged, and subsequently a strand issue has been identified for one of those samples) use the option --flip-subset, for example
plink --file data --flip list.txt --flip-subset mylist.txt --recode

where mylist.txt is a text file containing the individuals (family ID, individual ID) to be flipped.

HINT When merging two datasets, it is clearly very important that the two sets of SNPs are concordant in terms of positive or negative strand. Whereas some mismatches will be easy to spot as more than two alleles will be observed in the merged dataset, other instances will not be so easy to spot, i.e. for A/T and C/G SNPs.

Using LD to identify incorrect strand assignment in a subset of the sample

If cases and controls have been genotyped separately and then the data merged, it is always possible that strand has been incorrectly or incompletely assigned to each SNP, meaning that the merged data may contain a number of SNPs for which the allele coding differs between cases and controls (or between any other grouping, such as collection site, etc).

If the two mis-matched groups correspond to cases and controls exactly, then rare SNPs will show a very strong association with disease (e.g. 5% MAF in cases, 95% in controls) and be easy to spot as potential problems. More common SNPs could show intermediate levels of association that might be easier to confuse with a real signal.

A simple approach to detect some proportion of such SNPs uses differential patterns of LD in cases versus controls: the command --flip-scan will query each SNP, and calculate the signed correlation between it and a set of nearby SNPs in cases and controls separately (of course, with the --pheno command, case and control status can be set to represent any binary split of the sample).

For each index SNP, PLINK identifies other SNPs in which the absolute value of the genotypic correlation is above some threshold. For these SNP pairs, it counts the number of times the signed correlation is different in sign between cases and controls (a negative LD pair) versus the same (a positive LD pair). For example, the command
plink --bfile mydata --flip-scan

produces the output file
     plink.flipscan
with the fields
     CHR     Chromosome
     SNP     SNP identifier for index SNP
     BP      Base-pair position
     A1      Minor allele code
     A2      Major allele code
     F       Allele frequency (A1 allele)
     POS     Number of positive LD matches
     R_POS   Average correlation of these 
     NEG     Number of negative LD matches
     R_NEG   Average correlation of these
     NEGSNPS The SNPs showing negative correlation
For example, the majority of this file should show SNPs have a NEG value of 0; the value of POS will be zero or greater, depending on the extent of LD. For example:
    CHR         SNP       BP  A1  A2      F  POS  R_POS  NEG  R_NEG  NEGSNPS
      1   rs9439462  1452629   T   C      0    0     NA    0     NA (NONE)
      1   rs1987191  1457348   C   T      0    0     NA    0     NA (NONE)
      1   rs3766180  1468016   C   T  0.285    2  0.893    0     NA (NONE)
However, occasionally one might observe different patterns of results. Of particular interest is when one SNP shows a large number of NEG SNPs. For example, here we show rs2240344 and nearby SNPs, all of which have at least one NEG SNP (lines truncated)

   CHR          SNP         BP  A1  A2       F  POS    R_POS  NEG    R_NEG  NEGSNPS
    14   rs12434442   72158039   T   C   0.249    5    0.515    1     0.46  rs2240344
    14    rs4899437   72190986   G   C   0.394    5    0.802    1    0.987  rs2240344
    14    rs2803980   72196284   G   A    0.41    5    0.808    1     0.95  rs2240344
    14    rs2240344   72197893   C   G   0.489    0       NA    7    0.807  rs12434442|rs4899437|...
    14    rs2286068   72198107   C   T   0.407    7    0.741    1    0.962  rs2240344
    14    rs7160830   72209491   T   C   0.414    6    0.801    1    0.922  rs2240344
    14   rs10129954   72220454   T   C   0.413    6    0.729    1     0.73  rs2240344
    14    rs7140455   72240734   T   C   0.469    4     0.72    1     0.64  rs2240344
This pattern of results quite clearly points to rs2240344 as being the odd man out: for 7 other SNPs, there is strong LD (r above 0.5) in either cases or controls, but with a SNP-SNP correlation in the other phenotype class that has the opposite direction. In contrast, there is not a single SNP for which both cases and controls have a consistent pattern of LD. For the nearby SNPs, all of which have only 1 NEG SNP, it is with rs2240344. So, in this particular case, it would suggest that stand is flipped in either cases or controls.

To display the specific sets of correlations in cases and controls for each SNP, add the option
     --flip-scan-verbose
which generates a file
     plink.flipscan.verbose
which lists for any SNP with at least one NEG pair of LD values, the correlations between the index SNP and the other flanking SNPs, showing the correlation in cases (R_A) and controls (R_U):
 CHR_INDX   SNP_INDX     BP_INDX  A1_INDX    SNP_PAIR    BP_PAIR  A1_PAIR     R_A      R_U
     14    rs2240344    72197893     C     rs12434442   72158039     T     -0.504    0.416
     14    rs2240344    72197893     C      rs4899437   72190986     G     -0.99     0.983
     14    rs2240344    72197893     C      rs2803980   72196284     G     -0.969    0.931
     14    rs2240344    72197893     C      rs2286068   72198107     C     -0.971    0.952
     14    rs2240344    72197893     C      rs7160830   72209491     T     -0.935    0.91
     14    rs2240344    72197893     C     rs10129954   72220454     T     -0.782    0.679
     14    rs2240344    72197893     C      rs7140455   72240734     T     -0.671    0.609
Here we see a clear pattern in which the correlation is similar between cases and controls in magnitude but has the opposite direction, strongly suggestive of a strand flip problem for this C/G SNP. In this case, the allele frequency turns out to be quite different between cases and controls (60% versus 40%) but the LD approach would have clearly detected this particular SNP being flipped in either cases or controls even if the true allele frequency were exactly 50%. This latter class of SNP would not cause problems of spurious association in single SNP analysis, but it could cause severe problems in haplotype and imputation analysis.

Naturally, if a SNP does not show strong LD with nearby SNPs, then this approach will not be able to resolve strand issues. Also, if more than one SNP in a region shows strand flips, or if there is a higher level of mis-coding alleles in general, then this approach may indicate that there are problems (many NEG scores above 0) but it might be less clear how to remedy them.

To know which to resolve (cases or controls) one would need to look at the frequency in other panels, or even the correlations, e.g. in HapMap. Ideally, one would only need to do this for a small number of SNPs if any. The --flip and --flip-subset commands described above can then be used to flip the appropriate genotypes.

Finally, the default threshold for counting can be changed by the following command:
     --flip-scan-threshold 0.8
The default is set at 0.5 (i.e. the pair needs to have a correlation of 0.5 or greater in either cases or controls). The number of flanking SNPs with are considered for each index SNP can be modified with the commands
     --ld-window 10
to set the number of SNPs considered upstream and downstream; the maximum physical distance away from the index SNP (1Mb by default) is specified in kb with the command:
     --ld-window-kb 500

Merge two filesets

To merge two PED/MAP files:
plink --file data1 --merge data2.ped data2.map --recode --out merge

The --merge option must be followed by 2 arguments: the name of the second PED file and the name of the second MAP file. A --recode (or --make-bed, etc) option is necessary to output the newly merged file; in this case, --out option will create the files merge-recode.ped and merge-recode.map.

The --merge option can also be used with binary PED files, either as input or output, but not as the second file: i.e.
plink --bfile data1 --merge data2.ped data2.map --make-bed --out merge

will create merge.bed, merge.fam and merge.bim, as the --make-bed option was used instead of the --recode option. Likewise, the data1.* files point to a binary PED file set.

If the second fileset (data2.*) were in binary format, then you must use --bmerge instead of --merge
plink --bfile data1 --bmerge data2.bed data2.bim data2.fam --make-bed --out merge

which takes 3 parameters (the names of the BED, BIM and FAM files, in that order).

The two filesets can either overlap completely, partially, or not at all both in terms of markers and individuals. Imputed genotypes will be set to missing (i.e. if SNP_B is not measured in the first file, but it is in the second, then any individuals in the first file who are not also present in the second file will be set to missing for SNP_B.

By default, any existing genotype data (i.e. in data1.ped) will not be over-written by data in the second file (data2.ped). By specifying a --merge-mode this default behavior can be changed. The modes are:
     1    Consensus call (default)
     2    Only overwrite calls which are missing in original PED file
     3    Only overwrite calls which are not missing in new PED file
     4    Never overwrite
     5    Always overwrite mode
     6    Report all mismatching calls (diff mode -- do not merge)
     7    Report mismatching non-missing calls (diff mode -- do not merge)
The default (mode 1) behaviour is to call the merged genotype as missing if the original and new files contain different, non-missing calls; otherwise: i.e.
                                Merge mode
    data1.ped ,  data2.ped  ->  1    2    3    4    5    
    ---------    ---------      -----------------------
     0/0      ,   0/0       ->  0/0  0/0  0/0  0/0  0/0
     0/0      ,   A/A       ->  A/A  A/A  A/A  0/0  A/A
     A/A      ,   0/0       ->  A/A  A/A  A/A  A/A  0/0
     A/A      ,   A/T       ->  0/0  A/A  A/T  A/A  A/T
Modes 6 and 7 effectively provide a means for comparing two PED files -- no merging is performed in these cases; rather, a list of mismatching SNPs is written to the file
     plink.diff
They should also report the concordance rate in the LOG file, based on all SNPs that feature in both sets.

A warning will be given if the chromosome and/or physical position differ between the two MAP files.

NOTE Alleles must be exactly coded to match: that is, PLINK will not assume that a {1,2,3,4} SNP coding maps onto a {A,C,G,T} coding. You can use the --allele1234 and --alleleACGT commands prior to merging to convert datasets and then merge these consistently coded files (you cannot convert and merge on the fly, i.e. simply do putting --allele1234 on the command line along with --merge will not work: you need to use --allele1234 and --make-bed first).

Merge multiple filesets

To merge more than two standard and/or binary filesets, it is often more convenient to specify a single file that contains a list of PED/MAP and/or BED/BIM/FAM files and use the --merge-list option. Consider, for an extreme example, the case where each fileset contains only a single SNP, and that there are thousands of these files -- this option would help build a single fileset, in this case.

For example, consider we had 4 PED/MAP filesets (labelled fA.* through fD.*) and 4 binary filesets, labelled fE.* through fH.*). Then using the command
plink --file fA --merge-list allfiles.txt --make-bed --out mynewdata

would create the binary fileset
     mynewdata.bed
     mynewdata.bim
     mynewdata.fam
(alternatively, the --recode option could have been used instead of --make-bed to generate a standard ASCII PED/MAP fileset). In this case, the file allfiles.txt was a list of the to-be-merged files, one set per row:
     fB.ped fB.map
     fC.ped fC.map
     fD.ped fD.map
     fE.bed fE.bim fE.fam
     fF.bed fF.bim fF.fam
     fG.bed fG.bim fG.fam
     fH.bed fH.bim fH.fam

Important Each fileset must be on a line by itself: lines with two files are interpreted as PED/MAP filesets; lines with three files are interpreted as binary BED/BIM/FAM filesets. The files on a line must always be in this order (PED then MAP; BED then BIM then FAM)

Note In this case the first of the 8 files must be the starting file, i.e. associated with --file on the command line; this file only contains the 8-1 remaining files therefore. The final mynewdata.* files will contain information from all 8 files.

The --merge-mode option can also be used with the --merge-list option, as described above: however, it is not possible to specify the "diff" features (i.e. modes 6 and 7).

Extract a subset of SNPs: command line options

There are multiple ways to extract just specific SNPs for analysis; this section describes options that use the command-line directly; the next section describes other methods that read a file containing the information.
Based on a single chromosome (--chr)
To analyse only a specific chromosome use
plink --file data --chr 6

Based on a range of SNPs (--from and --to)
To select a specific range of markers (that must all fall on the same chromosome) use, for example:
plink --bfile mydata --from rs273744 --to rs89883

Based on single SNP (and window) (--snp and --window)
Alternatively, you can specify a single SNP and, optionally, also ask for all SNPs in the surrounding region, with the --window option:
plink --bfile mydata --snp rs652423 --window 20

which extracts only SNPs within +/- 20kb of rs652423.
Based on multiple SNPs and ranges (--snps)
Alternatively, the newer --snps command is more flexible but slower than the previously described --snp and --from/--to commands. The --snps command will accept a comma-delimited list of SNPs, including ranges based on physical position. For example,
plink --bfile mydata --snps rs273744-rs89883,rs12345-rs67890,rs999,rs222

selects the same range as above (rs273744 to rs89883) but also the separate range rs273744 to rs89883 as well as the two individual SNPs rs999 and rs222. Note that SNPs need not be on the same chromosome; also, a range can span multiple chromosomes (the range is defined based on chromosome code order in that case, as well as physical position, i.e. a range from a SNP on chromosome 4 to one on chromosome 6 includes all SNPs on chromosome 5). No spaces are allowed between SNP names or ranges, i.e. it is
     --snps rs1111-rs2222,rs3333,rs4444
and not
     --snps rs1111 - rs2222, rs3333 ,rs4444
Hint As mentioned above, unlike other methods mentioned above, --snps will load in all the data before extracting what it needs, whereas --snp only loads in what it needs, as so is a much faster way to extract a region from a very large dataset: as a result, if you really do want only a single SNP or a single range, use --snp (with --window) or some variant of the from/--to commands.
Based on physical position (--from-kb, etc)
One can also select regions based on a window defined in terms of physical distance rather than SNP ID, using the command: e.g.
plink --bfile mydata --chr 2 --from-kb 5000 --to-kb 10000

to select all SNPs within this 5000kb region on chromosome 2 (when using --from-kb and --to-kb you always need to specify the chromosome with the --chr option).

HINT Two alternate forms of the --from-kb command are --from-bp and --from-mb that take a parameter in terms of base-pair position or megabase position, instead of kilobase (to be used with the corresponding --to-bp and --to-mb options).
Based on a random sampling (--thin)
To keep only a random 20% of SNPs, for example, add the flag
     --thin 0.2

 

All the above options can be used either with standard pedigree files (i.e. using --ped or --file) or with binary format pedigree (BED) files (i.e. using --bfile). One must combine this option with the desired analytic (e.g. --assoc), summary statistic (e.g. --freq) or data-generation (e.g. --make-bed) option.

Extract a subset of SNPs: file-list options

To extract only a subset of SNPs, it is possible to specify a list of required SNPs and make a new file, or perform an analysis on this subset, by using the command
plink --file data --extract mysnps.txt

where the file is just a list of SNPs, one per line, e.g.
     snp005
     snp008
     snp101
Alternatively, you can use the command --range to modify the behavior of --extract and --exclude. If the --range flag is added, then instead of a list of SNPs, PLINK will expect a list of chromosomal ranges to be given instead, one per line.
plink --file data --extract myrange.txt --range

All SNPs within that range will then be excluded or extracted. The format of myrange.txt should be, one range per line, whitespace-separated:
     CHR     Chromosome code (1-22, X, Y, XY, MT, 0)
     BP1     Start of range, physical position in base units
     BP2     End of range, as above
     LABEL   Name of range/gene
For example,
     2 30000000 35000000  R1
     2 60000000 62000000  R2
     X 10000000 20000000  R3
would extract/exclude all SNPs in these three regions (5Mb and 2Mb on chromosome 2 and 10Mb on chromosome X).

Based on an attribute file (--attrib)
See below
Based on a set file (--gene)
Finally, if a SET file is also specified, you can use the --gene option to extract all SNPs in that gene/region. For example, if the SET file genes.set contains two genes:
     GENE1
     rs123456
     rs10912
     rs66222
     END

     GENE2
     rs929292
     rs288222
     rs110191
     END
then
plink --file mydata --set genes.set --gene GENE2 --recode

would, for example, create a new dataset with only the 3 SNPs in GENE2. One must combine these options with the desired analytic (e.g. --assoc), summary statistic (e.g. --freq) or data-generation (e.g. --make-bed) option.

Remove a subset of SNPs

To re-write the PED/MAP files, but with certain SNPs excluded, use the option
plink --file data --exclude mysnps.txt

where the file mysnps.txt is, as for the --extract command, just a list of SNPs, one per line. As described above, the --range command can modify the behaviour of --exclude in the same manner as for --extract.

One must combine this option with the desired analytic (e.g. --assoc), summary statistic (e.g. --freq) or data-generation (e.g. --make-bed) option.

NOTE Another way of removing SNPs is to make the physical position negative in the MAP file (this can not be done for binary filesets (e.g. the *.bim file).

Make missing a specific set of genotypes

To blank out a specific set of genotypes, use the following commands, e.g.
	--zero-cluster test.zero  --within test.clst
in conjunction with other data analysis, file generation or summary statistic commands, where the file test.zero is a list of SNPs and clusters, and test.clust is a standard cluster file.

If the original PED file is
     1  1 0 0 1 1   A A  C C  A A 
     2  1 0 0 1 1   C C  A A  C C 
     3  1 0 0 1 1   A C  A A  A C 
     4  1 0 0 1 1   A A  C C  A A 
     5  1 0 0 1 1   C C  A A  C C 
     6  1 0 0 1 1   A C  A A  A C 
     1b 1 0 0 1 1   A A  C C  A A 
     2b 1 0 0 1 1   C C  A A  C C 
     3b 1 0 0 1 1   A C  A A  A C 
     4b 1 0 0 1 1   A A  C C  A A 
     5b 1 0 0 1 1   C C  A A  C C 
     6b 1 0 0 1 1   A C  A A  A C 
and the MAP file is
     1 snp1 0 1000
     1 snp2 0 2000
     1 snp3 0 3000
and the list of SNPs/clusters to zero out in test.zero is
     snp2   C1
     snp3   C1
     snp1   C2
and the cluster file test.clst is
     1b 1 C1
     2b 1 C1
     3b 1 C1
     4b 1 C1
     5b 1 C1
     6b 1 C1
     2  1 C2
     3  1 C2
then the command
plink --file test --zero-cluster test.zero --within test.clst --recode

results in a new PED file, plink.ped,
     1  1 0 0 1  1  A A C C A A
     2  1 0 0 1  1  0 0 A A C C
     3  1 0 0 1  1  0 0 A A A C
     4  1 0 0 1  1  A A C C A A
     5  1 0 0 1  1  C C A A C C
     6  1 0 0 1  1  A C A A A C
     1b 1 0 0 1  1  A A 0 0 0 0
     2b 1 0 0 1  1  C C 0 0 0 0
     3b 1 0 0 1  1  A C 0 0 0 0
     4b 1 0 0 1  1  A A 0 0 0 0
     5b 1 0 0 1  1  C C 0 0 0 0
     6b 1 0 0 1  1  A C 0 0 0 0
i.e. with the appropriate genotypes zeroed out.

HINT See the section on handling obligatory missing genotype data, which can often be useful in this context.

Extract a subset of individuals

To keep only certain individuals in a file, use the option:
plink --file data --keep mylist.txt

where the file mylist.txt is, as for the --remove command, just a list of Family ID / Individual ID pairs, one set per line, i.e. one person per line. (fields can occur after the 2nd column but they will be ignored -- i.e. you could use a FAM file as the parameter of the --keep command, or have comments in the file. For example
   F101   1 
   F1001  2_B
   F3033  1_A  Drop this individual because of consent issues   
   F4442  22
would be fine.

One must combine this option with the desired analytic (e.g. --assoc), summary statistic (e.g. --freq) or data-generation (e.g. --make-bed) option.

Remove a subset of individuals

To remove certain individuals from a file
plink --file data --remove mylist.txt

where the file mylist.txt is, as for the --keep command, just a list of Family ID / Individual ID pairs, one set per line, i.e. one person per line (although, as for --keep, fields after the 2nd column are allowed but they will be ignored).

One must combine this option with the desired analytic (e.g. --assoc), summary statistic (e.g. --freq) or data-generation (e.g. --make-bed) option.

Filter out a subset of individuals

Whereas the options to keep or remove individuals are based on files containing lists, it is also possible to specify a filter to include only certain individuals based on phenotype, sex or some other variable.

The basic form of the command is --filter which takes two arguments, a filename and a value to filter on, for example:
plink --file data --filter myfile.raw 1 --freq

implies a file myfile.raw exists which has a similar format to phenotype and cluster files: that is, the first two columns are family and individual IDs; the third column is expected to be a numeric value (although the file can have more than 3 columns), and only individuals who have a value of 1 for this would be included in any subsequent analysis or file generation procedure. e.g. if myfile.raw were
     F1  I1   2
     F2  I1   7
     F3  I1   1
     F3  I2   1
     F3  I3   3
then only two individuals (F3 I1 and F3 I2) would be included based on this filter for the calculation of allele frequencies. The filter can be any integer numeric value.

As with --pheno and --within, you can specify an offset to read the filter from a column other than the first after the obligatory ID columns. Use the --mfilter option for this. For example, if you have a binary fileset, and so the FAM file contains phenotype as the sixth column, then you could specify
plink --bfile data --filter data.fam 2 --mfilter 4

to select cases only; i.e. cases have the value 2, and this is the 4th variable in the file (i.e. the first two columns are ignored, as these are the ID columns).

Because filtering on cases or controls, or on sex, or on position within the family, will be common operations, there are some shortcut options that can be used instead of --filter. These are
     --filter-cases
     --filter-controls
     --filter-males
     --filter-females
     --filter-founders
     --filter-nonfounders
These flags can be used in any circumstances, e.g. to make a file of control founders,
plink --bfile data --filter-controls --filter-founders --make-bed --out newfile

or to analyse only males
plink --bfile data --assoc --filter-males

IMPORTANT Take care when using these with options to merge filesets: the merging occurs before these filters.

Attribute filters for markers and individuals

One can define an attribute file for SNPs (or for individuals, see below) that is simply a list of user-defined attributes for SNPs. For example, this might be a file
     snps.txt
which contains
     rs0001   exonic 
     rs0007   candidate
     rs0010   failed exonic 
     rs0012   nssnp
These codes can be whatever you like, as is appropriate for your study; a SNP can have multiple, white-space delimited attributes. Not all SNPs need appear in this file; SNPs not in the dataset are allowed to appear (they are just ignored); the order does not need to be the same. Each SNP should only be listed once however. A SNP can be listed by itself without any attributes (for example, to ensure it is not excluded when filtering to exclude SNPs with a certain attribute, see below).

To filter SNPs on these, use the command (combined with some other data generation or analysis option)
     --attrib snps.txt exonic
for example, to extract only exonic SNPs (rs0001 and rs0010 in this example, assuming they've been coded this way).

To exclude SNPs that match the attribute, preface the attribute with a minus sign on the command line, e.g.
     --attrib snps.txt -failed
to extract only non-failed SNPs. Finally, multiple filters can be combined in a comma-delimited list
     --attrib snps.txt exonic,-failed
would select exonic SNPs that did not fail. If a SNP does not feature in the attribute file, it will always be excluded.

NOTE Within match type, multiple matches are treated as logical ORs; between positive and negative matches as AND. For example, matching on A,B,-C,-D implies individuals with ( A or B ) and not ( C or D )

This approach works similarly for individuals, except the command is now --attrib-indiv, e.g.
     --attrib-indiv inddat.txt  sample1,fullinfo
and the attribute file starts with family ID and individual ID before listing any attributes, e.g.
     F1  1   sample2 
     F2  1   sample1 
     F3  1   sample2 fullinfo 
     ...

Create a SET file based on a list of ranges

Given a list of ranges in the following format (4 columns per row; no header file)
 
     Chromosome
     Start base-pair position 
     End base-pair position
     Set/range/gene name
then the command
plink --file mydata --make-set gene.list --write-set

will generate the file
     plink.set
in the standard set file format. The command --make-set-border takes a single integer argument, allowing for a certain kb window before and after the gene to be included, e.g. for 20kb upstream and downstream:
plink --file mydata --make-set gene.list --make-set-border 20 --write-set

HINT The --make-set command doesn't necessarily have to be used with --write-set. Rather, it can be used anywhere that --set can be used, to make sets on the fly. Similar, --set and --write-set can be combined, e.g. to create a new, filtered set file.

Options for --make-set
To collapse all ranges into a single set (i.e. to generate one set that corresponds to all SNPs in a gene, from a list of gene co-ordinates, for example), use
     --make-set-collapse-all SETNAME
along with --make-set, where SETNAME is any single word that you use to name to newly created set. To make a set file of all SNPs not in the specified ranges, add
     --make-set-complement-all SETNAME
Optionally, the range file can contain a fifth column, to specify groups of ranges. Sets can be constructed which collapse over these groups. That is, the input for --make-set is now
     Chromosome
     Start base-pair position 
     End base-pair position
     Set/range/gene name
     Group label
e.g.
  1   10001   20003  GENE1  PWAY-A
  8   80001   99995  GENE2  PWAY-A
  12  1001    10001  GENE3  PWAY-B
  5   110001 127362  GENE4  PWAY-B
  ...
Normally, the fifth column will just be ignored, unless the command
      --make-set-collapse-group
is added, which creates sets of SNPs that correspond to each group (i.e. PWAY-A, PWAY-B, etc, in this example) rather than each gene/region (i.e. GENE1, etc). The command
     --make-set-complement-group
works in a similar manner, except now forming sets of all SNPs not in the given group of ranges.

HINT See the resources page for pre-compiled RefSeq gene-lists that can be used here.

Tabulate set membership for all SNPs

It is possible to create a table that maps SNPs to sets, given a --set file has been specified, with the --set-table command, e.g.
./plink --bfile mydata --set mydata.set --set-table

which generates a file
     plink.set.table
which contains the fields
     SNP                SNP identifier
     CHR                Chromosome code
     BP                 Base-pair physical position
     First set name     Membership of first set
     Second set name    Membership of second set
     ...
For each row, a series of 0s and 1s indicate whether or not each SNP in the dataset is in a given SET. This format can be useful for subsequent analyses (i.e. it can be easily lined up with other result files, e.g. from --assoc).

SNP-based quality scores

PLINK supports quality scores for SNPs and, described in the next section, genotypes. These can be used to filter on user-defined thresholds. The command --qual-scores indicates the file containing the scores. Scores are assumed to be numbers between 0 and 1, a higher number representing better quality. The threshold at which SNPs are selected can be set with the command --qual-threshold. For example,
./plink --bfile mydata --qual-scores myscores.txt --qual-threshold 0.8 --make-bed --out qc-data

where myscores.txt is a text file of SNPs and scores, e.g.
     rs10001 0.87
     rs10002 0.46
     rs10003 1.00
     ...
will remove SNPs with scores less than 0.8. The additional flag --qual-max-threshold can be used to specify a maximum threshold also (i.e. to select low-quality SNPs only). Not all SNPs need be in the file (the SNP is left in, in this case; the order can be different, it can contain SNPs not in the data).

Genotype-based quality scores

Quality scores for each genotype, rather than each SNP, can also be applied to PLINK datasets, using the --qual-geno-scores command, e.g.
./plink --bfile mydata --qual-geno-scores gqual.txt --qual-geno-threshold 0.99 --assoc

(with a similar --qual-geno-max-threshold command as well).

The file containing the genotype quality scores should have the following format:
  Q FID IID SNPID score
e.g.
 Q fam1 ind1 rs10001 0.873
 Q fam1 ind1 rs10002 0.998
 ...
Not all genotypes need be in this file. Rather than have a very large file, one could only list genotype scores that are below some threshold, for example, assuming most genotypes are of very good quality. Genotypes not in the this file will be untouched. This format is designed to accept wildcards, as follows. Every item should start with a Q character, to allow PLINK to check the correctness of the file format. Consider this example file,
     Q  A 1  rs1234 0.986
     Q  B 1  rs1234 0.923
     Q  A 1  rs5678 0.323
     Q  B 1  rs5678 0.97
that lists two genotypes for people with FID/IID A/1 and B/1 for SNPs rs1234 and rs5678. If a score if below threshold, it is set to missing in the data. The order of this file is arbitrary; not all individuals/SNPs need appear.

PLINK accepts wildcards in this file, to allow for different data formats to be specified. With a person wild-card, PLINK expects all quality scores for that SNP, in order as in the FAM or PED file, e.g.
     Q * rs1234 0.986 0.923
     Q * rs5678 0.323 0.97
With a SNP wildcard, PLINK exects all SNPs for a given person:
     Q A 1 * 0.986 0.323
     Q B 1 * 0.923 0.97
All these formats can be mixed together in a single file. These can be combined (in which case, PLINK expects all individuals for the first SNP, all for the second SNP, etc)
     Q * * 0.986 0.923 0.323 0.97

WARNING This option is recently added in beta-stage of development. Currently, a wild card looks to the current data to get the list of individuals and SNPs to loop over. This could cause a problem if the file has been filtered, etc. The next release will include commands to specify the order of individuals and SNPs, e.g.
  --qual-people-list mysamples.lst
where mysamples.lst is a file with 2 columns (FID/IID), and
  --qual-geno-snp-list mysnp.lst
where mysnp.lst is list of SNPs. This way if somebody is in the quality score file but they have been removed from the actual genotype dataset (or added), then this can be handled properly without needing to change the whole quality score file.
 

This document last modified Saturday, 10-Oct-2009 11:28:11 EDT