统计遗传学：第八章，基因型数据质控
大家好，我是飞哥，本章节是理论+实操，干货满满，这里我将书中的数据用代码进行了实现，你可以下载相关的数据，用我整理好的代码进行操作，666！
本书其它章节介绍
第八章：基因型数据处理和质控
注意：
本书的数据，我已经下载整理好了，下载本书的电子版pdf+数据+代码，链接：书籍及配套代码领取--统计遗传分析导论
主要内容
本章节包括：
- 了解如何使用命令行
- 打开并使用PLINK二进制文件
- 将PLINK文件重新编码为其他格式
- 了解数据管理的基础，以选择特定标记或个体子样本的信息
- 获取等位基因频率、表型、，和缺失值
- 合并不同的基因文件
- 将表型与PLINK文件相关联
- 在个体、标记和全基因组关联研究水平上理解和执行质量控制程序
简介：使用遗传数据
主要应用的软件是R语言和plink软件。
上一章向读者介绍了不同类型的基因组数据，本章的目的是为那些不习惯在命令行环境中工作并且从未使用过计算机程序PLINK的人提供如何使用遗传数据的温和介绍。在下一节中，我们将更详细地描述命令行和PLINK。在概述了使用PLINK的基本知识（例如调用PLINK、打开文件和导入数据）之后，我们将描述基本的数据管理。这包括选择个体和标记以及合并不同的遗传文件。然后，我们说明了如何产生描述性统计数据，包括等位基因频率、表型和缺失值。最后，我们概述了通过个体和标记对遗传数据的质量控制，以及全基因组关联研究的质量控制，并进行了简要总结。读者可以选择在自己的计算机上积极地完成本章中的练习，我们将在下一节中介绍如何做到这一点。
PLINK软件入门
命令行
我们在第7章介绍了PLINK[1.2]，它是用于QC和GWAS分析的最流行的开源自由软件程序之一。它由Shaun Purcell及其同事开发，便于多种类型的数据处理和使用。有关如何安装PLINK的详细说明，请参阅本书附录1。还有许多在线教程（http://zzz.bwh.harvardedu/plink/tutorial.shtml，说明和plink文档格式。本章简要介绍了您将需要的基本知识和主要命令；它绝不是详尽无遗的。我们强烈建议您直接参考可用的大量在线文档和资料。具体取决于您的操作系统（如Windows、macOS、Linux/Unix），某些命令和进程可能不同。在本书中，我们使用了macOS和Linux系统的示例，Windows用户可以在附录1和方框8.1中找到更多详细信息。对于希望积极阅读本章教程的读者，我们建议您首先进行以下四个步骤：
- 第一步：新建一个文件夹，用于存储数据
- 第二步：下载安装plink软件
- 第三步：下载示例数据
- 第四步：解压文件
注意，这里我已经将需要的示例数据下载了，目录如下：
image-20220713112342136
第八章的数据如下：
image-20220713112401358
在您的计算机上安装PLINK后，您可以在命令行中键入我们在本书中显示的所有命令。命令行是直接向计算机操作系统键入命令的界面。我们需要使用命令行，因为使用PLINK需要一个等待命令的活动shell。 考虑这一点的一种方法是，命令行是与计算机程序交互的一种方式。在每个连续的文本行（即命令行）中键入一个命令，shell充当命令语言解释器或处理行的程序。shell就像一个语言解释器，它接受文本命令，然后将其转换为适当的操作系统函数。那些使用Unix或熟悉20世纪70年代和80年代MS-DOS或Apple DOS界面的人将已经习惯于在这种环境中工作。由于引入了许多具有点击、菜单驱动操作的图形用户界面（GUI），一些人可能已经失去了命令行界面的艺术。方框8.1比较了操作系统之间的一些命令行差异。 本书中介绍的ost软件使用命令行，因此没有图形用户界面。然而，在不同操作系统之间执行命令的方式存在一些差异。下面列出了Windows、Mac OS和Linux之间的一些差异。由于Mac OS和Linux都基于Unix，因此这两个操作系统之间的大多数命令没有差异。
常见的命令，在不同系统下的代码如下：
image-20220713112502484
附录1详细讨论了Windows用户的命令行。简而言之，Windows 10用户可以下载并安装Ubuntu，这将有效地让他们像在Linux操作系统中一样执行这些命令。对于Mac用户，要访问终端或Unix命令提示符，首先需要打开Utilities文件夹。要执行此操作，请单击位于屏幕底部基座的Finder lcon。在Finder窗口的左侧面板上查找应用程序，然后在Applications窗口内向下滚动，直到找到Utilities，然后单击它以打开。在右侧面板上，您将找到端子。当vou双击它时，它应该如图8.1所示。 对于那些在阅读本章的同时进行练习的人，我们建议您现在就准备好命令行提示符。命令行以光标前的提示符开始（$或>）。在大多数情况下，当前目录的提示会在提示之前显示。与R一样，使用PLINK时，当前目录很重要。默认情况下，PLINK会将数据和结果文件加载并保存到该目录中。要更改目录，请使用：cd。我们在方框8.2中概述了普林克的基本规则。
从终端打开plink
如图8.2所示，PLINK命令由几个参数组成。虽然参数的顺序是任意的，但典型的命令是从调用PLINK开始的。 然而，在执行此操作之前，您需要确保您位于存储PLINK的正确目录中。PLINK将安装在安装时指定的目录中（见附录1和方框8.1）。要查找您正在使用的当前目录，请键入；pwd（即打印工作目录）。如果PLINK不在此目录中，则需要在命令行中键入vou存储PLINK的目录行。例如，这可以是：/usr/local/bin/plink。如果需要更改目录，请使用cd命令，例如：cd/Users/yourname/plink。 在这些示例中，我们展示了在光标后运行命令。当您进行更高级的分析时，您将从脚本运行PLINK，就像在R或其他程序中一样。 如图8.2所示，PLINK是通过键入来调用的/在Mac或Linux计算机上点击并点击。exe（在旧版本的Windows计算机上）-假设PLINK位于您当前的工作目录中（请参见方框8.2）。接下来是遗传输入文件的名称，然后是我们要对这些文件执行的操作，最后是指定输出文件的名称。注意，这些命令总是以两个破折号开头
plink用户指南：
1 PLINK是一个命令行程序，因此它需要在该环境中运行（即DOS窗口或文本中定义的Unix终端）。 2 将文件放在与相同的目录中/plink（或plink.exe）fle，以便从一个目录中进行所有分析。 3 所有结果都写入具有各种不同扩展名的文件。 4 请记住始终检查日志文件，它是控制台输出的注释、警告和错误（请参阅方框8.3）。 5 PLINK没有累积内存。每次运行都会像加载新数据一样加载数据，以前的所有过滤器等都会丢失。 6 拼写和准确的命令语法极为重要。这包括您将经常使用的双破折号或减号（--）。 7 您不能将所有选项相互组合。plink并不总是警告你，但例如，基本单倍型测试不能有协变量。 8 目录路径由“斜杠”符号分隔（/在Unix和OSx中），而windows使用“反斜杠”符号（）。在Unix和macOS中，反斜杠用于打断不同的代码行。（有关操作系统之间的更多差异，请参见方框8.1。）9 9 PLINK是为人类遗传数据而设计的（我们在本书中介绍过），并假设基因组有23条染色体，其中chr23为X染色体。对于那些关注其他物种的读者，请确保您指定了这一点（例如，39条染色体的狗）。 10 值得您花时间和精力查阅在线提供的优秀文档http://www.cog-genomics.org/plink2/.
计算机可读但非人类可读。我们所有的可执行程序都存储在二进制文件中，大多以数字数据文件的形式。 我们调用这个例子——bfile，它由三个相互链接的文件组成：example.bim, example.fam, example.bed，--bfile（而不是--file）指定输入数据为二进制格式。接下来是您希望采取的选项。在我们的例子中，-freq指示PLINK计算等位基因频率。还有许多其他选项，例如用于关联分析的--recode或--assoc。请注意，在单个命令行中可以使用多个选项，并且顺序不相关，因为PLINK实现了默认顺序。 最后，您需要指定一个输出文件名，在下面的示例中是——out -example-freg。数据fle的后缀将立即添加到PLINK中。 在本例中，命令-freg将生成带有后缀的输出。如果未指定输出文件的名称，则默认名称为plink（小写）。新用户应该注意，PLINK会在没有询问的情况下覆盖具有相同名称的文件，因此在使用此命令时保持警惕很重要。推荐的做法是始终重命名输出文件以避免混淆。因此，您不使用默认值，而是遵循一致的文件命名约定。保持文件名简短且有意义，不要使用空格，如果使用不同的版本号，则以可以轻松检索最新版本的方式进行排序（2019v5example.txt或example v5 2019.txt）。使用有意义的目录名也很重要，稍后您将能够找到（并记住）这些目录名。建议对链接文件使用类似的名称，例如hapmap。ped和hapmap。地图此外，尽量避免使用诸如*：/<>|“？[]；=+&f$之类的字符，因为它们可能很难找到或打开，避免使用以后难以区分的常用词（例如版本、草稿），并尽量避免文件名和路径中不必要的重复和冗余。
在终端下运行脚本
因为我们只运行简单的命令来让读者熟悉如何使用PLINK，所以我们通常只在命令行中键入所有命令。然而，一旦你们进行了更高级的分析，你们就会想要编写你们的脚本，然后在终端中运行它。 要在终端中运行脚本，需要键入：sh/path/To/file，然后按enter键。 例如，键入：
sh hello_world.sh
操作plink文件
数据中的文件：
hapmap ceu数据PLINK可用于导入不同格式的数据。第7章已经详细讨论了最常用的格式之一是更紧凑的二进制文件，由三个同名文件和后缀组成。bed, bim和fam。我们使用的示例数据是hapmap ceu数据，您应该已经将其下载到指定的文件夹中。CEU是指具有北方和西方血统的犹他州居民，是HapMap中1l人口之一。hapmap ceu数据的文件为：hapmap ceu。
- bed包含基因型（在文本编辑器中不可读）hapmap ceu的压缩信息。
- fam包含关于个体单倍型ceu的信息
- bim包含关于SNP的信息
将二进制数据变为文本数据
将二进制数据变为文本数据，变为ped和map数据
我们在前一章中注意到，PLINK有自己的遗传数据格式，不同于Excel或Stata等统计软件的经典矩形结构。 通过将数据存储在不同的文件中，我们可以获得样本和基因变体的信息。任何数据分析的第一步都是通过生成一些基本的汇总统计数据来了解数据。然而，hapmap ceu数据中的上述三个链接文件是无法读取的二进制格式。可以使用选项将二进制文件转换为人类可读的文件集——使用下面的命令重新编码。请注意，只有当工作目录是/User/YourName/Chapter8目录（或您之前创建的目录，并且PLINK二进制文件与您的操作系统兼容）时，此命令才起作用。旧系统上的Windows用户可能会收到警告输出，请参阅附录1。通过在PLINK命令中包含选项--bfile，我们调用了一个二进制文件。我们特别说“调用”而不是打开，因为它不是一个可以打开和观察的可读文件。
命令：
plink --bfile hapmap-ceu --recode --out hapmap-ceu
上面的命令显示，唯一的区别是plink.exe。直接打印在终端上的一系列输出告诉我们PLINK在做什么，并提供有关文件的一些基本描述信息。2下面我们显示打印到终端屏幕的输出（分为两部分）。输出的第一行显示版本号vl。90b6.7和软件发布日期（2018年12月2日）、PLINK版本（本例中为PLINK 1.90b6.7）和版权。然后，我们看到日志文件保存在默认的“plink.log”中，正如我们前面提到的。但还要注意，由于我们使用--out hapmap ceu命令更改了输出文件的名称，因此现在将成为输出文件。
日志：
$ plink --bfile hapmap-ceu --recode --out hapmap-ceu
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-ceu.log.
Options in effect:
--bfile hapmap-ceu
--out hapmap-ceu
--recode
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--recode ped to hapmap-ceu.ped + hapmap-ceu.map ... done.
Warning: 59 het. haploid genotypes present (see hapmap-ceu.hh ); many commands
treat these as missing.
输出的第二部分报告了有关文件中标记和个体数量的重要信息。在下面的示例中，我们看到PLINK从中加载了2239392个变体。fam文件中60名个体（30名男性和30名女性）的bim文件。我们还看到有60名创始人和0名非创始人出席。创始人是指数据集中没有父母的个人。PLINK只使用创始人计算等位基因频率，不包括非创始人（例如兄弟姐妹对数据集）。PLINK还可能报告一系列注释、警告和错误，当它检测到可能有错误或在某些方面不标准的潜在问题时，但它不会停止PLINK命令的执行。有关警告的更详细讨论和解释，请参阅方框8.3，并在下面的输出中注明。另一个重要的统计数据是基因分型率（如下所示为0.992022），这是个体可用标记的平均比例（具有非缺失数据）[3]。在本章后面，我们将了解如何使用这些信息来评估数据的质量。
常见的plink报错：
普林克有三种类型的错误和警告，其严重程度从“注意”的轻微警告到“警告”，再到更致命的“错误”注释只是简单地告诉你一些可能有用的信息，但它并不表示错误或错误。上述输出中的注释表明不存在表型，对于本例而言，这不是问题，因为我们没有执行关联分析或加载表型数据。为此，我们需要使用--pheno命令。一个警告提醒您可能有错误，但对您的分析来说并非致命。 警告是：警告：59 het。存在单倍体基因型（见hapmapceu.hh）；许多命令将其视为缺失，这通常与X染色体假常染色体区域中的男性杂合呼叫有关。假常染色体区域（PAR1，PAR2）是X和Y染色体上的同源序列。强烈建议您特别注意警告，因为在PLINK 2.0中，它可以被视为实际错误。建议检查中的变体。hh文件，如果它们都接近X染色体的开始或结束，请使用命令--split-X来解决这个问题。错误是导致PLINK终止的严重问题。这可能包括诸如“无输入数据集”、“内存不足”等消息。请尝试使用--内存和/或--并行标志或“所有人已删除”请参阅PLINK网站和丰富的资源，以解释这些问题以及如何更详细地处理这些问题。
他利用了这个优势。ped和。映射文件而不是二进制文件是我们可以通过在文本编辑器中打开它们来直观地检查这些文件。这可以帮助我们理解数据的结构，但需要更多内存，并且是一种计算效率较低的存储方法。我们建议您在处理大型文件时使用二进制数据（例如，在处理数千个个体的全基因组数据时）。
导入其它格式的数据
选项--make bed可用于将其他格式的数据转换为PLINK二进制文件。这将创建一个更紧凑的数据版本，不仅可以节省空间，还可以加快分析速度。例如，以下命令可用于转换。vcf文件转换为PLINK二进制文件。A.vcf文件是指1000基因组项目文本变量调用格式，其中包含变量信息。包括样本ID和基因型调用文本文件。我们将在QC一节中更详细地定义和描述基因型调用。请注意，PLINK的当前版本仅适用于基因型数据（见第7章）。输入的变异自动分配给最可能的等位基因。如前一章所述。PLINK 2.0（在撰写本书时正在开发）将有效存储插补质量信息。
代码：
plink --vcf ALL.chr21.vcf.gz --make-bed --out test_vcf
日志：
$ plink --vcf ALL.chr21.vcf.gz --make-bed --out test_vcf
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 test_vcf.log.
Options in effect:
--make-bed
--out test_vcf
--vcf ALL.chr21.vcf.gz
15236 MB RAM detected; reserving 7618 MB for main workspace.
--vcf: test_vcf-temporary.bed + test_vcf-temporary.bim + test_vcf-temporary.fam
written.
1105538 variants loaded from .bim file.
2504 people (0 males, 0 females, 2504 ambiguous) loaded from .fam.
Ambiguous sex IDs written to test_vcf.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 2504 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.999787.
1105538 variants and 2504 people pass filters and QC.
Note: No phenotypes present.
--make-bed to test_vcf.bed + test_vcf.bim + test_vcf.fam ... done.
上面数据中，包括1105523个SNP，2504个样本数据。
PLINK可以读取各种格式的数据。我们在这里提供了一个命令列表，可用于从不同软件导入数据，包括来自23andMe的个人数据。我们仅将其作为示例。
plink还支持其它格式的数据：
- Oxford format
plink --gen oxford__example.gen --sample oxford example.sample --make-bed --out oxford_example - Bgen format:
plink --bgen bgen_example --make-bed --out bgen_example - 23andme format:
plink --23file 23andme_example --make-bed --out 23andme_example
数据管理
plink还可用于管理和转换数据。例如，我们可能希望将分析限制在个体子集或某些标记上。如果我们需要合并数据集进行分析，PLINK可以用于确保报告的遗传变异等位基因以相同的方式编码。我们在另一章中讨论了协调以不同方式编码的SNP的方法。当使用来自不同基因分型平台的数据时，这一点尤为重要。在本节中，我们首先描述如何选择个体和标记，然后如何合并不同的遗传文件。
选择个体和位点
选择个体
可用于选择或排除某些个人。这可以通过为PLINK提供一个文件来实现，该文件包含将在分析中包括（使用--keep）或排除（使用--remove）的个人ID。该文件必须是一个以空格/制表符分隔的文本文件，第一列中有族ID，第二列中有族ID。
--keep选项可用于从样本中选择个体。--remove选项执行相反的操作，并从分析中排除文件中列出的个人。
代码：
plink --bfile hapmap-ceu --keep list.txt --make-bed --out selectedIndividuals
日志：
$ plink --bfile hapmap-ceu --keep list.txt --make-bed --out selectedIndividuals
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 selectedIndividuals.log.
Options in effect:
--bfile hapmap-ceu
--keep list.txt
--make-bed
--out selectedIndividuals
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
--keep: 60 people remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--make-bed to selectedIndividuals.bed + selectedIndividuals.bim +
selectedIndividuals.fam ... done.
Warning: 59 het. haploid genotypes present (see selectedIndividuals.hh ); many
commands treat these as missing.
类似地，我们可以通过使用--keep fam和--remove fam选项来选择或排除整个族。当您使用--make bed选项时，缺失率和等位基因频率的阈值过滤器会自动设置为不排除任何人。您可以使用-mind、-geno和-maf手动指定这些过滤器来排除人。这也是可以应用命令-extract/-exclude和-keep/-remove的时候。 例如，如果您想创建一个新文件，其中只包含高基因型至少95%个完整的个体，您可以添加--mind 0.05命令。
plink --bfile hapmap-ceu --make-bed --mind 0.05 --out highgeno
$ plink --bfile hapmap-ceu --make-bed --mind 0.05 --out highgeno
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 highgeno.log.
Options in effect:
--bfile hapmap-ceu
--make-bed
--mind 0.05
--out highgeno
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
0 people removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--make-bed to highgeno.bed + highgeno.bim + highgeno.fam ... done.
Warning: 59 het. haploid genotypes present (see highgeno.hh ); many commands
上面的代码，会生成三个前缀为highgeno的文件，注意，不要尝试打开bed文件，因为它是二进制文件。
选择位点
可能是您对一个SNP、一组SNP或特定区域感兴趣。PLINK还可以用于将特定变体的基因型信息提取到单独的较小文件中。例如，如果你有1000个基因组数据，你可能对单核苷酸多态性列表感兴趣，例如单核苷酸多态性3的列表。我们可以通过使用--SNPs选项指定要选择的snp的名称，或者通过提供一个文件（使用选项--include）来实现这一点，该文件包含我们要包含在新文件中的变体的标记名。选项--exclude可用于从文件中删除某些变体。以下示例说明了如何选择单个变体的基因型。变体rs9930506是FTO基因中的一个单核苷酸多态性，一些研究表明该基因与BMI和体重增加有关。可以使用以下命令选择此变体。在上一节中，我们演示了如何从direct to consumer company 23andMe读取数据。如果您有这些数据，并且想知道您携带了多少风险等位基因rs9930596，您可以导入您自己的个人数据，并使用以下命令进行检查。rs9930596是FTO基因的一种遗传变异，G等位基因与肥胖风险增加有关。
代码：
plink --file hapmap-ceu --snps rs9930506 --make-bed --out rs9930506sample
日志：
$ plink --file hapmap-ceu --snps rs9930506 --make-bed --out rs9930506sample
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 rs9930506sample.log.
Options in effect:
--file hapmap-ceu
--make-bed
--out rs9930506sample
--snps rs9930506
15236 MB RAM detected; reserving 7618 MB for main workspace.
.ped scan complete (for binary autoconversion).
Performing single-pass .bed write (2239392 variants, 60 people).
--file: rs9930506sample-temporary.bed + rs9930506sample-temporary.bim +
rs9930506sample-temporary.fam written.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
--snps: 1 variant remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
1 variant and 60 people pass filters and QC.
Note: No phenotypes present.
--make-bed to rs9930506sample.bed + rs9930506sample.bim + rs9930506sample.fam
... done.
如果查看bim文件，内容如下：
$ cat rs9930506sample.bim
16 rs9930506 0 52387966 G A
合并数据
合并基因型数据
在这种类型的分析中，我们经常处理多个文件。例如，通常情况下，基因组数据是通过染色体存储的，以避免在互联网上难以传输且需要大量计算的巨大文件。因此，通常情况下，你有22个独立的常染色体文件，性染色体的文件，也许还有一个线粒体DNA的fle，并希望将所有内容合并到一个文件中。或者，您可能有来自不同子组的文件，例如来自不同基因型中心的文件。在其他情况下，可能需要合并来自不同研究的文件以创建单个文件。合并基因文件需要相当小心。在一个文件中测量的变异可能不会在另一个文件中测量，并且可能具有不同的等位基因或碱基对位置。PLINK中的命令-bmerge确保这些问题得到正确解决。默认情况下，命令--bmerge将这些类型的不匹配设置为missing。然而，可以使用选项——合并模式来制定不同的规范。如果需要同时合并多个文件，就像合并染色体特定文件一样，请使用包含不同基因型文件名的文件和选项——merge-1ist。下面的示例合并了两个不同的基因文件。 但请注意，它仅作为说明性示例提供（因此不执行）。
比如现在有两个plink的二进制文件，HapMap_founder和HapMap__nonfounders进行合并，用--bmerge:
plink --bfile HapMap_founders --bmerge HapMap__nonfounders -make-bed --out merged file
如果有多个plink文件，用--merge-list
合并表型数据
遗传分析的主要目的是研究基因型和表型之间的关系。表型可以在fam文件。要将表型“附加”到遗传文件，我们可以在PLINK中使用--pheno命令，如上面的示例所示。在这里，我们使用两个新文件：1kg_EU_qc和BMI_pheno.txt。首先看看fam文件，可以看到第六列表型数据都是-9，都是缺失。
$ head 1kg_EU_qc.fam
0 HG00096 0 0 1 -9
0 HG00097 0 0 2 -9
0 HG00099 0 0 2 -9
0 HG00100 0 0 2 -9
0 HG00101 0 0 1 -9
0 HG00102 0 0 2 -9
0 HG00103 0 0 1 -9
0 HG00104 0 0 2 -9
0 HG00106 0 0 2 -9
0 HG00108 0 0 1 -9
再看一下bim文件：
$ head 1kg_EU_qc.bim
1 rs1048488 0 760912 C T
1 rs3115850 0 761147 T C
1 rs2519031 0 793947 G A
1 rs4970383 0 838555 A C
1 rs4475691 0 846808 T C
1 rs1806509 0 853954 C A
1 rs7537756 0 854250 G A
1 rs28576697 0 870645 C T
1 rs7523549 0 879317 T C
1 rs3748592 0 880238 A G
查看一下表型数据：
$ head BMI_pheno.txt
FID IID BMI
0 HG00096 25.022827
0 HG00097 24.853638
0 HG00099 23.689295
0 HG00100 27.016203
0 HG00101 21.461624
0 HG00102 20.673635
0 HG00103 25.71508
0 HG00104 25.252243
0 HG00106 22.765049
合并命令：
plink --bfile 1kg_EU_qc --pheno BMI_pheno.txt --make-bed --out 1kg_EU_BMI
日志：
$ plink --bfile 1kg_EU_qc --pheno BMI_pheno.txt --make-bed --out 1kg_EU_BMI
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 1kg_EU_BMI.log.
Options in effect:
--bfile 1kg_EU_qc
--make-bed
--out 1kg_EU_BMI
--pheno BMI_pheno.txt
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
379 people (178 males, 201 females) loaded from .fam.
379 phenotype values present after --pheno.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 379 founders and 0 nonfounders present.
Calculating allele frequencies... done.
851065 variants and 379 people pass filters and QC.
Phenotype data is quantitative.
--make-bed to 1kg_EU_BMI.bed + 1kg_EU_BMI.bim + 1kg_EU_BMI.fam ... done.
查看更新后的fam文件，可以看到已经有内容了。
$ head 1kg_EU_BMI.fam
0 HG00096 0 0 1 25.0228
0 HG00097 0 0 2 24.8536
0 HG00099 0 0 2 23.6893
0 HG00100 0 0 2 27.0162
0 HG00101 0 0 1 21.4616
0 HG00102 0 0 2 20.6736
0 HG00103 0 0 1 25.7151
0 HG00104 0 0 2 25.2522
0 HG00106 0 0 2 22.765
0 HG00108 0 0 1 23.069
描述性统计
与任何类型的分析一样，了解数据并进行初始描述性分析很重要。PLINK还可以用于获取有关用于分析的数据或基因型变体的信息。
基因频率
基因频率可以通过PLINK中的--freg命令计算。输出文件（后缀为.frq）包含关于基因型的等位基因和次要等位基因频率（MAF）以及每个SNP的等位基因代码的信息。
命令：
plink --bfile hapmap-ceu --freq --out Allele_Frequency
日志：
$ plink --bfile hapmap-ceu --freq --out Allele_Frequency
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 Allele_Frequency.log.
Options in effect:
--bfile hapmap-ceu
--freq
--out Allele_Frequency
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
--freq: Allele frequencies (founders only) written to Allele_Frequency.frq .
Warning: 59 het. haploid genotypes present (see Allele_Frequency.hh ); many
commands treat these as missing.
查看结果，结果文件是Allele_Frequency.frq
$ head Allele_Frequency.frq
CHR SNP A1 A2 MAF NCHROBS
1 rs12565286 C G 0.0678 118
1 rs12138618 A G 0.05833 120
1 rs3094315 G A 0.1552 116
1 rs3131968 A G 0.125 120
1 rs12562034 A G 0.09167 120
1 rs2905035 A G 0.125 120
1 rs12124819 G A 0.3417 120
1 rs2980319 A T 0.125 120
1 rs4040617 G A 0.125 120
在此之前，您将看到以下列：CHR（染色体数或性染色体代码）；SNP（变体名称，大多数SNP的rsID）：A1（等位基因1，通常是次要等位基因[即，频率较低）；A2（等位基因2，通常是主要等位基因）、MAF（等位基因1频率）和NCHROB（等位基因观察数），也可以进行多种其他分析，例如使用-in选项通过分类变量进行分层。
缺失统计
样本和基因型缺失统计
检查缺失值对于评估数据质量很重要。有两种类型的缺失值需要考虑。第一种是个人。对于每个个体，我们计算缺失的所有变体的比例。其次，基因变异也可能有缺失值。对于每个遗传变异，我们计算已进行基因分型的个体比例和缺失值的比例。PLINK用于导出此信息的命令缺失，该命令生成两个输出文件：a。个人缺失信息的imiss文件和SNP缺失的.lmiss文件。
命令：
plink --bfile hapmap-ceu --missing --out missing_data
日志：
$ plink --bfile hapmap-ceu --missing --out missing_data
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 missing_data.log.
Options in effect:
--bfile hapmap-ceu
--missing
--out missing_data
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
--missing: Sample missing data report written to missing_data.imiss, and
variant-based missing data report written to missing_data.lmiss.
Warning: 59 het. haploid genotypes present (see missing_data.hh ); many
commands treat these as missing.
查看样本缺失的文件：
然后，您将看到以下输出。在单个缺失fle（IMIS）中，它的结构是每行代表一个个体。这些列涉及：FID（家族ID）、IID（家族ID内）、MISS PHENO（缺失表型的是/否指标）、N\u MISS（缺失基因型调用数）、N\u GENO（潜在有效调用数）和F\u MISS（缺失调用率）。
$ head missing_data.imiss
FID IID MISS_PHENO N_MISS N_GENO F_MISS
1334 NA12144 Y 15077 2239392 0.006733
1334 NA12145 Y 19791 2239392 0.008838
1334 NA12146 Y 13981 2239392 0.006243
1334 NA12239 Y 14072 2239392 0.006284
1340 NA06994 Y 16080 2239392 0.007181
1340 NA07000 Y 26113 2239392 0.01166
1340 NA07022 Y 17467 2239392 0.0078
1340 NA07056 Y 12133 2239392 0.005418
1341 NA07034 Y 20425 2239392 0.009121
查看位点缺失的文件：
在这里，您将看到这些列由CHR（染色体代码）、SNP（变体标识符）、N MISS（缺失基因型调用的数量，不包括强制性缺失）、N GENO（潜在值调用的数量）和F MISS（缺失调用率）表示。仅在族内数据中存在的列是CLST（簇标识符）和N-CLST（簇大小）。
$ head missing_data.lmiss
CHR SNP N_MISS N_GENO F_MISS
1 rs12565286 1 60 0.01667
1 rs12138618 0 60 0
1 rs3094315 2 60 0.03333
1 rs3131968 0 60 0
1 rs12562034 0 60 0
1 rs2905035 0 60 0
1 rs12124819 0 60 0
1 rs2980319 0 60 0
1 rs4040617 0 60 0
数据筛选
- --filter-controls，提取control文件
- --filter-males，提取male文件
- --filter-females，提取female文件
- --fileter-nonfounders，无亲本信息的个体？？？
数据质控
在第4章中，我们已经简要介绍了有关GWASs的内容，质量控制（QC）是处理遗传数据的标准程序。不良的研究设计和基因型调用错误会在遗传研究中引入系统性偏见，尤其是在关联研究中，破坏结果的有效性。遗传研究应该有一个严格的质量控制方案，该方案已经在研究的初始设计中计划好了（见第4章，我们在那里讨论了GWAS研究设计）。为了确保OC协议独立于结果，该协议通常在外部存储库（如开放科学框架）中预先注册(https://osfio)，另外。 来自大型GWAS联盟的研究人员通常旨在确保所有OC程序由独立的中心执行，尽管这是相当可变的。质量控制协议的建立是为了最大限度地降低关联性普遍膨胀的风险，以及在给定研究中，当与疾病没有真正相关的变异与疾病显著相关时，出现假阳性错误的风险。数据质量差可能与多种问题有关，例如表型测量差或不一致，这会降低统计能力（即使用统计程序检测真实关联的能力）。良好的质量控制方案确保数据在研究中具有可比性，并可用于后续分析。我们现在描述的三种主要QC类型是：（1）每个QC，（2）每个标记QC，（3）全基因组关联荟萃分析QC。我们将简要介绍主要程序以及如何在PLINK中实现它们。
样本质控
样本缺失质控
--mind
第一步是确保样本中的个体拥有高质量的数据。 全基因组数据的每个QC包括设置过滤器，从样本中删除可能因数据质量低而在分析中引入偏差的个体。 每个QC一般包括五个步骤，即个体识别：1。DNA质量差（低呼叫率和缺失基因型，区分如下）；2、常染色体杂合度高，表明可能存在样本污染或杂合度低，这可能是由于近亲繁殖所致；3、性别信息不一致；（四）重复或者相关的；和5。来自不同的祖先群体。
识别基因型质量差的个体这类分析中的一个挑战是缺少来自基因型缺失率高（即低基因型调用率）的个体的数据，以及在很大一部分个体中缺失的SNP。在这里，我们讨论了影响基因型调用率和基因型准确性的DNA质量的巨大差异。有时需要从样本中移除低质量的DNA样本。当我们观察到缺失基因型比例较高的个体时，就有可能检测到这一点。通常情况下，缺失基因型超过3-7个的个体会被从分析中删除，选择mind和缺失截止值的规格，例如，对于5个缺失率，0.05。这可以在PLINK中指定如下：
代码：
plink --bfile 1kg_hm3 --mind 0.05 --make-bed --out 1kg_hm3_mind005
日志：
$ plink --bfile 1kg_hm3 --mind 0.05 --make-bed --out 1kg_hm3_mind005
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 1kg_hm3_mind005.log.
Options in effect:
--bfile 1kg_hm3
--make-bed
--mind 0.05
--out 1kg_hm3_mind005
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
0 people removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
851065 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--make-bed to 1kg_hm3_mind005.bed + 1kg_hm3_mind005.bim + 1kg_hm3_mind005.fam
... done.
样本杂合质控
--het
识别常染色体杂合性杂合性是指携带特定SNP的两个不同等位基因（见第1章第1.5节）。 个体的杂合率是杂合基因型的比例。 个体内高水平的杂合性可能表明样本质量低，而低水平的杂合性可能是由于近亲繁殖。因此，我们在分析中排除了杂合性水平高和低的个体。原因可能是近亲繁殖，但也有样品污染。PLINK中计算杂合度的基本命令描述如下（选项：--het）。
日志：
plink --bfile 1kg_hm3 --het --out 1kg_hm3_het
代码：
$ plink --bfile 1kg_hm3 --het --out 1kg_hm3_het
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 1kg_hm3_het.log.
Options in effect:
--bfile 1kg_hm3
--het
--out 1kg_hm3_het
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
851065 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--het: 851065 variants scanned, report written to 1kg_hm3_het.het .
创建wo文件：1kg hm3 het。het和。日志每个样本的杂合度计算为杂合子基因型调用数与非缺失调用总数的比率。杂合度统计可以使用标准软件（如Excel或R）进行检查。平均杂合度值非常高或很低的异常值可以使用单独PLINK命令中的--exclude选项过滤掉。规则通常是删除偏离特定样本杂合度均值±3个标准差的个体[4]。图8.3显示了使用前一个PLINK命令的输出的杂合度统计分布。红线表示OC阈值与平均值的±3个标准偏差。杂合度较高或较低的观察值可以从连续分析中删除。要生成此直方图，请在RStudio中使用以下命令（如果尚未安装R和RStudio，请参阅附录1以获取有关如何安装R和RStudio的信息）：
用R语言作图：
heterogeneity_stats<-read.table("1kg__hm3_het.het",
header=T)
attach(heterogeneity _ stats)
het<-(N.NM.-O.HOM.)/N.NM.
min<-mean(het)-3*sd(het)
max<-mean(het)+3*gd(het)
hist(het,col="lightblue",xlim=c(0.25,0.37),
main="Histogram of Heterogeneity")
abline(v = c(min,max),col="red",1wd=3,lty=2)
image-20220713145043478
性别检测错误
对性别信息不一致的个体的识别表明，男性只有一个X染色体拷贝，因此性染色体中的任何标记都不能杂合。因此，使用来自X染色体的基因型数据可以检查与已确定性别的不一致性。通常，当基因型调用算法检测到X染色体标记的男性杂合子时，它将该基因型称为缺失。 获得样本的正确性别信息非常重要，尤其是在按性别分层进行分析时。要计算此值，请运行以下命令，包括--check sex选项。
男性的X染色体近亲繁殖系数应大于0.8，而雌性的X染色体近亲繁殖系数应小于0.2。PLINK将估计的近亲繁殖系数与确定的性别进行比较，并将值超出这些阈值的个体作为可能的不匹配信号。不一致的性别信息可能是实验室错误造成的，如果你的样本中有很多人有这种信息，你应该更详细地检查数据。文件中报告了X染色体近亲繁殖系数F。在上一个PLINK命令中获得的sexcheck。
代码：
plink --bfile hapmap-ceu --check-sex --out hapmap_sexcheck
日志：
$ plink --bfile hapmap-ceu --check-sex --out hapmap_sexcheck
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_sexcheck.log.
Options in effect:
--bfile hapmap-ceu
--check-sex
--out hapmap_sexcheck
15236 MB RAM detected; reserving 7618 MB for main workspace.
2239392 variants loaded from .bim file.
60 people (30 males, 30 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 60 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992022.
2239392 variants and 60 people pass filters and QC.
Note: No phenotypes present.
--check-sex: 56526 Xchr and 0 Ychr variant(s) scanned, no problems detected.
Report written to hapmap_sexcheck.sexcheck .
Warning: 59 het. haploid genotypes present (see hapmap_sexcheck.hh ); many
commands treat these as missing.
查看结果：
$ head hapmap_sexcheck.sexcheck
FID IID PEDSEX SNPSEX STATUS F
1334 NA12144 1 1 OK 0.9999
1334 NA12145 2 2 OK -0.06528
1334 NA12146 1 1 OK 0.9999
1334 NA12239 2 2 OK 0.05498
1340 NA06994 1 1 OK 0.9999
1340 NA07000 2 2 OK -0.1001
1340 NA07022 1 1 OK 0.9998
1340 NA07056 2 2 OK -0.03786
1341 NA07034 1 1 OK 0.9999
重复样本检测
重复或相关个体的识别重要的是检查个体的无意重复和隐性关联，或者换句话说，vour样本中的个体可能是近亲的可能性。为了识别重复和相关个体，可以根据基因型SNP（仅常染色体）上共享的等位基因的平均比例，为每对个体计算一个称为状态识别（IBS）的度量。通过使用全基因组SNP及其衍生的IBS，可以估计一对个体的共同祖先程度（identity by generation.IBD）。IBS等于1的成对个体要么是单卵双生子，要么是数据中的重复个体。兄弟姐妹的预期IBD为0.5，二级亲属为0.25。由于基因分型错误和种群结构，这些理论值经常存在差异。通常，IBS>0.1875的成对个体从分析中删除。基于线性混合模型（如BOLT-LMM）的最新方法不要求我们删除相关个体，因为它们在计算中考虑了相关性。第9章说明了如何计算相关性以及如何基于IBS统计筛选个体。
群体分层
不同祖先个体的识别：如第3章所述，群体分层引起的混淆是关联研究中偏见的主要来源。因此，通常按大陆血统进行分层分析。识别高水平祖先群体的一种常见做法是对大量变体进行主成分分析。下一章第9.4节将对此进行更详细的论述。
位点质控
第二组质量控制分析侧重于变体的数据质量。我们依次采取了几个步骤来消除可能在研究中引入偏见的低质量变体。每个标记QC包括五个步骤：207处理遗传数据，第一部分1。排除低呼叫率SNP 2。在某些类型的分析中去除具有极低等位基因频率（罕见变体）的SNP 3。识别和排除与HardyWeinberg平衡（HWE）4极度偏离的变体。在病例对照研究中，排除组间呼叫率极为不同的SNP 5。在插补单核苷酸多态性的情况下，从插补质量低的变体研究中排除，重要的是去除可能影响分析的次优标记，从而增加假阳性的数量。重要的是要注意哪些标记被排除在外，以及与此相关的潜在LD方面。下面，我们对执行每标记QC分析所需的步骤进行了最低限度的解释。
缺失质控
--geno
低呼叫率SNPs低呼叫率是指基因型缺失率高的情况，或者换句话说，SNPs在很大一部分个体中缺失。在这种情况下，“调用”用于表示对一个唯一SNP的估计。通常，呼叫率低于95的标记被排除在分析之外。 这可以通过PLINK命令中的选项--geno和缺失截止值的指定来指定，例如，对于5缺失率，为0.05。
$ plink --bfile 1kg_hm3 --geno 0.05 --make-bed --out 1kg_hm3_geno
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 1kg_hm3_geno.log.
Options in effect:
--bfile 1kg_hm3
--geno 0.05
--make-bed
--out 1kg_hm3_geno
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
0 variants removed due to missing genotype data (--geno).
851065 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--make-bed to 1kg_hm3_geno.bed + 1kg_hm3_geno.bim + 1kg_hm3_geno.fam ... done.
它会生成三个文件。
次等位基因频率
第1章（第1.3.1节）回顾了具有低次要等位基因频率的NP，次要等位基因频率（MAF）是一个位点上第二常见的等位基因在给定人群中出现的频率。例如，在HapMap项目中。MAF为5？的单核苷酸多态性？。05）或更高。具有低MAF的SNP是罕见的。因此，分析人员只包括高于某个MAF阈值的SNP。排除阈值低的MAF有两个原因。首先，在低MAF的情况下，缺乏检测任何真正SNP特征关联的能力。其次，这些SNP通常更容易出现基因分型错误。由于基因分型错误和检测关联的能力较低，低频SNP更容易产生偏差，因此通常从研究中删除。我们可以根据MAF阈值1，选择MAF 0.01排除SNP。
$ plink --bfile 1kg_hm3 --maf 0.01 --make-bed --out 1kg_hm3_maf
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 1kg_hm3_maf.log.
Options in effect:
--bfile 1kg_hm3
--maf 0.01
--make-bed
--out 1kg_hm3_maf
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
4581 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
846484 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--make-bed to 1kg_hm3_maf.bed + 1kg_hm3_maf.bim + 1kg_hm3_maf.fam ... done.
创建wo文件，1kg hm3 het。het和。日志MAF阈值的选择在很大程度上取决于样本量。如果样本量较大，可以使用较低的MAF阈值。虽然没有固定的规则，但大样本被认为是一个类似N=100.000、MAF阈值0.01的值，而中等样本N=10.000。MAF阈值0.05；低于此值将被视为一个小样本。通常情况下，MAF小于1-2的SNP被排除在研究之外，但样本量小的研究可能会应用更高的阈值。另请注意，这1-2也随着插补的进步而减少。
哈温平衡质控
与第3章（第3.5节）中详细描述的哈代-温伯格平衡（HWE）的偏差，HWE假设人口无限大，没有选择、突变或迁移，如果不违反任何条件，则基因型和等位基因频率在世代中是恒定的。因此，偏离HWE可以表明基因分型错误、进化以及对等位基因和基因型频率之间关系的担忧。违反HWE表明基因型频率与预期显著不同，观察到的频率不应显著不同。例如，如果等位基因A的频率为0.20，等位基因T的频率为0.80，则基因型AT的预期频率为20.20.8=0.32。在GWASs中，通常假设与HWE的偏差是基因分型错误的结果。HWE测试的PLINK命令是--HWE，后跟违反显著性阈值的规范。
$ plink --bfile 1kg_hm3 --hwe 0.00001 --make-bed --out 1kg_hm3_hwe
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 1kg_hm3_hwe.log.
Options in effect:
--bfile 1kg_hm3
--hwe 0.00001
--make-bed
--out 1kg_hm3_hwe
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
--hwe: 0 variants removed due to Hardy-Weinberg exact test.
851065 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--make-bed to 1kg_hm3_hwe.bed + 1kg_hm3_hwe.bim + 1kg_hm3_hwe.fam ... done.
他的作品有四部新的FLE，。日志bim，。fam和。床排除因你是否具有二元或数量（即连续）特征而不同。对于二元性状，规则通常是HWE p值在病例中<10-10，在对照组中<10-6。对于数量性状，建议HWE p值<10-6[4]。文献中采用了不同的p值阈值，不同的参考文献具有不同的水平。
所有质控条件合并
结合不同的质量控制过滤器以删除所有失败的SNP多个质量控制过滤器，我们可以同时应用之前在个体和标记水平上涵盖的命令。文件个人失败。txt包括所有未通过个体水平QC的个体（例如，极端杂合性或相关个体）。
$ plink --bfile 1kg_hm3 --mind 0.03 --geno 0.05 --maf 0.01 --hwe 0.00001 --exclude individuals_failQC.txt --make-bed --out 1kg_hm3_QC
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 1kg_hm3_QC.log.
Options in effect:
--bfile 1kg_hm3
--exclude individuals_failQC.txt
--geno 0.05
--hwe 0.00001
--maf 0.01
--make-bed
--mind 0.03
--out 1kg_hm3_QC
15236 MB RAM detected; reserving 7618 MB for main workspace.
851065 variants loaded from .bim file.
1092 people (525 males, 567 females) loaded from .fam.
--exclude: 851065 variants remaining.
0 people removed due to missing genotype data (--mind).
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1083 founders and 9 nonfounders present.
Calculating allele frequencies... done.
0 variants removed due to missing genotype data (--geno).
--hwe: 0 variants removed due to Hardy-Weinberg exact test.
4581 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
846484 variants and 1092 people pass filters and QC.
Note: No phenotypes present.
--make-bed to 1kg_hm3_QC.bed + 1kg_hm3_QC.bim + 1kg_hm3_QC.fam ... done.
GWAS分析Meta分析数据质控
在我们关于GWASs的章节中进行了相当详细的描述，大多数研究都是基于关联结果的元分析，该分析使用了数百万个基因变体组合了多个数据集。为了避免II型错误（误报），GWAS研究人员花费大量精力在最终荟萃分析QC阶段识别归因于低质量数据的异常值。筛选出对一系列检查不稳健的变体，不包括在分析中。在这一阶段应用的典型过滤器去除了等位基因频率低、插补质量低、等位基因频率与参考样本相差很大或由特定研究驱动的关联结果无法在其他地方复制的变体。在某些情况下，这可能是之前的一些步骤的重复（例如删除次要等位基因频率较低的SNP），这些步骤通常在荟萃分析之前筛选出来。 然而，有几个原因可以解释为什么需要在荟萃分析阶段重复这些步骤。首先，在单个研究中没有显示任何问题的变体在多队列研究中汇总时可能会有额外的问题，正如我们在GWAS一章中指出的那样，这些研究可能非常复杂，并且使用来自世界各地数十个个体遗传学研究的数据。例如，2016年关于人类生殖行为的GWAS基于对63个不同队列的性别特异性分析的三种表型分析[5]。这意味着对数百个FLE进行了广泛检查，其中包含数百万个基因变体的关联结果。该研究中的QC由两个独立的中心（位于英国牛津和荷兰鹿特丹）进行，该中心分析数据并比较HLTER和诊断工具，以确保报告的结果基于最佳质量变体。这样一个大型项目的质量控制通常需要几个月的时间。在本节中，我们总结了GWASs质量控制程序的主要步骤和诊断测试。 尽管研究人员不想放弃“真实”的结果，但GWAS分析师通常倾向于保守的方法，以确保低质量的标记被排除在分析之外。这些是可能的过滤器，可以应用于全基因组关联结果的质量控制。
1、GWAS的第一步是使用公共引用协调文件中标记的碱基对位置。例如，美国国家生物技术信息中心d 37。 2、还需要评估关联结果的信息，尤其是来自多个文件的信息。从样本中删除效应等位基因或替代等位基因信息缺失的标记，或效应估计值、标准误差或p值不可信的标记。 3、非双等位基因（与两个等位基因相关）或单态（无变异或与多态性相反）的标记被排除在最终结果列表之外。 4、罕见变异的结果通常有问题，可能会影响结果，MAF低于1的标记通常会被删除。然而，该阈值可能取决于研究的样本量。一些小型研究往往会从MAF小于5的变体中得出有问题的结果？虽然大型生物库可能不会对不太常见的变体产生问题。 5、插补质量会影响关联结果的质量。插补质量低于0.7的分析标记下降很常见。
诊断检查在此之后，通常会对应用上述过滤器后剩余的SNP进行额外的诊断检查。一个非常常见的图是等位基因频率图（见图8.4）。这检查了变体（1）是否以相同的方式编码，等位基因频率和链方向是否存在错误，以及（2）在研究中是否具有相似的等位基因频率。通常会绘制一张图，其中每个研究中的等位基因频率与参考文献中的等位基因频率相对应
image-20220713151250900
数量充足，通常为1000个基因组。图8.4显示了一个等位基因频率图，该图在我们的大型GWAS人类生殖行为联合会的OC中使用[5]。该图是一个散点图，y轴上的预期初生年龄等位基因频率（AFB）和x轴上参考联盟（1000个基因组）的等位基因频率，因为我们需要绘制约1000万个点，所以我们只绘制与参考样本偏离的点。特别是，这里我们只绘制等位基因频率差异为0.2的点（因此是对角线）。在这种情况下，可以看到一些变体的等位基因频率低于参考样本。 第二个诊断检查对于检查结果的异质性非常有用，即所谓的森林图（见图8.5）。森林图是元分析中常用的诊断工具。在该图中，取自我们上面讨论的同一研究，绘制了给定变体（rs2777888）的不同研究的关联结果及其在研究中使用的所有数据源的置信区间。这是结果异质性的可视化表示。该图验证了所有结果具有相同的作用方向。较大的研究应具有较小的置信区间，给出更精确的估计。例如，在这种情况下，来自deCODE genetics和23andMe的大样本包含在我们的荟萃分析中。框的大小还表示数据的样本大小。林地通常与异质性测试相结合。 其他QC诊断工具包括报告的p值与Z分数的p值的图（也称为PZ图）和分位数-分位数（Q-Q）图，其中观察到的p值与我们在无关联情况下观察到的理论分位数绘制。O-O图可能揭示未知的人口分层。我们在第12章更详细地讨论Q-Q图（见第12.1.3节）。同样值得注意的是，我们期望在一个具有良好遗传特性的GWAS的O-O图中出现插入。当考虑基因型数据中的错误和偏差时，需要注意的另一个方面是批次效应。批次效应是指由于实验室条件、试剂批次或其他差异而产生的测量误差。有关这一主题的优秀综述，请参阅Leek及其同事2010年的一项研究[6]。
结论
本章基本介绍了如何开始使用PLINK处理遗传数据。我们的目标是，您现在已经掌握了如何使用命令行、调用和打开PLINK以及导入数据处理此类数据的基本知识。数据管理在这些类型的分析中非常重要，因此，我们试图概述一些关于如何选择个体和标记、合并不同的遗传文件或添加表型的基础知识。您现在还应该能够检查描述性统计，包括等位基因频率和缺失值过滤。我们总结了这类数据分析中最重要的步骤之一，即单个标记和GWAS级质控过程的质控。在第9章中，我们将向您展示如何进行更高级的分析，包括关联分析、GWA、独立SNP的识别、，计算PCA，并使用另一个称为GCTA（全基因组复杂性状分析）的程序[7]。
在本章中，我们研究了欧洲祖先群体。请访问以下网站，您将在那里找到其他祖先群体：http://z.bwh哈佛。edu/plink/res.shtml。 下载中文和日文样本（CHB和JPT），并进行以下分析：1。将文件重新编码到地图和。ped文件。 2、计算等位基因频率和缺失值。 3.使用以下标准进行QC分析：a.调用ate>959%b.次要等位基因频率>3%c.Hardy-Weinberg平衡测试<0.0001
plink学习资料
Online PLINK Tutorial: http://zz.bwh.harvard.edu/plink/tutorial.shtml.
PLINK webpage: http:// www.cog-genomics.org/plink2/.
YouTube PLINK tutorial by Broad Institute: https://wwwyoutube.com/watch?v=ppBJqCMSBYk.