GWAS 学习笔记 | 从理论到实践
原文：A Tutorial on Conducting Genome‐wide Association Studies: Quality Control and Statistical Analysis.[1]
1. 引言
在过去的二十年中，越来越多的人开始研究遗传风险因素对人类行为的影响，进行这些研究所需的技术和分析工具也变得越来越唾手可得。
全基因组关联研究（GWAS）的目的是鉴定单核苷酸多态性（SNPs），等位基因频率随表型特征值的变化而系统地变化。鉴定这些与特征相关的 SNPs 可能揭示潜在的生物学机制。
名词解释
•Clumping ：识别并选择每个 LD block 中最显著的 SNP（即 p 值最低）以进行进一步分析。这样可以减少 SNP 之间的相关性，同时保留具有最强统计证据的 SNP。•Co-heritability ：疾病之间遗传关系的量度。 基于 SNP 的共遗传性是 SNP 对疾病效应之间协方差的比例。•Gene ：DNA 中编码分子（例如，蛋白质）的核苷酸序列。•Heterozygosity 杂合性：对于特定 SNP 的两种不同等位基因的携带。个体的杂合率是杂合基因型的比例。个体内高水平的杂合性可能表明样品质量低，而低水平的杂合性可能是近亲繁殖所致。•Individual-level missingness 个体级别的缺失：特定个体缺少的 SNP 数量。较高的缺失表明 DNA 可能质量较差或存在技术问题。•Linkage disequilibrium（LD）连锁不平衡：给定种群中同一染色体上不同基因座等位基因之间非随机关联的一种度量。当其等位基因的关联频率高于随机分类下的预期频率时，SNP 位于 LD 中。LD 涉及 SNP 之间的模式。•Minor allele frequency （MAF）次要等位基因频率：这是在特定位置上出现频率最低的等位基因的频率。大多数研究的 power 不足以检测表型与 MAF 低的 SNP 的关联，因此需要过滤这些 SNP。•Population stratification 群体分层：研究中是否存在多个亚人群（例如，具有不同种族背景的个人）。由于等位基因频率在亚群之间可能不同，因此群体分层可能导致假阳性关联和/或掩盖真实关联。筷子基因就是一个很好的例子，由于群体分层的现象而导致得到 SNP 可以用来解释用筷子吃饭的习惯的结论[2]。•Pruning ：选择近似连锁平衡的标记子集的方法。在 PLINK 中，此方法使用染色体特定窗口（区域）内 SNP 之间的 LD 强度，并根据用户指定的 LD 阈值仅选择不相关的SNP。与 clumping 相比，这种方法过滤时不考虑 SNP 的 p 值。•Relatedness ：个体之间在遗传上有多强的关联性。常规的 GWAS 研究假定所有受试者都是无关的（即，没有任何一个个体比二级亲属更接近）。若数据集中包括亲属关系，不进行适当校正的话，可能对 SNP 效应大小的标准误的估计导致一定偏差。目前已有用于分析家族数据的特定工具。•Sex discrepancy 性别差异：表型信息中的性别与根据基因型确定的性别之间的差异。这个差异可以验证实验室中的样品是否混淆。注意，只有在计算了性染色体（X和Y）上的 SNP 时，才能进行此测试。•Single nucleotide polymorphism （SNP）单核苷酸多态性：在基因组中特定位置发生的单核苷酸（即A，C，G或T）变异。SNP 通常以两种不同的形式存在（例如，A与T）。这些不同的形式称为等位基因。包含两个等位基因的 SNP 有三种不同的基因型（例如，AA，AT和TT）。•SNP‐heritability：这是分析中一定集合内 SNP 解释的性状的表型变异分数。•SNP‐level missingness：这是样本中缺少特定 SNP信息的个体数量。具有高度缺失的SNP 可能导致误差。•摘要统计数据：这些是进行 GWAS 后得到的结果，包括染色体位置， SNP 的位置，SNP（rs）标识符，MAF，效应大小（比值/β），标准误 和 p 值的信息。GWAS的摘要统计信息通常可以在研究人员之间免费访问或共享。•哈迪-温伯格定律：指在理想状态下，各等位基因的频率在遗传中是稳定不变的，即保持着基因平衡。该定律运用在生物学、生态学、遗传学。条件：①种群足够大；②种群个体间随机交配；③没有突变；④没有选择；⑤没有迁移；⑥没有遗传漂变。[3]违反 HWE 法则表明基因型频率与预期值显着不同（例如，如果等位基因A 的频率= 0.20，等位基因T的频率= 0.80；基因型AT的预期频率为2 * 0.2 * 0.8 = 0.32），并且观察到的频率不应有显着差异。在 GWAS 中，通常假设与 HWE 的差异是基因分型错误的结果。病例中的 HWE 阈值通常不如对照组中的阈值严格，因为在病例中违反HWE法则可表明。
迄今为止，GWAS 已成功地揭示了许多与精神病风险有关的 SNPs ，包括精神分裂症，自闭症谱系障碍等等。这些结果表明，精神病性状受许多常见和罕见的 SNP 的影响。同时，HapMap 项目和千人基因组计划也推动了人们对人类基因组遗传结构的不断了解，促进了 GWAS 的发展。
由于 GWAS 结果显示单个 SNP 的作用量很小，因此精神病学领域的研究人员开始对多 SNP 的共同作用效应产生了兴趣。多基因风险评分（ polygenic risk score，PRS）分析将多个 SNP 的效应大小合并为一个可用于预测疾病风险的分数。根据每个人所携带的风险突变数量可计算相应 PRS 分值，并由一个独立的大规模 GWAS 得出的 SNP 效应大小进行加权。这样，该分数可以说明每个个体对于特定表型的总遗传风险，可用于临床预测或筛查。对于精神病特征，PRS 也与 case-control 状态显著相关。然而也有研究表明，它判别准确性还不足以应用于临床。
即使最近的 GWAS 研究已鉴定出许多与表型性状显著相关的 SNP，但社会科学家和临床医生对遗传学领域的贡献仍不断增进我们对已识别与特定行为，认知或神经相关的风险 SNP 的进一步理解。但这个前提是，我们需要对遗传数据执行严格的 QC 和统计分析，以避免由于多种潜在的混淆因素（例如群体分层）而导致的假阳性结果。此外，我们需要对遗传力计算有相当的了解，才能避免进行 power 不足的研究。有关如何进行 power 分析的更多信息，可参阅本组的另一篇教程[4]。
教程大纲如下：
•对数据进行质量控制（QC），并使用适当的方法排除种族异质性。•控制潜在混杂因素的同时，对 SNP 与表型性状进行关联性检验。•进行 PRS 分析。
教程所用到的 R 和 Unix 脚本以及示例数据下载：https://github.com/MareesAT/GWA_tutorial/ 。
关于基因型填充的内容未包括在本文的范围之内，可参见 van Leeuwen 及其同事 2015 年的文章[5]。
2. 软件
主要用到 PLINK 1.07 进行 QC 和统计分析，可从 http://zzz.bwh.harvard.edu/plink 下载。PLINK 1.9 beta 包含相同的参数功能，但速度更快 https://www.cog-genomics.org/plink/1.9/ ，某些步骤会用到 R。
数据格式
PLINK 既可以读取文本格式的文件也可以读取二进制文件。读取大型文本文件时，建议使用二进制文件。
文本格式的 PLINK 数据包含两个文件：
•包含相关个体及其基因型的信息（*.ped）；•包含遗传标记的信息（*.map）。
二进制 PLINK 数据由三个文件组成：
•包含样本标识符（ID）和基因型（*.bed）的二进制文件；•包含相关个体信息文件（*.fam）；•遗传标记文件（ *.bim）。
例如，我们要进行一项躁郁症研究，*.bed 文件将包含所有患者和健康对照组的基因分型结果；*.fam 文件将包含与受试者相关的数据（与研究中其他参与者的家庭关系，性别和临床诊断信息）；*.bim 文件将包含 SNP 物理位置的信息。使用协变量进行分析通常需要第四个文件，其中包含每个个体的协变量值。
' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E "null")
图1 各种常用 PLINK 文件概述
基本的 PLINK 命令
图2 PLINK 命令行的结构
•‐‐file {your_file} 指定文本文件•‐‐bfile {your_file} 指定为二进制文件。
然后，可以添加所有其他的必需选项，例如，‐‐assoc 选项以执行关联分析，如图2 所示；此特定选项将告诉 PLINK 对目标表型的每个 SNP 进行卡方检验。
•‐‐out {outfile} 为输出文件前缀（后缀将根据需要由PLINK添加）。
更多 PLINK 选项参见 http://zzz.bwh.harvard.edu/plink/ 。
3. 遗传数据的质控
QC 是任何 GWAS 研究都应进行的重要步骤。如果不进行全面的质量控制，GWAS 将无法生成可靠的结果。造成这些错误的原因可能有很多种，例如，DNA 样品质量差，DNA 与芯片的杂交差，基因型探针的性能差以及样品混合或污染等原因。
3.1 使用 HapMap 数据产生模拟数据集
为了用真实的数据来说明所有分析步骤，我们使用 HapMap 项目（ http://hapmap.ncbi.nlm.nih ）的公开数据集得到一个模拟数据集（
N
= 207）。在本教程中，为了创建种族同质的数据集，模拟数据中仅包含北欧和西欧（CEU）的犹他州居民。由于数据的样本量相对较小，因此模拟数据中的遗传效应量会比复杂性状的遗传研究中通常观察到的值更大。注意，真正检测复杂性状的遗传危险因素需要更大的样本量。
3.2 质控步骤概述
1：SNP 和个体水平的缺失
•--geno：排除大部分受试者中缺失的 SNP。•‐‐mind：排除基因型缺失率高的个体。
我们建议首先基于宽松的阈值（0.2;> 20％）过滤 SNP 和个体，过滤掉丢失程度很高的 SNP 和个体。然后再应用更严格的过滤参数（0.02）。注意，应先执行 SNP 过滤再进行个体水平过滤。
# R 和 Unix 脚本以及示例数据下载：https://github.com/MareesAT/GWA_tutorial/# 解压 1_QC_GWAS.zip 并进入目录# 检查每个人和每个 SNP 的缺失情况，并绘制直方图plink --bfile HapMap_3_r3_1 --missing# output: plink.imiss and plink.lmiss，这些文件分别为每个人丢失的 SNP 的比例和每个 SNP 丢失的比例。# 可视化缺失结果Rscript --no-save hist_miss.R# 删除缺少高度缺失的 SNP 和样本# 在示例数据中，以下两个 QC 命令不会删除任何 SNP 和样本。# 建议从非严格阈值开始进行质量控制# 删除缺失比例 >0.2 的SNPplink --bfile HapMap_3_r3_1 --geno 0.2 --make-bed --out HapMap_3_r3_2# 删除缺失比例 >0.2 的样本plink --bfile HapMap_3_r3_2 --mind 0.2 --make-bed --out HapMap_3_r3_3# 删除缺失比例 >0.02 的SNPplink --bfile HapMap_3_r3_3 --geno 0.02 --make-bed --out HapMap_3_r3_4# 删除缺失比例 >0.02 的样本plink --bfile HapMap_3_r3_4 --mind 0.02 --make-bed --out HapMap_3_r3_5
2：性别差异
•‐‐check‐sex：根据 X 染色体杂合/纯合率检查数据集中记录的性别与真实性别之间的差异。
如果许多样本存在此类差异，则应仔细检查数据。男性 X 染色体纯合度估计值 > 0.8，女性 <0.2 。
# 检查性别差异# 女性的受试者的 F 值小于 0.2# 不满足这些要求的个体会被 PLINK 标记为 "PROBLEM"plink --bfile HapMap_3_r3_5 --check-sex# 可视化性别检查结果Rscript --no-save gender_check.R#这些检查表明，有一名女性存在性别差异，F值# 以下两个脚本可用于处理具有性别差异的个人# 注意，请使用以下两个选项之一来生成 bfile hapmap_r2# 1) 删除有性别差异的样本grep "PROBLEM" plink.sexcheck| awk '{print$1,$2}'> sex_discrepancy.txt# 此命令生成状态为 PROBLEM 的个人列表。plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMa# 此命令删除状态为 PROBLEM 的个人列表# 2) impute-sex.# plink --bfile HapMap_3_r3_5 --impute-sex --make-bed --out HapMap_3_r3_6# 这会将基于基因型信息的性别推算到您的数据集中
3：次要等位基因频率（MAF）
•‐‐maf：仅输出高于 MAF 阈值的 SNP。
MAF 低的 SNP 很罕见，因此缺乏检测 SNP 表型关联的 power 。这些 SNP 也更容易出现基因分型错误。MAF 阈值应取决于样本量，较大的样本可以使用较低的 MAF 阈值。对于大样本（N = 100,000）与中等样本（N = 10000），通常将 0.01 和 0.05 用作 MAF 阈值。
# 生成仅包含常染色体 SNP 的 bfile 文件，并删除较低次要等位基因频率（MAF）的SNP# 仅选择常染色体SNP（即从 1 号到 22 号染色体）awk '{ if ($1 >= 1 && $1 <= 22) print $2 }' HapMap_3_r3_6.bim > snp_1_22.txtplink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7# 生成 MAF 分布图。plink --bfile HapMap_3_r3_7 --freq --out MAF_checkRscript --no-save MAF_check.R# 删除 MAF 频率较低的 SNPplink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8# 剩下 1073226 SNP#根据样本量的不同，常规 GWAS 的常规 MAF 阈值在 0.01 或 0.05 之间
4：哈迪-温伯格平衡（HWE）
•‐‐hwe：过滤偏离 Hardy-Weinberg 平衡的标记。
这是基因分型错误的常见指标，它也可能指示进化选择。对于二元性状，建议：HWE p 值 <1e-10（病例）和 <1e-6（对照）。不太严格的病例阈值可避免丢弃选择中的疾病相关 SNP 。对于数量性状，建议 HWE p 值<1e-6。
# 删除不符合 Hardy-Weinberg 平衡（HWE）的SNP# 查看所有 SNP 的 HWE p值分布plink --bfile HapMap_3_r3_8 --hardy# 筛选 HWE p值低于 0.00001 的 SNP，作为下一步 Rscript 的输入，这可以放大可视化严重偏离的 SNPawk '{ if ($9 <0.00001) print $0 }' plink.hwe>plinkzoomhwe.hweRscript --no-save hwe.R# 默认情况下，plink 中的 --hwe 选项仅过滤 control# 因此，我们使用两个步骤，首先对 control 使用严格的 HWE 阈值，然后对 case 数据使用较不严格的阈值plink --bfile HapMap_3_r3_8 --hwe 1e-6 --make-bed --out HapMap_hwe_filter_step1# case 的 HWE 阈值仅过滤掉与 HWE 极为不同的 SNP# 第二个 HWE 步骤仅针对 case，因为在 control 中，所有 HWE p值 <hwe 1e-6 的 SNP 已被删除plink --bfile HapMap_hwe_filter_step1 --hwe 1e-10 --hwe-all --make-bed --out HapMap_3_r3_9
5：杂合率
•排除杂合率高或低的个体。
此类偏差可能表明样品受到污染或近亲繁殖。我们建议删除样本杂合率均值偏离 ±3 SD 的个体。
# 生成受试者杂合率分布图# 删除杂合率偏离均值超过 3 sd 的个体# 对一组非高度相关的 SNP 执行杂合性检查。# 因此，要生成非高度相关的 SNP 列表，我们排除高反转区域（inversion.txt [高LD区域]），并使用 --indep-pairwise 命令修剪 SNP。# 参数“ 50 5 0.2”分别代表：窗口大小，每一步要移动窗口的 SNP 数量以及同时在所有其他 SNP 上进行回归的 SNP 的多重相关系数plink --bfile HapMap_3_r3_9 --exclude inversion.txt --range --indep-pairwise 50 5 0.2 --out indepSNP# 注意，请勿删除文件 indepSNP.prune.in，我们将在后续步骤中用到这个文件plink --bfile HapMap_3_r3_9 --extract indepSNP.prune.in --het --out R_check# 此文件包含您修剪后的数据集# 杂合率分布图Rscript --no-save check_heterozygosity_rate.R# 以下代码生成了杂合率偏离均值超过 3 sd 的个体样本列表# 对于数据处理，我们建议使用 UNIX。但是，当执行统计计算时，R 可能更方便，因此在此步骤中使用Rscript：Rscript --no-save heterozygosity_outliers_list.R# 上面命令的输出：fail-het-qc.txt# 当使用我们的示例数据时，此列表包含 2 个个体（即，有两个个体的杂合率偏离均值超过3个SD）# 通过删除文件中的所有引号并仅选择前两列，使该文件与 PLINK 兼容sed 's/"// g' fail-het-qc.txt | awk '{print$1, $2}'> het_fail_ind.txt# 删除杂合率异常值plink --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10
6：相关性
•‐‐genome：计算两两样本之间的血缘同一性。使用独立的 SNP（pruning）进行此分析，并将其限制为仅常染色体。•‐‐min：设置阈值并输出关联性高于所选阈值的样本列表。可以检测到与 pihat> 0.2 相关的个体（即，二级亲戚）。
这些相关性可能会干扰关联分析。如果你有一个基于家庭的数据集（例如，父母后代），则不需要删除相关对，但统计分析应考虑家庭之间的相关性。但是，对于基于人口的样本，我们建议使用 pihat 阈值为 0.2。
# 检查您分析的数据集个体之间的关联性# 假设样本是随机的，我们将排除本教程中 pihat 阈值 0.2 以上的所有个体# 检查 pihat> 0.2 的样本之间的关系plink --bfile HapMap_3_r3_10 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2# 已知 HapMap 数据集包含父子关系# 以下命令将使用 z 值专门显示这些父子关系awk '{ if ($8 >0.9) print $0 }' pihat_min0.2.genome>zoom_pihat.genome# 可视化评估关系类型Rscript --no-save Relatedness.R# 生成的图在 Hapmap 数据中显示了大量的相关个体（解释图； PO = 亲子后代，UN = 无关个体），这也在意料之中，因为这个数据集是就是这样构造的。# 通常，应使用特定的基于家庭的方法来分析基于家庭的数据。在本教程中，出于教程的目的，我们将随机总体样本中的相关性视为神秘的相关性。# 在本教程中，我们将从数据集中删除所有“关联性”。# 为了证明大多数相关性是由于父母后代所致，我们仅包括创建者（数据集中没有父母的个体）。plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11# 现在，我们将再次寻找 pihat> 0.2 的个体plink --bfile HapMap_3_r3_11 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2_in_founders# 文件'pihat_min0.2_in_founders.genome'显示，在排除所有non-founders 之后，HapMap 数据中仅剩下 1 对 pihat 大于 0.2 的两个人# 根据 Z 值，他们可能是全同胞或双合子双胞胎。值得注意的是，在 HapMap 数据中他们的家庭身份（FID）并不相同# 对于 pihat> 0.2 的每对“相关”的样本，我们建议删除 call rate 最低的那个人plink --bfile HapMap_3_r3_11 --missing# 使用 UNIX 文本编辑器（例如vim）检查“相关对”中哪个 call rate 最高。# 生成 Pihat 大于 0.2 的样本的 FID 和 IID 列表，以检查谁的 call rate 较低# 在我们的数据集中，个体 13291 NA07045 的 call rate 较低vi 0.2_low_call_rate_pihat.txti13291 NA07045# 在键盘上按esc！:x# 按键盘上的 Enter# 如果有多个“相关”配对样本，可以使用相同的方法扩展上面生成的列表# 删除 pihat> 0.2 的“相关”配对样本中 call rate 最低的个体plink --bfile HapMap_3_r3_11 --remove 0.2_low_call_rate_pihat.txt --make-bed --out HapMap_3_r3_12
7：群体分层
•‐‐genome：计算两两样本之间的血缘同一性。使用独立的 SNP（pruning）进行此分析，并将其限制为仅常染色体。•‐‐cluster ‐‐mds‐plot k：根据 IBS ，生成数据的 k 维子结构。k 是维度（通常为10）。这是 QC 的重要步骤，由多个程序组成。此步骤将在“控制群体分层”一节中详细描述。
通常，如果样本包含多个种族群体（例如，非洲人，亚洲人和欧洲人），建议分别对每个种族群体进行关联测试，并使用适当的方法，例如荟萃分析，以合并结果。如果你的样本仅来自单个种族，则可以通过下面介绍的方法校正群体分层。
## 下载千人基因组数据 ### 这个来自千人基因组的文件包含来自不同种族背景的 629 个人的遗传数据# 注意，此文件非常大（> 60 GB）wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz# 建议使用 aspera 下载: ascp -QT -l 300m -P33001 -L- -k1 -i asperaweb_id_dsa.openssh fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz ./# 将 vcf 转换为 Plink 格式plink --vcf ALL.2of4intersection.20100804.genotypes.vcf.gz --make-bed --out ALL.2of4intersection.20100804.genotypes# 值得注意的是，文件 “ALL.2of4intersection.20100804.genotypes.bim” 包含不带 rs 标识符的 SNP，这些 SNP 用 . 表示。这也可以在文件 “ALL.2of4intersection.20100804.genotypes.vcf.gz” 中观察到。要查看该文件，请使用以下命令：zmore ALL.2of4intersection.20100804.genotypes.vcf.gz# 千人基因组数据中缺少的 rs 标识符对于本教程来说不是什么问题# 但是，为了便于实践，我们将为缺少 rs-identifier 的 SNP 分配唯一的标识符（即，带有 . 的SNP）plink --bfile ALL.2of4intersection.20100804.genotypes --set-missing-var-ids @:#[b37]\$1,\$2 --make-bed --out ALL.2of4intersection.20100804.genotypes_no_missing_IDs# 缺失的标识符替换成如下格式：1:11508[b37]A,G## 对千人基因组数据进行 QC# 根据缺失的基因型数据删除突变plink --bfile ALL.2of4intersection.20100804.genotypes_no_missing_IDs --geno 0.2 --allow-no-sex --make-bed --out 1kG_MDS# 根据缺失的基因型数据删除个体plink --bfile 1kG_MDS --mind 0.2 --allow-no-sex --make-bed --out 1kG_MDS2# 根据缺失的基因型数据删除突变plink --bfile 1kG_MDS2 --geno 0.02 --allow-no-sex --make-bed --out 1kG_MDS3# 根据缺失的基因型数据删除个体plink --bfile 1kG_MDS3 --mind 0.02 --allow-no-sex --make-bed --out 1kG_MDS4# 基于 MAF 删除突变plink --bfile 1kG_MDS4 --maf 0.05 --allow-no-sex --make-bed --out 1kG_MDS5# 从千人基因组数据集中提取 HapMap 数据集中存在的突变awk '{print$2}' HapMap_3_r3_12.bim > HapMap_SNPs.txtplink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6# 从 HapMap 数据集中提取千人基因组数据集中存在的突变awk '{print$2}' 1kG_MDS6.bim > 1kG_MDS6_SNPs.txtplink --bfile HapMap_3_r3_12 --extract 1kG_MDS6_SNPs.txt --recode --make-bed --out HapMap_MDS# 现在，两个数据集包含完全相同的突变信息## 数据集必须是相同的版本。更改千人基因组数据版本awk '{print$2,$4}' HapMap_MDS.map > buildhapmap.txt# buildhapmap.txt 每行包含一个 SNP-id 和物理位置plink --bfile 1kG_MDS6 --update-map buildhapmap.txt --make-bed --out 1kG_MDS7# 1kG_MDS7 和 HapMap_MDS 现在具有相同的版本## 合并 HapMap 和千人基因组数据集# 在将千人基因组数据与 HapMap 数据合并之前，我们要确保文件是可合并的，为此，我们要执行 3 个步骤：#1）确保 HapMap 和千人基因组数据集中所使用的参考基因组相同。#2）解决链的问题。#3）删除经过前两个步骤仍在数据集中有所不同的 SNP。# 以下步骤在命令方面可能是非常技术性的，但我们只是比较两个数据集并确保它们对应#1）设置参考基因组awk '{print$2,$5}' 1kG_MDS7.bim > 1kg_ref-list.txtplink --bfile HapMap_MDS --reference-allele 1kg_ref-list.txt --make-bed --out HapMap-adj# 1kG_MDS7 和 HapMap-adj 对所有 SNP 都具有相同的参考基因组# 此命令将生成一些警告，无法分配 A1 等位基因#2）解决链问题# 检查潜在的链的问题awk '{print$2,$5,$6}' 1kG_MDS7.bim > 1kGMDS7_tmpawk '{print$2,$5,$6}' HapMap-adj.bim > HapMap-adj_tmpsort 1kGMDS7_tmp HapMap-adj_tmp |uniq -u > all_differences.txt# 这两个文件中有 1624 处不同，其中一些可能是由于链问题所致## 翻转 SNP 以解决链问题# 打印 SNP 标识符并删除重复项awk '{print$1}' all_differences.txt | sort -u > flip_list.txt# 生成包含 812 个 SNP 的文件。这些是两个文件之间不对应的SNP# 翻转 812 个非对应的SNPplink --bfile HapMap-adj --flip flip_list.txt --reference-allele 1kg_ref-list.txt --make-bed --out corrected_hapmap# 检查翻转后仍存在问题的SNPawk '{print$2,$5,$6}' corrected_hapmap.bim > corrected_hapmap_tmpsort 1kGMDS7_tmp corrected_hapmap_tmp |uniq -u > uncorresponding_SNPs.txt# 此文件表明他们之间存在 84 个差异#3）从 HapMap 和千人基因组中删除有问题的 SNPawk '{print$1}' uncorresponding_SNPs.txt | sort -u > SNPs_for_exlusion.txt# 上面的命令生成了包含 42 个SNP的列表，这些列表在翻转和设置参考基因组后导致 HapMap 和千人基因组数据集之间出现 84 个差异# 从两个数据集中删除 42 个有问题的SNPplink --bfile corrected_hapmap --exclude SNPs_for_exlusion.txt --make-bed --out HapMap_MDS2plink --bfile 1kG_MDS7 --exclude SNPs_for_exlusion.txt --make-bed --out 1kG_MDS8# 将 HapMap 与千人基因组数据合并plink --bfile HapMap_MDS2 --bmerge 1kG_MDS8.bed 1kG_MDS8.bim 1kG_MDS8.fam --allow-no-sex --make-bed --out MDS_merge2# 注意，我们应该完全了解 HapMap 和千人基因组数据集之间的样本重叠。但对于本教程而言，这不重要。## 对由千人基因组数据锚定的HapMap-CEU数据执行MDS。# 使用一组修剪后的SNPplink --bfile MDS_merge2 --extract indepSNP.prune.in --genome --out MDS_merge2plink --bfile MDS_merge2 --read-genome MDS_merge2.genome --cluster --mds-plot 10 --out MDS_merge2### MDS-plot# 下载包含千人基因组数据集的种群信息的文件wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/20100804.ALL.panel# 文件 20100804.ALL.panel 包含千人基因组的个体的种群代码# 将人口代码转换为超级人口代码（即AFR，AMR，ASN和EUR）awk '{print$1,$1,$2}' 20100804.ALL.panel > race_1kG.txtsed 's/JPT/ASN/g' race_1kG.txt>race_1kG2.txtsed 's/ASW/AFR/g' race_1kG2.txt>race_1kG3.txtsed 's/CEU/EUR/g' race_1kG3.txt>race_1kG4.txtsed 's/CHB/ASN/g' race_1kG4.txt>race_1kG5.txtsed 's/CHD/ASN/g' race_1kG5.txt>race_1kG6.txtsed 's/YRI/AFR/g' race_1kG6.txt>race_1kG7.txtsed 's/LWK/AFR/g' race_1kG7.txt>race_1kG8.txtsed 's/TSI/EUR/g' race_1kG8.txt>race_1kG9.txtsed 's/MXL/AMR/g' race_1kG9.txt>race_1kG10.txtsed 's/GBR/EUR/g' race_1kG10.txt>race_1kG11.txtsed 's/FIN/EUR/g' race_1kG11.txt>race_1kG12.txtsed 's/CHS/ASN/g' race_1kG12.txt>race_1kG13.txtsed 's/PUR/AMR/g' race_1kG13.txt>race_1kG14.txt# 创建您自己的数据的 racefileawk '{print$1,$2,"OWN"}' HapMap_MDS.fam>racefile_own.txt# 连接 racefilecat race_1kG14.txt racefile_own.txt | sed -e '1i\FID IID race' > racefile.txt# 生成群体分层图Rscript MDS_merged.R# 输出文件 MDS.pdf 证明我们的“own”数据属于千人基因组数据的欧洲组。因此，我们不必删除该个体。# 但是，出于教程目的，我们在下面提供了脚本来过滤出人口分层离群值。请执行以下脚本，以便为下一个步骤生成适当的文件。## 排除种族离群值# 在HapMap数据中选择低于阈值的个人。cut-off 不是固定的阈值，而是必须根据前两个维度的可视化确定。为了排除种族离群值，需要围绕关注的人群群设置阈值。awk '{ if ($4 <-0.04 && $5 >0.03) print $1,$2 }' MDS_merge2.mds > EUR_MDS_merge2# 在 HapMap 数据中提取这些人plink --bfile HapMap_3_r3_12 --keep EUR_MDS_merge2 --make-bed --out HapMap_3_r3_13# 注意，由于我们的 HapMap 数据确实没有包含任何种族离群值，因此在此步骤中未删除任何个体。但是，如果我们的数据将包括我们设定的阈值以外的个人，则这些个人将被删除。## 基于MDS创建协变量# 仅对 HapMap 数据执行 MDS，而不会出现种族异常。随后，在第三篇教程的关联分析中，将 10 个 MDS 维的值用作协变量。plink --bfile HapMap_3_r3_13 --extract indepSNP.prune.in --genome --out HapMap_3_r3_13plink --bfile HapMap_3_r3_13 --read-genome HapMap_3_r3_13.genome --cluster --mds-plot 10 --out HapMap_3_r3_13_mds# 将 .mds 文件的格式更改为 plink 协变量文件awk '{print$1, $2, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13}' HapMap_3_r3_13_mds.mds > covar_mds.txt#在下一部分教程中，我们将执行全基因组关联分析，将 covar_mds.txt 中的值用作协变量，以调整剩余的群体分层。
4. 控制群体分层
如先前所述，群体分层是 GWAS 系统误差的重要来源。研究表明，在单个种族中也可能存在一些群体分层的现象。因此，检测和控制群体分层是至关重要的 QC 步骤。
有多种方法可以校正群体分层。在本教程中，我们将介绍 PLINK 中包含的一种方法：多维缩放（MDS）。该方法计算个体之间共享的全基因组等位基因平均比例，以生成每个个体遗传变异的数量指标（组分）。我们可以绘制各个组成部分的分数来探索这些样本是否存在在遗传上彼此相较预期更相似。例如，在一项包括来自亚洲和欧洲的受试者的研究中，MDS 分析表明，亚洲人彼此之间在遗传上比欧洲人更相似。
图3 举例说明了这种分析。基于 MDS 分析的异常个体应在进一步的分析中剔除。在排除了这些个体之后，必须进行新的 MDS 分析，并且将其主成分用作关联分析的协变量，以校正群体内的群体分层。需要包括多少成分取决于人口结构和样本量，但是精神遗传学界普遍接受最多包含 10 个成分。
图3 1KG 相对于 own 数据的 MDS 图
5. 二元/数量性状的关联统计检验
经过 QC 和 MDS 成分计算后，数据已准备就绪，可用于后续关联检验。根据预期的性状或疾病的遗传模型以及所研究的表型性状的性质，可以选择适当的统计检验。在随附的教程中，我们提供了适用于二元特征（例如，酒精依赖患者与健康对照）或数量性状（例如，每周消耗的酒精饮料数量）进行关联分析的脚本。
PLINK 提供了一种自由度（1 df）的等位基因检验，其中特征值，或二元性状的对数，根据风险等位基因数目（ minor allele [a] vs. major allele [A] ）的函数线性增加或减少。此外，还提供非加性检验，例如基因型关联检验（ 2 df: aa vs. Aa vs AA ），显性基因作用检验（ 1 df: [aa & Aa] vs AA ）和隐性检验基因作用检验（ 1 df: aa vs [Aa & AA] ）。但是，非加性检验并未得到广泛应用，因为在实践中检测非加性的统计 power 很低。通过在 PLINK 中使用基于R的插件，可进行一些更复杂的分析（例如 Cox 回归分析和治愈模型）。
5.1 二元性状计算
在 PLINK 中，SNP 与二元性状的关联（ 1 = 不受影响，2 = 影响; 0 和 -9 表示缺失；这些为 PLINK 中的默认选项，可更改 ）可通过 ‐‐assoc 或 ‐‐logistic 进行分析。
•‐‐assoc 参数进行卡方关联检验，不允许引入协变量。•‐‐logistic，将执行 logistic 回归分析，该分析允许包含协变量。‐‐logistic 比 ‐‐assoc 更加灵活，但也会增加计算时间。
5.2 数量性状计算
在 PLINK 中，可使用 ‐‐assoc 和 ‐‐linear 来计算 SNP 与数量性状的关联。
•当 PLINK 检测到数量性状（ 即1、2、0或其他值 ）时，‐‐assoc 将通过执行 Student's t 检验的渐近版本以比较两个均值来，不允许引入协变量。•‐‐linear 则以每个单独的 SNP 作为预测变量执行线性回归分析。与 ‐‐logistic参数类似，‐‐linear 可以引入协变量，但比 ‐‐assoc 要慢一些。
# 对于二元特征# assocplink --bfile HapMap_3_r3_13 --assoc --out assoc_results# 注意，--assoc 选项不允许校正诸如主成分（PC）/ MDS成分之类的协变量，这使其不太适合进行关联分析。# logistic# 在此逻辑分析中，我们将使用 10 个主成分作为协变量。我们使用从上一教程中计算得出的 MDS 成分：covar_mds.txtplink --bfile HapMap_3_r3_13 --covar covar_mds.txt --logistic --hide-covar --out logistic_results# 注意，我们使用 --hide-covar 参数仅输出文件中SNP的累加结果# 删除 NA 值，因为可能导致之后的作图报错。awk '!/'NA'/' logistic_results.assoc.logistic > logistic_results.assoc_2.logistic# 从这些 GWAS 分析获得的结果将在最后一步中显示。这还将显示数据集是否包含任何全基因组范围内的显著SNP。#注意，如果采用定量结果度量，则应将 --logistic 选项替换为 --linear。也可以将 --assoc 选项用于定量结果度量（如前所述，该选项不允许使用协变量）。
5.3 多重检验校正
基因分型芯片可以同时进行多达 400 万个标记的基因分型，从而产生大量的统计检验，产生多重检验负担。基因型填充也可能会进一步增加检验关联的数量。各种模拟表明，无论研究的实际 SNP 密度如何，用于我们广泛使用的全基因组显著性阈值 5×10-8 足以控制整个基因组中独立 SNP 的数量。在测试非洲人口时，由于这些个体之间更大的遗传多样性（可能接近 1.0×10-8 ），因此需要更严格的阈值。
Bonferroni 校正，Benjamini-Hochberg 错误发现率（FDR）和置换检验是三种广泛使用的确定全基因组显著性的方法。Bonferroni 校正旨在控制至少出现一个假阳性结果的可能性，它使用公式 0.05/n 计算矫正后的 p 值，其中 n 是检验的 SNP 数。然而，如前所述，由于连锁不平衡（LD），许多 SNP 是相关的，并非独立。因此，这种方法通常过于保守，导致假阴性结果的比例增加。
假设 SNP 之间是独立的，FDR 控制 FDR 值低于固定阈值的所有信号中假阳性的预期比例。该方法不如 Bonferroni 校正保守。注意，控制 FDR 并不具有统计显著性的概念；这仅是一种最小化假阳性预期比例的方法。此外，这种方法由于 SNP 自身的局限性，因此 p 值不是独立的，而也这是 FDR 方法的前提假设。PLINK提供了 ‐‐adjust 参数，可使我们轻松进行 Bonferroni 和 FDR 校正，。
最后，我们还可以用置换的方法来进行多重检验校正。我们将结果度量标签随机排列多次（比如1,000–1,000,000）次，从而有效消除了结果度量与基因型之间的任何真实关联。对于所有排列的数据集，然后进行统计检验。这提供了无关联的零假设下检验统计量和 p 值的经验分布。随后将从观察数据中获得的原始测试统计量或 p 值与 p 值的经验分布进行比较，以确定根据经验调整的 p 值。要使用此方法，需要使用两个 PLINK 参数 ‐‐assoc 和 ‐‐mperm 组合以生成两个 p 值：EMP1，经验性p 值（未校正），以及 EMP2，经多次检验校正经验性 p 值。此过程的计算量很大，尤其是在需要许多次置换的情况下。
# 多重检验# 除了传统的全基因组范围显著性阈值 5.0E-8 以外，还有多种方法可以进行多重检验，下面我们介绍几个。# adjustplink --bfile HapMap_3_r3_13 -assoc --adjust --out adjusted_assoc_results# 该文件提供 Bonferroni 矫正后的 p-value, 包括 FDR 和其他值。## Permutation# 这是一个计算密集型步骤# 减少计算时间，我们仅对 22 号染色体的 SNP 的子集执行此测试# EMP2 列为校正的 p值# 生成 SNP 的子集awk '{ if ($4 >= 21595000 && $4 <= 21605000) print $2 }' HapMap_3_r3_13.bim > subset_snp_chr_22.txt# 根据上述步骤中生成的 SNP 子集过滤 bfileplink --bfile HapMap_3_r3_13 --extract subset_snp_chr_22.txt --make-bed --out HapMap_subset_for_perm# Perform 1000000 perrmutations.plink --bfile HapMap_subset_for_perm --assoc --mperm 1000000 --out subset_1M_perm_result# 从最低到最高 p 值对数据进行排序sort -gk 4 subset_1M_perm_result.assoc.mperm > sorted_subset.txt#检查结果head sorted_subset.txt# 生成曼哈顿和QQ图#R 版本 > = 3.0.0。#如果您将 .assoc 文件或 assoc.logistic 文件的名称更改了，Rscripts 也要修改相应名称，否则这些脚本将无法运行。# 以下 Rscripts 需要 qqman R包Rscript --no-save Manhattan_plot.RRscript --no-save QQ_plot.R## 祝贺您成功完成了 GWAS 分析！！# 如果您还对学习如何进行多基因风险评分（PRS）分析感兴趣，可参考下一部分教程。# PRS 教程独立于之前的教程
6. PRS分析
6.1 计算 PRS
单一变异关联分析一直是 GWAS 中的主要方法，但这一般需要非常大的样本量才可以进行。相比之下，PRS 分析并非仅识别单个 SNP，而是将整个基因组中的遗传风险汇总到一个多基因评分中（ 简化示例请参见图4 ）。在这种方法中，需要大量的发现样本才能确定每个 SNP 的权重。随后，在样本量较小的独立目标样品中，基于遗传 DNA 图谱和这些权重计算多基因得分。根据经验，需要大约 2,000 名受试者的目标样本才有足够的 power 来识别所解释的很大一部分方差。此外，发现数据集和目标样本应具有相同数量的样本，直到目标样本包括 2,000 个个体为止。如果有更多样本可用，这些样本应该包括在发现样本中，以最大程度地估计效应量。虽然 PRS 不足以在个体水平上预测疾病风险，但它已成功用于显示性状内和性状之间的重要关联。例如，对精神分裂症的 PRS 分析首次发现，根据常见 SNP（来自发现样本）的影响估算出的发展为精神分裂症的遗传风险的总体量度显示出与疾病风险显著相关，同时在独立（目标）样本中也与精神分裂症风险相关。尽管可用的样本量太小无法检测出全基因组范围内的重要SNP，但通过这种分析仍发现了显著的相关性。此外，用于精神分裂症的 GWAS（发现样本）已被用于预测具有多种表型的目标样本中的遗传风险，例如双相情感障碍等等。
图4 将三个单核苷酸多态性（SNP）汇总为一个多基因风险评分（PRS）的例子。权重可以是 beta 或比值比的对数，具体取决于分析的是数量性状还是二元性状。
为了进行 PRS 分析，从发现集 GWAS 获得特定表型的权重（数量性状的 beta 值和二元性状的比值比的对数）。在目标样本中，将根据他或她携带的风险等位基因数量的加权总和乘以特定性状权重，为每个个体计算 PRS。对于许多复杂的性状，SNP 效应大小已经公开（例如，参阅 https://www.med.unc.edu/pgc/downloads）。
尽管原则上所有常见的 SNP 都可以用于 PRS 分析中，但习惯上先将 GWAS 结果 clump（ 参见 clumping ），然后再计算风险评分。p 值的阈值通常用于删除几乎没有或没有关联证据的 SNP（ 例如，仅保留 p 值<0.5 或 <0.1 的SNP 。通常，我们会执行多个 PRS 分析，选取多个阈值。
6.2 进行多基因风险预测分析
一旦为目标样品中的所有个体计算了 PRS，就可以将这些分数用于（逻辑）回归分析中，以预测任何预期与目标性状表现出遗传相似的性状。预测精度可以通过回归分析的 (pseudo‐)R2 表示。重要的是，在回归分析中至少包含一些 MDS 成分作为协变量，以校正群体分层。为了估计有多少变化可以由 PRS 所解释，可以把包括协变量（例如，MDS组分）模型的 R2 和包括协变量+PRS 的模型 R2 进行比较。若由于 PRS 而导致 R2 增加就说明遗传风险因素增加了预测的精度。
PRS 的预测准确性主要取决于所分析特征的（共）遗传性，SNP 数量和discovery 样本的大小。目标样本的大小仅影响 R2 的可靠性，如果感兴趣特征的（共）遗传性以及使用的 discovery 样本的大小足够大，即使是样本量为几千的目标样本也可以获得一个显著的 R2。PRS 计算 power 的代码参考：https://sites.google.com/site/fdudbridge/software 。
PRSice[6] 是一个进行 PRS 分析的程序。它会进行 clumping，设定 p 值阈值，MDS 成分，并绘制图表。下面的教程将介绍用 PRSice 进行 PRS 分析，参考 https://github.com/MareesAT/GWA_tutorial/ （4_PRS.doc）。其他可进行 PRS 分析的软件，例如 PLINK（‐‐score）和 LDpred 。
PRSice-2 安装
wget https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_linux.nightly.zipunzip PRSice_linux.nightly.zip# OSX 64-bit: https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_mac.nightly.zip# Windows 32-bit: https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_win32.nightly.zip# Windows 64-bit: https://github.com/choishingwan/PRSice/releases/download/2.2.11/PRSice_win64.nightly.zipunzip PRSice_linux.nightly.zip# 下载需要安装的 R 包Rscript PRSice.R --dir .
使用 PRSice 进行 PRS 分析
这里用 PRSice 自带的示例数据为例。
•二元性状
Rscript PRSice.R --dir . \--prsice ./PRSice \--base TOY_BASE_GWAS.assoc \--target TOY_TARGET_DATA \--thread 1 \--stat OR \--binary-target T
•数量性状
Rscript PRSice.R --dir . \--prsice ./PRSice \--base TOY_BASE_GWAS.assoc \--target TOY_TARGET_DATA \--thread 1 \--stat BETA \--beta \--binary-target F
base 参数是指具有基本 base 样本（也称为发现样本或训练样本）的统计信息文件。 这些汇总统计信息包含每个遗传变异的效应值和 p 值。target 参数是指 target 样本二进制 plink 格式的基因型数据前缀，也支持 BGEN 格式。 base 样本和 target 样本也称为验证样本和测试样本，target 样本应独立于 base 样本。
为简单起见，我们在此分析中未包括主成分以及其他协变量，但在进行您自己的分析时，我们强烈建议添加协变量进行分析。
解释结果
默认情况下，PRSice 会输出两个图和几个文本文件。第一个图是PRSice_BARPLOT_<date>.png （图S1）。此图显示了不同阈值得到的关联结果对应的 R2 分布。 如图 S1 所示，使用 p 值 0.4463 的模型在目标样品中的 p 值最高，为 4.7*10-18 。 但通常样本量较小的多基因风险评分分析的预测值相对较低（ R2 约为5％）。
图S1 PRSice 条形图
第二个图是 PRSice_HIGH-RES_PLOT_<date> .png（图S2）显示了许多不同的 p值阈值，以及对应 R2 的 p 值（黑点），绿线为趋势线。y 轴最高点代表最优解。
图S2 PRSice 高分辨率图
这两个图均表明，许多影响 base 样本性状的 SNP 可用于预测 target 样品中的性状。 注意，这两个性状可以相同也可以不同。 如果使用相同的性状，则预测值与该性状的遗传性（以及 base 样本的样本量）有关。 如果分析不同的性状，则预测值与两个性状之间的遗传重叠有关。无论哪种方式，多基因风险评分分析通常都表明，宽松 p 值阈值的模型通常比更严格阈值的模型有更好的预测能力，这表明许多统计学上无关紧要的 SNP 在多基因性状上具有预测价值。
文本文件则详细记录了每个 p 值的 R2。
# head PRSice.prsice# 不同阈值对应的 SNP 数量，SNP 解释度，回归系数Set Threshold R2 P Coefficient Standard.Error Num_SNPBase 0.00025005 0.0133696 0.000119198 -0.197266 0.0512675 2Base 0.00030005 0.00824473 0.00155473 -0.225204 0.0711709 3Base 0.00040005 0.0089725 0.000519956 -0.350267 0.100934 5Base 0.00045005 0.0101339 0.000174268 -0.445497 0.118683 6Base 0.00065005 0.00532975 0.00605014 -0.402003 0.146446 8Base 0.00070005 0.00876654 0.000367565 -0.549246 0.154181 9Base 0.00080005 0.00233607 0.0618813 -0.369219 0.197745 13Base 0.00085005 0.00153157 0.130151 -0.342923 0.226575 15Base 0.00095005 0.000124324 0.666585 -0.100725 0.233789 16# head PRSice.summary# 最佳阈值的分析结果Phenotype Set Threshold PRS.R2 Full.R2 Null.R2 Prevalence Coefficient Standard.Error P Num_SNP- Base 0.4463 0.0520082 0.0520082 0 - 86.288 9.96334 4.69493e-18 36759# head PRSice.best# 每个个体的分值FID IID In_Regression PRSCAS_1 CAS_1 Yes -0.00599501328CAS_2 CAS_2 Yes -0.00631017938CAS_3 CAS_3 Yes -0.00227495325CAS_4 CAS_4 Yes -0.00204360007CAS_5 CAS_5 Yes -0.000830676955CAS_6 CAS_6 Yes -0.00224943517CAS_7 CAS_7 Yes -0.000687589983CAS_8 CAS_8 Yes -0.00413102565CAS_9 CAS_9 Yes 0.00256661049
7.写在最后
本文对遗传分析背后的理论（例如 GWAS 和 PRS），必要的 QC 步骤以及软件和方法的进行了基本介绍以及实践操作。
基于 GWAS 研究，已有很多公共资源使我们受益良多。例如，从 GTEx 我们可以获得 SNP 与基因表达之间的关联信息。Ensembl 和 FUMA 是功能注释的常用工具。此外，我们还可以免费获得对所研究的精神病性状或疾病的遗传结构产生重要见解的研究方法。例如，GCTA 和 LD 评分回归分析已用于估计基于 SNP 的遗传力。正是这些基于表型性状与基因中的多个 SNP 之间的关联以及通路/基因集的分析方法，增加了我们对精神疾病生物学通路的认识。详细讨论所有可用的 GWAS 下游工具和资源超出了本文的范围。有关 GWAS 后续分析的内容，建议参考 Reed 及其同事的文章。[7]
引用链接
[1] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6001694/[2] https://www.nature.com/articles/4000662[3] https://baike.baidu.com/item/%E5%93%88%E8%BF%AA-%E6%B8%A9%E4%BC%AF%E6%A0%BC%E5%AE%9A%E5%BE%8B/6265394?fr=aladdin[4] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6878611/[5] https://www.ncbi.nlm.nih.gov/pubmed/26226460[6] http://www.prsice.info/step_by_step/[7] https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5019244/
友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
•生信技能树的2019年终总结，你的生物信息学成长宝藏•2020学习主旋律，B站74小时免费教学视频为你领路•全国巡讲全球听（买一得五），你的生物信息学入门课