使用Plink对CNV做GWAS分析（一）

Introduction

随着高通量测序技术的发展，这些年关于CNV（拷贝数变异）的研究也逐渐多了起来。以往CNV主要通过PCR或者CGH芯片等技术进行查找，现在由于二代测序的价格飞速下降，通过将二代测序产生的reads 回帖（map）到参考基因组上来检测CNV的流程和软件逐渐建立并且流行了起来。本文主要介绍得到CNV后如何使用Plink进行下游的GWAS（全基因组关联分析）。

软件版本

Plink v1.07
Plink最新的版本是1.9(16 May 2016)，但是由于我在Plink官网上查到的CNV关联分析流程中几个参数在新版中被移除了（？！！！）。目前还没搞明白在新版中是怎么做的，所以就先用旧的Plink做一版。

原始数据

我拿到的CNV数据是用我们实验室自己开发的软件CNVcaller，从BAM文件中得到的。软件暂时还没发表，等发表了我再来安利一发。
此外，还需要你要与之关联的表型数据（phenotype），数量性状或者质量性状等等都可以。

使用Plink进行分析

1.创建CNV list

把你的CNV按照Plink的要求进行整理。格式如下：

FID	IID	CHR	BP1	BP2	TYPE	SCORE	SITE
P1	P1	4	71338469	71459318	1	27	0
P1	P1	5	31250352	32213542	1	0	0
P1	P1	7	53205351	53481230	3	0	0

文件总共有8列，每一列代表的意思分别如下：

 FID     Family ID
 IID     Individual ID
 CHR     Chromosome
 BP1     Start position (base-pair)
 BP2     End position (base-pair)
 TYPE    Type of variant, e.g. 0,1 or 3,4 copies
 SCORE   Confidence score associated with variant 
 SITES   Number of probes in the variant

注：如果是自然群体，没有家系关系的话，FID和IID起一样的名字就可以了。而SCORE和SITES两项不会直接用于关联分析，所以不知道就直接全部写0就OK了。

The SCORE and SITES values are not used in any direct way, except potentially as variates to filter segments on, as described below. That is, the values of these do not fundamentally impact the way analysis is performed by PLINK itself (they might alter the meaning of the results of course, e.g. if including low-confidence calls into the analysis!). In other words, if whatever software was used to generate the CNV calls does not supply some conceptually similar values, it is okay to simply put dummy codes (e.g. all 0) in these two fields.

此部分详见这里

2.创建FAM file

.fam文件储存的是样本的性状、性别等信息。
这个文件没有header，每行一个样本，六列信息，分别是：

Family ID ('FID')
Within-family ID ('IID'; cannot be '0')
Within-family ID of father ('0' if father isn't in dataset)
Within-family ID of mother ('0' if mother isn't in dataset)
Sex code ('1' = male, '2' = female, '0' = unknown)
Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

示例如下：

P1	P1	0	0	1	20
P2	P2	0	0	1	30

注1：如果没有family ID, parental ID, sex, and/or phenotype columns也行，使用下面的参数就OK了（适用于.fam和.ped文件）。
--no-fid
--no-parents
--no-sex
--no-pheno

注2：第六列的表型值，如果写的不是1、2、-9、0中的任意一个，则会被判定为质量性状。

此部分详见这里

3.生成MAP文件

这一步是用之前制作的CNV list，使用Plink自动生成.map文件。CNV的map文件和plink通用的map文件一样。
第一列是Chromosome code。
第二列是Variant identifier，是ID加上CNV的起始或者终止位置。
第三列是Position in morgans or centimorgans，在CNV的map文件中都用0来替代。
第四列是该标记在染色体上的位置Base-pair coordinate。
命令如下：

plink --cnv-list plink.cnv --cnv-make-map

其中，plink.cnv就是第一步整理出来的CNV list。运行此命令后会自动生成plink.cnv.map

4.关联分析

有了前面三个文件，我们就可以开始做关联分析了。此部分以数量性状为例。
命令如下：

plink --map plink.cnv.map --fam mydata.fam --cnv-list mydata.cnv --mperm 50000 --noweb

如果你前面那三个文件名前缀统一的话（a.cnv, a.cnv.map, a.fam），你也可以用这条命令：

plink --cfile a

之后你会得到一堆文件，其中.summary后缀的文件包含六列，含义分别如下：

 CHR    Chromosome code
 SNP    Dummy label for map position
 BP     Physical position (base-pairs)
 NCNV   Number of individiuals with a CNV here
 M1     QT mean in individuls with a CNV here
 M0     QT mean in individuals without a CNV here

打开.mperm文件，里面包含置换检验的P值：EMP1和EMP2。

EMP1    Empirical p-value, per region
EMP2    Empirical p-value, corrected for all tests

注1：--noweb是说Skipping web check。--mperm是置换检验。
注2：QT mean应该是数量性状的均值的意思。
注3：质量性状默认使用单侧检验，数量性状默认使用双侧检验。如果想使用单侧检验，添加参数：--cnv-test-1sided

结语

至此，使用Plink对CNV做GWAS分析的第一部分就讲完了。我们在.mperm中可以得到想要的结果，包括CNV和对应的P值。
上面的分析在plink软件中是属于Rare copy number variant (CNV)
分析，而其中的好多参数在plink 1.9中被移除了，而且在这里我也没找到如何导入协变量以消除群体结构的影响，GWAS的模型也不可选。
因此，下一章我们将学习Common copy number polymorphism (CNP) data的GWAS分析。