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

  1. 哈迪-温伯格平衡
  2. 哈迪-温伯格平衡精确检验检验原理
  3. 使用PLINK进行HWE检验
  4. 参考

哈迪-温伯格平衡

回顾: 哈迪温伯格平衡 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时的问题:
  1. 对应位点rsID不存在,可能是新变异等等原因,通常可以以chr:pos:ref:alt的形式替代。但还有个问题就是Alt allele不存在,比如rs123456 对应chr1:123456的 T>C,A 而你手里的数据是chr1:123456的 T>G, 那问题来了,这应不应该给他们相同的rsID?仅凭位点和类型来说应该给,或许下个版本的dbsnp会加上这个变异,但其实我也没有明确的答案(欢迎评论区讨论),不过实际操作中我会倾向于保守一点,用chr:pos:ref:alt 而不是rsID来表示。
  2. 如上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的唯一性特点,是不是有点违背初衷,本末倒置了。
  3. 还有一个问题就是手头数据里的变异是否已经标准化? 未标准化的变异的chrpos是不准确的,进行左对齐与节俭原则的标准化后可能产生位移,用未标准化chrpos匹配时可能会错位匹配到其他相邻的位点上。比如手头的变异可能是 chr1:123456:AA:AT ,标准化后则是chr1:123457:A:T,向后移了一位,如果你看过1000genome的原始数据就会发现这样的情况大量存在,所以应当注意(参考:GWASLab:变异的标准化 Variant Normalization
  4. 0起点还是1起点的参考系问题,处理数据时应该注意,这里不做过多赘述。(GWASLab:LiftOver 基因组坐标变换 与 01坐标系统

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

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

解决办法

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

GWAS Sumstats Harmonization GWAS数据的协调统一

背景介绍

在进行meta分析之前,我们首先要对gwas的sumstats进行预处理,这一步看似简单,但却是高质量meta分析中必不可少的一步,本文将简单介绍预处理过程中需要注意的事项,并简单介绍一款实用的预处理软件。

不同GWAS sumstats通常格式千差万别,单单统一格式就已经让人烦恼,而更令人万念俱灰的是,不同sumstats的参考的基因组版本并不相同,allele的正负链也不统一,SNPID有的用rsID(rsID的版本又是个问题~_~),有的又用chrpos…

针对以上种种的不统一,我们都要逐一进行处理,要注意事项主要包括:

  1. 参考基因组的版本:hg19还是hg38 或其他
  2. 表示变异等位的正负链是否统一
  3. 如果是a1,a2这种表示形式,则需要明确effect allele与non-efffect allele
  4. 对于palindromic SNP (A/T, 或是G/C如何处理?)
  5. 表示效应量的是beta还是OR

Pipeline

预处理简要的pipeline:

  1. 确定参考所用的VCF文件, 最好包含相应群体的allele frequency(AF)数据
  2. 对于基因组版本不一致的sumstats首先liftover转换成参考vcf的基因组版本
  3. 对sumstats进行整理并统一beta与OR的使用,
    1. 数量表型用beta与SE
    2. 二分类表型用OR与OR_95L,OR_95U
  4. 对不同群体的sumstats,使用相应的AF,利用harmoniser工具分别进行harmonise
    1. 设定依据AF判断正负链的阈值
  5. 对于处理后的gwas sumstats进行统一的质控 QC
    1. INFO
    2. MAF等

Genetics-Sumstats-Harmoniser工具使用教程

下载与安装

harmonise可以使用的工具

该工具由opentargets项目组成员开发

github链接:https://github.com/opentargets/genetics-sumstat-harmoniser

# Download 从github下载
git clone <https://github.com/opentargets/sumstat_harmoniser.git>
cd sumstat_harmoniser

# Install dependencies into isolated environment 创建环境
conda env create -n sumstat_harmoniser --file environment.yaml

# Activate environment 激活环境
source activate sumstat_harmoniser

# Run tests 检验是否安装成功
cd tests
./run_tests.py

环境要求 Requirements

参考数据 Reference data

可以使用的参考数据:

以global biobank meta-analysis initiative为例,其参考数据使用了gnomAD的数据:

https://gnomad.broadinstitute.org/downloads

示例:

INFO字段中包含了不同族裔中的allele frequency的信息:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
1	10108	rs62651026	CAACCCT	C	46514.3	RF	AF_eas=0;AF_nfe=0;AF_fin=0.0227273
1	10109	rs376007522	AACCCT	A	89837.3	RF	AF_eas=0;AF_nfe=0.0512821;AF_fin=0.0588235

使用方法

输入选项

需要两个文件:一是gwas的sumstats,二是参考的vcf文件

Input files:
  --sumstats <file>     GWAS summary statistics file
  --vcf <file>          Reference VCF file. Use # as chromosome wildcard.

其中gwas sumstats还需要指定对应的列名:

Input column names:
  --chrom_col <str>     Chromosome column 染色体
  --pos_col <str>       Position column 碱基位置
  --effAl_col <str>     Effect allele column 效应等位
  --otherAl_col <str>   Other allele column 非效应等位
  --beta_col <str>      beta column
  --or_col <str>        Odds ratio column 
  --or_col_lower <str>  Odds ratio lower CI column
  --or_col_upper <str>  Odds ratio upper CI column
  --eaf_col <str>       Effect allele frequency column
  --rsid_col <str>      rsID column in the summary stat file

输出选项

输出3个文件,Harmonise后的文件,Harmonise过程中的统计文件,以及forward/reverse/palindromic variants的数量的文件,通过以下选项指定:

Output files:
  --hm_sumstats <file>  Harmonised sumstat output file (use 'gz' extension to
                        gzip)
  --hm_statfile <file>  Statistics from harmonisation process output file.
                        Should only be used in conjunction with --hm_sumstats.
  --strand_counts <file>
                        Output file showing number of variants that are
                        forward/reverse/palindromic

其他选项

Other args:
  --only_chrom <str>    Only process this chromosome
  --in_sep <str>        Input file column separator (default: tab)
  --out_sep <str>       Output file column separator (default: tab)
  --na_rep_in <str>     How NA are represented in the input file (default: "")
  --na_rep_out <str>    How to represent NA values in output (default: "")
  --chrom_map <str> [<str> ...]
                        Map summary stat chromosome names, e.g. `--chrom_map
                        23=X 24=Y`

Harmonisation mode 模式选择:

可以选择以下四种模式:[infer|forward|reverse|drop]

对于 palindromic variants 的处理方式为:

(a) infer strand from effect-allele freq, 根据AF推测

(b) assume forward strand, 假设全部都为正链

(c) assume reverse strand, 假设全部为负链

(d) drop palindromic variants 舍弃所有palindromic variants

AF字段名与阈值

如果选择通过AF推测正负链(–palin_mode infer),那么还需要制定vcf文件中af的字段名与推测的阈值 –af_vcf_field <str> VCF info field containing alt allele freq (default:AF_NFE)

指定VCF文件中INFO里包含AF信息的字段

–infer_maf_threshold <float>Max MAF that will be used to infer palindrome strand (default: 0.42)

推测palindrome strand的最大MAF的阈值

使用例

如果输入的sumstats为以下格式:

CHROMOSOME      POSITION        MARKERNAME      NEA     EA      EAF     INFO    N       BETA    SE      P       OR      OR_95L  OR_95U

那么harmonise时,可以使用如下的代码:

sumstat_harmoniser \
 --sumstats ${sumstats} \
 --vcf gnomad.vcf.gz \    #以gnomad为参考
 --eaf_col EAF \
 --af_vcf_field AF_eas \
 --chrom_col CHROMOSOME \
 --pos_col POSITION \
 --effAl_col EA \
 --otherAl_col NEA \
 --beta_col BETA \
 --or_col OR \
 --or_col_lower OR_95L \
 --or_col_upper OR_95U \
 --palin_mode infer \
 --hm_sumstats ${out}.hm \
 --hm_statfile ${out}.stats

harmonise后的结果为原文件前加上以下harmonise后的列

hm_varid hm_rsid hm_chrom hm_pos hm_other_allele hm_effect_allele hm_beta hm_OR hm_OR_lowerCI hm_OR_upperCI hm_eaf hm_code

参考

https://gnomad.broadinstitute.org/downloads

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

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

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


UMAP 简介

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

UMAP的理论基础基于流形理论(manifold theory)与拓扑分析

主要基于以下假设:

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

总体来看,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)

注意:

  1. UMAP只接受完整的数据,对于有缺失的基因型数据要事先对缺失的SNP进行填补以得到完整的数据,填补方法可参考:GWASLab:Eagle2单倍型定相工具 Haplotype phasing
  2. 对于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与染色体位置信息的匹配或转换是生物信息学研究中必不可少的,但却又是十分繁琐的一项步骤,很多同学都在纠结这个问题,本文将总结常用的转换方法,以供参考:

转换前需要考虑的问题:

  1. 变异的标准化 :Variant Normalization 变异的标准化
  2. 参考基因组的版本 :人类参考基因组 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/

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

Variant Normalization 变异的标准化

背景

VCF文件是一种灵活的,可以表示包括SNPs, indels 和CNV等多种变异类型的文件格式,但是VCF文件中变异的表示方式并不是唯一的,对同一变异的表示方式上的混淆会对后续分析造成影响。所以我们需要对这样的变异进行标准化,使其拥有唯一的(unique)表示方式。

变异标准化的定义

变异表示的标准化包含两个部分:

1,节俭原则 2,左对齐。

这两部分分别对应变异的长度与位置。


节俭原则 Parsimony

在表示一个变异(variant )时,节俭原则(parsimony )是指,在不使任何等位长度为0的情况下,用最少的核苷酸数来表示这个变异。

所以节俭原则意味着两点:

1 尽可能的用最少的核苷酸数来表示

2 但等位(allele)的长度至少为1不能为0

当一个变异的两个等位最左侧的核苷酸相同时,且删去这个相同的核苷酸后不会造成等位长度为零时,我们称这个变异在左侧是有冗余核苷酸的。

我们来看下面这个例子,前三组多核苷酸多态(Multi Nucleotide Polymorphism ,MNP)都有冗余的核苷酸,最后一组蓝色的变异的则负荷上述的节俭原则。

左对齐 Left alignment

左对齐的意思是,将变异起始位置向左移动到不能在移动为止。这是一个与插入和缺失变异相关的概念,目的是明确变异的位置(相对于明确长度的节俭原则)。

左对齐的定义如下:

当且仅当在保持其等位长度一定的情况下,无法再向左移动起始位置时,这个变异是左对齐的。

我们来看下面这个例子(一个未标准化的短重复序列):

红色:替代等位为空,不符合节俭原则

绿色:没有左对齐,但节俭

橙色:左对齐但不右节俭

蓝色:左对齐但不左节俭

紫色:左对齐,而且节俭 √

如何进行变异标准化:

目前有多种软件支持对VCF文件进行标准化,我们只用提供参考基因组以及输入文件就可以轻松的进行标准化:

例如,使用常用的bcftools进行标准化,命令如下所示:

bcftools norm -f ref.fa in.vcf -O z > out.vcf.gz

参考:

https://genome.sph.umich.edu/wiki/Variant_Normalization#Left_alignment

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 方法的异同如下所示:

根据保留主要用途
PruningLD的R2处于LD的一对SNP中MAF最高的PCA
ClumpingLD的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/

LiftOver 基因组坐标变换

本文内容:

  1. liftover 简介与网页版工具
  2. liftover 命令行工具使用方法
  3. 基于0的坐标系统与基于1的坐标系统之间的转换

人类参考基因组不同版本之间基因坐标不同,进行研究时我们需要统一基因组坐标系,从一个版本的坐标系向另一个坐标系转换的过程就称之为liftover,UCSC,ensembl等为我们提供了便捷的工具,本文以 UCSC 的liftOver工具为例:

首先UCSC liftOver工具提供了即开即用的网页版:http://genome.ucsc.edu/cgi-bin/hgLiftOver

选择目标物种,新旧基因组版本,粘贴或上传原文件就可以开始liftover。


但更多时候我们需要使用命令行的工具,以下,本文主要介绍liftOver的命令行版本,下载地址:http://hgdownload.soe.ucsc.edu/downloads.html

liftover的url: http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver

下载完成后添加到环境路径中,并通过chmod +x 添加执行权限,然后在命令行中输入liftOver验证是否安装成功:(如果成功会有如下界面)

liftOver主要语法如下:

liftOver <输入文件> <chain文件> <输出文件> <unmapped文件>

input/output 可以使用bed格式文件chain file 则需要我们根据转化前后的数据库在ucsc上下载,无法转换的条目则输出到unmapped文件里。

这里以hg19->hg38为例,进行liftover,首先下载hg19tohg38的chain文件,下载地址为: http://hgdownload.soe.ucsc.edu/downloads.html#source_downloads

点击hg19的liftover files后下载对应文件:

下载完成后不需解压即可开始使用:

liftOver input.bed hg19ToHg38.over.chain output.bed unmapped.txt

input.bed

就转换成了: output.bed


在进行LiftOver时,有一点需要我们注意,那就是文件中变异位点起始的编号,是基于0还是基于1的。

因为 UCSC 使用基于0的坐标系统,而 Ensembl 等使用基于1的坐标系统 ,不同工具切换时应该注意这一不同。

除此以外,一些文件格式是基于1的(GFF, SAM , VCF),而另一些是基于0的(BED,BAM),不同文件间转换时也需要注意。

基于0和基于1的坐标系统示可以理解为:

基于1的坐标系:对核苷酸直接编号。

基于0的坐标系:对两个核苷酸之间的间隙编号。

enter image description here

在表示单核苷酸或多个核苷酸变异时,

基于1的坐标系:直接使用变异位点的编号。

基于0的坐标系:变异两边的位置作为起止。

enter image description here

表示insert或deletion时,

基于1的坐标系:Deletion直接使用相应位点编号,insertion则是插入位置两边的核苷酸编号。

基于0的坐标系: Deletion 是插入位置两边的间隙编号表示, Insertions 则直接由插入间隙的编号表示。

基于1的坐标系 与 基于0的坐标系 互相转换时,伪代码如下:

从 基于0的坐标系 向 基于1的坐标系 转换:

if (type=SNV){start=start+1; end=end;}
if (type=DEL){start=start+1; end=end;}
if (type=INS){start=start; end=end+1;}

从 基于1的坐标系 向 基于0的坐标系 转换:

if (type=SNV){start=start-1; end=end;}
if (type=DEL){start=start-1; end=end;}
if (type=INS){start=start; end=end-1;}

参考:

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

https://genome.ucsc.edu/cgi-bin/hgLiftOver

http://hgdownload.soe.ucsc.edu/downloads.html#liftover

PLINK 1.9 / 2 基本使用方法

更新版(英文): https://cloufield.github.io/GWASTutorial/04_Data_QC/

Plink 是全基因组关联分析中最为常用的软件,其主要用途是对于原始数据的QC质控(relatedness, population structure等),数据格式转换(ped/map, bed/bim/fam, VCF, bgen等),基础数据的计算与统计(MAF,fst,grm等),以及基于线性回归和逻辑斯蒂回归的关联分析(linear/ logistic)等。 本文主要讲解Plink的基本使用方法以及常用功能。

Plink目前常用的有两个版本,plink1.9以及plink2。 Plink2主要针对日益增长的数据量,进行了多线程优化,大幅提升了处理效率。

plink链接:

plink1.9:https://www.cog-genomics.org/plink/

plink2: https://www.cog-genomics.org/plink/2.0/

首先介绍plink1.9中最为常用的数据格式:

1.ped,pedigree file。ped文件一般包含6+2N列,第一至六列分别为 1. Family ID 2. Individual ID 3.Mother ID 4.Father ID 5.Sex 6.Phenotype。第六列以后为各个SNP的等位基因,两列一组,可以使用具体的碱基,也可以使用拷贝数(0,1)。

2. map,与ped文件相伴随的文件,主要包含ped文件中SNP的位置信息。一般包含4列。分别是1. 染色体号 2.SNP ID 3.遗传图距(单位为摩根或厘摩,通常分析不需要这一列,使用哑值(dummy value) 0 填充) 4.碱基对坐标。每行一个SNP,顺序与ped文件中的SNP相对应。

因为纯文本格式占用大量储存空间,实际操作中尽量使用二进制格式,一组ped/map文件可转换成一组bed/bim/fam文件。

3.bed ,二进制格式,存储基因型,可以想象成ped文件中除去前6列,剩下基因型数据组成的矩阵

4.bim ,纯文本,存储SNP索引 (map文件 + 两列allele信息 )

5.fam ,纯文本,存储pedigree (ped文件的前六列)

6.phenotype/covariates file (optional) ,表型与协变量文件,纯文本,该文件非必须,但表型与协变量通常使用单独的纯文本文件提供(为了准备与使用上的便捷)。该文件通常包含1.family ID 2.Individual ID, 第三列及以后为表型或协变量(常用的如年龄,性别,主成分等)。


接下来简单介绍plink的命令行语法

主要可分成三部分 1.输入 2.操作 3.输出

1.输入主要是上述的ped/map 或 bed/bim/fam文件

#ped/map: –file test 如为一组文件可只输入前缀

–bfile test

plink [输入文件] [操作] [输出文件]
plink --file input_prefix --assoc --out output_prefix

plink 格式的转换

参考文献:

1.plink官方文档

2.A tutorial on conducting genome-wide association studies: Quality control and statistical analysis