表示疾病的相关名词辨别 Distinguishing synonyms of disease

疾病（一般） Illness

强调长期，持续的状态，非急性的症状

疾病（一般）Sickness

表示疾病时，与illness可互换

另外可表示恶心，想吐，晕车晕船等

可用于比喻，对某事感到失望，恶心

疾病（一般）Ailment

疾病，小恙，强调的是程度较轻

疾病（一般）Disease

相对明确的病因，病理变化或临床诊断等

障碍，失调，紊乱 Disorder

通常由身体或心理某种功能的失常导致（病因通常不明确）

例如心理障碍等

Anxiety disorder 焦虑症
Obsessive-compulsive disorder 强迫症

症状 Symptom

强调病人的主观感受

征候 Sign

（客观表现）医学征象，征候、体征、病征等

指征，适应证 Indication

指某项操作或治疗能够得以实施的特征 (适应证)

→ 禁忌证 Contraindication

注意这里中文用的是“证”而不是“症”，医师决定某项操作可否实施不仅取决于症状，还取决于其他证据，所以用证。

综合征 Syndrome

一组相关联的症状和体征共同出现的疾病或功能障碍

Metabolic syndrome

疾病（状态）Condition

泛指所有疾病，强调的是状态

外伤，损伤 Injury

外界创伤因素引起皮肉、筋骨、脏腑等组织结构破坏，及其局部和全身反应疾病的统称。

残疾，失能 Disability

由于某种机能损伤导致的，身体能力受限或缺失

病态（一般） Morbidity

描述范围广泛的术语，一般指出现症状或疾病的状态

英文： the state of being symptomatic or unhealthy for a disease or condition （Hernandez, J. B. R., & Kim, P. Y. (2022). ）

常用于 Morbidity rate （流行病学），意为患病率，或罹患率。

共病，合并症 Comorbidity

与原发疾病同时存在且相互独立的一种或多种疾病或临床状态。

并发症 Complication

一种疾病在发展过程中引起另一种疾病或症状的发生，后者即前者的并发症。

使用Token向Pypi上传Python模块包

创建Token

使用Pypi账号登陆后，在Account settings中的API tokens里，点击Add API Token，就会进入如下所示的界面

输入Token的名称，选择范围后点击Create Token

拷贝保存Token

出于安全原因Token创建后只显示一次，要及时拷贝保存。

使用时有两种方法

（1）使用时用户名为 ‘__token__’, 密码为Token

twine upload -u "__token__" -p "<API_token>"

（2）在HOME中创建如下的文件($HOME/.pypirc)，并粘贴Token

[distutils]
  index-servers =
    pypi
    PROJECT_NAME

[pypi]
  username = __token__
  password = # either a user-scoped token or a project-scoped token you want to set as the default
[PROJECT_NAME]
  repository = https://upload.pypi.org/legacy/
  username = __token__
  password = # a project token

参考

https://pypi.org/help/#apitoken

https://twine.readthedocs.io/en/stable

GWAS中的有效样本容量 GWAS Effective sample size

GWAS中的有效样本容量对于后续分析十分重要，错误或不准确的估计会影响对GWAS结果的进一步解读。这里简单介绍病例对照GWAS中有效样本容量的定义，以及与一些容易混淆的概念的区分。

病例对照GWAS中有效样本容量

病例对照GWAS中有效样本容量(Effective sample size, ESS)通常定义为在病例与对照比例平衡 (Case: Control = 1:1; 也就是病例，对照各占比50%) 的研究里，得到等效检验效能的样本量。这样定义使得不同研究间的样本量可以进行比较。

注意: 有效样本容在不同的前后文背景下有更一般化的定义，也就是在简单随机样本的情况下，能够达到与目标样本某项数值同等化精度的样本量大小，例如在问卷调查中，Markov chain Monte Carlo中，时间序列分析中等。本文只针对病例对照GWAS中的Effective sample size进行介绍。

假设v为病例所占总体的比例，N为总样本量，即病例加对照，Neff为有效样本量，对于某个SNP其效应量beta的se可以写为

则有效样本容量可以通过下式计算

在其他文献中也时常看到其等价的表示方式

N_eff = 4 * Ncase * Ncontrol / (Ncase + Ncontrol)

或者

N_eff = 4/(1/ncase + 1/ncontrol)

同时还有其他的计算方法（假设的case, control比例不同）

https://www.nature.com/articles/nprot.2014.071

N_eff = 2/(1/ncase + 1/ncontrol)

荟萃分析时有效样本量的计算

需要注意的问题 : 单个研究的有效样本量的加和不等于整体prevalence计算得到的有效样本量

两种推荐的计算方法:

简单方法：单独的研究有效样本量的简单加和
没有单独的研究有效样本量数据时可基于荟萃后的sumstats近似估计： 4/(2pq x SE^2). 其中pq为 MAF*（1-MAF）

详细推导和分析可以参考

Grotzinger, A. D., de la Fuente, J., Privé, F., Nivard, M. G., & Tucker-Drob, E. M. (2023). Pervasive downward bias in estimates of liability-scale heritability in genome-wide association study meta-analysis: a simple solution. Biological psychiatry, 93(1), 29-36.

与有效群体大小区别

有效群体大小 (Effective population size, Ne) 群体遗传学中一个重要的概念，其描述的是等效的理想化Wright–Fisher population群体的大小，决定了由genetic drift导致的群体构成发生的变化。其定义中的等效是对于某种遗传数值而言，可以指 allele variance (variance effective population size) 或 inbreeding coefficient (inbreeding effective population size)。概念上与有效样本量有类似之处，但具体所指对象不同，注意不要混淆。

估计数量性状的线性混合模型GWAS中的有效样本量

目前大多GWAS都采用了线性混合模型，其优点是允许纳入存在亲缘关系的个体，但会影响有效样本量的大小，对于这部分的讨论可以参考

Ziyatdinov, A., Kim, J., Prokopenko, D., Privé, F., Laporte, F., Loh, P. R., … & Aschard, H. (2021). Estimating the effective sample size in association studies of quantitative traits. G3, 11(6), jkab057.

参考

https://github.com/GenomicSEM/GenomicSEM/wiki/2.1-Calculating-Sum-of-Effective-Sample-Size-and-Preparing-GWAS-Summary-Statistics

https://sphweb.bumc.bu.edu/otlt/mph-modules/bs/bs704_power/bs704_power_print.html

多基因风险分数PRS临床应用的争论

背景

近期一篇发表在BMJ Medicine的有关PRS在疾病筛查，预测以及风险分层的文章引发了激烈讨论。

文章：Hingorani, A. D., Gratton, J., Finan, C., Schmidt, A. F., Patel, R., Sofat, R., … & Wald, N. J. (2023). Performance of polygenic risk scores in screening, prediction, and risk stratification: secondary analysis of data in the Polygenic Score Catalog. BMJ medicine, 2(1).

Hingorani等人的文章主要研究对象为PGS Catalog中已发表PGS的性能指标估计值(hazard ratio, odds ratio等), 并将其转化为临床检验常用的指标 DR5 （detection rate for a 5% false positive rate）进行二次分析。这篇文章得出了负面的结论，PRS在疾病筛查，预测以及风险分层等方面表现较差，对PRS的强调与其在健康体系中的实际效果不成正比。

该文章发表于2023年10月17日，随即遭到了曼彻斯特大学D Gareth Evans教授的激烈质疑，“他们关于不要使用PRS的论断就跟说你不应该用任何风险因子评估疾病风险一样。 一个有两倍PRS的女性，如果没有乳腺癌家族史，第一次怀孕过早，和过晚的初潮，其风险也只有平均值。与其他风险因素合并使用可以显著降低文中提到的假阳性以及发现有真的高风险的个体。” 同一天谢菲尔德大学的Harry Hill也对此提出了相似的质疑。几天后塔尔图大学的Padrik等人则举出乳腺癌的例子进行质疑，他们认为评估PRS对临床的潜在影响，需要在各个疾病的特殊临床背景下进行分析。

几周后本文的第一作者Hingorani对上述质疑一一进行了回应，主要举例数据并引用其他研究说明，即使结合其他风险因子，PRS对于筛查性能的提升是很小的，上述几位的质疑对该文章的主要结论没有影响。

上个月月底，Samuel A. Lambert等一众该领域知名学者（包括PGS Catalog的主要贡献者）联名发文反驳该文的观点，强调PRS既不是诊断性检查，也不是单独的风险因子，并举例说明PRS的应用与其成本效益，指出该文有缺陷的建模与不完整的理论前提。同时致信杂志编辑 “他们的结论基于了不完整的前提，没有考虑临床背景。我们认为他们的文章并没有推动PGS潜在临床应用的讨论，并且会误导临床工作者与大众” 。

目前编辑还没有回复，上面的回复原文可以在 https://bmjmedicine.bmj.com/content/2/1/e000554.responses 找到。有兴趣的同学可以看一看原文。

一些个人理解

（主要基于 Lewis, C. M., & Vassos, E. (2020). Polygenic risk scores: from research tools to clinical instruments. Genome medicine, 12(1), 1-11.）

PRS反映的是个体遗传因素相关的疾病风险，对于个体，定义上PRS是独立风险位点的风险等位拷贝数的加权之和。这个模型简单实用，但还有很多不完美之处，例如这个计算方法通常只考虑加性遗传结构，没有考虑可能的基因间相互作用或基因环境相互作用等等，反过来说这些也都是未来PRS方法研究可能的方向。

对于PRS可能的临床应用（疾病预测，风险分层等），应当考虑到PRS的性能会受到多方面影响，例如疾病的多因子性质，遗传信号测量中可能出现的不准确或错误的问题等等。对于复杂疾病，PRS单独使用时性能性能不会太高，解读时，应当认为PRS为传统风险因子预测模型的补充，而非替代。

对于PRS的理解要避免基因决定论的错误印象。如果将PRS单独用于复杂疾病预测就会类似于刻舟求剑，完全忽略了周围不断变化的因素的影响。概念上来讲，个体的遗传易感性（Genetic liability）是固定的，但其所引起的风险却是动态变化的，这个变化依赖于变化的因素例如年龄，环境暴露，家族史，个人病史等等。一个最简单的例子就是，假如一个人有很高的酗酒的遗传风险，但如果他因为一些原因从来没有见过酒，那这个遗传风险自然就无从谈起。

对于PRS的实际应用，还存在诸多需要解决的问题，但相关研究的出发点需要基于对PRS的科学理解，才能避免偏差。

参考

（批判性参考）Hingorani, A. D., Gratton, J., Finan, C., Schmidt, A. F., Patel, R., Sofat, R., … & Wald, N. J. (2023). Performance of polygenic risk scores in screening, prediction, and risk stratification: secondary analysis of data in the Polygenic Score Catalog. BMJ medicine, 2(1). https://bmjmedicine.bmj.com/content/2/1/e000554

Lewis, C. M., & Vassos, E. (2020). Polygenic risk scores: from research tools to clinical instruments. Genome medicine, 12(1), 1-11.

scDRS 单细胞疾病相关分数

基因组学与单细胞RNA测序的结合

目前复杂疾病研究方法热点之一就是多组学多方法的结合，近来多种新的结合方法中，scDRS是比较有代表性的将基因组学与单细胞RNA测序的结合的方法，类似的方法还有sc-linker等，这类方法也属于一个新的交叉领域单细胞遗传学（Single-cell Genetics）。

传统的基因组学与单细胞测序的结合的方法，例如MAGMA或是LDSC-SEG本质上还是基于细胞特异表达的基因集而进行的富集检测，而近来以scDRS为代表的方法则利用scRNA-seq的表达矩阵进一步深入至单个细胞的层面，能够或得更高的分辨率，这对解析疾病异质性，疾病关联的细胞亚群等方面可以发挥巨大作用（就是让精准医疗更精准）。

scDRS的全称是 single-cell disease relevance score, 单细胞疾病相关分数，正如其名称类似于多基因分风险分数的PRS，该方法结构上或多或少也类似PRS的构建和计算方法，不过这里的个体是一个一个的细胞。PRS用于评估个体疾病的风险，而scDRS则评估细胞是否高表达疾病相关的基因。注意这里的R，一个是risk风险（有方向），一个是relevance相关（没有方向），概念上的差异。

scDRS的方法概况

首先从GWAS结果构建疾病基因集：使用MAGMA和目标疾病的GWAS sumstats进行基因水平的关联检验，得到每个基因与疾病关联的Z分数，选定前1000个基因作为假定的疾病基因（这个步骤不是方法的重点，除了MAGMA也可以使用其他类似方法构建疾病基因集）
然后计算单个细胞的疾病分数：scDRS会对每个细胞进行计算，量化假定疾病基因的整体表达。为了最大化检验效能，会根据MAGMA所得Z分数进行加权，同时根据每个基因在单细胞测序中特异的技术性噪音进行进行加权。
最后scDRS对所有基因集和细胞，标准化其原始的疾病分数以及原始的对照分数。然后基于所有基因集和细胞的标准化后分数的经验分布，计算每个细胞的P值。
利用得到的分数可以进行多种下游分析，分析包括（1）单个细胞层面的关联性检验（2）细胞类型关联性检验以及（3）基因关联性检验等，这些分析的P值均通过MC检验（Monte Carlo 蒙特卡洛检验）获得。

用例（下游分析）

分析与表型相关的细胞类型

相比于传统方法LDSC或MAGMA，scDRS可以检验细胞类型的相关性，还可以检验同一细胞类型内的异质性。

作者列举了 22个表型与19种细胞的关联结果的热力图，Y轴为表型，X轴为细胞类型，方框表示显著相关，×表示细胞类型内的异质性，颜色深浅表示显著关联细胞的比例。

发现存在异质性的细胞亚群

作者研究了自身免疫性疾病中T细胞的异质性，11个T细胞群中（a），与IBD相关的T细胞构成了4个新的群（b）。

基因-分数的关联分析

scDRS检验基因是否与GWAS由来的基因集中的基因共表达。

相比于MAGMA，scDRS能更准确地识别出疾病相关的基因

典型的分析流程

只需要GWAS的Sumstats和单细胞RNA测序的数据，好消息是这两个都可以很容易从公开数据库中获得。

具体流程官方文档以及大牛博客已经写得很详细了，

分析代码可以参考 https://martinjzhang.github.io/scDRS/

以及 https://zhuanlan.zhihu.com/p/592128325

注意的点

scDRS的假设是基因集里的基因与疾病有关，与疾病有关的基因会在疾病相关的细胞群体（可以是疾病细胞或健康细胞）里高表达，与基因的方向无关。 scDRS并不是假设基因集里的基因会在疾病细胞里高表达。注意这里概念上细节的差别。（https://github.com/martinjzhang/scDRS/issues/42）
基本数据要求：为了能够得到足够的检验效能，GWAS的heritability z-score最好大于5,或样本量大于10万。（https://martinjzhang.github.io/scDRS/faq.html#which-gwas-and-scrna-seq-data-to-use）
Seurat格式的数据需要转换为scanpy使用的h5ad格式（表达矩阵不能有负值，Seurat的scaled.data里表达矩阵会有负值，转换时要注意）(https://github.com/martinjzhang/scDRS/issues/44)；转换可以使用 SeuratDisk https://mojaveazure.github.io/seurat-disk/articles/convert-anndata.html
为了增加检验效能，单细胞RNA测序可以事先进行 imputation （https://github.com/martinjzhang/scDRS/issues/32）
-adj-prop 可以调整细胞类别的比例 (当某些细胞种类比例过高的时候) （https://github.com/martinjzhang/scDRS/issues/32）
显著细胞的数太少：正常现象，仍然可以进行group分析（检验效能通常更高）（https://martinjzhang.github.io/scDRS/faq.html#scdrs-detected-few-significant-cells-fdr-0-2）

参考

Zhang, M. J., Hou, K., Dey, K. K., Sakaue, S., Jagadeesh, K. A., Weinand, K., … & Price, A. L. (2022). Polygenic enrichment distinguishes disease associations of individual cells in single-cell RNA-seq data. Nature genetics, 54(10), 1572-1580.

群体遗传学中种族使用上的区分 Race/Ethnicity/Ancestry

英文中的常用的表示种族的词语包括Race，Ethnicity以及Ancestry。但在中文中通常都翻译成种族。本文就这些词的使用区分做简单介绍与讨论。

种族概念的区分

首先介绍Population，这是一个最为广义的词语，可以用于表示任何一群体，可大可小。通常含义基于上下文，没有明确区分。例如，在中国人群中的全基因组关联分析就可以说成 GWAS in a Chinese population.

Race, 种族(人种)，是一个由社会构建的区分系统，但该系统基于对内在的生物学特征或差异错误的认知，典型的例子便是物理特征（诸如肤色）以及社会文化的特征。举例，种族歧视应当被消除。

Ethnicity, 种族 (民族)，是一个表示某一群体的社会政治概念，通常有相连的地理位置，基于共同的遗产或相似文化，例如语言，宗教信仰等。举例，中国有五十六个民族，这个民族就是Ethnic group，汉族Han Chinese 就是一个Ethnic group。Ethnicity与Ancestry容易混淆的点在于，多数情况下Ethnicity所表示的群体通常情况下也会有共同的家系或是遗传继承，但有一些地区Ethnicity表示的仅为社会文化实体而没有遗传学基础。

Ancestry, 种族 (族裔/祖先)，是一个更为复杂的概念，包括了生物学以及社会学的成分。在西方，这个词通常反应群体的社会文化以及所来自大陆的起源，而在东方，以及南半球，这个词通常反映家系或是遗传继承。多数情况下，ancestry是群体遗传学文章中更应当使用的词语。举例，使用频率较高的有 European ancestry， East Asian ancestry， South Asian ancestry等等。

举一个例子来综合上述概念，某研究组收集了中国人群的基因数据用于GWAS研究，那这个群体泛称就可以是一个中国人群体 a Chinese population，其中有汉族和傣族，这里的族就是ethinc group（Han Chinese 和 Chinese Dai），而整个群体在群体遗传学上则都属于East Asian Ancestry。

群体遗传学领域使用上的区分

一个核心上的区别点就在于是否主观与客观， race以及ethinitity存在主观成分，而ancestry则为客观描述性的词语，反映基因组中的某些固定特征。在生物学或遗传学文章中，单纯描述遗传学意义的种族时应使用客观性的词语，即ancestry。

群体遗传学中跨种族跨群体的英文使用

简单来说应当使用 cross-population, cross-ancestry, multi-population 或 multi-ancestry 而不是 trans-ethnic

原因

trans有多种含义，应当使用更准确且而不引起歧义的cross或者multi
ethnic包含社会学成分，存在易变的主观成分，应当使用ancestry，或更广义的population

基于一些历史原因，早期的文章常常混用，早期的文章中例如

Brown, B. C., Ye, C. J., Price, A. L., & Zaitlen, N. (2016). Transethnic genetic-correlation estimates from summary statistics. The American Journal of Human Genetics, 99(1), 76-88.

中使用了，Transethnic，但其含义应为cross-ancestry，比较合适的用例如

Momin, M. M., Shin, J., Lee, S., Truong, B., Benyamin, B., & Lee, S. H. (2023). A method for an unbiased estimate of cross-ancestry genetic correlation using individual-level data. Nature Communications, 14(1), 722.

参考

Kachuri, L., Chatterjee, N., Hirbo, J. et al. Principles and methods for transferring polygenic risk scores across global populations. Nat Rev Genet (2023). https://doi.org/10.1038/s41576-023-00637-2

Kamariza, M., Crawford, L., Jones, D., & Finucane, H. (2021). Misuse of the term ‘trans-ethnic’in genomics research. Nature Genetics, 53(11), 1520-1521.

GWAS中的赢家诅咒与其校正 Winner’s curse correction

GWAS中的赢家诅咒 Winner’s curse

GWAS中的赢家诅咒是指遗传效应的大小由于GWAS中的筛选过程（通过全基因组显著阈值筛选lead SNP）而被系统性地过高估计。

赢家诅咒本用来指代在拍卖中类似的现象。即使一件拍卖品对所有买家来说都有相同的价值（出价是无偏的），最后拍得物品的赢家很可能过高估计了拍卖偏的内在价值。类比于GWAS，lead SNP即为赢家，而它的效应量可能过高估计了真实的遗传效应。

赢家诅咒的校正 WC correction

假设观察到的 $\beta_{Observed}$ 的近似分布为:

$\beta_{Observed} \sim N(\beta_{True},\sigma^2)$

$\beta_{Observed}$ 的一个例子

$c$ : 显著性阈值对应的Z分数

上面的式子等价于

${{\beta_{Observed} - \beta_{True}}\over{\sigma}} \sim N(0,1)$

${{\beta_{Observed} - \beta_{True}}\over{\sigma}}$ 的一个例子

在通过阈值筛选的情况下， $\beta_{Observed}$ 的近似抽样分布（实际上为一个截断正态分布 truncated normal distribution）为：

$f(x,\beta_{True}) ={{1}\over{\sigma}} {{\phi({{{x - \beta_{True}}\over{\sigma}}})} \over {\Phi({{{\beta_{True}}\over{\sigma}}-c}) + \Phi({{{-\beta_{True}}\over{\sigma}}-c})}}$

其中

$|{{x}\over{\sigma}}|\geq c$

$\phi(x)$ : 标准正态分布的概率密度函数
$\Phi(x)$ : 标准正态分布的累积分布函数

从以上的近似抽样分布可以得到，筛选出来的SNP的效应量的期望分布为：

$E(\beta_{Observed}; \beta_{True}) = \beta_{True} + \sigma {{\phi({{{\beta_{True}}\over{\sigma}}-c}) - \phi({{{-\beta_{True}}\over{\sigma}}-c})} \over {\Phi({{{\beta_{True}}\over{\sigma}}-c}) + \Phi({{{-\beta_{True}}\over{\sigma}}-c})}}$

$\beta_{Observed}$ is biased.
偏差的大小由 $\beta_{True}$ , SE $\sigma$ , 以及用于筛选的显著性阈值决定.

公式推导可以参考 Ghosh, A., Zou, F., & Wright, F. A. (2008). Estimating odds ratios in genome scans: an approximate conditional likelihood approach. The American Journal of Human Genetics, 82(5), 1064-1074. 中的Appendix A

用这个式子便可以对效应量进行赢家诅咒的校正。

winnerscurse R包

可以使用这个R包进行赢家诅咒的校正。

https://amandaforde.github.io/winnerscurse/articles/winners_curse_methods.html

参考

Bazerman, M. H., & Samuelson, W. F. (1983). I won the auction but don’t want the prize. Journal of conflict resolution, 27(4), 618-634.
Göring, H. H., Terwilliger, J. D., & Blangero, J. (2001). Large upward bias in estimation of locus-specific effects from genomewide scans. The American Journal of Human Genetics, 69(6), 1357-1369.

Zhong, H., & Prentice, R. L. (2008). Bias-reduced estimators and confidence intervals for odds ratios in genome-wide association studies. Biostatistics, 9(4), 621-634.
Ghosh, A., Zou, F., & Wright, F. A. (2008). Estimating odds ratios in genome scans: an approximate conditional likelihood approach. The American Journal of Human Genetics, 82(5), 1064-1074.

Also see reference: https://amandaforde.github.io/winnerscurse/articles/winners_curse_methods.html

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

哈迪-温伯格平衡

回顾：哈迪温伯格平衡 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

GWAS检验效能 Power analysis for GWAS

第一类错误，第二类错误以及检验效能

该表列举了零假设 $H_0$ 与统计学检验结果（是否拒绝原假设 $H_0$ ）之间的关系

	H0 为真	H0 为假
不拒绝原假设	真阴性 : $1 - \alpha$	第二类错误 (伪阴性) : $\beta$
拒绝原假设	第一类错误 (伪阳性) : $\alpha$	真阳性 : $1 - \beta$

$\alpha$ : 显著性水平

根据定义，检验效能（ statistical power ）指某检验正确地拒绝零假设的概率，也就是上表中的真阳性（ True positive）。

$Power = Pr ( Reject\ | H_0\ is\ False) = 1 - \beta$

影响检验效能的因素 Factors affecting power

总的样本量 Total sample size
病例与对照的比例 Case and control ratio
变异的效应量大小 Effect size of the variant
风险等位的频率 Risk allele frequency
显著性阈值 Significance threshold

非中心参数 Non-centrality parameter

非中心参数：非中心参数（Non-centrality parameter; NCP）用于描述零假设 $H_0$ 与备择假设 $H_1$ 之间差异的程度。

考虑如下的线性模型：

$y = \mu +\beta x + \epsilon$

误差项的方差为：

$\sigma^2 = Var(y) - Var(x)\beta^2$

通常情况下单个SNP所能解释的表型的方差是极其有限的，所以我们可以近似地认为

$\sigma^2 \thickapprox Var(y)$

在哈迪温伯格平衡下，有

$Var(x) = 2f(1-f)$

$f$ : 该变异的等位频率（allele frequency）

自由度为1的 $\chi^2$ 分布的非中心参数NCP则为

$\lambda = ({{\beta}\over{SE_{\beta}}})^2$

数量表型的检验效能

$\lambda = ({{\beta}\over{SE_{\beta}}})^2 \thickapprox N \times {{Var(x)\beta^2}\over{\sigma^2}} \thickapprox N \times {{2f(1-f) \beta^2 }\over {Var(y)}}$

显著性阈值: $C = CDF_{\chi^2}^{-1}(1 - \alpha,df=1)$

$CDF_{\chi^2}^{-1}(x)$ : $\chi^2$ 分布的累积分布函数的反函数

$Power = Pr(\lambda > C ) = CDF_{\chi^2}(C, ncp = \lambda,df=1)$

$CDF_{\chi^2}(x, ncp= \lambda)$ : 非中心参数NCP为 $\lambda$ 的 $\chi^2$ 分布的累积分布函数

病例对照表型的检验效能 Power for large-scale case-control genome-wide association studies

设

$P_{case}$ : 在病例中风险等位的频率 Risk allele frequency in cases
$N_{case}$ : 病例的样本量 Number of cases. The total allele count for cases is then $2N_{case}$ .
$P_{control}$ : 在对照中风险等位的频率 Risk allele frequency in controls
$N_{control}$ : 对照的样本量 Number of control. The total allele count for control is then $2N_{control}$ .

这种情况下零假设为 : $P_{case} = P_{control}$ ，即风险等位的频率在病例中与对照中是一样的。

检验两个正态分布的比例的不同时，所用的统计量为

$z = {{P_{case} - P_{control}}\over {\sqrt{ {{P_{case}(1 - P_{case})}\over{2N_{case}}} + {{P_{control}(1 - P_{control})}\over{2N_{control}}} }}}$

显著性阈值: $C = \Phi^{-1}(1 - \alpha / 2 )$

$Power = Pr(|Z|>C) = 1 - \Phi(-C-z) + \Phi(C-z)$

计算GWAS统计效能的网页工具 GAS power calculator

GAS power calculator工具实现了上述的计算方法，可以通过网页工具，指定参数后进行计算。

GAS power calculator

示例：

参考

https://cloufield.github.io/GWASTutorial/20_power_analysis/
Skol, A. D., Scott, L. J., Abecasis, G. R., & Boehnke, M. (2006). Joint analysis is more efficient than replication-based analysis for two-stage genome-wide association studies. Nature genetics, 38(2), 209-213.
Johnson, J. L., & Abecasis, G. R. (2017). GAS Power Calculator: web-based power calculator for genetic association studies. BioRxiv, 164343.
Sham, P. C., & Purcell, S. M. (2014). Statistical power and significance testing in large-scale genetic studies. Nature Reviews Genetics, 15(5), 335-346.

参数扩展 Parameter Expansion

本文简单介绍linux命令行中的参数扩展的部分内容。在实际生物信息学分析中，参数扩展是十分实用的小技巧。可以用来快速处理大量的重命名文件，获取文件路径，获取文件类型，替换文件文件夹的路径等等

参数扩展可以用于对字符串的操作，主要用途包括以下：

字符串大小写替换
字符串替换
字符串删减
字符串切片
字符串长度

Case modification 大小写转换

${PARAMETER^} ：第一个字母大写
${PARAMETER^^} ：全部大写
${PARAMETER,} ：第一个字母小写
${PARAMETER,,} ：全部小写
${PARAMETER~} ：第一个字母大小写反转
${PARAMETER~~} ：全部大小写反转

# 转换大写
$ X="abc"
$ echo ${X^}
Abc
$ echo ${X^^}
ABC

# 转换小写
$ Y="ABC"
$ echo ${Y,}
aBC
$ echo ${Y,,}
abc

# 翻转大小写
$ echo ${X~}
Abc
$ echo ${X~~}
ABC
$ echo ${Y~~}
abc
$ echo ${Y~}
aBC

Substring replace 字符串替换

${PARAMETER/PATTERN/STRING} ：将第一个PATTERN 替换为STRING
${PARAMETER//PATTERN/STRING} ：将所有PATTERN 替换为STRING
${PARAMETER/PATTERN} ：删除第一个PATTERN （将第一个PATTERN 替换为nullstring）
${PARAMETER//PATTERN} ：删除将所有PATTERN

# 从字符串开头开始匹配
% 从字符串末尾开始匹配


MYSTRING="xxxxxxxxxx"
echo ${MYSTRING/x/y}
yxxxxxxxxx

echo ${MYSTRING//x/y}
yyyyyyyyyy

echo ${MYSTRING/#x/y} 
yxxxxxxxxx
echo ${MYSTRING/%x/y} 
xxxxxxxxxy

echo ${MYSTRING//x/}
# is equivalent to
echo ${MYSTRING//x}

# change suffix for a file 
MYSTRING="abc.txt"
echo ${MYSTRING/txt/md}
abc.md

# change prefix for a file 
echo ${MYSTRING/abc/edf}
edf.txt

Substring removal 字符串删减

${PARAMETER#PATTERN} ：删除第一个PATTERN
${PARAMETER##PATTERN} ：删除所有PATTERN
${PARAMETER%PATTERN} ：删除第一个PATTERN （从末尾开始）
${PARAMETER%%PATTERN} ：删除所有PATTERN （从末尾开始）

$ MYSTRING="This is a simple sentence with some simple words."
$ echo ${MYSTRING#*simple}
sentence with some simple words.

$ echo ${MYSTRING##*simple}
words.

$ echo ${MYSTRING%simple*}
This is a simple sentence with some

$ echo ${MYSTRING%%simple*}
This is a

$ Get name without extension
MYFILEPATH="/home/work/example.txt"
$ echo ${MYFILEPATH%.*}
/home/work/example

# Get extension
$ echo ${MYFILEPATH##*.}
txt

# Get directory name
echo ${MYFILEPATH%/*}
/home/work

# Get filename
$ echo ${MYFILEPATH##*/}
example.txt

String slicing 字符串切片

${PARAMETER:OFFSET} : OFFSET指从起始位置的偏倚
${PARAMETER:OFFSET:LENGTH} : OFFSET指从起始位置的偏倚，可以为空，LENGTH 则指偏移后截取的长度，负数则指截取到到末位前的长度。

a=myfile.txt.gz

echo $a
myfile.txt.gz

echo ${a::-3}
myfile.txt

echo ${a:2}
file.txt.gz

echo ${a:2:-3}
file.txt

String length 字符串长度

${#PARAMETER} ：获取变量所指字符串的长度

$ MYSTRING="abcd"
$ echo ${#MYSTRING} 
4

Indirection 间接指向

${!PARAMETER} `间接指向变量所指字符串所表示点变量

$ A=1
$ B="A"
$ echo $B
A
$ echo ${!B}
1

参考

https://wiki.bash-hackers.org/syntax/pe

https://www.gnu.org/software/bash/manual/html_node/Shell-Parameter-Expansion.html