MGISEQ-2000、HiSeq 2000 與NovaSeq 6000平臺全基因組重測序數據的比較分析

李偉寧，劉 剛，周 榮，唐中林，劉劍鋒*，孫飛舟*

（1.中國農業大學動物科學技術學院，北京 100193；2.全國畜牧總站畜禽遺傳資源保存利用中心，北京 100193；3.中國農業科學院北京畜牧獸醫研究所，北京 100193；4.中國農業科學院農業基因組研究所，廣東深圳 518120）

全基因組重測序（以下簡稱“重測序”）被廣泛用于變異檢測[1]、遺傳成分鑒定[2]和多態性分析[3]等研究，現已成為疾病預測及診斷[4]、動植物分子育種[5-6]等領域最常用的分析方法之一。目前最常用的二代測序平臺是Illumina（美國）測序平臺，其在基因測序市場所占份額近75%。HiSeq 2000 是Illumina 在2010 年發布的一款測序儀，其將人類基因組測序費用降至1 萬美元以下，大量生物公司和科研機構均采購了該測序儀，目前它仍是市場上的主流測序儀。NovaSeq 6000是Illumina 在2017 年發布的一款被譽為里程碑式產品的測序儀，可以搭配4 種不同的流動槽靈活地開展不同通量要求的測序任務，有望將基因組的測序費用進一步降至100 美元。國內的測序平臺研發仍在起步階段，2014 年6 月深圳華大智造科技股份有限公司（簡稱“華大智造”，英文名稱為MGI）推出了BGISEQ-1000 和BGISEQ-100 2 個二代測序平臺，是國家食品藥品監督管理總局首次批準注冊的第二代基因測序診斷產品，隨后幾年華大智造也陸續發布了多款適用于不同場景的二代測序儀。MGISEQ-2000[7]是華大智造在2017 年9 月推出的一款產品，該平臺采用了該公司自主研發的CoolMPS[8]高通量測序試劑和DNA 納米球測序技術，可在測序過程中實現高準確性、低重復序列率和低標簽跳躍率，其憑借卓越的性能表現及超高性價比在眾多測序平臺中脫穎而出。

相關研究比較了Illumina 與華大智造二代測序平臺的性能表現。Huang 等[9]研究發現BGISEQ-500的比對質量要好于HiSeq 2500，且二者的SNP 檢測準確性均在99% 以上。Korostin 等[10]對比分析了MGISEQ-2000 和HiSeq 2500 的測序數據，發現在原始數據質量、變異檢測方面二者表現相似，但MGISEQ-2000 的比對質量要優于HiSeq 2500。但上述研究的樣本量較少，平臺也較為單一。MGISEQ-2000和NovaSeq 6000 作為華大智造和Illumina 同期發布的2 款測序儀在性能表現上是否有明顯差異，NovaSeq 6000 與Illumina 早期發布的HiSeq 2000 相比性能又是否有明顯提升，還需要通過實際數據進一步驗證。本研究基于MGISEQ-2000、HiSeq 2000 和NovaSeq 6000 3個平臺對同一批樣品進行重測序分析，比較不同平臺的性能表現和測序穩定性，為研究者在選擇不同測序平臺時提供參考。

1 材料與方法

1.1 實驗動物及樣品 本研究實驗動物為9 頭公豬（4頭馬身豬、5 頭大河豬），采集所有實驗動物的耳組織，浸沒于75%的酒精，置于-80℃冰箱保存備用。采用天根生物科技（北京）有限公司生產的組織基因組DNA提取試劑盒，嚴格按照產品說明書提取實驗樣本耳組織的全基因組DNA，采用酶標儀測定待實驗樣本的基因組DNA 濃度與純度，檢測合格后進行后續文庫構建。

1.2 全基因組重測序及SNP 芯片分型 將9 頭公豬的DNA 樣本各自均分為3 份，按照MGISEQ-2000、HiSeq 2000 和NovaSeq 6000 3 個測序平臺的標準建庫流程分別對每份樣本進行建庫，在3 個平臺上均采用paired-end、PE150（即雙端150 bp 讀長）對樣本進行全基因組重測序，測序深度大于20X。全基因組SNP芯片在基因分型上具有很高的準確性[11]，通常作為評價測序數據SNP 檢測準確性的金標準[12-13]。本研究用豬50K 芯片（參考基因組版本為Sscrofa10.2）對所有樣本進行了SNP 分型，作為評價重測序數據SNP 檢測準確性的依據。

1.3 序列比對及變異檢測 本研究分析步驟及所用軟件如圖1 所示，括號中為該步驟所用軟件。3 個平臺測序獲得的reads 即為原始數據，將各平臺各樣本的兩個文庫數據合并后進行后續分析。原始數據質量由fastp[14]統計，質控和adapter 接頭去除由Trim Galore[15]（0.6.1）執行，所用參數為“--stringency 3 --length 20 -e 0.1”，序列比對由BWA[16]（0.7.17）的mem 算法處理，所用參數為“-t 6 -M–R "@RG ID:id LB:id PL:ILLUMINA SM:id"”（其中id 為自定義的樣本編號），參考基因組版本為Sscrofa11.1。比對后的sam文件用GATK[17]（4.0.12.0）的ReorderSam 排序后由Samtools[18]（1.9）轉為二進制格式bam 文件，隨后用GATK 的SortSam、MarkDuplicates 依次進行排序和重復reads 標記工作，重復reads 只進行標記而不剔除（--REMOVE_DUPLICATES false），隨后用GATK的BaseRecalibrator 獲取堿基質量校正的校準表文件，“--known-sites”所用dbsnp 庫版本為150，另一參數為“--bqsr-baq-gap-open-penalty 30”，最后用GATK的ApplyBQSR 利用上述校準表文件對堿基質量進行校正，獲得最終的bam 文件。根據此bam 文件評價比對質量，Samtools 的flagstat 模塊用于統計雙端reads 比對率，插入片段由Picard[19]的CollectInsertSizeMetrics統計，平均比對深度和位點覆蓋超過10X、20X 的比例由Mosdepth[20]輸出。以上所列參數外的其余參數均為相應軟件的默認參數。

圖1 數據分析流程及所用軟件

將最終獲得的bam 文件用GATK 的Haplotype Caller 進行個體水平的變異檢測，使用參數“-ERC GVCF”得到各平臺各樣本的gvcf 文件，然后分別將各平臺各樣本的gvcf 文件用GATK 的CombineGVCFs合并，獲得3 個平臺各包含9 個樣本的vcf 文件，再用GATK 的GenotypeGVCFs 分別基于3 個平臺各自的vcf 文件進行群體水平的變異檢測，得到3 個平臺各自的單個vcf 文件，最后用GATK 對檢測得到的SNP 和INDEL 執行過濾操作。SNP 的過濾參數為“QD<2.0||M Q<40.0||FS>60.0||SOR>3.0||MQRankSum<-12.5||ReadPos RankSum<-8.0”，INDEL 的過濾參數為“QD<2.0||F S>200.0||SOR>10.0||MQRankSum<-12.5||ReadPosRank Sum<-8.0”。使用VCFtools[21]（0.1.16）統計各平臺位點數及三者共有位點所占比例，同時用SnpSift[22]對SNP 位點進行注釋（本研究中所使用的dbsnp 庫版本均為150），然后統計各平臺檢測的SNP 位點中被dbsnp庫收錄的位點所占比例。

不同平臺結果之間的差異顯著性采用SPSS 20.0 的配對樣本t檢驗進行檢驗，統計檢驗的顯著水平（雙側）設為P<0.05，2 個變量之間的相關系數由Excel 2016中的CORREL 函數計算得出（皮爾遜相關系數）。

1.4 SNP 分型及準確性評價 將50K 芯片（Sscrofa10.2版本）位點坐標在UCSC 數據庫轉為Sscrofa11.1 版本。Sscrofa11.1 和Sscrofa10.2 中的一些DNA 序列片段是反向互補的關系，找出50K 芯片中位于這些序列上的位點，根據堿基互補原則將其轉換，使得測序數據檢測SNP 與基因芯片SNP 處于同一條鏈。隨后指定dbsnp庫位點的參考堿基（ref）和突變堿基（alt）為基因芯片位點的ref 及alt，用plink[23]（1.90）將plink 格式的芯片位點文件轉為vcf 格式文件。將9 個樣本分型后只存在缺失（./.）和野生純合（0/0）2 種類型的位點進行剔除。此處的“野生純合”是指2 個等位基因均與參考基因組堿基一致的位點，這些位點在重測序數據的變異檢測過程中未被判定為SNP。各個測序平臺vcf 文件中的ref 和alt 與基因芯片vcf 文件中一致，可以通過直接比較位點基因型是否一致來判斷測序數據SNP 判型是否正確，即測序平臺檢測SNP 和基因芯片SNP 二者在某個位點上均為0/0（或0/1、1/0、1/1 三者之一）時，認為該測序平臺在該位點上判型正確，判型一致的位點數與位點總數的比值作為SNP 檢測準確性的評價指標。

2 結果

2.1 原始數據質量 各平臺重測序后的原始數據統計結果見表1。本研究中3 個平臺測序后的原始數據量介于61.9～83.9 Gbp 之間，平均測序數據量均在70 Gbp 以上，符合測序深度20X 以上的要求，且3 個平臺之間的原始數據量無顯著差異。GC 含量介于42%～46% 之間，與參考基因組的42% 相近，表明測序過程中出現序列偏向的可能性較低。NovaSeq 6000 的Q30 以上reads所占比例為91.71%，略高于HiSeq 2000（91.46%）（P>0.05），且二者均高于MGISEQ-2000（86.39%）（P<0.01）。HiSeq 2000 和NovaSeq 6000 的重復reads比例分別為17.17% 和14.57%，高于MGISEQ-2000（0.51%）（P<0.05）。

表1 原始reads 數據量及質量（平均值±標準差）

圖2 為3 個平臺測序數據中不同位置堿基的質量值分布。HiSeq 2000 和NovaSeq 6000 平臺之間的堿基質量值分布的差異較小，而MGISEQ-2000 平臺不同位置的堿基質量稍低于其他2 個平臺，且波動范圍較大，3 個平臺reads 不同位置的堿基質量值均在30 以上。雖然3 個平臺的原始測序數據質量存在一定差異，但都達到后續分析的要求。

圖2 reads 不同位置上的堿基質量分布

2.2 比對質量 MGISEQ-2000 的平均雙端比對率為96.20%，高于HiSeq 2000（95.49%）和NovaSeq 6000（95.37%）（P<0.01），后兩者間的差異不顯著。從圖3 可以看出，MGISEQ-2000 和HiSeq 2000 的結果一致性較高（相關系數r=0.94），而NovaSeq 6000 各個樣本之間的比對率差異較大，測序穩定性較差。圖4 為統計的插入片段長度。3 個平臺的平均插入片段長度介于322～382 bp 之間，且較為集中，能在一定程度上反映出建庫質量較好。

圖3 雙端reads 比對到參考基因組的比例

圖4 平均插入片段長度

3 個平臺的平均比對深度和位點覆蓋深度超過20X 的位點所占比例見圖5。根據比對結果，3 個平臺的平均比對深度均超過20X，達到送測要求。MGISEQ-2000 的平均比對深度為27.35X，高于HiSeq 2000（23.10X）和NovaSeq 6000（24.44X）（P<0.05），而Illumina 2 個平臺的平均比對深度無顯著差異。MGISEQ-2000、HiSeq 2000 和NovaSeq 6000 覆蓋深度在10X 以上的位點比例均在99%以上，分別為99.56%、99.45%和99.59%。3 個平臺在覆蓋深度超過20X 的位點比例的結果差異較大，MGISEQ-2000 為82.78%，高于其他2 個平臺（P<0.05），而NovaSeq 6000（73.11%）高于HiSeq 2000（65.11%）（P<0.05）。

圖5 平均比對深度和覆蓋>20X 位點比例

2.3 SNP 變異檢測 進行群體水平的變異檢測后獲得3個vcf 文件（每個平臺1 個），各平臺的變異檢測情況見表2。3 個平臺所得到的SNP 位點數相似，與Kang等[24]用13 頭豬的樣本（10 個品種）在HiSeq 2000 平臺檢測得到的結果2 812 萬相當。在SNP 和INDEL 數量上，MGISEQ-2000 多于其他2 個平臺。3 個平臺的Ti/Tv 均值均為2.40，與Lee 等[25]的研究結果相似。3個平臺共有SNP 位點數為27 359 678 個，在各個平臺位點總數的占比均達到95% 以上，檢出位點一致性較高，另外3 個平臺所檢測SNP 中dbsnp 庫收錄位點比例均達到80%以上。

表2 不同平臺測序數據變異檢測結果統計

2.4 SNP 判型準確性 將參考基因組為Sscrofa10.2 的50K 芯片轉為Sscrofa11.1 版本的過程中，有1 677 個位點未在Sscrofa11.1 中匹配到，同時將缺失和野生純合位點剔除后剩余35 871 個位點。3 個平臺各個樣本的SNP 判型準確性見圖6。MGISEQ-2000、HiSeq 2000和NovaSeq 6000 檢出50K 芯片中SNP 位點的比例分別為97.50%、97.43% 和97.40%，MGISEQ-2000 檢出的50K 芯片位點數高于其他2 個平臺，而Illumina 的2個平臺之間結果相近。以基因芯片的判型結果為參考標準，MGISEQ-2000 和HiSeq 2000 的平均準確性均達到97.21%，且二者各樣本的準確性高度一致（r=0.94），NovaSeq 6000 的S1～S7 樣本的準確性與其他2 個平臺相似，但S8 和S9 的準確性較低，分別為77.06% 和76.85%。NovaSeq 6000 的S8 和S9 2 個樣本判型與芯片不一致的位點中，2 個樣本同時出現錯判的位點約占50%。這些位點中，MGISEQ-2000 和HiSeq 2000與芯片基因型一致而NovaSeq 6000 與芯片不一致的位點占90% 以上（S8 為90.20%，S9 為90.94%），MGISEQ-2000、HiSeq 2000 和芯片陣列中判型為純合位點而NovaSeq 6000 中判型為雜合位點的位點占87%以上（S8 為87.78%，S9 為90.29%）。

圖6 3 個平臺各個樣本的判型準確性

為了分析S8 和S9 樣本位點判型錯誤的原因，本研究選擇了S8 樣本的一個SNP 位點（chr1:252 645）進行了分析。該位點在基因芯片、MGISEQ-2000 和HiSeq 2000 中基因分型均為T/T，而NovaSeq 6000中基因分型為T/G。在NovaSeq 6000 中共有26 條reads 比對到覆蓋該位點的區域，其中21 條在該位點的堿基為A/T（正/ 反鏈），其余5 條堿基為G/C。在NovaSeq 6000 的S8 個體的最終比對文件中找出該位點堿基為G/C 的reads（共5 條，read1～5），同時選擇了該位點堿基為T/A 的reads 作為參考（共2 條，read6～7），用BLAST[26]軟件將以上得到的7 條reads比對到參考基因組（Sscrofa11.1）。7 條reads 均比對到了1 號染色體，且完全匹配的堿基在148 個及以上，比對結果可信且與BWA 軟件一致，可排除軟件比對算法不同帶來的差異。圖7 展示了該位點（chr1:252 645）前后30 bp 的比對情況，其中橫線表示read 在該位點的堿基與第一行中的參考序列對應位置的堿基一致，同時顯示了所有reads 在1 號染色體252 645 位置上的堿基。S8 在該位點出現G 等位基因即由read 1～5 引起，可判斷該判型錯誤出現在原始reads 上，測序錯誤造成了SNP 的分型錯誤。

圖7 chr1:252 645 的局部比對情況

3 討 論

3.1 原始數據質量 Q30 以上reads 比例不僅受樣本類型、文庫質量、插入片段長度等因素影響，還與測序試劑和光信號采集過程等因素有關。雖然本研究中MGISEQ-2000 的Q30 以上reads 比例低于HiSeq 2000和NovaSeq 6000，但其已遠遠超過了該平臺宣傳手冊上Q30>75% 的性能參數[7]。Korostin 等[10]的研究中華大智造的MGISEQ-2000 的Q30 以上reads 比例同樣低 于Illumina 的HiSeq 2500 平 臺。MGISEQ-2000 測序重復率顯著低于Illumina 平臺的原因可能是其采用了CoolMPS 測序試劑套裝和序列片段線性擴增的建庫方式。本研究只對重復reads 進行了標記而未將其刪除，相關研究發現測序分析過程中是否去除重復的reads 對后續分析影響不大[27]。雖然3 個測序平臺在原始測序數據上表現出一定的差異，但三者均達到了測序要求，可以進行下游的數據分析。

3.2 序列比對質量 從表1 中可以看到，NovaSeq 6000在原始數據量、Q20 及Q30 以上reads 比例上均高于其他2 個平臺，但其在雙端比對率上卻低于后兩者。另外，在圖3 中可以看到，NovaSeq 6000 各個樣本之間的雙端比對率差異較大，而其他2 個平臺差異較小且不同樣本之間變化一致，表明NovaSeq 6000 的測序穩定性不 如MGISEQ-2000 和HiSeq 2000。NovaSeq 6000 在平均比對深度和覆蓋深度大于20X 的位點比例上高于HiSeq 2000，但低于MGISEQ-2000（P<0.05），但應注意3 個平臺測序數據量的差異給比對質量評價帶來的影響。在本研究中，MGISEQ-2000 雖然在原始數據量上少于其他2 個平臺，但其比對質量卻好于Illumina 的HiSeq 2000 和NovaSeq 6000 平臺。以上結果說明原始測序數據的數據量和Q20、Q30 等不能直接反映比對質量，所以在選擇測序服務時可以對序列比對環節的質量做出要求，進一步保證測序數據質量。

3.3 變異檢測結果 在變異檢測上，本研究主要分析了SNP 和INDEL 2 種變異類型。HiSeq 2000 和NovaSeq 6000 2 個平臺不僅在變異檢測數目上非常接近，其共有位點和dbsnp 庫收錄位點比例也均高于MGISEQ-2000，表明Illumina 不同時期發布的測序平臺仍能保持較高的一致性。而MGISEQ-2000 變異檢測數目高于其他2 個平臺，這可能與其平均比對深度和覆蓋深度在20X 以上的位點所占比例高于其他2 個平臺有關。總體上看，MGISEQ-2000、HiSeq 2000 和NovaSeq 6000 檢測出的SNP 和INDEL 變異數目較為接近。

3.4 SNP 判型準確性 在與基因芯片位點比較的過程中，3 個平臺都存在未檢出SNP 芯片中所有位點的情況，原因可能是測序深度不夠，一些雜合位點未被檢出，相關研究發現，重測序數據檢測出所有雜合位點要求測序深度在33X 以上[12]。實驗所用的豬為中國地方品種（馬身豬和大河豬），而比對使用的參考基因組（Sscrofa11.1）為杜洛克品種，可能一些地方豬種的特異位點不能被檢測到，但通過分析發現實驗中馬身豬和大河豬兩者在雙端比對率、平均比對深度和覆蓋>20X 的位點比例上無顯著差異。在基因分型結果上，MGISEQ-2000與HiSeq 2000 位點基因型都與芯片位點高度吻合，而NovaSeq 6000 除S8 和S9 外的樣本也與芯片結果一致。本研究中剔除了SNP 芯片中只由野生型純合和缺失2種類型組合的位點，這些位點因在群體內基因型一致，在變異檢測中未被判定為SNP，所以其實際上與芯片位點基因型是一致的，這使得本研究中的重測序數據檢測的SNP 與基因芯片的共有位點數和判型準確性偏低。部分重測序檢測SNP 位點的基因型與SNP 芯片不一致，表現形式為測序檢出的SNP 為缺失或多等位基因，可能的原因是該SNP 位點周圍存在INDEL 等變異，導致該位點上的堿基移位或缺失。NovaSeq 6000 平臺S8 和S9 樣本的原始數據量均在72 Gbp 以上，達到Q30 以上的reads 比例也均在91%以上，且兩者的雙端比對率分別為71%和64%，與S1～S7 樣本類似，所以在原始數據質量和比對質量上不能發現這2 個樣本存在問題。從NovaSeq 6000 的S8 樣本的局部比對結果中可以看出，判型錯誤的原因在于測序過程中出現了堿基錯判，其原因可能是DNA 延伸時連接了錯誤的堿基或者是光信號采集中出現了誤讀，可以考慮進一步采用PCR 和Sanger 法測序驗證這些位點的準確性。

4 結 論

MGISEQ-2000 在重復reads 比例和比對質量方面均優于HiSeq-2000 和NovaSeq-6000，在SNP 變異檢測的準確性上與HiSeq-2000 相當且高于NovaSeq-6000。NovaSeq 6000 在原始數據和序列比對上優于HiSeq 2000，而在SNP 檢測準確性上低于HiSeq 2000，且存在測序上的系統性誤差。綜合而言，HiSeq-2000 的測序質量與近幾年推出的二代測序相比未表現出明顯差距，而MGISEQ-2000 平臺重測序表現性能穩定、質量可靠，在實際應用上有明顯的優勢和應用價值。