plink的GWAS数据处理作业流程.docx

Data management
Generate binary fileset
--make-bed
--make-bed creates a new PLINK 1 binary fileset, after applying sample/variant filters and other operations below. For example,
plink --file text_fileset --maf 0.05 --make-bed --out binary_fileset
does the following:
1. Autogenerate binary_fileset-temporary.bed + .bim + .fam. (The MAF filter has not yet been applied at this stage. See the Order of operations page for more details.)
2. Read binary_fileset-temporary.bed + .bim + .fam. Calculate MAFs. Remove all variants with MAF < 0.05 from the current analysis.
3. Generate binary_fileset.bed + .bim + .fam. Any samples/variants removed from the current analysis are also not present in this fileset. (This is the --make-bed step.)
4. Delete binary_fileset-temporary.bed + .bim + .fam.
In contrast, the fileset left behind by --keep-autoconv is just the result of step 1.
--make-just-bim
--make-just-fam
--make-just-bim is a variant of --make-bed which only generates a .bim file, and --make-just-fam plays the same role for .fam files. Unlike most other PLINK commands, these do not require the main input to include a .bed file (though you won't have access to many filtering flags when using these in no-.bed mode).
Use these cautiously. It is very easy to desynchronize your binary genotype data and your .bim/.fam indexes if you use these commands improperly. If you have any doubt, stick with --make-bed.
Generate text fileset
--recode <01 | 12> <23 | A | A-transpose | AD | beagle | beagle-nomap | bimbam | bimbam-1chr | compound-genotypes | fastphase | fastphase-1chr | HV | HV-1chr | lgen | lgen-ref | list | oxford | rlist | structure | transpose | vcf | vcf-fid | vcf-iid> <tab | tabx | spacex | bgz | gen-gz> <include-alt> <omit-nonmale-y>
--recode-allele [filename]
--recode creates a new text fileset, after applying sample/variant filters and other operations. By default, the fileset includes a .ped and a .map file, readable with --file.
· The '12' modifier causes A1 (usually minor) alleles to be coded as '1' and A2 alleles to be coded as '2', while '01' maps A1→0 and A2→1. (PLINK forces you to combine '01' with --{output-}missing-genotypewhen this is necessary to prevent missing genotypes from becoming indistinguishable from A1 calls.)
· The '23' modifier causes a 23andMe-formatted file to be generated. This can only be used on a single sample's data (a one-line --keep file may come in handy here). There is currently no special handling of the XY pseudo-autosomal region.
· The 'AD' modifier causes an additive (0/1/2) + dominant (het = 1, otherwise 0) component file, suitable for loading from R, to be generated. 'A' is the same, except without the dominance component.
o By default, A1 alleles are counted; this can be customized with --recode-allele. --recode-allele's input file should have variant IDs in the first column and allele IDs in the second.
o By default, the header line for .raw files only names the counted alleles. To include the alternate allele codes as well, add the 'include-alt' modifier.
o Haploid additive components are 0/2-valued instead of 0/1-valued, to maintain a consistent scale on the X chromosome.
See also --R.
· The 'A-transpose' modifier causes a variant-major additive component file to be generated. This can also be used with --recode-allele.
· The 'beagle' modifier causes unphased per-autosome .dat and .map files, readable by BEAGLE 3.3 and earlier, to be generated, while 'beagle-nomap' generates a single .dat file (no chromosome splitting occurs in this case).
· The 'bimbam' modifier causes a BIMBAM-formatted fileset to be generated. If your input data only contains one chromosome, you can use 'bimbam-1chr' instead to write a two-column .pos.txt file.
· If all allele codes are single-character, you can use the 'compound-genotypes' modifier to omit the space between each pair of allele codes in a single genotype call when generating a .ped + .map fileset. You will need to use the --compound-genotypes flag to load this data in PLINK 1.07, but it's not needed for PLINK 1.9.
· The 'fastphase' modifier causes per-chromosome fastPHASE files to be generated. If your input data only contains one chromosome, you can use 'fastphase-1chr' instead to exclude the chromosome number from the file extension.
· The 'HV' modifier causes a Haploview-format .ped + .info fileset to be generated per chromosome. 'HV-1chr' is analogous to 'fastphase-1chr'.
· The 'lgen' modifier causes a long-format fileset, loadable with --lfile, to be generated. 'lgen-ref' is equivalent to PLINK 1.07 --recode-lgen --with-reference.
· The 'list' modifier causes a genotype-based list to be generated. This does not produce a .fam or .map file.
· The 'oxford' modifier causes a Oxford-format .gen + .sample fileset to be generated. If you also include the 'gen-gz' modifier, the .gen file is gzipped.
· The 'rlist' modifier causes a rare-genotype fileset to be generated (similar to --list's output, but with .fam and .map files, and without homozygous major genotypes).
· With the 'list' and 'rlist' formats, the 'omit-nonmale-y' modifier causes nonmale genotypes to be omitted on the Y chromosome.
· The 'structure' modifier causes a Structure-format file to be generated.
· The 'transpose' modifier causes a transposed text fileset, loadable with --tfile, to be generated.
· The 'vcf', 'vcf-fid', and 'vcf-iid' modifiers result in production of a VCFv4.2 file. 'vcf-fid' and 'vcf-iid' cause family IDs and within-family IDs respectively to be used for the sample IDs in the last header row, while 'vcf' merges both IDs and puts an underscore between them (in this case, a warning will be given if an ID already contains an underscore).
If the 'bgz' modifier is added, the VCF file is block-gzipped. (Gzipping of other --recode output files is not currently supported.)
The A2 allele is saved as the reference and normally flagged as not based on a real reference genome ('PR' INFO field value). When it is important for reference alleles to be correct, you'll usually also want to include --a2-allele and --real-ref-alleles in your command.
· The 'tab' modifier makes the output mostly tab-delimited instead of mostly space-delimited when the format permits both delimiters. 'tabx' and 'spacex' force all tabs and all spaces, respectively. (See this page for guidelines on swapping tabs/spaces in other contexts.)
For example,
plink --bfile binary_fileset --recode --out new_text_fileset
generates new_text_fileset.ped and new_text_fileset.map from the data inbinary_fileset.bed + .bim + .fam, while
plink --bfile binary_fileset --recode vcf-iid --out new_vcf
generates new_vcf.vcf from the same data, removing family IDs in the process.
Irregular output coding
--output-chr [MT code]
Normally, autosomal/sex/mitochondrial chromosome codes in PLINK output files are numeric, e.g. '23' for human X. --output-chr lets you specify a different coding scheme by providing the desired human mitochondrial code; supported options are '26' (default), 'M', 'MT', '0M', 'chr26', 'chrM', and 'chrMT'. (PLINK 1.9 correctly interprets all of these encodings in input files.)
--output-missing-genotype [char]
--output-missing-phenotype [string]
--output-missing-genotype allows you to change the character (normally the --missing-genotype value) used to represent missing genotypes in PLINK output files, while --output-missing-phenotype changes the string (normally the --missing-phenotype value) representing missing phenotypes.
Note that these flags do not affect --{b}merge/--merge-list or the autoconverters, since they generate files that may be reloaded during the same run. Add --make-bed if you want to change missing genotype/phenotype coding when performing those operations.
Set blocks of genotype calls to missing
--zero-cluster [filename]
If clusters have been defined, --zero-cluster takes a file with variant IDs in the first column and cluster IDs in the second, and sets all the corresponding genotype calls to missing. See the PLINK 1.07 documentation for an example.
This flag must now be used with --make-bed and no other output commands (since PLINK no longer keeps the entire genotype matrix in memory).
Heterozygous haploid errors
--set-hh-missing
Normally, heterozygous haploid and nonmale Y chromosome genotype calls are logged to plink.hh and treated as missing by all analysis commands, but left undisturbed by --make-bed and --recode (since, once gender and/or chromosome code errors have been fixed, the calls are often valid). If you actually want --make-bed/--recode to erase this information, use --set-hh-missing. (The scope of this flag is a bit wider than for PLINK 1.07, since commands like --list and --recode-rlist which previously did not respect --set-hh-missing have been consolidated under --recode.)
Note that the most common source of heterozygous haploid errors is imported data which doesn't follow PLINK's convention for representing the X chromosome pseudo-autosomal region. This should be addressed with --split-x below, not --set-hh-missing.
--set-mixed-mt-missing
Mitochondrial DNA is subject to heteroplasmy, so PLINK 1.9 permits 'heterozygous' genotypes and treats MT more like a diploid than a haploid chromosome. However, some analytical methods don't use mixed MT genotype calls, and instead assume that no 'heterozygous' MT calls exist. The --set-mixed-mt-missing flag can be used with --make-bed/--recode to export a dataset with mixed MT calls erased.
X chromosome pseudo-autosomal region
--split-x [last bp position of head] [first bp position of tail] <no-fail>
--split-x [build code] <no-fail>
--merge-x <no-fail>
PLINK prefers to represent the X chromosome's pseudo-autosomal region as a separate 'XY' chromosome (numeric code 25 in humans); this removes the need for special handling of male X heterozygous calls. However, this convention has not been widely adopted, and as a consequence, heterozygous haploid 'errors' are commonplace when PLINK 1.07 is used to handle X chromosome data. The new --split-x and --merge-x flags address this problem.
Given a dataset with no preexisting XY region, --split-x takes the base-pair position boundaries of the pseudo-autosomal region, and changes the chromosome codes of all variants in the region to XY. As (typo-resistant) shorthand, you can use one of the following build codes:
· 'b36'/'hg18': NCBI build 36/UCSC human genome 18, boundaries 2709521 and
· 'b37'/'hg19': GRCh37/UCSC human genome 19, boundaries 2699520 and
· 'b38'/'hg38': GRCh38/UCSC human genome 38, boundaries 2781479 and
By default, PLINK errors out if no variants would be affected by the split. This behavior may break data conversion scripts which are intended to work on e.g. VCF files regardless of whether or not they contain pseudo-autosomal region data; use the 'no-fail' modifier to force PLINK to always proceed in this case.
Conversely, in preparation for data export, --merge-x changes chromosome codes of all XY variants back to X (and 'no-fail' has the same effect). Both of these flags must be used with --make-bed and no other output commands.
Mendel errors
--set-me-missing
In combination with --make-bed, --set-me-missing scans the dataset for Mendel errors and sets implicated genotypes (as defined in the --mendel table) to missing.
· --mendel-duos causes samples with only one parent in the dataset to be checked, while --mendel-multigen causes (great-)ngrandparental data to be referenced when a parental genotype is missing.
· It is no longer necessary to combine this with e.g. '--me 1 1' to prevent the Mendel error scan from being skipped.
· Results may differ slightly from PLINK 1.07 when overlapping trios are present, since genotypes are no longer set to missing before scanning is complete.
Fill in missing calls
--fill-missing-a2
It can be useful to fill in all missing calls in a dataset, e.g. in preparation for using an algorithm which cannot handle them, or as a 'decompression' step when all variants not included in a fileset can be assumed to be homozygous reference matches and there are no explicit missing calls that still need to be preserved.
For the first scenario, a sophisticated imputation program such as BEAGLE or IMPUTE2 should normally be used, and --fill-missing-a2 would be an information-destroying operation bordering on malpractice. However, sometimes the accuracy of the filled-in calls isn't important for whatever reason, or you're dealing with the second scenario. In those cases you can use the --fill-missing-a2 flag (in combination with --make-bed and no other output commands) to simply replace all missing calls with homozygous A2 calls. When used in combination with --zero-cluster/--set-hh-missing/--set-me-missing, this always acts last.
You may want to combine this with --a2-allele below.
Update variant information
--set-missing-var-ids [template string]
--new-id-max-allele-len [n]
--missing-var-code [missing ID string]
Whole-exome and whole-genome sequencing results frequently contain variants which have not been assigned standard IDs. If you don't want to throw out all of that data, you'll usually want to assign them chromosome-and-position-based IDs.
--set-missing-var-ids provides one way to do this. The parameter taken by these flags is a special template string, with a '@' where the chromosome code should go, and a '#' where the base-pair position belongs. (Exactly one @ and one # must be present.) For example, given a .bim file starting with
chr1 . 0 10583 A G
chr1 . 0 886817 C T
chr1 . 0 886817 CATTTT C
chrMT . 0 64 T C
'--set-missing-var-ids @:#[b37]' would name the first variant 'chr1:10583[b37]', the second variant 'chr1:886817[b37]'... and then error out when naming the third variant, since it would be given the same name as the second variant. (Note that this position overlap is actually present in 1000 Genomes Project phase 1 data.)
To maintain unique IDs in this situation, you can include '$1' and '$2' in your template string as well; these refer to the first and second allele names in ASCII-sort order. So, if we're using a bash shell, we can try again with
--set-missing-var-ids @:#[b37]\$1,\$2
which would name the first variant 'chr1:10583[b37]A,G', the second variant 'chr1:886817[b37]C,T', the third variant 'chr1:886817[b37]C,CATTTT', and the fourth variant 'chrMT:64[b37]C,T'. Note the extra backslashes: they are necessary in bash because '$' is a reserved character there.
You may still get a small number of duplicate ID errors when using '$1' and '$2'. If indels are involved, it is likely that the ambiguity cannot be resolved by PLINK 1 at all, because it matters which allele is the reference allele1. Instead, you must e.g. use a shell script to manually name variants in your original VCF file; see this blog post by Giulio Genovese for a detailed discussion. We apologize for the inconvenience; PLINK 2.0 will extend --set-missing-var-ids