全基因組數(shù)據(jù)分析軟件PLINK在統(tǒng)計(jì)遺傳學(xué)教學(xué)中的應(yīng)用

呂洪超，張瑞杰，姜永帥，段煉，李晉

哈爾濱醫(yī)科大學(xué)

全基因組數(shù)據(jù)分析軟件PLINK在統(tǒng)計(jì)遺傳學(xué)教學(xué)中的應(yīng)用

呂洪超，張瑞杰，姜永帥，段煉，李晉

哈爾濱醫(yī)科大學(xué)

統(tǒng)計(jì)遺傳學(xué)是一門理論實(shí)踐相結(jié)合的課程，是生物信息學(xué)專業(yè)重要的必須課程。我們?cè)谥v授理論知識(shí)的同時(shí)，注重培養(yǎng)學(xué)生的實(shí)踐操作技能，引入了PLINK等當(dāng)前國(guó)際流行的統(tǒng)計(jì)分析軟件。通過指導(dǎo)學(xué)生利用PLINK軟件對(duì)一套類風(fēng)濕性關(guān)節(jié)炎的真實(shí)案例的分析，即促進(jìn)了學(xué)生理論聯(lián)系實(shí)際，又激發(fā)了學(xué)習(xí)的積極主動(dòng)性，為將來順利開展科研訓(xùn)練打下基礎(chǔ)。

統(tǒng)計(jì)遺傳學(xué)；PLINK軟件；關(guān)聯(lián)分析

人類基因組計(jì)劃和人類基因組單體型圖計(jì)劃完成以后，隨著基因分型技術(shù)的發(fā)展與成熟，統(tǒng)計(jì)遺傳學(xué)方法已成為研究多基因復(fù)雜疾病易感性關(guān)系和基因定位的主要技術(shù)手段。在生物信息學(xué)領(lǐng)域，針對(duì)現(xiàn)代遺傳學(xué)的研究尤其是分子遺傳學(xué)研究逐漸成一門新興學(xué)科《統(tǒng)計(jì)遺傳學(xué)》。這是一門綜合運(yùn)用遺傳學(xué)與數(shù)學(xué)的理論和方法，歸納整合群體遺傳學(xué)、遺傳流行病學(xué)、數(shù)量遺傳學(xué)、生態(tài)遺傳學(xué)和分子遺傳學(xué)等分支學(xué)科內(nèi)容，是遺傳學(xué)、統(tǒng)計(jì)學(xué)、信息學(xué)高度交叉的學(xué)科。目前，經(jīng)過多年的教學(xué)探索，教研室已經(jīng)形成了比較完善的統(tǒng)計(jì)遺傳教學(xué)內(nèi)容和教學(xué)手段。現(xiàn)行的教學(xué)內(nèi)容包含了基因頻率與基因型頻率、遺傳連鎖分析與關(guān)聯(lián)分析、群體遺傳結(jié)構(gòu)與分化分析、QTL基因定位、單體型分析等內(nèi)容，闡明了遺傳變異規(guī)律及方法論，用于指導(dǎo)復(fù)雜疾病的分子分型與基因定位、單體型分析等研究實(shí)踐［1］。

統(tǒng)計(jì)遺傳學(xué)作為一門側(cè)重實(shí)際應(yīng)用的學(xué)科，實(shí)驗(yàn)教學(xué)是統(tǒng)計(jì)遺傳學(xué)課程教學(xué)的重要組成部分，也是統(tǒng)計(jì)遺傳學(xué)教學(xué)成功的重要保證。我們?cè)O(shè)置了12學(xué)時(shí)的實(shí)驗(yàn)課，結(jié)合當(dāng)前國(guó)際上最流行的統(tǒng)計(jì)遺傳學(xué)軟件，把理論課中學(xué)習(xí)到的大部分方法在計(jì)算機(jī)上得以實(shí)現(xiàn)。教學(xué)過程中涉及到的相關(guān)軟件很多，比如遺傳關(guān)聯(lián)研究的PLINK［2］，單倍體關(guān)聯(lián)研究的Haploview［3］，循證醫(yī)學(xué)相關(guān)的Revman，STATA。相對(duì)而言，遺傳關(guān)聯(lián)研究軟件PLINK應(yīng)用廣泛，功能強(qiáng)大，在統(tǒng)計(jì)遺傳學(xué)應(yīng)用中占據(jù)重要的地位，因此，筆者將重點(diǎn)介紹其在本科統(tǒng)計(jì)遺傳學(xué)教學(xué)中的應(yīng)用。

1　PLINK軟件及其功能簡(jiǎn)介

PLINK是一個(gè)免費(fèi)、開源的全基因組關(guān)聯(lián)分析工具集合，是由人類遺傳研究中心（CHGR），馬薩諸塞州總醫(yī)院（MGH），哈佛大學(xué)和麻省理工學(xué)院的Broad研究所等機(jī)構(gòu)科研人員所開發(fā)。PLINK主要針對(duì)基因型/表型數(shù)據(jù)的分析，近年內(nèi)容不斷豐富，包含了拷貝數(shù)變異（CNVs）分析、Meta分析等等。此外，軟件除了使用命令行分析之外，基于JAVA開發(fā)了圖形用戶界面gPLINK，并與單體型可視化分析軟件Haploview建立了關(guān)聯(lián)，對(duì)于后續(xù)可視化展示，結(jié)果的注釋和存儲(chǔ)都起到了有力的支撐作用。

PLINK功能強(qiáng)大，包含的功能分析模塊眾多，是統(tǒng)計(jì)遺傳學(xué)領(lǐng)域常用軟件。主要的功能模塊包括：數(shù)據(jù)處理，質(zhì)量控制的基本統(tǒng)計(jì)，群體分層分析，單位點(diǎn)的基本關(guān)聯(lián)分析（病例-對(duì)照數(shù)據(jù)），家系數(shù)據(jù)的傳遞不平衡檢驗(yàn)，多點(diǎn)連鎖分析，單倍體關(guān)聯(lián)分析，拷貝數(shù)變異分析，Meta分析等等。

2　PLINK在統(tǒng)計(jì)遺傳教學(xué)中的應(yīng)用

《統(tǒng)計(jì)遺傳學(xué)》是哈爾濱醫(yī)科大學(xué)校級(jí)精品課程，經(jīng)過教研室教師多年來教學(xué)實(shí)踐的不斷改進(jìn)，緊密契合當(dāng)前學(xué)科研究的前沿，形成了完整的內(nèi)容體系。課程的大部分理論內(nèi)容都可以通過PLINK軟件分析模塊實(shí)現(xiàn)，因此，我們把PLINK軟件作為本課程案例教學(xué)主要的分析實(shí)現(xiàn)軟件。隨著大數(shù)據(jù)時(shí)代基因分型技術(shù)的發(fā)展，高通量測(cè)序成本的降低，大量真實(shí)測(cè)序數(shù)據(jù)的獲取，都為我們開展統(tǒng)計(jì)遺傳學(xué)實(shí)踐教學(xué)提供了有力的支撐。我們獲得了威康信托基金會(huì)病例對(duì)照協(xié)會(huì)（WTCCC）和北美類風(fēng)濕性關(guān)節(jié)炎協(xié)會(huì)（NARAC）等國(guó)際組織機(jī)構(gòu)提供的多套遺傳變異數(shù)據(jù)，為同學(xué)初步利用真實(shí)數(shù)據(jù)開展科研實(shí)踐活動(dòng)提供了鍛煉的機(jī)會(huì)。結(jié)合課堂的理論內(nèi)容，通過同學(xué)對(duì)這些案例以討論或者研討的方式，分析問題和解決問題。這種教學(xué)方式易于為學(xué)生接受，能很好地克服統(tǒng)計(jì)遺傳學(xué)理論學(xué)習(xí)中的困難。下面我們以GAW16類風(fēng)濕性關(guān)節(jié)炎數(shù)據(jù)為例，對(duì)PLINK軟件在教學(xué)過程中的使用作一簡(jiǎn)單介紹。

案例：類風(fēng)濕性關(guān)節(jié)炎（RA）是是一種慢性炎癥性自身免疫疾病。我們獲取了一套NARAC提供的全基因組測(cè)序數(shù)據(jù)，總樣本數(shù)為2062個(gè)病例對(duì)照數(shù)據(jù)（868為cases、1194為controls），包含了22條常染色體的531689個(gè)SNPs［4］。針對(duì)本套數(shù)據(jù)，我們結(jié)合教學(xué)理論知識(shí)，指導(dǎo)學(xué)生完成部分的數(shù)據(jù)分析工作。現(xiàn)僅截取號(hào)染色體的一段數(shù)據(jù)為例，對(duì)關(guān)聯(lián)分析過程中的一些關(guān)鍵內(nèi)容作一簡(jiǎn)單介紹。同時(shí)，考慮到本專業(yè)課學(xué)生已經(jīng)擁有較好的編程能力，擁有C++，JAVA，R，MATLAB等課程學(xué)習(xí)經(jīng)歷，所以PLINK教學(xué)實(shí)踐采用的是命令行分析教學(xué)。

1）數(shù)據(jù)格式類型

PLINK軟件輸入數(shù)據(jù)類型主要為L(zhǎng)inkage格式純文本文件，包含：“*.ped”文件，數(shù)據(jù)文件的前六項(xiàng)固定為家系、個(gè)體、父親、母親、性別和患病狀況的編碼，此后各項(xiàng)為擬分析的單核苷酸多態(tài)性（SNPs）的基因型編碼，各變量之間以空格間隔（見下表1）；“*.map”文件，按SNP所屬染色體編碼、SNP名稱、遺傳距離、物理距離的順序排列。

命令語句：plink--file RA（假設(shè)RA為文件名，包含RA.ped和RA.map）。

2）數(shù)據(jù)質(zhì)量控制

這是對(duì)于數(shù)據(jù)的過濾過程，一般考慮的因素是基因型缺失率，哈迪溫波格平衡，統(tǒng)計(jì)結(jié)果從中得到滿足條件的數(shù)據(jù)，命令如下：

plink--file RA--mind 0.05--geno 0.05--maf 0.01--hwe0.001。

上述命令行表示，提取數(shù)據(jù)中滿足以下條件的SNP數(shù)據(jù)：個(gè)體基因型缺失率小于0.05，位點(diǎn)基因型缺失率小于0.05，位點(diǎn)次等位基因頻率大于0.05，哈迪溫波格平衡檢驗(yàn)大于0.001。分析結(jié)果會(huì)產(chǎn)生一個(gè)名稱為plink的純文本文件。

表1　RA病例對(duì)照部分?jǐn)?shù)據(jù)（部分）

表2　PLINK軟件關(guān)聯(lián)分析結(jié)果（部分）

3）基本的關(guān)聯(lián)分析

病例對(duì)照數(shù)據(jù)的關(guān)聯(lián)分析對(duì)應(yīng)的遺傳模型有多種，包括等位基因模型，顯性效應(yīng)模型，隱形效應(yīng)模型，加性模型等。假如我們研究的位點(diǎn)有C和A兩個(gè)等位基因，其中CC表示野生基因型，CA表示雜合型突變基因型，AA表示純合型突變基因型，則相應(yīng)的等位基因模型：A VS C；隱性模型：AA VS（CA+CC）；顯性模型：（CA+ AA）VS CC；加性模型：AA VS CC。比如我們執(zhí)行等位基因模型分析，命令為：plink--file RA-assoc。

運(yùn)行將生成結(jié)果文件“plink.assoc”，部分結(jié)果如表2。

上述結(jié)果可以看出rs17806289卡方檢驗(yàn)的p值為5.87E-10，小于0.001，具有統(tǒng)計(jì)學(xué)意義，并且優(yōu)勢(shì)比OR=0.43，顯示該位點(diǎn)等位基因A可能是一個(gè)保護(hù)因素。

當(dāng)然，如果考慮某一位點(diǎn)SNP等位基因頻數(shù)不滿足列聯(lián)表檢驗(yàn)的要求（頻數(shù)小于5），可以采用Fisher精確檢驗(yàn)的方法，選擇PLINK命令參數(shù)“--fisher”即可。更進(jìn)一步，如果考慮分析其他遺傳模型，可以選擇PLINK命令參數(shù)“--model”。為鞏固理論知識(shí)，加強(qiáng)學(xué)生學(xué)習(xí)的興趣，這兩部分內(nèi)容均以課堂作業(yè)形式交給學(xué)生以討論組的形式完成。很好的激發(fā)了學(xué)生的學(xué)習(xí)熱情，加強(qiáng)了動(dòng)手能力。

3　理論實(shí)踐結(jié)合的思考

《統(tǒng)計(jì)遺傳學(xué)》是一門理論與實(shí)踐結(jié)合非常緊密的課程，也是一門與當(dāng)前科學(xué)研究前沿緊密結(jié)合的課程。針對(duì)理論教學(xué)中抽象的統(tǒng)計(jì)原理，復(fù)雜的推導(dǎo)以及繁瑣的運(yùn)算，我們通過實(shí)驗(yàn)案例教學(xué)，一方面，加深學(xué)生對(duì)課堂學(xué)習(xí)的印象，能熟練掌握軟件分析工具，增強(qiáng)動(dòng)手操作能力；另一方面，我們采用真實(shí)案例數(shù)據(jù)，在很大程度上也很好的調(diào)動(dòng)了學(xué)生學(xué)習(xí)的積極性，促使他們主動(dòng)地查閱相關(guān)文獻(xiàn)資料，深化學(xué)習(xí)的內(nèi)容。此外，我們鼓勵(lì)學(xué)生參與教師的科研活動(dòng)，運(yùn)用所學(xué)統(tǒng)計(jì)遺傳學(xué)知識(shí)解決實(shí)際問題。在今后教學(xué)中，我們將繼續(xù)進(jìn)行教學(xué)方法的創(chuàng)新與改進(jìn)，引導(dǎo)學(xué)生自主的利用所學(xué)知識(shí)解決實(shí)際問題，為將來開展科研訓(xùn)練打下堅(jiān)實(shí)的基礎(chǔ)。

［1］張明明，張瑞杰，李晉，呂洪超.統(tǒng)計(jì)遺傳學(xué)精品課建設(shè)探索與實(shí)踐［J］.黑龍江科技信息，2010（24）∶185-185.

［2］Purcell S1，Neale B，Todd-Brown K，et al.PLINK∶a tool set for whole-genome association and population-based linkage analyses［J］.Am J Hum Genet.2007 Sep；81（3）∶559-75.

［3］Barrett JC，F(xiàn)ry B，Maller J，Daly MJ.Haploview∶analysis and visualization of LD and haplotype maps［J］.Bioinformatics.2005 Jan 15；21（2）∶263-5.

［4］Cupples LA，Beyene J，Bickeb?ller H，et al.Genetic Analysis Workshop 16∶Strategies for genome-wide association study analyses. BMC Proc.2009 Dec 15；3 Suppl 7∶S1.