GWAS QC – 性别检查 与 阈值选择 Sex check

性别错误 (Sex discrepancy):有时数据录入会有错误,我们应基于X染色体的杂合性,确定是否有性别录入错误的个体。PLINK默认对于男性,X染色体的纯合性(homozygosity)估计值应高于0.8,女性应低于0.2,但实际情况更为复杂,本文简要对性别检查 与 阈值选择进行讲解。

本文使用数据来源:https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608

PLINK中提供了check-sex选项,通过计算X染色体的纯合性来判断样本性别是否有错误:—check-sex [女性阈值] [男性阈值] 默认的参数为 0.2 0.8,但该标准过于严格,并不适用与实际情况,下面会讲解。首先介绍基本操作:

#首先进行LD-pruning
plink \
        --bfile HapMap_3_r3_1 \
        --indep-pairwise 50 5 0.2 \
        --out HapMap_3_r3_1

#估计X染色体的纯合性
plink \
        --bfile HapMap_3_r3_1 \
        --extract HapMap_3_r3_1.prune.in \
	--check-sex \
        --out HapMap_3_r3_1

注意:

  1. 要提前确认X染色体上PAR区域已经被分离 (可以使用—split-x 选项来进行分离)
  2. 默认的阈值是,男性>0.8, 女性<0.2

log文件如下:

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 HapMap_3_r3_1.log.
Options in effect:
  --bfile HapMap_3_r3_1
  --check-sex
  --extract HapMap_3_r3_1.prune.in
  --out HapMap_3_r3_1

191875 MB RAM detected; reserving 95937 MB for main workspace.
Allocated 3037 MB successfully, after larger attempt(s) failed.
1457897 variants loaded from .bim file.
165 people (80 males, 85 females) loaded from .fam.
112 phenotype values loaded from .fam.
--extract: 166240 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 53 nonfounders present.
Calculating allele frequencies... done.
Warning: 35 het. haploid genotypes present (see HapMap_3_r3_1.hh ); many
commands treat these as missing.
Total genotyping rate is 0.996691.
166240 variants and 165 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
are missing.)
--check-sex: 3993 Xchr and 0 Ychr variant(s) scanned, 1 problem detected.
Report written to HapMap_3_r3_1.sexcheck .

–check-sex: 3993 Xchr and 0 Ychr variant(s) scanned, 1 problem detected.从日志文件中可以看出我们检测出了1个不一致的样本。上述代码运行结果如下:

head HapMap_3_r3_1.sexcheck

FID       IID       PEDSEX       SNPSEX       STATUS            F
   1328   NA06989            2            2           OK      0.02151
   1377   NA11891            1            1           OK            1
   1349   NA11843            1            1           OK            1
   1330   NA12341            2            2           OK     0.009813
   1444   NA12739            1            1           OK            1
   1344   NA10850            2            2           OK       0.0533
   1328   NA06984            1            1           OK       0.9989
   1463   NA12877            1            1           OK            1
   1418   NA12275            2            2           OK      -0.0747

前两列为FID与IID,三四列为PED文件中的性别,与通过X染色体估计的性别,第五列表示是否一致,不一致时会表示为PROBLEM,第六列是估计的F系数。可以提取出不一致的样本:

cat HapMap_3_r3_1.sexcheck | grep PROBLEM
1349   NA10854            2            1      PROBLEM       0.9933

#将有问题的样本FID与IID输出到sex_error_to_remove.id文件中
cat HapMap_3_r3_1.sexcheck | grep PROBLEM | awk '{print $1,$2}'> sex_error_to_remove.id

后续PLINK分析可以使用 —remove sex_error_to_remove.id 来剔除性别错误的样本.以上结果的F系数分布如图所示:

左边为女性,右边为男性

实际情况:上面的例子是一种很理想的分布,可以看到样本被清晰地分开,但实际情况要复杂很多,例如在千人基因组项目Phase1的样本中,一些女性样本的F系数也大于0.8,即使进行LD-pruning后,女性样本最大的F系数也超过了0.66。所以一般来说,我们会进行两次—check-sex,第一次用于观察分布,在看到一个右侧紧凑的峰(男性),与左侧平缓的峰后(女性),根据其分界线来决定判定阈值,(真实的数据的分布大多如下图所示)

这里常用的阈值为 -check-sex 0.7 0.8

#如果使用真实数据,通过目视决定阈值后,更改参数,重新判定性别
plink \
        --bfile HapMap_3_r3_1 \
        --extract HapMap_3_r3_1.prune.in \
	   --check-sex 0.7 0.8 \
        --out HapMap_3_r3_1

参考:

https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608

Basic statistics

发表评论

Fill in your details below or click an icon to log in:

WordPress.com 徽标

您正在使用您的 WordPress.com 账号评论。 注销 /  更改 )

Twitter picture

您正在使用您的 Twitter 账号评论。 注销 /  更改 )

Facebook photo

您正在使用您的 Facebook 账号评论。 注销 /  更改 )

Connecting to %s