Introduction

隨著高通量測(cè)序技術(shù)的發(fā)展，這些年關(guān)于CNV（拷貝數(shù)變異）的研究也逐漸多了起來。以往CNV主要通過PCR或者CGH芯片等技術(shù)進(jìn)行查找，現(xiàn)在由于二代測(cè)序的價(jià)格飛速下降，通過將二代測(cè)序產(chǎn)生的reads 回帖（map）到參考基因組上來檢測(cè)CNV的流程和軟件逐漸建立并且流行了起來。本文主要介紹得到CNV后如何使用Plink進(jìn)行下游的GWAS（全基因組關(guān)聯(lián)分析）。

軟件版本

Plink v1.07
Plink最新的版本是1.9(16 May 2016)，但是由于我在Plink官網(wǎng)上查到的CNV關(guān)聯(lián)分析流程中幾個(gè)參數(shù)在新版中被移除了（？！?。。Ｄ壳斑€沒搞明白在新版中是怎么做的，所以就先用舊的Plink做一版。

原始數(shù)據(jù)

我拿到的CNV數(shù)據(jù)是用我們實(shí)驗(yàn)室自己開發(fā)的軟件CNVcaller，從BAM文件中得到的。軟件暫時(shí)還沒發(fā)表，等發(fā)表了我再來安利一發(fā)。
此外，還需要你要與之關(guān)聯(lián)的表型數(shù)據(jù)（phenotype），數(shù)量性狀或者質(zhì)量性狀等等都可以。

使用Plink進(jìn)行分析

1.創(chuàng)建CNV list

把你的CNV按照Plink的要求進(jìn)行整理。格式如下：

文件總共有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

注：如果是自然群體，沒有家系關(guān)系的話，F(xiàn)ID和IID起一樣的名字就可以了。而SCORE和SITES兩項(xiàng)不會(huì)直接用于關(guān)聯(lián)分析，所以不知道就直接全部寫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.創(chuàng)建FAM file

.fam文件儲(chǔ)存的是樣本的性狀、性別等信息。
這個(gè)文件沒有header，每行一個(gè)樣本，六列信息，分別是：

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)

示例如下：

注1：如果沒有family ID, parental ID, sex, and/or phenotype columns也行，使用下面的參數(shù)就OK了（適用于.fam和.ped文件）。
--no-fid
--no-parents
--no-sex
--no-pheno

注2：第六列的表型值，如果寫的不是1、2、-9、0中的任意一個(gè)，則會(huì)被判定為質(zhì)量性狀。

此部分詳見這里

3.生成MAP文件

這一步是用之前制作的CNV list，使用Plink自動(dòng)生成.map文件。CNV的map文件和plink通用的map文件一樣。
第一列是Chromosome code。
第二列是Variant identifier，是ID加上CNV的起始或者終止位置。
第三列是Position in morgans or centimorgans，在CNV的map文件中都用0來替代。
第四列是該標(biāo)記在染色體上的位置Base-pair coordinate。
命令如下：

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

其中，plink.cnv就是第一步整理出來的CNV list。運(yùn)行此命令后會(huì)自動(dòng)生成plink.cnv.map

4.關(guān)聯(lián)分析

有了前面三個(gè)文件，我們就可以開始做關(guān)聯(lián)分析了。此部分以數(shù)量性狀為例。
命令如下：

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

如果你前面那三個(gè)文件名前綴統(tǒng)一的話（a.cnv, a.cnv.map, a.fam），你也可以用這條命令：

plink --cfile a

之后你會(huì)得到一堆文件，其中.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文件，里面包含置換檢驗(yàn)的P值：EMP1和EMP2。

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

注1：--noweb是說Skipping web check。--mperm是置換檢驗(yàn)。
注2：QT mean應(yīng)該是數(shù)量性狀的均值的意思。
注3：質(zhì)量性狀默認(rèn)使用單側(cè)檢驗(yàn)，數(shù)量性狀默認(rèn)使用雙側(cè)檢驗(yàn)。如果想使用單側(cè)檢驗(yàn)，添加參數(shù)：--cnv-test-1sided

結(jié)語(yǔ)

至此，使用Plink對(duì)CNV做GWAS分析的第一部分就講完了。我們?cè)?mperm中可以得到想要的結(jié)果，包括CNV和對(duì)應(yīng)的P值。
上面的分析在plink軟件中是屬于Rare copy number variant (CNV)
分析，而其中的好多參數(shù)在plink 1.9中被移除了，而且在這里我也沒找到如何導(dǎo)入?yún)f(xié)變量以消除群體結(jié)構(gòu)的影響，GWAS的模型也不可選。
因此，下一章我們將學(xué)習(xí)Common copy number polymorphism (CNP) data的GWAS分析。