多倍體同源區段二代測序生物信息學分析關鍵參數優化

武建楠，陳肯，王歡，龐鉑實，周宇荀，肖君華，李凱*

（1.東華大學化學化工與生物工程學院，上海 201620；2.上海農林職業技術學院，上海 201699）

單核苷酸多態性（Single nucleotide polymorphism,SNP） 是由單堿基突變引起的DNA 序列多態性， 作為新一代遺傳標記具有數量豐富、分布廣泛和覆蓋度高等特點[1]。 與其他分子標記相比，SNP 標記在基因分型過程可以做到高度自動化、高通量檢測[2]，已成為包括玉米在內的二倍體植物研究中的首選標記類型[3]。 2010 年滑峰等[4]發現CYP2A12基因的同源序列對SNP 分析產生嚴重影響。異源多倍體SNP 開發與鑒定的困境主要在于靶區段常存在大量同源序列， 存在SNP、部分同源序列變異（Homoeologous sequence variant,HSV）[5]和旁系同源序列變異（Paralogous sequence variant,PSV）[6]等多種變異類型。 在多倍體遺傳研究中，80%以上假陽性SNP 是HSV 和PSV 而不能用于開發遺傳標記[7]。

多倍體物種SNP 分型的常用技術包括Sanger 測序、SNP 芯片、 競爭性等位基因特異性PCR（Kompetitive allele-specific polymerase chain reaction，KASP）等。然而，它們依賴特異性擴增或特異性雜交， 而多倍體物種亞基因組間高度同源， 靶向特異性大大降低。 第二代測序（Next generation sequencing，NGS） 技術的發展使得測序費用大幅度降低，測序效率極大提高。 目前已有研究者嘗試采用生物信息學方法區分多倍體SNP 和HSV（也稱部分同源SNP），即通過不同的過濾方法[8]識別并消除部分同源SNP。 該方法在不同倍性的物種如花生（異源四倍體）、小麥（異源六倍體）和草莓（異源八倍體）的模擬數據中均表現出很高的準確性[9]，有望得到廣泛應用。

本研究的對象為異源四倍體陸地棉（Gossypium hirsutumL.），是由亞洲棉（G. arboreum，AA）和雷蒙德氏棉（G. raimondii，DD）雜交加倍而來[10]。在A12 染色體上隨機選取了3 個靶標區段，并對136 個陸地棉樣本的DNA 用多重PCR 靶向擴增子測序[11]。 將參考基因組作為關鍵參數進行調整優化，對其分型結果進行比較，研究同源區段對多倍體SNP 鑒定與生物信息學分析的影響。 優化生物信息學方案后獲得了準確的分型信息，提高了SNP 分型的準確度，為棉花等多倍體作物全基因組SNP 基因型分析以及定制育種打包檢測奠定了基礎。

1 材料與方法

1.1 植物材料與DNA 提取

參考序列來源于陸地棉參考基因組Ghir-NAU （https://cottonfgd.org/about/download.html）。在136 個不同陸地棉樣本中，采用匡猛等[12]的方法快速提取DNA，DNA 的OD260/OD280比值為1.8～2.0， 通過瓊脂糖凝膠電泳檢測DNA質量。

1.2 文庫制備與測序

三個靶向擴增子區段均位于陸地棉A12 染色體（表1）。 利用本實驗室前期設計的靶向測序的多重PCR 引物系統TASPrimer[13]，對潛在變異位點設計特異性PCR 引物（表2），由上海百力格生物技術有限公司合成引物。PCR 反應的產物大小為215～237 bp，引物長度為19～30 bp，退火溫度為57～63 ℃，GC 含量為28%～63%， 并設計了136 對含有索引序列和通用序列的條形碼引物來區分不同樣本， 使用Illumina 公司的P5、P7 建庫引物和多重PCR 法[14]建庫。將不同樣本混合之后，使用瓊脂糖凝膠對混合文庫進行電泳分離，切膠回收純化文庫。 建庫產物送金唯智生物（蘇州金唯智生物科技有限公司）進行高通量測序（測序儀為Illumina X-10），上機前用Agilent 2100 Bioanalyzer 對該測序文庫進行質量控制。

表1 靶標區段及潛在SNP 變異位點基本信息Table 1 Basic information of target segment and potential SNP variant sites

1.3 生物信息分析流程

基于已標注SNP 位點的參考基因組、原始數據和棉花基因組數據庫CottonFGD[15]，對關鍵參數即參考基因組（Reference）調整優化進行生物信息學分析， 流程如圖1 所示。 原始數據含136個樣本， 用FastQC 軟件進行質量控制，用Cutadapt 軟件切去接頭序列，最終得到有效數據樣本130 個。

圖1 生物信息分析流程Fig. 1 Bioinformatics analysis process

采用BWA 軟件[16]與參考基因組進行比對得到SAM 文件，用自定義腳本對匹配的讀長（Read）數進行統計生成比對數據表，并用SAMtools[17]轉化為BAM 文件。 使用SAMtools mpileup 命令對目標SNP 位點進行變異分析，并分別計算比對到區段1 、區段2 和區段3 的每個樣本對應讀長的r值，生成散點圖。r＝Nv/Nt，其中Nv是變異讀長數，Nt是比對到參考基因組的讀長總數。r值是判斷雜合性的指標[18]，若用R 表示參考堿基、A 表示變異堿基，那么：當0＜r≤0.2 時，表示某一樣本中該位點為純合基因型且與參考序列一致即RR；當0.8≤r＜1 時，表示該位點為純合變異基因型AA；r為其他數值時， 表示此樣本該位點為雜合基因型RA。

再用GATK[19]對區段中的所有位點進行變異分析，最后生成報表檢驗基因型。 GATK 報表為vcf 格式，其中0/0 或1/1 均表示純合，分別與參考序列或變異序列一致；0/1 表示雜合。

在棉花數據庫（https://cottonfgd.org/）對靶標區段進行比對，確定高相似度的同源序列。 將靶標區段與同源區段分別作為參考序列進行優化調整分析。

2 結果與分析

2.1 基于靶標區段作為參考序列的分析

樣本總讀長數為11 031 472， 比對到3 個靶標區段的讀長相對接近， 共約占總測序讀長的70%，3 個擴增子覆蓋率為100%。

表1 中3 個區段5 個潛在SNP 位點的有效測序深度均值分別為8 400×、8 534×、8 403×、7 716×、12 793×。 利用SAMtools 對3 個靶標區段所有目標SNP 位點進行變異分析， 結果顯示：區段1 和區段2 的目標SNP 位點的r值接近1或0，說明基因型純合；區段3 的143 位點的r值大多接近0.5，表明基因型多為雜合（圖2）。

圖2 3 個靶標區段目標SNP 位點等位基因的讀長所占比例Fig. 2 Proportion of the target SNP locus alleles reads in the three target segments

用GATK 進行全片段SNP 檢出（表3），結果表明：區段1 的3 個潛在SNP 位點和區段2 的1個潛在SNP 位點在130 個樣本中100%為純合子；而區段3 的潛在143 位點在130 個樣本中有近99%為雜合，只在1 個樣本中為純合；此外，在區段3 上還檢出了5 個新的變異位點（位點48、78、174、180、182）且均為雜合子。

表3 靶標區段作為參考序列的變異檢測與分型結果Table 3 Target segment as a reference sequence for variance detection and typing results

2.2 關鍵參數優化與分析

根據棉花數據庫BLAST 結果， 靶標區段3 的同源序列位于D12 染色體47 602 597～47 602 836 bp， 兩者相似性達96.28%（圖3），BWA 軟件比對時無法區分區段3 及其同源序列（同源區段3）， 因此考慮將區段3 和同源區段3同時作為參考序列。

圖3 區段3 與同源區段3 的BLAST 比對情況Fig. 3 BLAST comparison of segment 3 and homologous 3

首先， 將同源區段3 單獨作為參考序列，其潛在SNP 位點分析結果與區段3 作為參考序列時的分析結果一致。然后，將區段3 與同源區段3同時作為參考序列分析，比對到區段3 與同源區段3 的讀長總數為3 462 670，其中約48%比對到區段3，約52%比對到其同源區段。

隨機選取1 個樣本，從以區段3 作為參考序列分析的SAM 文件中， 隨機選取比對到區段3的部分讀長，再從以區段3 與其同源區段3 同時作為參考序列分析的SAM 文件中提取相同數量、 相同序列號的讀長進行比對和比較可以看到，只有靶標區段作為參考序列時得到4 個雜合位點（圖4A）；在參考序列中加入同源區段3 后（圖4B），不同亞基因組的讀長分開，能夠直觀判斷出SNP 和亞基因組內的多態性即HSV。

圖4 不同參考序列部分對比結果比較Fig. 4 Comparison of mapping reads with different reference sequences

利用SAMtools 對區段3 和同源區段3 中目標SNP 位點進行變異分析，區段3 的143 位點r值接近1 即基因型為純合TT， 與GATK 報表中143 位點基因型一致（表4）；而其同源序列對應位點141 位點的r值均接近0 即基因型為純合AA。

比較使用區段3、同源區段3、區段3 加同源區段3 作為參考序列的GATK 結果：在靶標區段3 單獨作為參考序列時， 區段3 的143 位點在大部分材料中被鑒定為雜合TA， 在個別材料中被鑒定為AA， 而把靶標區段與其同源區段同時作為參考序列時，143 位點被正確地鑒定為AA 或TT； 區段3 的另外3 個SNP 其實是亞基因組間的HSV；48 位點SNP 和78 位點SNP 是2 個新的SNP（表4）。

表4 靶標區段和同源區段同時作為參考序列時變異檢測與分型結果Table 4 Variation detection and typing results with the target and segment and homologous segment as reference sequences

3 討論

陸地棉A 亞基因組與D 亞基因組具有極高的同源性[20]，存在大量的同源區段干擾SNP 分型。正如Kaur 等[21]所述，在異源多倍體植物中，同源區段間SNP、HSV 和PSV 都有可能存在，比對時由于高度相似的同源位點映射錯誤，SNP 準確分型較為困難， 而且會增加SNP 的假陽性率，其中關注最多的是區分SNP 與HSV[22]。 PSV（本研究未涉及）大多出現在基因家族成員中，由于橫向同源基因的干擾， 使開發SNP 難度大大增加，如陸地棉GhNAC家族[23]。

多倍體植物的分子育種需要精確分型，但常規的技術手段難以區分同源區段[7,24]。基于Sanger測序成本高、通量低，且只能由套峰判斷雜合性，而目前最流行的高通量基因分型平臺是基于雜交的SNP 陣列和各種使用NGS 的基因分型[25]。特定技術方案有可能解決同源干擾問題，如李杏瑜等[26]利用花生的近緣二倍體和高分辨率溶解曲線（High resolution melting，HRM）方法區分SNP 和HSV。Oliver 等[27]利用生物信息學分析從燕麥（異源六倍體）127 000 個重疊群中篩除HSV，獲得了9 448 個候選SNP。 Byers 等[28]嘗試在棉花中僅擴增帶有SNP 的亞基因組， 雖然成功率低于預期，仍在陸地棉中鑒定出11 834 個SNP。 總之，對不同多倍體的基因分型有不同的最佳策略[29]。 但這些方法或是基于雜交而缺少特異性，或是只適用于沒有參考序列的多倍體， 因此需要開發新策略，以高通量方式區分多倍體中部分同源變異與真正等位基因SNP。結合多重PCR 靶向測序[11]與生物信息學方案[14,30]，本研究通過對陸地棉樣本的常規分析與參考序列調整優化，證明了多倍體中同源區段的存在會影響生物信息學分析與SNP 基因型鑒定，表明多倍體不同亞基因組SNP分型需要個性化的生物信息學方案。

本研究對靶標區段3 潛在SNP 位點常規分析的分型結果為雜合， 但從陸地棉遺傳學角度，出現幾乎全為雜合子的基因型分布是不合理的；結合多倍體特性，我們猜測并證明了是同源區段干擾導致的這一結果。 由此，本文提出把同源區段與靶標區段同時作為參考序列的方法，區分測序數據中來自不同亞基因組的讀長，不僅能得到正確基因型，還能排除HSV 造成的假陽性SNP。區段3 和同源區段3 的分析結果也正好與Clevenger 等[8]提出的“如果亞基因組相差3%左右，則僅允許差異小于3%的reads 比對到正確的亞基因組上”這一結論吻合。

4 結論

通過生物信息學方案優化可有效區分不同亞基因組讀長并準確鑒定SNP 的基因型， 提高SNP 分型的準確度， 有助于多倍體品種的檢測，同時在多倍體育種方面可以提供DNA 水平上的信息。

致謝：

感謝華中農業大學林忠旭教授提供寶貴的植物材料。