• 數(shù)據(jù)

• 1 二進(jìn)制文件

• 2. plink二進(jìn)制文件變?yōu)槲谋疚募╬ed和map）

• 3. plink將vcf轉(zhuǎn)化為plink文件

• 4. 提取樣本和SNP

• 4.1 提取樣本

• 4.2 提取SNP

• 5. plink和表型數(shù)據(jù)合并

• 6. 數(shù)據(jù)匯總

• 6.1 次等位基因頻率（maf）

• 6.2 缺失

• 6.3 哈溫檢測(cè)

• 6.4 雜合率檢測(cè)

• 7. 數(shù)據(jù)質(zhì)控qc

• 7.1 樣本質(zhì)控--mind

• 7.2 位點(diǎn)質(zhì)控--geno

• 7.3 maf質(zhì)控--maf

• 7.4 哈溫質(zhì)控 --hwe

數(shù)據(jù)

《統(tǒng)計(jì)遺傳學(xué)》中的章節(jié)介紹，有關(guān)代碼實(shí)操部分，單獨(dú)列出來(lái)，進(jìn)行展示。

我已經(jīng)下載整理好了，下載本書(shū)的電子版pdf+數(shù)據(jù)+代碼，鏈接：書(shū)籍及配套代碼領(lǐng)取--統(tǒng)計(jì)遺傳分析導(dǎo)論

1 二進(jìn)制文件

文件中包括二進(jìn)制的三個(gè)文件：

2. plink二進(jìn)制文件變?yōu)槲谋疚募╬ed和map）

命令

 plink --bfile hapmap-ceu --recode --out hapmap-ceu

• --bfile 是指定二進(jìn)制plink的前綴名稱(chēng)
• --recode是生成文本ped和map的二進(jìn)制文件
• --out是指定輸出的結(jié)果文件前綴名稱(chēng)

日志

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hapmap-ceu.log.
Options in effect:
  --bfile hapmap-ceu
  --out hapmap-ceu
  --recode

1031523 MB RAM detected; reserving 515761 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Warning: 59 het. haploid genotypes present (see hapmap-ceu.hh ); many commands
treat these as missing.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--recode ped to hapmap-ceu.ped + hapmap-ceu.map ... done.

上面的信息有：

• 共有2239392個(gè)SNP位點(diǎn)
• 共有60個(gè)個(gè)體，其中30個(gè)males，30個(gè)female
• 結(jié)果生成hapmap-ceu.ped和hapmap-ceu.map文件

3. plink將vcf轉(zhuǎn)化為plink文件

文件：ALL.chr21.vcf.gz

plink將vcf文件變?yōu)閜link的二進(jìn)制文件（bed和bim和fam）。

命令：

plink --vcf ALL.chr21.vcf.gz --make-bed --out test_vcf

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to test_vcf.log.
Options in effect:
  --make-bed
  --out test_vcf
  --vcf ALL.chr21.vcf.gz

1031523 MB RAM detected; reserving 515761 MB for main workspace.
--vcf: test_vcf-temporary.bed + test_vcf-temporary.bim + test_vcf-temporary.fam
written.
1105538 variants loaded from .bim file.
2504 people (0 males, 0 females, 2504 ambiguous) loaded from .fam.
Ambiguous sex IDs written to test_vcf.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 2504 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.999787.
1105538 variants and 2504 people pass filters and QC.
Note: No phenotypes present.
--make-bed to test_vcf.bed + test_vcf.bim + test_vcf.fam ... done.

4. 提取樣本和SNP

• 可以使用--keep 和 --remove 進(jìn)行樣本ID的提取或者刪除
• 可以使用--extract和 --exclude進(jìn)行SNP的提取或者刪除

樣本ID的示例：兩列FID和IID，沒(méi)有行頭：

SNP的ID的示例：一列SNP名稱(chēng)，沒(méi)有行頭：

4.1 提取樣本

代碼：

 plink --bfile hapmap-ceu --keep list.txt --make-bed --out selectedIndividuals

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to selectedIndividuals.log.
Options in effect:
  --bfile hapmap-ceu
  --keep list.txt
  --make-bed
  --out selectedIndividuals

1031523 MB RAM detected; reserving 515761 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
--keep: 60 people remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Warning: 59 het. haploid genotypes present (see selectedIndividuals.hh ); many
commands treat these as missing.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--make-bed to selectedIndividuals.bed + selectedIndividuals.bim +
selectedIndividuals.fam ... done.

4.2 提取SNP

代碼：

plink --bfile hapmap-ceu --extract snp_name.txt --make-bed --out selectedSNP

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to selectedSNP.log.
Options in effect:
  --bfile hapmap-ceu
  --extract snp_name.txt
  --make-bed
  --out selectedSNP

1031523 MB RAM detected; reserving 515761 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
--extract: 20 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.995833.
20 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--make-bed to selectedSNP.bed + selectedSNP.bim + selectedSNP.fam ... done.

也可以提取單個(gè)SNP：

 plink --file hapmap-ceu --snps rs9930506 --make-bed --out rs9930506sample

結(jié)果：

 rs9930506sample.bed + rs9930506sample.bim + rs9930506sample.fam

5. plink和表型數(shù)據(jù)合并

如果想要把表型數(shù)據(jù)和基因型數(shù)據(jù)合并，需要整理的表型格式：FID，IID，y三列。

FID IID BMI
0 HG00096 25.022827
0 HG00097 24.853638
0 HG00099 23.689295
0 HG00100 27.016203
0 HG00101 21.461624
0 HG00102 20.673635
0 HG00103 25.71508
0 HG00104 25.252243
0 HG00106 22.765049

看一下原始的fam文件：

0 HG00096 0 0 1 -9
0 HG00097 0 0 2 -9
0 HG00099 0 0 2 -9
0 HG00100 0 0 2 -9
0 HG00101 0 0 1 -9
0 HG00102 0 0 2 -9
0 HG00103 0 0 1 -9
0 HG00104 0 0 2 -9
0 HG00106 0 0 2 -9
0 HG00108 0 0 1 -9

合并命令：

 plink --bfile 1kg_EU_qc --pheno BMI_pheno.txt --make-bed --out 1kg_EU_BMI

生成文件：

1kg_EU_BMI.bed + 1kg_EU_BMI.bim + 1kg_EU_BMI.fam

看一下新的fam文件：

0 HG00096 0 0 1 25.0228
0 HG00097 0 0 2 24.8536
0 HG00099 0 0 2 23.6893
0 HG00100 0 0 2 27.0162
0 HG00101 0 0 1 21.4616
0 HG00102 0 0 2 20.6736
0 HG00103 0 0 1 25.7151
0 HG00104 0 0 2 25.2522
0 HG00106 0 0 2 22.765
0 HG00108 0 0 1 23.069

可以看到，已經(jīng)合并成功了。

6. 數(shù)據(jù)匯總

6.1 次等位基因頻率（maf）

查看基因頻率的統(tǒng)計(jì)結(jié)果，用--freq

命令：

plink --bfile hapmap-ceu --freq --out Allele_Frequency

日志：

PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to Allele_Frequency.log.
Options in effect:
  --bfile hapmap-ceu
  --freq
  --out Allele_Frequency

15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
--freq: Allele frequencies (founders only) written to Allele_Frequency.frq .
Warning: 59 het. haploid genotypes present (see Allele_Frequency.hh ); many
commands treat these as missing.

查看結(jié)果，結(jié)果文件是Allele_Frequency.frq

$ head Allele_Frequency.frq
 CHR          SNP   A1   A2          MAF  NCHROBS
   1   rs12565286    C    G       0.0678      118
   1   rs12138618    A    G      0.05833      120
   1    rs3094315    G    A       0.1552      116
   1    rs3131968    A    G        0.125      120
   1   rs12562034    A    G      0.09167      120
   1    rs2905035    A    G        0.125      120
   1   rs12124819    G    A       0.3417      120
   1    rs2980319    A    T        0.125      120
   1    rs4040617    G    A        0.125      120

如果對(duì)其進(jìn)行質(zhì)控，用--maf 0.01，會(huì)去掉maf小于0.01的位點(diǎn)。

6.2 缺失

缺失包括樣本缺失率統(tǒng)計(jì)和位點(diǎn)缺失率統(tǒng)計(jì)。

命令：

plink --bfile hapmap-ceu --missing --out missing_data

日志：

$ plink --bfile hapmap-ceu --missing --out missing_data
PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to missing_data.log.
Options in effect:
  --bfile hapmap-ceu
  --missing
  --out missing_data

15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
--missing: Sample missing data report written to missing_data.imiss, and
variant-based missing data report written to missing_data.lmiss.
Warning: 59 het. haploid genotypes present (see missing_data.hh ); many
commands treat these as missing.

查看樣本缺失的文件：

$ head missing_data.imiss
 FID       IID MISS_PHENO   N_MISS   N_GENO   F_MISS
1334   NA12144          Y    15077  2239392 0.006733
1334   NA12145          Y    19791  2239392 0.008838
1334   NA12146          Y    13981  2239392 0.006243
1334   NA12239          Y    14072  2239392 0.006284
1340   NA06994          Y    16080  2239392 0.007181
1340   NA07000          Y    26113  2239392  0.01166
1340   NA07022          Y    17467  2239392   0.0078
1340   NA07056          Y    12133  2239392 0.005418
1341   NA07034          Y    20425  2239392 0.009121

查看位點(diǎn)缺失的文件：

$ head missing_data.lmiss
 CHR          SNP   N_MISS   N_GENO   F_MISS
   1   rs12565286        1       60  0.01667
   1   rs12138618        0       60        0
   1    rs3094315        2       60  0.03333
   1    rs3131968        0       60        0
   1   rs12562034        0       60        0
   1    rs2905035        0       60        0
   1   rs12124819        0       60        0
   1    rs2980319        0       60        0
   1    rs4040617        0       60        0

6.3 哈溫檢測(cè)

--hardy

代碼：

plink --bfile 1kg_EU_qc --hardy

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to plink.log.
Options in effect:
  --bfile 1kg_EU_qc
  --hardy

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
--hardy: Writing Hardy-Weinberg report (founders only) to plink.hwe ... done.

結(jié)果：

$ head plink.hwe
 CHR         SNP     TEST   A1   A2                 GENO   O(HET)   E(HET)            P
   1   rs1048488  ALL(NP)    C    T            9/136/234   0.3588   0.3238      0.03894
   1   rs3115850  ALL(NP)    T    C            8/135/236   0.3562    0.319      0.02412
   1   rs2519031  ALL(NP)    G    A             0/11/368  0.02902   0.0286            1
   1   rs4970383  ALL(NP)    A    C           22/149/208   0.3931   0.3796       0.5881
   1   rs4475691  ALL(NP)    T    C           13/118/248   0.3113   0.3078            1
   1   rs1806509  ALL(NP)    C    A           67/171/141   0.4512   0.4809       0.2402
   1   rs7537756  ALL(NP)    G    A           14/124/241   0.3272   0.3206       0.8725
   1  rs28576697  ALL(NP)    C    T           37/161/181   0.4248   0.4278       0.9045
   1   rs7523549  ALL(NP)    T    C             0/28/351  0.07388  0.07115            1
[dengfei@localhost 08_chapter]$ wc -l plink.hwe
851066 plink.hwe

6.4 雜合率檢測(cè)

--het

代碼：

plink --bfile 1kg_EU_qc --het

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to plink.log.
Options in effect:
  --bfile 1kg_EU_qc
  --het

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
851065 variants and 379 people pass filters and QC.
Note: No phenotypes present.
--het: 851065 variants scanned, report written to plink.het .

結(jié)果：

$ head plink.het
 FID       IID       O(HOM)       E(HOM)        N(NM)            F
   0   HG00096       582990    5.838e+05       851065    -0.003127
   0   HG00097       582246    5.838e+05       851065    -0.005911
   0   HG00099       582456    5.838e+05       851065    -0.005125
   0   HG00100       582700    5.838e+05       851065    -0.004212
   0   HG00101       581527    5.838e+05       851065    -0.008602
   0   HG00102       585010    5.838e+05       851065     0.004432
   0   HG00103       583984    5.838e+05       851065    0.0005925
   0   HG00104       586189    5.838e+05       851065     0.008844
   0   HG00106       583868    5.838e+05       851065    0.0001584
$ wc -l plink.het
380 plink.het

7. 數(shù)據(jù)質(zhì)控qc

7.1 樣本質(zhì)控--mind

代碼：

 plink --bfile 1kg_EU_qc --mind 0.1 --recode --out re1

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re1.log.
Options in effect:
  --bfile 1kg_EU_qc
  --mind 0.1
  --out re1
  --recode

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
0 people removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
851065 variants and 379 people pass filters and QC.
Note: No phenotypes present.
--recode ped to re1.ped + re1.map ... done.

結(jié)果：

$ wc -l re1.*
        24 re1.log
    851065 re1.map
       379 re1.ped
    851468 總用量

7.2 位點(diǎn)質(zhì)控--geno

代碼：

plink --bfile 1kg_EU_qc --geno 0.1 --recode --out re1

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re1.log.
Options in effect:
  --bfile 1kg_EU_qc
  --geno 0.1
  --out re1
  --recode

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
0 variants removed due to missing genotype data (--geno).
851065 variants and 379 people pass filters and QC.
Note: No phenotypes present.
--recode ped to re1.ped + re1.map ... done.

結(jié)果：

$ wc -l re1*
        24 re1.log
    851065 re1.map
       379 re1.ped
    851468 總用量

7.3 maf質(zhì)控--maf

代碼：

plink --bfile 1kg_EU_qc --maf 0.01 --recode --out re1

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re1.log.
Options in effect:
  --bfile 1kg_EU_qc
  --maf 0.01
  --out re1
  --recode

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
17767 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
833298 variants and 379 people pass filters and QC.
Note: No phenotypes present.
--recode ped to re1.ped + re1.map ... done.

結(jié)果：

$ wc -l re1*
        25 re1.log
    833298 re1.map
       379 re1.ped
    833702 總用量

7.4 哈溫質(zhì)控 --hwe

代碼：

plink --bfile 1kg_EU_qc --hwe 1e-5 --recode --out re1

日志：

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re1.log.
Options in effect:
  --bfile 1kg_EU_qc
  --hwe 1e-5
  --out re1
  --recode

1031523 MB RAM detected; reserving 515761 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
--hwe: 25 variants removed due to Hardy-Weinberg exact test.
851040 variants and 379 people pass filters and QC.
Note: No phenotypes present.
--recode ped to re1.ped + re1.map ... done.

結(jié)果：

$ wc -l re1*
        24 re1.log
    851040 re1.map
       379 re1.ped
    851443 總用量