QC – GWASLab – GWAS实验室

哈迪-温伯格平衡

回顾：哈迪温伯格平衡 Hardy– Weinberg equilibrium

哈迪-温伯格平衡精确检验检验原理

假设有N个无亲缘关系的样本（对应有2N个等位）

在哈迪温伯格平衡下，在N个样本的群体中观察到有 $n_{AB}$ 个样本为AB基因型的精确概率为：

$P(N_{AB} = n_{AB} | N, n_A) = {{2^{n_{AB}}}N!\over{n_{AA}!n_{AB}!n_{BB}!}} \times {{n_A!n_B!}\over{n_A!n_B!}}$

计算哈迪温伯格平衡精确检验的统计量时，我们需要把概率小于观察到的概率（ $n_{AB}$ 个样本为AB基因型）的情况的概率进行加和，如下所示：

$P_{HWE} = \sum_{n^{*}_{AB}} I[P( N_{AB} = n_{AB}|N, n_A)$

$\geqq P(N_{AB} = n^{*}_{AB} | N, n_A)] \times P(N_{AB} = n^{*}_{AB} | N, n_A)$

$I(x)$ 为一个指示函数. 如果x为真, $I(x) = 1$ ; 否则, $I(x) = 0$ .

实际使用软件计算时，通常会采用一些近似方法来避免大量的计算，可以参考PLINK中的HWE的算法。

使用PLINK进行HWE检验

PLINK提供了计算哈迪温伯格平衡精确检验的统计量--hardy以及基于统计量进行过滤--hwe的选项：

plink \
    --bfile ${genotypeFile} \
    --hardy \
    --out plink_results

输出结果如下， P列即为哈迪温伯格平衡精确检验的结果：

$ head plink_results.hwe
 CHR              SNP     TEST   A1   A2                 GENO   O(HET)   E(HET)            P 
   1      1:13273:G:C  ALL(NP)    C    G             1/61/442    0.121   0.1172       0.7113
   1      1:14599:T:A  ALL(NP)    A    T             1/88/415   0.1746   0.1626       0.1625
   1      1:14604:A:G  ALL(NP)    G    A             1/88/415   0.1746   0.1626       0.1625
   1      1:14930:A:G  ALL(NP)    G    A             4/409/91   0.8115   0.4851    1.679e-61
   1      1:69897:T:C  ALL(NP)    T    C            7/111/386   0.2202   0.2173            1
   1      1:86331:A:G  ALL(NP)    G    A             0/88/416   0.1746   0.1594      0.02387
   1      1:91581:G:A  ALL(NP)    A    G          137/228/139   0.4524      0.5      0.03271
   1     1:122872:T:G  ALL(NP)    G    T            1/259/244   0.5139   0.3838     8.04e-19
   1     1:135163:C:T  ALL(NP)    T    C             1/91/412   0.1806   0.1675       0.1066

或者可以通过--hwe 1e-6 直接过滤掉P小于1e-6的SNP

plink \
    --bfile ${genotypeFile} \
    --hwe 1e-6 \
    --out plink_results

参考

https://www.cog-genomics.org/plink/1.9/dev#exact

https://www.cog-genomics.org/plink/1.9/basic_stats#hardy

Wigginton, J. E., Cutler, D. J., & Abecasis, G. R. (2005). A note on exact tests of Hardy-Weinberg equilibrium. The American Journal of Human Genetics, 76(5), 887-893. Link

rsID的介绍与chr:pos转换时的陷阱

很多小伙伴都觉得位置转换rsID是很麻烦的事情，有时会偷懒只用手头文件的chr:pos位置信息匹配rsID，但这样做带来的的问题却少有人讨论，本文将主要介绍什么是rsID，以及rsID在使用和转换中的一些常见问题。

本文内容
什么是rsID
主要优点
rsID可能表示的变异类型
（重点）rsID与chrpos转换时的常见错误
解决办法
参考

什么是rsID

rsID 就是 dbSNP的Reference SNP ID （缩写为rs 或者RefSNP），一个由dbSNP设定的，为了识别变异位点的一串数字编号。rsID设计上是非冗余的，也就是全局唯一的id，用户提交的变异会被归类整理注释，重复的变异会被整合。

主要优点

rsID无关参考基因组版本，不像chrpos会随版本变化而变化, rsID在不同版本间是一致的。对于群体遗传学或是精准医学的大规模的研究来说会更加方便，rsID提供了稳定的变异表示方法。（摘自官网，个人认为有时候rsID的转换带来的问题远超不转换的问题，有好有坏，但是传统上还是需要转换）

rsID表示变异的类型

rsID中的rs尽管是Reference SNP的首字母缩写，但实际上一些其他类型的变异也会被赋予rsID。（通常变异的长度小于50bp）

单核苷酸变异 Single nucleotide variation (SNV)
短多核苷酸变异 Short multi-nucleotide changes (MNV)
较小的短插入或删除 Small deletions or insertions (INDEL)
较小的短串联重复序列 Small STR repeats
逆转录转座子插入 retrotransposable element insertions

rsID是把双刃剑：仅凭chr：pos与rsID互相转换时的陷阱

仅使用chr:pos 转换 rsID时的问题：

对应位点rsID不存在，可能是新变异等等原因，通常可以以chr:pos:ref:alt的形式替代。但还有个问题就是Alt allele不存在，比如rs123456 对应chr1：123456的 T>C,A 而你手里的数据是chr1：123456的 T>G，那问题来了，这应不应该给他们相同的rsID？仅凭位点和类型来说应该给，或许下个版本的dbsnp会加上这个变异，但其实我也没有明确的答案（欢迎评论区讨论），不过实际操作中我会倾向于保守一点，用chr:pos:ref:alt 而不是rsID来表示。
如上rsID的介绍所述，rsID并不止只用来表示单一核苷酸的SNP，也会表示其他变异类型，这会导致同一位点有多个rsID表示的变异，最常见的就是某个位点同时有SNP和INDEL，仅凭chr:pos信息而不管allele的话会混淆并大量的错误匹配SNP与INDEL的rsID，后续功能分析会引起很大的不便，举个例子： rs123456 对应chr1：123456的 T>C ,而rs987654 同样对应chr1：123456这个位置，但是这个变异是个INDEL， T>TA，如果仅凭chr：pos匹配会混淆SNP与INDEL，虽然是同样的位置，但变异造成的影响会完全不同。解释时本应是rs987654这个INDEL造成的影响却错误地解释到rs123456这个SNP上，这种情况应该被避免。这么做破坏了rsID的唯一性特点，是不是有点违背初衷，本末倒置了。
还有一个问题就是手头数据里的变异是否已经标准化？ 未标准化的变异的chrpos是不准确的，进行左对齐与节俭原则的标准化后可能产生位移，用未标准化chrpos匹配时可能会错位匹配到其他相邻的位点上。比如手头的变异可能是 chr1:123456:AA:AT ,标准化后则是chr1:123457:A:T，向后移了一位，如果你看过1000genome的原始数据就会发现这样的情况大量存在，所以应当注意（参考：GWASLab：变异的标准化 Variant Normalization）
0起点还是1起点的参考系问题，处理数据时应该注意，这里不做过多赘述。（GWASLab：LiftOver 基因组坐标变换与 01坐标系统）

rsID 向 chr：pos 某参考基因组版本的位置转换时，会遇到的问题：

设计上rsID是唯一对应某个变异的，但实际上由于dbSNP版本的不同或其他原因，手头GWAS的sumstats里的rsID可能对应两个位置，而多个rsID又可能对应同一个位置上相同的变异
在对应参考基因组版本上的位置不存在等等

解决办法

rsID转换chrpos时要尽量明确原始数据的dbsnp版本，能确定版本的时候用对应版本，不能的时候要制定统一标准（为了研究的可重复性），转换时要使用统一的dbsnp的版本。

而chrpos转换rsID时，不贪多，不求快，老老实实用先确认标准化，然后利用注释的方法，也就是相应基因组版本的位置chr:pos以及 ref与alt全部与rsID全部匹配时才进行转换。

可以参考以下内容：

GWASLab：使用ANNOVAR对变异进行功能注释

GWASLab：SNP的rsID与位置信息的相互匹配 rsID/ chr:pos conversion

参考：

https://www.ncbi.nlm.nih.gov/snp/do

使用UMAP对基因组数据降维，对比PCA

本文内容:
UMAP简介
在群体遗传学中的应用（与PCA的不同）
UMAP使用方法
使用对示例千人基因组数据进行降维
参考

关键词：UMAP, PCA, t-SNE, PCA-UMAP, 基因组降维

UMAP 简介

UMAP（uniform manifold approximation and projection）是近年来新出现的一种相对灵活的非线性降维算法，目前在统计遗传学等领域也有了较为广泛的应用。

UMAP的理论基础基于流形理论（manifold theory）与拓扑分析，

主要基于以下假设：

存在一个数据均一分布的流形。
这个目标流形是局部相连的。
该算法的主要目标是保存此流形的拓扑结构。

总体来看，UMAP利用了局部流形近似，并拼接模糊单纯集合表示（local fuzzy simplicial set representation），以构建高维数据的拓扑表示。在给定低维数据时，UMAP会采取相似的手法构建一个等价的拓扑表示。最后UMAP会对低维空间中数据表示的布局进行最优化，以最小化高维和低维两种拓扑表示之间的交叉熵（cross-entropy）。

算法概览如下：

UMAP与其他降维方法对比：

快，准确：需要降维数据维度增加时，t-SNE计算时间呈指数型增长，umap则线性增长。

在群体遗传学中的应用：

在处理大样本基因组数据时，通常我们会对数据进行降维，以达到数据可视化，并发现存在亲缘关系的样本。

最常用的方法是主成分分析PCA：GWASLab：群体分层与主成分分析教程 Population structure & PCA

UMAP与PCA的不同：

但由于PCA投影找的是最大化方差的方向，通常会忽略掉其他方向的方差，也就是说PCA倾向于发现较大的人群结构，而忽略精细的结构，但UMAP的思路则与PCA不同，如第一节所述，其主要目的是保留其与邻接样本的拓扑结构，也就是更为精细的局部人群结构，而非整体结构。

目前UMAP已经被广泛应用于群体遗传学的研究之中。

下面简单介绍UMAP使用方法：

UMAP （是一个python包），使用方法:

通过pip或conda安装：

pip install umap-learn
#或者
conda install -c conda-forge umap-learn

使用文档详见：https://umap-learn.readthedocs.io/en/latest/parameters.html

一个最简单的示例：

#导入umap包
import umap   

#创建UMAP实例
reducer = umap.UMAP()   

#使用你的数据对该实例进行训练，得到嵌入后的结果
embedding = reducer.fit_transform(your_data)

UMAP() 的主要参数，可参考上述算法概览

n_neighbors ：邻接样本的数量，默认15
min_dist ：控制布局的参数，取值范围0-1，默认0.1
n_components ：维数，默认2

使用基因型数据进行实际演示：

教程改编与数据改编自：

Sakaue, S., Hirata, J., Kanai, M. et al. Dimensionality reduction reveals fine-scale structure in the Japanese population with consequences for polygenic risk prediction. Nat Commun 11, 1569 (2020). https://doi.org/10.1038/s41467-020-15194-z

代码与数据下载：

https://github.com/saorisakaue/Genotype-dimensionality-reduction

数据在umap_data中：1KG.selected.bed/bim/fam

使用软件：

PLINK，python（umap,pandas,matplotlib）

注意：

UMAP只接受完整的数据，对于有缺失的基因型数据要事先对缺失的SNP进行填补以得到完整的数据，填补方法可参考：GWASLab：Eagle2单倍型定相工具 Haplotype phasing
对于PLINK格式的数据，一般需要在LD-pruning 后，转换为 genotype matrix

输入文件预处理：

#如果数据有缺失请先进行填补

#首先进行LD-pruning
GENOTYPE="./data_umap/1KG.selected"

plink \
 --bfile ${GENOTYPE} \
 --indep-pairwise 50 5 0.2 \
 --maf 0.01 \
 --hwe 1E-6 \
 --out ${GENOTYPE}

#对二进制plink文件进行转换,转为0/1/2的文本文件
plink \
 --bfile ${GENOTYPE} \
 --extract ${GENOTYPE}.prune.in \
 --recode A \
 --out ${GENOTYPE}.pruned

# 2) 抽取出 genotype matrix
cat ${GENOTYPE}.pruned.raw | cut -d " " -f7- | awk 'NR>1{print}' > ${GENOTYPE}.pruned.geno.txt

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 ./data_umap/1KG.selected.log.
Options in effect:
  --bfile ./data_umap/1KG.selected
  --hwe 1E-6
  --indep-pairwise 50 5 0.2
  --maf 0.01
  --out ./data_umap/1KG.selected

191875 MB RAM detected; reserving 95937 MB for main workspace.
Allocated 3037 MB successfully, after larger attempt(s) failed.
120489 variants loaded from .bim file.
50 people (25 males, 25 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 50 founders and 0 nonfounders present.
Calculating allele frequencies... done.
--hwe: 14 variants removed due to Hardy-Weinberg exact test.
861 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
119614 variants and 50 people pass filters and QC.
Note: No phenotypes present.
Pruned 5127 variants from chromosome 1, leaving 4909.
Pruned 5171 variants from chromosome 2, leaving 4680.
Pruned 4052 variants from chromosome 3, leaving 3919.
Pruned 3398 variants from chromosome 4, leaving 3517.
Pruned 3594 variants from chromosome 5, leaving 3586.
Pruned 4346 variants from chromosome 6, leaving 3628.
Pruned 3298 variants from chromosome 7, leaving 3193.
Pruned 3421 variants from chromosome 8, leaving 3027.
Pruned 2854 variants from chromosome 9, leaving 2895.
Pruned 3445 variants from chromosome 10, leaving 3244.
Pruned 3235 variants from chromosome 11, leaving 3021.
Pruned 2936 variants from chromosome 12, leaving 3020.
Pruned 2281 variants from chromosome 13, leaving 2326.
Pruned 1916 variants from chromosome 14, leaving 2070.
Pruned 1875 variants from chromosome 15, leaving 1943.
Pruned 1848 variants from chromosome 16, leaving 2005.
Pruned 1589 variants from chromosome 17, leaving 1839.
Pruned 1744 variants from chromosome 18, leaving 1925.
Pruned 947 variants from chromosome 19, leaving 1208.
Pruned 1528 variants from chromosome 20, leaving 1620.
Pruned 874 variants from chromosome 21, leaving 887.
Pruned 759 variants from chromosome 22, leaving 914.
Pruning complete.  60238 of 119614 variants removed.
Marker lists written to ./data_umap/1KG.selected.prune.in and
./data_umap/1KG.selected.prune.out .
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 ./data_umap/1KG.selected.pruned.log.
Options in effect:
  --bfile ./data_umap/1KG.selected
  --extract ./data_umap/1KG.selected.prune.in
  --out ./data_umap/1KG.selected.pruned
  --recode A

191875 MB RAM detected; reserving 95937 MB for main workspace.
Allocated 3037 MB successfully, after larger attempt(s) failed.
120489 variants loaded from .bim file.
50 people (25 males, 25 females) loaded from .fam.
--extract: 59376 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 50 founders and 0 nonfounders present.
Calculating allele frequencies... done.
59376 variants and 50 people pass filters and QC.
Note: No phenotypes present.
--recode A to ./data_umap/1KG.selected.pruned.raw ... done.

接下来使用umap的默认参数进行降维处理：

import umap

prefix = sys.argv[1]
filename = prefix + ".pruned.geno.txt"
pre_data = pd.read_table(filename, delim_whitespace=True, header=None, low_memory=False)
data = pre_data.dropna(how='any', axis=1)

# UMAP
outfile = prefix + ".umap.txt"
d_embedding = umap.UMAP().fit_transform(data)
pd.DataFrame(d_embedding).to_csv(outfile, sep='\t', header=False, index=False)

查看d_embedding ，基因组数据已经降至2维，结果如下：

array([[ 29.74708  ,   3.0977843],
       [ 29.213293 ,   3.1712177],
       [  8.758959 ,  -1.9003359],
       [  8.377356 ,  -1.4676492],
       [ 29.458044 ,   3.3701472],
       [ 29.87717  ,   3.4556358],
       [ 29.3555   ,   3.8484392],
       [ 28.955818 ,   3.5634878],
       [ 29.105188 ,   2.8024642],
       [ 29.197195 ,   3.7704961],
       [ 29.347229 ,   4.2163496],
       [ 29.865662 ,   3.1569576],
       [  8.503107 ,  -1.3564333],
       [-24.176699 , -14.017474 ],
       [-23.141544 , -13.751009 ],
       [  8.204682 ,  -1.8065634],
       [ 29.94206  ,   2.8537552],
       [-23.501389 , -13.829134 ],
       [ 29.065151 ,   4.255709 ],
       [-23.86171  , -13.711441 ],
       [-23.297297 , -13.290715 ],
       [ 30.9667   , -10.477952 ],
       [-22.91213  , -13.484747 ],
       [-23.249172 , -14.22486  ],
       [-22.951435 , -13.851592 ],
       [ 30.514542 , -11.441013 ],
       [ 30.728407 , -11.229405 ],
       [ 31.13736  , -11.319468 ],
       [ 30.189419 , -10.865067 ],
       [ 30.20413  , -11.283942 ],
       [ 30.89431  , -11.590417 ],
       [ 29.568645 ,   2.9131896],
       [ 29.632095 ,   2.6519256],
       [-23.177608 , -14.538039 ],
       [-23.335829 , -13.95216  ],
       [  8.488959 ,  -2.0645304],
       [  8.731548 ,  -1.8205943],
       [  8.5880165,  -1.6314453],
       [-23.663996 , -13.476801 ],
       [-23.71588  , -14.627886 ],
       [-23.864864 , -14.236027 ],
       [-23.965855 , -14.178329 ],
       [-23.203787 , -14.796338 ],
       [-23.587828 , -14.443732 ],
       [-22.973204 , -14.444413 ],
       [ 30.451714 , -10.90268  ],
       [ 31.011879 , -11.050354 ],
       [ 30.863823 , -10.923892 ],
       [ 30.679998 , -11.574213 ],
       [ 30.77399  , -10.824773 ]], dtype=float32)

绘图后可以看到样本被清楚地分到了4个cluster：

可以自己试着改变UMAP（）的参数，看看有什么变化，这里不过多演示，（下图为UMAP(min_dist=1)的结果）

当然，也可以使用PCA后各个样本主成分PCs进行UMAP （PCA-UMAP）：

PCA方法：GWASLab：群体分层与主成分分析教程 Population structure & PCA

可以参考：Sakaue, S., Hirata, J., Kanai, M. et al. Dimensionality reduction reveals fine-scale structure in the Japanese population with consequences for polygenic risk prediction. Nat Commun 11, 1569

该文章对日本人群结构做了分析，并对比了各种降维方法的结果，强烈推荐。

参考

McInnes, L, Healy, J, UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction, ArXiv e-prints 1802.03426, 2018

Sakaue, S., Hirata, J., Kanai, M. et al. Dimensionality reduction reveals fine-scale structure in the Japanese population with consequences for polygenic risk prediction. Nat Commun 11, 1569

Diaz-Papkovich, A., Anderson-Trocmé, L. & Gravel, S. A review of UMAP in population genetics. J. Hum. Genet. 66, 85–91 (2021).

https://umap-learn.readthedocs.io

SNP的rsID与位置信息的相互匹配 rsID/ chr:pos conversion

本文内容：

1 个别SNP位点或rsID查询 
	dbsnp网页查询
2 少量转换可用的网页工具 （100,000 snp以内）
	SNP-nexus网页工具 
	vep注释工具网页版
3 大量转换使用的命令行工具 
	ANNOVAR命令行
	VEP命令行
4 参考

rsID与染色体位置信息的匹配或转换是生物信息学研究中必不可少的，但却又是十分繁琐的一项步骤，很多同学都在纠结这个问题，本文将总结常用的转换方法，以供参考：

转换前需要考虑的问题：

变异的标准化：Variant Normalization 变异的标准化
参考基因组的版本：人类参考基因组 Human reference genome hg19/hg38/GRCh37/GRCh38/b37/hs375d

1.个别SNP的转换：

直接使用dbsnp（详见：SNP数据库 – dbSNP）查询，查询个别SNP时方便快捷，可以直接搜rsID

也可以以chr:pos的形式搜rsID

2.少量转换可用的网页工具（100,000 snp以内）

2.1 SNPnexus : https://www.snp-nexus.org/v4/

首先输入用户信息，学术用途是免费的，使用自己的edu邮箱即可

之后选择assembly的版本

选择后可以通过多种方式提交自己要查询的SNP

2.1.1 这里我们通过上传文件的方式对rsID进行注释

输入文件 snp.txt 格式如下（其他格式的具体描述详见：https://www.snp-nexus.org/v4/guide/）

dbsnp   rs293794
dbsnp   rs1052133
dbsnp   rs3136820
dbsnp   rs2272615
dbsnp   rs2953993
dbsnp   rs1799782
dbsnp   rs25487
dbsnp   rs2248690
dbsnp   rs4918
dbsnp   rs1071592

然后是注释用的数据选择，各数据库版本见：https://www.snp-nexus.org/v4/guide/#data_source

因为我们只查询位置信息，就不勾选其他数据库，只使用默认的Ensembl

点击submit query后，稍作等待，结果就会显示出来，可以导出为VCF或txt文本格式。

2.1.2 位点转换为rsID时，输入文件如下：

格式为：< Type Name Position Allele1 Allele2 Strand > # Genomic position data for novel SNPs

chromosome	1	9262273	G	A	1
chromosome	1	45015006	G	A	1
chromosome	1	226982947	C	T	1
chromosome	2	189001450	G	A	1
chromosome	17	46010375	T	C	1

提交后即可查询，结果如下所示

2.2 VEP网页版

VEP注释工具使用详见（预留链接）

VEP网页版：https://asia.ensembl.org/info/docs/tools/vep/online/index.html

点击launch VEP

2.2.1 可以在input data处直接粘贴要查询的rsID，或是上传文件

点击输入框下方的example，可以查看可用的输入格式

这里以Variant identifiers 为例：

选择合适数据库，提交后可以看到查询状态

完成后点击View results

可查看基本信息，也可以导出

2.2.2 位点转换rsID时，使用如下Ensembl默认的输入格式即可：

3.大量转换时的命令行工具

3.1 可使用如上所述VEP工具的命令行版本

下载地址与文档见：https://asia.ensembl.org/info/docs/tools/vep/index.html

3.2 也可以使用ANNOVAR

安装与使用具体使用方法见：使用ANNOVAR 对Variants进行功能注释 Annotation POST-GWAS analysis

这里只简单介绍rsID与chr：pos互相转换的方法

3.2.1 rsID转chr：pos

使用ANNOVAR的convert2annovar.pl 程序与-format rsid选项即可注释位置信息，使用例如下所示，

（注意参考基因组的版本）

[kaiwang@biocluster ~/]$ cat example/snplist.txt 
rs74487784
rs41534544
rs4308095
rs12345678

[kaiwang@biocluster ~/]$ convert2annovar.pl -format rsid example/snplist.txt -dbsnpfile humandb/hg19_snp138.txt > snplist.avinput 
NOTICE: Scanning dbSNP file humandb/hg19_snp138.txt...
NOTICE: input file contains 4 rs identifiers, output file contains information for 4 rs identifiers
WARNING: 1 rs identifiers have multiple records (due to multiple mapping) and they are all written to output

[kaiwang@biocluster ~/]$ cat snplist.avinput 
chr2 186229004 186229004 C T rs4308095
chr7 6026775 6026775 T C rs41534544
chr7 6777183 6777183 G A rs41534544
chr9 3901666 3901666 T C rs12345678
chr22 24325095 24325095 A G rs74487784

3.2.2 chr：pos转rsID

位点信息转化为rsID时需要使用基于筛选的注释：

这里使用avsnp150，具体描述可见：https://annovar.openbioinformatics.org/en/latest/user-guide/download/

使用方法如下：

#首先下载注释用数据库：
annotate_variation.pl -downdb -webfrom annovar -buildver hg19 avsnp150 humandb/

#完成后即可进行转换
annotate_variation.pl ex1.avinput humandb/ -filter -build hg19 -dbtype avsnp150

输入文件使用上述的（类bed文件），ex1.avinput

chr2 186229004 186229004 C T
chr7 6026775 6026775 T C
chr7 6777183 6777183 G A
chr9 3901666 3901666 T C
chr22 24325095 24325095 A G

输出文件如下，ex1.avinput.hg19_avsnp150_dropped，

（注意某些SNP的rsID有合并等原因，版本不同rsID注释结果不一定相同）

avsnp150        rs4308095       chr2 186229004 186229004 C T
avsnp150        rs140195        chr22 24325095 24325095 A G
avsnp150        rs2228006       chr7 6026775 6026775 T C
avsnp150        rs766725467     chr7 6777183 6777183 G A
avsnp150        rs12345678      chr9 3901666 3901666 T C

注释rsID可能会遇到的问题：

https://annovar.openbioinformatics.org/en/latest/articles/dbSNP/

参考：

dbsnp：https://www.ncbi.nlm.nih.gov/snp/

annovar：https://annovar.openbioinformatics.org/en/latest/

vep：https://asia.ensembl.org/info/docs/tools/vep/index.html

snpnexus：https://www.snp-nexus.org/v4/

基因组控制与基因组膨胀系数λ Genomic control λGC

本文内容


什么是基因组控制 Genomic Control
基因组膨胀系数λGC （Genomic Inflation Factor） 的计算方法
GC的注意事项以及解决办法

什么是基因组控制 Genomic Control

基因组控制是矫正GWAS中因群体分层等原因而导致检验统计量膨胀的一种方法。

考虑一种简单病例对照研究，假设有N个样本，以及n个SNP，φ是样本中病例的比例(0<φ<1), 对于某个SNP，其在病例以及对照中的频率如下表所示：

为了检验关联性，我们会基于隐形，显性，以及加性遗传模型来计算自由度为1的卡方检验量。这里以加性遗传模型为例，对SNP进行趋势检验 ( Armitage’s trend test）。其计算公式为：

我们假设对于每个SNP，我们都计算了其相应的趋势检验的Y2l统计量，当这些SNP之间没有连锁不平衡LD，且不存在人群分层以及隐性关联时，检验量Y2服从卡方分布。

基因组控制（Genomic control，GC）模型假设是这些检验统计量会因为群体分层以及隐形关联等原因出现膨胀，膨胀系数为 λ （基因组膨胀系数 Genomic Inflation Factor）。同时GC模型假设这个膨胀系数在基因组上对于所有SNP是近似相等的。

2.基因组膨胀系数 Genomic Inflation Factor 的计算方法

基于此我们可以通过下式来估计λGC ，

取GWAS检验后所有卡方检验量的中间值，除以0.456，其中0.456为卡方检验量的期望值（卡方分布中，第50百分位数的卡方统计量，r语言中qchisq(0.5，1)对应的值）。之所以取中间值计算λGC是因为要避免异常值的影响。

λ越接近1，就表明不存在群体分层导致的统计量膨胀。

将GWAS检验后所有卡方统计量除以λ后重新计算p值得过程即为基因组控制 GC。

例如这个GWAS研究的QQ图，可以看到观测值有一个明显的系统性的抬升，这通常意味着样本中存在在群体分层，通过计算我们得到这个GWAS研究的基因组膨胀系数为 λ=1.17，

将原始的统计量除以1.17，重新计算p值后，可以看到之前的抬升得到有效控制。

3. GC的注意事项以及解决办法：

但要注意的是，GC假设基因组中只有少数的loci与表型相关，绝大多数被检验的SNP是与表型无关的，而目前的主流GWAS检验方法大多基于无限小模型（infinitesimal model）(详见：解释复杂疾病的四种主流模型 CDCV/RAME/infinitesimal/Broad-sense-heritability)，该模型假设所有SNP与表型都是有关的，只是效应量很小。这种情况下就不再适用GC。

其他解决人群分层等的办法包括：

严格的质控QC GWAS前的质控（QC） GWAS Quality control
通过LD score回归来判别统计量的膨胀是因为人群分层等外部影响，还是多基因遗传结构的影响。连锁不平衡分数回归 LD score regression
使用LMM模型（线性混合模型 – LMM Linear Mixed Model）并将PCA的主成分加入协变量。群体分层与主成分分析 Population structure & PCA

等等

参考：

Devlin, B. & Roeder, K. Genomic control for association studies. Biometrics 55, 997–1004 (1999).

Devlin, B. Genomic control, a new approach to genetic-based associationb studies. (2001).

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

注意：

要提前确认X染色体上PAR区域已经被分离（可以使用—split-x 选项来进行分离）
默认的阈值是，男性>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

为什么在PCA或估计GRM时要去除长LD区域 Remove long-LD region

本文内容

1.背景介绍
2.为什么在PCA等时要去除LongLD区域？
3.长LD区域起始位置的列表
4.使用PLINK去除长LD区域里的SNP：

1.背景介绍：

LD ：连锁不平衡 linkage disequilibrium LD

PCA ：群体分层与主成分分析 Population structure & PCA

2.为什么在QC时要去除LongLD区域？

人类基因组中存在若干长LD的区域，这些区域多位于染色体的着丝粒附近，还有一些位于HLA等区域。如下图所示：

这些区域跨度很长（长度超过2Mb），单次LD-pruning无法完全去除互相成LD的SNP，在进行诸如PCA，或是计算GRM，进行基于LMM模型的GWAS分析时，我们应当去除掉这些区域。

长LD区域的形成并不一定是因为选择，其他原因诸如倒位多态性（inversion polymorphism）也可能造成长LD区域的存在。在进行研究时，应当谨慎区分这些区域形成的原因。如果在计算模型中没有对这些长LD区域进行处理，就可能影响群体遗传结构中对于本地群体的估计，造成系统性的偏倚。

3.长LD区域起始位置的列表（hg38,hg19与hg18参考基因组版本）

hg38版本

Chr	Start	Stop
chr1	47761740	51761740
chr1	125169943	125170022
chr1	144106678	144106709
chr1	181955019	181955047
chr2	85919365	100517106
chr2	87416141	87416186
chr2	87417804	87417863
chr2	87418924	87418981
chr2	89917298	89917322
chr2	135275091	135275210
chr2	182427027	189427029
chr2	207609786	207609808
chr3	47483505	49987563
chr3	83368158	86868160
chr5	44464140	51168409
chr5	129636407	132636409
chr6	25391792	33424245
chr6	26726947	26726981
chr6	57788603	58453888
chr6	61109122	61357029
chr6	61424410	61424451
chr6	139637169	142137170
chr7	54964812	66897578
chr7	62182500	62277073
chr8	8105067	12105082
chr8	43025699	48924888
chr8	47303500	47317337
chr8	110918594	113918595
chr9	40365644	40365693
chr9	64198500	64200392
chr9	88958735	88959017
chr10	36671065	43184546
chr10	41693521	41885273
chr11	88127183	91127184
chr12	32955798	41319931
chr12	34639034	34639084
chr14	87391719	87391996
chr14	94658026	94658080
chr17	43159541	43159574
chr20	4031884	4032441
chr20	33948532	36438183
chr22	30060084	30060162
chr22	42980497	42980522

hg19版本

Chr	Start	Stop ID
1 48000000 52000000 1
2 86000000 100500000 2
2 134500000 138000000 3
2 183000000 190000000 4
3 47500000 50000000 5
3 83500000 87000000 6
3 89000000 97500000 7
5 44500000 50500000 8
5 98000000 100500000 9
5 129000000 132000000 10
5 135500000 138500000 11
6 25000000 35000000 12
6 57000000 64000000 13
6 140000000 142500000 14
7 55000000 66000000 15
8 7000000 13000000 16
8 43000000 50000000 17
8 112000000 115000000 18
10 37000000 43000000 19
11 46000000 57000000 20
11 87500000 90500000 21
12 33000000 40000000 22
12 109500000 112000000 23
20 32000000 34500000 24

hg18版本

Chr	Start	Stop	ID
1	48060567	52060567	hild1
2	85941853	100407914	hild
2	134382738	137882738	hild3
2	182882739	189882739	hild4
3	47500000	50000000	hild5
3	83500000	87000000	hild6
3	89000000	97500000	hild7
5	44500000	50500000	hild8
5	98000000	100500000	hild9
5	129000000	132000000	hild10
5	135500000	138500000	hild11
6	25500000	33500000	hild12
6	57000000	64000000	hild13
6	140000000	142500000	hild14
7	55193285	66193285	hild15
8	8000000	12000000	hild16
8	43000000	50000000	hild17
8	112000000	115000000	hild18
10	37000000	43000000	hild19
11	46000000	57000000	hild20
11	87500000	90500000	hild21
12	33000000	40000000	hild22
12	109521663	112021663	hild23
20	32000000	34500000	hild24
X	14150264	16650264	hild25
X	25650264	28650264	hild26
X	33150264	35650264	hild27
X	55133704	60500000	hild28
X	65133704	67633704	hild29
X	71633704	77580511	hild30
X	80080511	86080511	hild31
X	100580511	103080511	hild32
X	125602146	128102146	hild33
X	129102146	131602146	hild34

4.使用PLINK去除长LD区域里的SNP：

我们可以使用PLINK来去除长LD区域里的SNP，主要分为两步：

1.将上一节中的列表拷贝进high-ld.txt文件中（使用时记得去掉header），使用--make-set选项提取区域中的SNP

2.在分析时利用--exclude去除掉上一步所生成列表中的SNP

plink --file mydata --make-set high-ld.txt --write-set --out hild
plink --file mydata --exclude hild.set --recode --out mydatatrimmed

参考：

https://genome.sph.umich.edu/wiki/Regions_of_high_linkage_disequilibrium_(LD)

Price et al. (2008) Long-Range LD Can Confound Genome Scans in Admixed Populations. Am. J. Hum. Genet. 86, 127-147

更新：

20220905 修改表述错误，更新PCA链接，并增加hg38版本

SNP的LD剪枝与聚集 LD pruning & clumping

GWAS相关的研究中，很多时候我们需要从总的SNP数据中，基于SNP两两之间的LD，来抽取出一个不含互相关联SNP的子集，目前主要的两种方法分别是 LD pruning 与 clumping。

例如，

在进行主成分分析（PCA）时，我们需要事先对SNP进行LD pruning 以去除互相关联（处于LD）的SNP，以防止高LD区域过高的方差对结果的影响。

在计算风险分数PRS时，我们需要从显著的loci中选取具有代表性的SNP来计算分线分数，这时就需要进行clumping，基于LD的r2，以及GWAS所得到的p值，来筛选出这个LD区域中的代表SNP（重要性最高），这样我们可以获得更准确的风险分数。

LD pruning 与 clumping 方法的异同如下所示：

	根据	保留	主要用途
Pruning	LD的R2	处于LD的一对SNP中MAF最高的	PCA
Clumping	LD的R2 与 SNP的P值	处于 LD的一对SNP中P值最显著的	PRS

Pruning 与 clumping 的主要区别

具体算法上，可以简单理解为：

Pruning：选取第一个SNP，然后计算这个SNP与窗口区间里第二个，第三个，等等的r2，当检测到高的相关性时，就会从这一对SNP中去除MAF较低的那个，保留 MAF 高的，也就是说这个过程中可能会去除掉我们选的第一个SNP。完成后下一步就是选取下一个SNP，重复这个过程。

Clumping：首先会依据从GWAS得到的p值对SNP的重要性进行排序，然后选取排序后的第一个SNP，计算这个SNP与窗口区间里其他SNP的r2，当检测到高的相关性时，就会从这一对SNP中去除重要性低的那个，这个过程中我们选的第一个SNP一定会得到保留。完成后下一步就是选取 p值排序后的下一个 SNP，重复这个过程。

PLINK中提供了 Pruning 和 Clumping 的功能：

Pruning：

我们主要是用–indep-pairwise选项，也就是根据SNP两两之间的LD来pruning。

--indep-pairwise  <window size>['kb']  <step size (variant ct)>  <r^2 threshold>

例 --indep-pairwise 500 50 0.2
这三个参数代表的意思分别是： 窗口大小，每一步移动窗口的距离，以及判定关联的r2阈值

plink -bfile input --indep-pairwise 500 50 0.2 --out input_pruned

输出两个文件
input_pruned.prune.in    #pruning后保留的互不相关的SNP
input_pruned.prune.out  #去除掉的SNP

Clumping:

PLINK提供了多种参数选项，具体可以参考：https://www.cog-genomics.org/plink/1.9/postproc

参考：

https://www.cog-genomics.org/plink/1.9/postproc

https://www.cog-genomics.org/plink/1.9/ld

https://www.biostars.org/p/343818/

群体分层与主成分分析教程 Population structure & PCA

本文内容

原理简介
使用千人基因组项目数据与PLINK进行PCA
绘图
参考

1.原理简介

主成分分析（Principal Components Analysis, PCA）是一种常用的数据降维方法，在群体遗传学中被广泛用于识别并调整样本的群体分层问题。群体分层会导致GWAS研究中的虚假关联，考虑一个case-control研究，如下图所示，红色的群体在整体样本中占了case的大多数，那么一些对疾病并没有影响，但在红蓝两群体之间等位基因频率相差很大的genetic marker，就有可能会造成虚假关联（spurious association）。

假设我们有N个样本，M个SNP的数据，其对应的矩阵记为G，Gij则为第i个样本的第j个SNP，取值范围为0, 0.5, 1 ,对应有0，1，2个参考等位基因。

首先对G进行标准化（减去均值后除以标准差），得到的标准化后N x M矩阵记为Z，

其中pj为第j个SNP的等位基因频率。

于是我们可以计算得到遗传关系矩阵GRM

GRM中第i行第j列的项就表示个体i与个体j的遗传相似性的均值。

对GRM进行特征分解即可得到对应不同的群体结构特征的主成分PC。

其中V是n个特征向量Vn组成的n x n 的矩阵，D则对应由n个特征值λn组成的向量。

这个n个特征值由大到小，λ1>λ2>…>λn, 特征值的大小与相对应特征向量所能解释方差成正比（也就是，第一主成分PC1能解释最多的方差，以此类推）。

主成分可以被视为反映由于祖先不同而造成的样本遗传变异的连续方差轴。

对于某一主成分，有相似值的个体有着相似的对应祖先。

2 .下面介绍如何使用PLINK来进行PCA分析，并去除异常个体

本教程改编自： https://www.biostars.org/p/335605/

2.1软件与硬件要求

软件版本要求：

plink > v1.9

BCFtools (tested on v 1.3)

磁盘容量要求：

downloaded data (VCF.gz and tab-indices), ~ 15.5 GB
converted BCF files and their indices, ~14 GB
binary PLINK files, ~53 GB
pruned PLINK binary files, ~ <1 Gb

2.2 下载所用数据（本教程数据主要来自1000 Genomes Phase III imputed genotypes）

2.3 下载对应每条染色体的vcf文件

prefix="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr" ;

suffix=".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz" ;

for chr in {1..22}; do
    wget "${prefix}""${chr}""${suffix}" "${prefix}""${chr}""${suffix}".tbi ;
done

2.4 下载 1000 Genomes Phase III 的 ped文件，包含了个体性别，群体等基本信息

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20130606_sample_info/20130606_g1k.ped ;

2.5 下载参考基因组（GRCh37 / hg19）注：不同版本的区别详见

wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz ;
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai ;

gunzip human_g1k_v37.fasta.gz ;

2.6 将 1000 genome 的文件转化成 BCF 格式，这一步需要对原始vcf文件进行预处理：

将multi-allelic calls分离（很多位点有3个及以上等位基因），并依据参考基因组标准化indels (详见：Variant Normalization 变异的标准化)
将ID filed设为 CHROM:POS:REF:ALT 这个格式，确保唯一性（文件中有很多重复的rsID，却指向不同的位置）x ID -I +'%CHROM:%POS:%REF:%ALT' 先抹掉现有ID，然后重新将ID设为 CHROM:POS:REF:ALT
去掉重复的变异

for chr in {1..22}; do
    bcftools norm -m-any --check-ref w -f human_g1k_v37.fasta \
      ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz | \
      bcftools annotate -x ID -I +'%CHROM:%POS:%REF:%ALT' | \
        bcftools norm -Ob --rm-dup both \
          > ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf ;

    bcftools index ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf ;
done

2.7 将BCF转化成PLINK格式

for chr in {1..22}; do
    plink --noweb \
      --bcf ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.bcf \
      --keep-allele-order \
      --vcf-idspace-to _ \
      --const-fid \
      --allow-extra-chr 0 \
      --split-x b37 no-fail \
      --make-bed \
      --out ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;
done

2.8 对每条染色体 Prune variants

mkdir Pruned ;

for chr in {1..22}; do
    plink --noweb \
      --bfile ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes \
      --maf 0.10 --indep 50 5 1.5 \
      --out Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;

    plink --noweb \
      --bfile ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes \
      --extract Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.prune.in \
      --make-bed \
      --out Pruned/ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes ;
done

2.9 准备一个索引文件

find . -name "*.bim" | grep -e "Pruned" > ForMerge.list ;

sed -i 's/.bim//g' ForMerge.list ;

2.10 合并每条染色体的plink文件

如果有自己的文件需要合并，可以在这一步加进索引文件中

plink --merge-list ForMerge.list --out Merge ;

2.11 进行主成分分析 PCA

plink --bfile Merge --pca

3 . 绘图

需要之前下载好的ped文件，20130606_g1k.ped

以及群体的三个字母的编码对照表（如下所示，来自1000genomes.org的ftp）

#Populations

This file describes the population codes where assigned to samples collected for the 1000 Genomes project. These codes are used to organise the files in the data_collections' project data directories and can also be found in column 11 of many sequence index files.

There are also two tsv files, which contain the population codes and descriptions for both the sub and super populations that were used in phase 3 of the 1000 Genomes Project:

<ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/20131219.populations.tsv>
<ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/20131219.superpopulations.tsv>

###Populations and codes

         CHB	Han Chinese             Han Chinese in Beijing, China
         JPT	Japanese                Japanese in Tokyo, Japan
         CHS	Southern Han Chinese    Han Chinese South
         CDX	Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV	Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
         CHD	Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)
	
         CEU	CEPH                    Utah residents (CEPH) with Northern and Western European ancestry 
         TSI	Tuscan                  Toscani in Italia 
         GBR	British                 British in England and Scotland 
         FIN	Finnish                 Finnish in Finland 
         IBS	Spanish                 Iberian populations in Spain 
	
         YRI	Yoruba                  Yoruba in Ibadan, Nigeria
         LWK	Luhya                   Luhya in Webuye, Kenya
         GWD	Gambian                 Gambian in Western Division, The Gambia 
         MSL	Mende                   Mende in Sierra Leone
         ESN	Esan                    Esan in Nigeria
	
         ASW	African-American SW     African Ancestry in Southwest US  
         ACB	African-Caribbean       African Caribbean in Barbados
         MXL	Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR	Puerto Rican            Puerto Rican in Puerto Rico
         CLM	Colombian               Colombian in Medellin, Colombia
         PEL	Peruvian                Peruvian in Lima, Peru

         GIH	Gujarati                Gujarati Indian in Houston, TX
         PJL	Punjabi                 Punjabi in Lahore, Pakistan
         BEB	Bengali                 Bengali in Bangladesh
         STU	Sri Lankan              Sri Lankan Tamil in the UK
         ITU	Indian                  Indian Telugu in the UK

Should you have any queries, please contact info@1000genomes.org.

下面演示用PYTHON进行绘图：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

#读取文件
pca = pd.read_csv("1000_all.eigenvec"," ",header=None)
ped = pd.read_csv("20130606_g1k.ped","\t")
#重命名pca文件的列
pca = pca.rename(columns=dict([(1,"Individual ID")]+[(x,"PC"+str(x-1)) for x in range(2,22)]))

#join两个数据表
pcaped=pca.join(ped,on="Individual ID",how="inner")

#给个体加上super population的标签
pcaped["Superpopulation"]="unknown"
pcaped.loc[pcaped["Population"].isin(["JPT","CHB","CDX","CHS","JPT","KHV","CHD"]),"Superpopulation"]="EAS"
pcaped.loc[pcaped["Population"].isin(["BEB","GIH","ITU","PJL","STU"]),"Superpopulation"]="SAS"
pcaped.loc[pcaped["Population"].isin(["CLM","MXL","PEL","PUR"]),"Superpopulation"]="AMR"
pcaped.loc[pcaped["Population"].isin(["CEU","FIN","GBR","IBS","TSI"]),"Superpopulation"]="EUR"
pcaped.loc[pcaped["Population"].isin(["ACB","ASW","ESN","GWD","LWK","MSL","YRI"]),"Superpopulation"]="AFR"

#绘制PC1与PC2的散点图，以Super population分不同颜色
plt.figure(figsize=(15,15))
sns.scatterplot(data=pcaped,x="PC1",y="PC2",hue="Superpopulation",s=50)

#绘制PC1与PC2的散点图，以Population分不同颜色
sns.scatterplot(data=pcaped,x="PC1",y="PC2",hue="Population",s=50)

参考：

https://www.cog-genomics.org/plink/1.9/strat

Price, A. L. et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909 (2006).

https://www.biostars.org/p/335605/

GWAS前的质控（QC） GWAS Quality control

进行GWAS前，我们需要对数据进行严格的质控，以排除可能造成假阳性的因素。

GWAS的质控主要可以分为两部分：

SNP层面的质控 Marker level QC

SNP的丢失率（Missingness of SNPs）
MAF (Minor allele frequency)
哈迪温伯格平衡 Hardy– Weinberg equilibrium (HWE)

样本层面的质控 Sample-level QC

样本的丢失率 (Missingness of individuals)
性别错误 (Sex discrepancy) ：GWAS QC – 性别检查与阈值选择 Sex check
杂合性 (Heterozygosity)
关联 (Relatedness)
群体分层 (Population stratification)

SNP层面的质控：

SNP的丢失率（Missingness of SNPs）：我们需要排除call rate过低的snp，通常排除call rate <=98%的snp。
MAF (Minor allele frequency) : 对于GWAS来说，我们一般研究对象为common snps，通常排除maf<5% 或maf<1%的snp。如果是rare variants (即MAF<1%)，检验的power会不足。大样本量时可以加入rare vairants。
哈迪温伯格平衡 Hardy– Weinberg equilibrium (HWE) ：我们要排除掉严重偏离哈迪温伯格平衡的variants，这是常用的表示基因分型错误的指标，二分型性状（binary traits）对于case，排除的标准为HWE p值<1e −10，对于control ，HWE p值<1e−6。对于数量性状来说，通常以HWE p value <1e‐6为标准。

样本层面的质控 Sample-level QC

样本的丢失率 (Missingness of individuals)：与SNP的丢失率类似，我们要排除掉高基因型丢失率的个体，通常的标准同样为98%。
性别错误 (Sex discrepancy)：有时数据录入会有错误，我们应基于X染色体的杂合性，确定是否有性别录入错误的个体。通常对于男性，X染色体的纯合性（homozygosity）估计值应高于0.8，女性应低于0.2.
杂合性 (Heterozygosity)：排除掉有过高或过低杂合性的个体，过高或过低杂合性通常意味着样本污染或是近亲繁殖。通常排除掉在平均杂合率在三个标准差以外的个体。
亲缘（Relatedness）：基于所有样本两两之间的IBD来排除掉有亲缘关系的个体。通常的标准为需要排除二级亲属（second degree relatives），详见：预留链接。
群体分层（Population stratification）：要去除掉非目标群体的个体（ethinic outlier），主要通过PCA等方法，详见：群体分层与主成分分析 Population structure & PCA

参考：

Marees, A.T., de Kluiver, H., Stringer, S., Vorspan, F., Curis, E., Marie-Claire, C., and Derks, E.M. (2018). A tutorial on conducting genome-wide association studies: Quality control and statistical analysis. International Journal of Methods in Psychiatric Research 27.

标签： QC

哈迪-温伯格平衡精确检验 HWE

哈迪-温伯格平衡

哈迪-温伯格平衡精确检验检验原理

使用PLINK进行HWE检验

参考

rsID的介绍与chr:pos转换时的陷阱

什么是rsID

主要优点

rsID表示变异的类型

rsID是把双刃剑：仅凭chr：pos与rsID互相转换时的陷阱

解决办法

使用UMAP对基因组数据降维，对比PCA

SNP的rsID与位置信息的相互匹配 rsID/ chr:pos conversion

基因组控制与基因组膨胀系数λ Genomic control λGC

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

为什么在PCA或估计GRM时要去除长LD区域 Remove long-LD region

SNP的LD剪枝与聚集 LD pruning & clumping

群体分层与主成分分析教程 Population structure & PCA

GWAS前的质控（QC） GWAS Quality control