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

（1.中國(guó)農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，北京 100193；2.全國(guó)畜牧總站畜禽遺傳資源保存利用中心，北京 100193；3.中國(guó)農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193；4.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)基因組研究所，廣東深圳 518120）

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

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

1 材料與方法

1.1 實(shí)驗(yàn)動(dòng)物及樣品 本研究實(shí)驗(yàn)動(dòng)物為9 頭公豬（4頭馬身豬、5 頭大河豬），采集所有實(shí)驗(yàn)動(dòng)物的耳組織，浸沒于75%的酒精，置于-80℃冰箱保存?zhèn)溆谩２捎锰旄锟萍迹ū本┯邢薰旧a(chǎn)的組織基因組DNA提取試劑盒，嚴(yán)格按照產(chǎn)品說明書提取實(shí)驗(yàn)樣本耳組織的全基因組DNA，采用酶標(biāo)儀測(cè)定待實(shí)驗(yàn)樣本的基因組DNA 濃度與純度，檢測(cè)合格后進(jìn)行后續(xù)文庫(kù)構(gòu)建。

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

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

圖1 數(shù)據(jù)分析流程及所用軟件

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

不同平臺(tái)結(jié)果之間的差異顯著性采用SPSS 20.0 的配對(duì)樣本t檢驗(yàn)進(jìn)行檢驗(yàn)，統(tǒng)計(jì)檢驗(yàn)的顯著水平（雙側(cè)）設(shè)為P<0.05，2 個(gè)變量之間的相關(guān)系數(shù)由Excel 2016中的CORREL 函數(shù)計(jì)算得出（皮爾遜相關(guān)系數(shù)）。

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

2 結(jié)果

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

表1 原始reads 數(shù)據(jù)量及質(zhì)量（平均值±標(biāo)準(zhǔn)差）

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

圖2 reads 不同位置上的堿基質(zhì)量分布

2.2 比對(duì)質(zhì)量 MGISEQ-2000 的平均雙端比對(duì)率為96.20%，高于HiSeq 2000（95.49%）和NovaSeq 6000（95.37%）（P<0.01），后兩者間的差異不顯著。從圖3 可以看出，MGISEQ-2000 和HiSeq 2000 的結(jié)果一致性較高（相關(guān)系數(shù)r=0.94），而NovaSeq 6000 各個(gè)樣本之間的比對(duì)率差異較大，測(cè)序穩(wěn)定性較差。圖4 為統(tǒng)計(jì)的插入片段長(zhǎng)度。3 個(gè)平臺(tái)的平均插入片段長(zhǎng)度介于322～382 bp 之間，且較為集中，能在一定程度上反映出建庫(kù)質(zhì)量較好。

圖3 雙端reads 比對(duì)到參考基因組的比例

圖4 平均插入片段長(zhǎng)度

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

圖5 平均比對(duì)深度和覆蓋>20X 位點(diǎn)比例

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

表2 不同平臺(tái)測(cè)序數(shù)據(jù)變異檢測(cè)結(jié)果統(tǒng)計(jì)

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

圖6 3 個(gè)平臺(tái)各個(gè)樣本的判型準(zhǔn)確性

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

圖7 chr1:252 645 的局部比對(duì)情況

3 討 論

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

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

3.3 變異檢測(cè)結(jié)果 在變異檢測(cè)上，本研究主要分析了SNP 和INDEL 2 種變異類型。HiSeq 2000 和NovaSeq 6000 2 個(gè)平臺(tái)不僅在變異檢測(cè)數(shù)目上非常接近，其共有位點(diǎn)和dbsnp 庫(kù)收錄位點(diǎn)比例也均高于MGISEQ-2000，表明Illumina 不同時(shí)期發(fā)布的測(cè)序平臺(tái)仍能保持較高的一致性。而MGISEQ-2000 變異檢測(cè)數(shù)目高于其他2 個(gè)平臺(tái)，這可能與其平均比對(duì)深度和覆蓋深度在20X 以上的位點(diǎn)所占比例高于其他2 個(gè)平臺(tái)有關(guān)?？傮w上看，MGISEQ-2000、HiSeq 2000 和NovaSeq 6000 檢測(cè)出的SNP 和INDEL 變異數(shù)目較為接近。

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

4 結(jié) 論

MGISEQ-2000 在重復(fù)reads 比例和比對(duì)質(zhì)量方面均優(yōu)于HiSeq-2000 和NovaSeq-6000，在SNP 變異檢測(cè)的準(zhǔn)確性上與HiSeq-2000 相當(dāng)且高于NovaSeq-6000。NovaSeq 6000 在原始數(shù)據(jù)和序列比對(duì)上優(yōu)于HiSeq 2000，而在SNP 檢測(cè)準(zhǔn)確性上低于HiSeq 2000，且存在測(cè)序上的系統(tǒng)性誤差。綜合而言，HiSeq-2000 的測(cè)序質(zhì)量與近幾年推出的二代測(cè)序相比未表現(xiàn)出明顯差距，而MGISEQ-2000 平臺(tái)重測(cè)序表現(xiàn)性能穩(wěn)定、質(zhì)量可靠，在實(shí)際應(yīng)用上有明顯的優(yōu)勢(shì)和應(yīng)用價(jià)值。