甘秋云

(福州理工學(xué)院應(yīng)用科學(xué)與工程學(xué)院，福建 福州 350014)

1 引言

在基因組中，單核苷酸多態(tài)性SNP(Single Nucleotide Polymorphism)是由單個核苷酸的變異產(chǎn)生的DNA序列多態(tài)性，是生物群體中最為常見的一種可遺傳變異。研究發(fā)現(xiàn)，在人類基因組中，SNP數(shù)量可達(dá)300萬個甚至更多，幾乎每1 000個堿基對就會有1個SNP發(fā)生。SNP作為重要的分子標(biāo)記在遺傳研究中具有重要的研究價值，它為我們提供了一個強有力的工具，為進一步研究不同生物個體之間遺傳差異，分析影響不同個體性狀的生物學(xué)基礎(chǔ)、定位相關(guān)疾病及藥物反應(yīng)的基因提供了科學(xué)依據(jù)。

但是，隨著高通量測序技術(shù)的發(fā)展，DNA數(shù)據(jù)庫中的核酸序列呈幾何增長，如何將數(shù)據(jù)的積累向數(shù)據(jù)分析進行轉(zhuǎn)變，成為一門新的課題。海量的生物信息數(shù)據(jù)不僅促進了生物信息學(xué)的發(fā)展，同時也給數(shù)據(jù)挖掘提出了新的挑戰(zhàn)。因此，使用計算機工具，結(jié)合統(tǒng)計學(xué)和數(shù)學(xué)的研究方法，將計算機算法應(yīng)用于生物大數(shù)據(jù)分析，深入理解DNA序列的結(jié)構(gòu)和生物學(xué)功能，是生物信息學(xué)、生物統(tǒng)計學(xué)和計算機科學(xué)等多門學(xué)科的綜合運用，它體現(xiàn)了當(dāng)代科技的進步[1]。

本文通過分析ED(Euclidean Distance)算法和SNP-index算法，以擬南芥為例，通過DNA測序獲取序列讀段(reads)，基于 Linux平臺，對測序數(shù)據(jù)進行過濾、篩選，使用SOAPaligner短序列比對軟件對序列進行比對，最后通過算法實現(xiàn)，計算在不同算法下的SNP位點數(shù)量及各個SNP基因型比例。

2 SNP計算算法

當(dāng)今較為常用的檢測SNP軟件包括soapSNP[2]、GATK(Genome Analysis ToolKit)[3]、samtools(sequence alignment/map tools)[4]和bcftools(utilities for variant calling and manipulating VCFs and BCFs)[5]等。目前用于計算子代混池中SNP的主要有ED算法和SNP-index算法。

2.1 ED算法

ED算法在無親本的情況下即可定位SNP，具有很強的去除背景噪聲的能力。它的主要思想是通過計算不同混池中各個突變型的基因頻率距離，即ED值，利用距離差異來反映標(biāo)記與目標(biāo)區(qū)域的連鎖強度。ED值的計算公式如式(1)所示：

ED=(Amut-Awt)2+(Cmut-Cwt)2+

(Gmut-Gwt)2+(Tmut-Twt)2)1/2

(1)

其中，A，C，G，T為4種堿基符號，mut與wt分別表示突變型和野生型。根據(jù)序列比對結(jié)果(包括SNP位點集合測序深度)，計算當(dāng)前位點中各突變型堿基占測序讀段(reads)的比例，通過式(1)計算混池間的突變型基因頻率差異，即ED值。ED值計算方法示意圖如圖1所示。

Figure 1 Schematic diagram of ED algorithm calculating gene frequency distance

由于SNP位點集中可能存在SNP位點的偏差，出現(xiàn)假陽性位點，因此需要對計算所得的ED值進行擬合。按照ED值從小到大排序，將前1%的ED值作為篩選SNP位點的閾值[6]。

2.2 SNP-index算法

SNP-index算法以某一親本或參考基因組為參考序列，統(tǒng)計在某位點上含有突變基因型的reads占總reads的比例，該算法需要親本或參考序列。該算法的計算方法如式(2)和式(3)所示：

SNP-index(ab)=Mab/(Pab+Mab)

(2)

SNP-index(aa)=Maa/(Paa+Maa)

(3)

其中，Mab表示ab群體來自于突變型的深度；Pab表示ab群體來自于父本的深度；Maa表示aa群體來自于突變型的深度；Paa表示aa群體來自于父本的深度。圖2所示為SNP-index算法計算基因頻率示意圖。

Figure 2 Schematic diagram of SNP-index algorithm to calculate gene frequency

為了尋找混池之間基因型頻率的差異，需要通過計算混池間的基因型頻率差異進行標(biāo)記關(guān)聯(lián)分析，用Δ(SNP-index)表示，若Δ(SNP-index)值越接近1，說明標(biāo)記與目標(biāo)區(qū)域的關(guān)聯(lián)性越強，反之越弱。Δ(SNP-index)的計算方法如式(4)所示：

Δ(SNP-index)=

SNP-index(ab)-SNP-index(aa)

(4)

SNP-index算法需要親本的測序序列信息，這樣可以排除2個親本相對參考基因組共有的SNP，即去除背景噪聲。此外，親本檢測出來的SNP是和目標(biāo)性狀直接對應(yīng)的，這樣可以排除一些假陽性的SNP位點，即Δ(SNP-index)趨于1但是并非與目標(biāo)性狀有關(guān)聯(lián)的標(biāo)記位點[7]。

3 實驗過程

3.1 前期數(shù)據(jù)處理

(1)數(shù)據(jù)來源。

實驗前期主要通過對擬南芥F1代野生型與突變型進行雜交，在F2代群體中分別構(gòu)建突變型與野生型DNA池，并對混池樣品進行深度測序，最終獲得平均長度為76 bp的測序序列讀段(reads)共29 264 012條。

Figure 3 Format of sequence reads

Figure 4 Format of sequence alignment results

(2)過濾、篩選。

測序獲得的野生型與突變型2個樣本的測序序列讀段格式如圖3所示。從圖3中可知序列最后一列以0和1標(biāo)記，代表測序質(zhì)量。為了保證后期檢測SNP位點的reads可靠，實驗中剔除了測序質(zhì)量較低的reads，即測序結(jié)果中最后一列以0為標(biāo)記的讀段。

(3)比對。

在Linux系統(tǒng)下，使用SOAPaligner軟件進行序列比對。參考序列從擬南芥數(shù)據(jù)庫TAIR(http://www.arabidopsis.org/)中得到，將比對結(jié)果文件按染色體及染色體的ID號進行排序，根據(jù)比對結(jié)果文件提取有效數(shù)據(jù)進行算法實現(xiàn)。比對結(jié)果文件格式如圖4所示。

圖4所示的比對結(jié)果最后幾列信息中，“chr1”表示1號染色體；“6”表示reads第1個堿基比對到染色體上的位置；“C→33G4” 表示在參考序列染色體的第34位(0開始算第1位)有1個錯配堿基，即參考序列上的堿基為C，而比對序列對應(yīng)位置的堿基是G；“44M”表示匹配的數(shù)目為44個堿基數(shù)； “33C10”表示比對序列的前33個堿基是匹配的，第34個堿基是錯配的，參考序列上的堿基應(yīng)為C，隨后又有10個堿基是匹配的。

(4)比對結(jié)果處理。

將比對的結(jié)果文件作為SNP查找的文件來源。因此，根據(jù)比對結(jié)果，主要獲取堿基匹配信息。本文實驗中主要獲取染色體號(每條染色體號對應(yīng)一個結(jié)果文件)、堿基在染色體上的位置、參考堿基和堿基配對信息(包括正確配對和錯配)。例如(3)中比對結(jié)果中獲取得到的數(shù)據(jù)為“Chr1 39 C G”，表示在1號染色體的第39個(reads第1個堿基比對到染色體上的位置是6，錯配堿基在reads上的位置為33，最終定位到染色體上的位置為39)位點上參考堿基為C，錯配堿基為G。最后，統(tǒng)計同一條染色體同一位置上的堿基配對情況，同時記錄當(dāng)前位置的覆蓋深度(這里主要指短序列比對到參考序列的數(shù)量，在本文實驗中以染色體同一個位置的記錄數(shù)量作為覆蓋深度)。

圖5a是1號染色體在第39個位點處的堿基配對情況，第4列為參考序列堿基，第5列為配對堿基。對配對結(jié)果文件進行統(tǒng)計，最終獲取當(dāng)前位點的堿基配對數(shù)目及錯配堿基數(shù)目，結(jié)果如圖5b所示，其中第3列表示參考堿基，第4列表示從圖5a結(jié)果中統(tǒng)計得到的與參考堿基配對上的堿基數(shù)為10，即覆蓋深度為10×，第5列和第6列表示該位點上的錯配堿基為G，數(shù)目為4。

Figure 5 Base pairing on chromosomal sites and statistical results

3.2 實現(xiàn)過程

為了更好地比較2種算法的SNP計算結(jié)果，實驗中將位點堿基配對結(jié)果分別使用ED和SNP-index 2種算法進行實現(xiàn)。在ED算法中，分別計算野生型和突變型的4種堿基的基因頻率，根據(jù)式(1)計算ED值。在SNP-index算法中，根據(jù)式(2)和式(3)計算對應(yīng)位點SNP-index值。圖6為計算結(jié)果。

Figure 6 SNP calculation results

圖6中，第1列為染色體號;第2列表示當(dāng)前計算的ED值或SNP-index結(jié)果;第3列和第4列用于判斷潛在的SNP位點；若2列的堿基型不相同，則認(rèn)為當(dāng)前位點為潛在的SNP位點;第5列為當(dāng)前測序深度;最后1列表示染色體位點位置。由于可能存在測序錯誤和比對錯誤等其他因素，導(dǎo)致一些不可靠的多態(tài)性位點存在，實驗中對測序深度值進行了篩選，主要過濾了深度小于或等于2或大于或等于50的SNP位點，以便提高SNP檢測的準(zhǔn)確性。

由于某些SNP位點可能存在偏差，設(shè)定一個閾值作為篩選可靠SNP位點的標(biāo)準(zhǔn)。根據(jù)計算結(jié)果，將前1%的ED結(jié)果作為篩選SNP位點的閾值，統(tǒng)計在ED算法下各個染色體的SNP位點數(shù)量，并輸出對應(yīng)位點堿基符號的差異。在實現(xiàn)SNP-index算法前，需要計算對應(yīng)位點的SNP- index值，實驗中將基因頻率低于0.3的位點進行剔除。但是，由于測序深度有限以及測序的隨機性帶來的波動，對單個SNP計算SNP-index值往往不穩(wěn)定。為了降低不可靠的SNP位點帶來的偏差，本文實驗主要計算某個區(qū)間內(nèi)的SNP-index平均值，實驗中選用1 Mb的窗口為一個區(qū)間，以10 kb為步長進行滑動，采用滑動窗口的方式統(tǒng)計某個窗口區(qū)間中所有SNP的SNP-index平均值作為當(dāng)前區(qū)間的SNP-index結(jié)果[7]，同時計算各個染色體的SNP位點數(shù)量，并輸出對應(yīng)位點堿基符號的差異。

4 實驗結(jié)果分析

4.1 平均覆蓋深度

通過對不同染色體同一位點的覆蓋深度進行統(tǒng)計，分別計算2種算法下，1～5號染色體的覆蓋深度及全基因組平均覆蓋深度。從表1中可得，ED算法得到的1～5號染色體的平均覆蓋深度分別為8.6×,12.3×,10.4×,12.5×和11.6×，全基因組平均覆蓋深度為9×。SNP-index算法得到的1～5號染色體的平均覆蓋深度分別為10.2×,14.5×,13.4×,11.3×和10.8×，全基因組平均覆蓋深度為12×。圖7是2種算法得到5條染色體平均覆蓋深度圖。

Table 1 Average coverage depth of arabidopsis thaliana genome with different algorithms

Figure 7 Average coverage depth of arabidopsis thaliana genome with different algorithms

4.2 SNP位點統(tǒng)計

由于測序中存在許多不確定的因素可能造成測序錯誤的發(fā)生，最終導(dǎo)致檢測結(jié)果出現(xiàn)虛假SNP位點。為保證SNP位點數(shù)據(jù)的可靠性，在前期的數(shù)據(jù)處理中對測序讀段進行過濾、篩選，同時通過不同算法分別對樣本序列進行SNP位點檢測，最終取得全基因組范圍內(nèi)的潛在SNP位點。表2是分別基于ED算法和SNP-index算法，計算擬南芥1～5號染色體SNP位點的結(jié)果。

從表2中可以發(fā)現(xiàn)，2種算法得到的擬南芥各個染色體對應(yīng)的SNP位點接近。其中，5號染色體SNP數(shù)目最多，平均達(dá)到70 000個，其余4條染色體的SNP位點均為50 000左右。但是，ED算法得到的SNP數(shù)目相對于SNP-index算法的更多，分布更廣。通過ED算法得到的擬南芥全基因組SNP分布密度約為1.68 kb-1，通過SNP-index算法得到的擬南芥全基因組SNP分布密度約為1.62 kb-1。不論哪種算法，從結(jié)果中可以發(fā)現(xiàn)平均每1 kb即有1～2個SNP發(fā)生，圖8是2種算法下擬南芥全基因組SNP數(shù)目對比圖。

Table 2 Statistics of SNP number and distribution density of arabidopsis thaliana genome with different algorithms

Figure 8 Comparison of SNP numbers in arabidopsis thaliana genome with different algorithms

4.3 SNP基因型統(tǒng)計

常見的SNP為二等位多態(tài)性，主要包括轉(zhuǎn)換和顛換2種類型。正常堿基的配對為A-T和C-G，由于錯配導(dǎo)致的堿基錯誤可分為6種類型，分別為A/G、A/T、A/C、C/T、C/G和G/T，根據(jù)SNP篩選結(jié)果，統(tǒng)計分析這6種SNP類型在全基因組范圍內(nèi)所占的百分比。表3和表4分別是2種算法下每條染色體的SNP基因型比例，表5是不同算法下擬南芥全基因組SNP基因型的統(tǒng)計結(jié)果。

從表3和表4中可以發(fā)現(xiàn)，同一種算法下，不同染色體對應(yīng)的相同的SNP基因型比例大致相同。從表5中可以發(fā)現(xiàn),基于不同的2種算法計算得到的同一種SNP類型的發(fā)生比例無明顯差異。全基因組中，6種SNP基因型分布比例各不相同，

Table 3 SNP genotype proportion per chromosome based on ED algorithm

Table 4 SNP genotype proportion per chromosome based on SNP-index algorithm

Table 5 Distribution proportion of SNP genotypes obtained by different algorithms

A-G和C-T發(fā)生概率最大，平均比例高達(dá)28%以上，說明SNP類型中轉(zhuǎn)換發(fā)生的頻率大于顛換。在轉(zhuǎn)換中，C-T的概率略高于A-G，以C-T轉(zhuǎn)換為主，其原因可能是由于CpG的C堿基是甲基化的，最容易發(fā)生突變，容易自發(fā)脫氨基形成胸腺嘧啶T，因此CpG也變?yōu)橥蛔儫狳c[8]。其中，C/G的基因型是所有SNP類型中比例最低的，約為8.0%，這可能是因為G/C與生物DNA的許多結(jié)構(gòu)特征有密切的聯(lián)系，在基因組中的含量最重，因此基因組中C-G比例偏低[9]。雖然不同算法下得到的SNP數(shù)量存在差異，但通過實驗結(jié)果得到的SNP類型的數(shù)據(jù)統(tǒng)計結(jié)果與現(xiàn)已知的生物SNP型研究結(jié)果相近，表明當(dāng)前計算分析得到的擬南芥測序數(shù)據(jù)SNP位點檢測結(jié)果是可靠的。圖9是不同算法下擬南芥全基因組SNP基因型分布比例圖。

Figure 9 Distribution of SNP genotypes in arabidopsis thaliana genome with different algorithms

5 結(jié)束語

通過測序擬南芥野生型和突變型樣本DNA，分析統(tǒng)計基于ED算法和SNP-index算法得到的擬南芥潛在的SNP位點，結(jié)果發(fā)現(xiàn)ED算法得到的SNP位點數(shù)量更多，分布更廣，相對分布密度大于SNP-index算法的，但是2種算法得到的SNP分析結(jié)果接近。ED算法無需親本，因此當(dāng)無親本時，可以選擇ED算法進行計算。若有親本序列時，2種算法都可以選擇，但是以SNP-index算法結(jié)果為準(zhǔn)。通過本文的前期數(shù)據(jù)分析，以SNP位點為分子標(biāo)記，為后期的突變基因的定位提供了可靠的數(shù)據(jù)基礎(chǔ)和理論依據(jù)。

隨著時代的發(fā)展及人類基因組計劃的逐步實施，海量的生物學(xué)數(shù)據(jù)促進了生物信息學(xué)的發(fā)展，同時也給數(shù)據(jù)挖掘工作提出了新的課題和挑戰(zhàn)。通過計算機對生物測序數(shù)據(jù)SNP進行檢測分析，是一種新型、快速、簡便和低廉的方法，為揭示不同個體之間遺傳差異，分析影響不同的生物性狀、疾病及藥物反應(yīng)的相關(guān)基因定位提供了科學(xué)依據(jù)和技術(shù)支持。