黃文功，浦月霞，劉艷偉，蘇紅梅，韋博尤，安春梅，譚福洋，王平陽(yáng)，閉立輝

2020年廣西蠶業(yè)技術(shù)推廣站育成抗血液型膿病新品種——桂蠶8 號(hào)［1］，4 個(gè)雜交親本“錦”（9MN）“繡”（芙H）“壯”（7MH）“麗”（87H）均為連續(xù)添食BmNPV病原篩選固定而形成，抗血液型膿病特征明顯，綜合性狀優(yōu)良。為研究4個(gè)親本基因組水平的特征，本研究抽樣并委托深圳華大基因股份有限公司開(kāi)展全基因組重測(cè)序，將測(cè)序文庫(kù)比對(duì)到家蠶參考基因組上，評(píng)估比對(duì)質(zhì)量，開(kāi)展單核苷酸多態(tài)性（SNP）、插入缺失（Indel）、結(jié)構(gòu)性變異（SV）、拷貝數(shù)變異（CNV）等分析，為下一步挖掘優(yōu)異基因做好前期工作。

1 分析方法

1.1 全基因組重測(cè)序建庫(kù)測(cè)序和原始數(shù)據(jù)過(guò)濾方法

采用BGISEQ 平臺(tái)構(gòu)建DNA 小片段文庫(kù)，片段讀取長(zhǎng)度為150 bp，每個(gè)樣本的采集數(shù)據(jù)量為20 G，用SOAPnuke［2］對(duì)數(shù)據(jù)進(jìn)行接頭、低質(zhì)量reads、雜質(zhì)數(shù)據(jù)和去N過(guò)濾后，最終得到高質(zhì)量的可用數(shù)據(jù)CleanData，過(guò)濾參數(shù)為“filter-n0.1-l20-q0.5-Q2-G”。

1.2 比對(duì)方法

應(yīng)用BWA［3-4］的“mem-t 5-M”算法，將CleanData比對(duì)到西南大學(xué)研發(fā)的參考基因組SilkDB 3.0（https：//silkdb.bioinfotoolkits.net）上，生成比對(duì)結(jié)果文件sam 文件，用samtools 的sort 工具將sam 文件轉(zhuǎn)變?yōu)榕判蚝骲am 文件，用qualimap分析比對(duì)質(zhì)量，再用samtools軟件過(guò)濾比對(duì)質(zhì)量小于30的reads。

1.3 SNP和InDel檢測(cè)方法

應(yīng)用GATK 軟件開(kāi)展SNPs 和InDels 檢測(cè)，用Annovar軟件對(duì)SNPs和InDels進(jìn)行注釋和統(tǒng)計(jì)［5］，檢測(cè)的過(guò)程如下。

（1）利用Picard的Mark Duplicate工具標(biāo)記比對(duì)結(jié)果中的duplication數(shù)據(jù)，避免PCR-duplication對(duì)后續(xù)變異檢測(cè)結(jié)果的影響。

（2）使用GATK的HaplotypeCaller工具進(jìn)行檢測(cè)獲得所有位點(diǎn)的GVCF文件，使用GenotypeGVCFs工具進(jìn)行變異檢測(cè)SNP和Indel變異位點(diǎn)，再進(jìn)行過(guò)濾。

（3）使用bcftools［6］拆分得到單個(gè)樣品的變異結(jié)果。

（4）基于基因組的gtf 文件分別對(duì)每個(gè)樣品的SNP和Indel位點(diǎn)進(jìn)行Annovar［5］注釋。

1.4 SV檢測(cè)方法

應(yīng)用BreakDancer［7］檢測(cè)SV，檢測(cè)參數(shù)為“-q 35”，過(guò)濾參數(shù)為“-m 100 -x 1000000 -s 30 -d 5”。

1.5 CNV檢測(cè)方法

使用Control-FREEC（v11.5）［8］軟件檢測(cè)樣品里面的拷貝數(shù)變異（CNV）區(qū)域。

2 結(jié)果與分析

2.1 重測(cè)序數(shù)據(jù)產(chǎn)出與質(zhì)量

單個(gè)品種測(cè)序數(shù)據(jù)經(jīng)過(guò)濾后的Clean reads 約2億條，堿基數(shù)量近300億bp，GC含量在37.90%左右，比參考基因組GC含量值（38.32%）低；Q20在95.73%～96.48%之間，Q30在90.61%～92.20%之間，重測(cè)序數(shù)據(jù)質(zhì)量可靠。

2.2 比對(duì)結(jié)果分析

選用的參考基因組（SilkDB 3.0）大小為454 710 009 bp，統(tǒng)計(jì)比對(duì)到參考基因組上各染色體區(qū)域的覆蓋度和測(cè)序深度分布情況，桂蠶8號(hào)4個(gè)親本reads比對(duì)率均接近99%，平均測(cè)序深度介于53.25×至66.65×之間，覆蓋深度從1×至51×區(qū)域的百分比介于96.09%～62.82%之間，參考基因組被均勻覆蓋，隨機(jī)性良好，詳見(jiàn)表1、圖1。

表1 比對(duì)質(zhì)量情況

圖1 參考基因組的覆蓋分布圖

2.3 單核苷酸多態(tài)性（SNP）分析

單核苷酸多態(tài)性（Single Nucleotide Polymorphisms，SNPs）是由單個(gè)核苷酸改變引起的染色體基因組的多樣性，包括單個(gè)堿基的轉(zhuǎn)換（transition，Ts）和顛換（transversion，Tv）。本研究采用GATK 軟件進(jìn)行多樣本call 變異，得到所有樣品的SNP 基因型信息，經(jīng)過(guò)濾后挑選每個(gè)樣品的位點(diǎn)信息，去掉缺失和與參考基因組一致的位點(diǎn)后統(tǒng)計(jì)每個(gè)樣品的SNP 總數(shù)，包括雜合SNP、純合SNP、轉(zhuǎn)換（Ts）、顛換（Tv）的數(shù)量，4 個(gè)雜交親本SNP 總數(shù)均超過(guò)600 萬(wàn)，純合SNP比例均為98.09%，雜合SNP占比較低，轉(zhuǎn)換（Ts）與顛換（Tv）的比值均為1.27左右，詳見(jiàn)表2。

表2 SNP數(shù)量及類型信息統(tǒng)計(jì)

對(duì)SNP變異類型進(jìn)行匯總比較，4個(gè)品種同類型堿基變異豐度基本一致，均為A 與G、C 與T 之間堿基轉(zhuǎn)換（Ts）發(fā)生位點(diǎn)較高，A 與T、A 與C、G 與C、G與T之間堿基顛換（Tv）較低，詳見(jiàn)圖2。

圖2 各種類型堿基突變豐度柱形圖

將得到的每個(gè)變異位點(diǎn)所在的基因信息注釋后，統(tǒng)計(jì)染色體基因的上下游區(qū)域、外顯子、內(nèi)含子、基因間區(qū)中變異位點(diǎn)的數(shù)量，同時(shí)判斷位于外顯子區(qū)域的SNP 突變是否會(huì)導(dǎo)致蛋白質(zhì)序列發(fā)生變化，即突變是同義突變還是非同義突變（見(jiàn)表3）。4個(gè)親本SNP 在上下游區(qū)域、外顯子、內(nèi)含子、基因間區(qū)等區(qū)域數(shù)據(jù)規(guī)律基本一致，大多數(shù)SNP 突變位于基因間區(qū)和內(nèi)含子區(qū)域，位于這2個(gè)區(qū)域的突變大多數(shù)屬于無(wú)義突變，一般不會(huì)對(duì)基因表達(dá)產(chǎn)生影響。位于基因上游區(qū)域的SNP位點(diǎn)處于16萬(wàn)個(gè)左右。位于終止密碼子區(qū)域的SNP 突變位點(diǎn)低于1 000 個(gè)，會(huì)導(dǎo)致形成終止密碼子的變異位點(diǎn)數(shù)量比引起終止密碼子消失SNP 多。在外顯子區(qū)域，引進(jìn)同義突變的SNP 均超過(guò)10 萬(wàn)個(gè)，引起非同義突變的有均超過(guò)3.7 萬(wàn)個(gè)。在可變剪切位點(diǎn)、基因下游等位置亦有大量SNP分布。

表3 SNP突變注釋信息統(tǒng)計(jì)表

2.4 插入缺失（InDel）分析

InDel 是指在DNA 上發(fā)生的核苷酸或核苷酸序列的插入和缺失，采用GATK 檢測(cè)InDel 位點(diǎn)信息，先得到所有樣品的基因型信息（genotype），再?gòu)闹刑暨x每個(gè)樣品的InDel 信息，經(jīng)過(guò)濾后統(tǒng)計(jì)樣品的InDel 總數(shù)、雜合Indel、純合InDel 的數(shù)量，去掉缺失和與參考基因組一致的位點(diǎn)后，得到單個(gè)樣品的InDel 位點(diǎn)，然后統(tǒng)計(jì)每個(gè)樣品的InDel 總數(shù)，雜合Indel、純合InDel 的數(shù)量。每個(gè)親本的Indel 總數(shù)約140萬(wàn)，其中插入位點(diǎn)為66萬(wàn)～68萬(wàn)，缺失位點(diǎn)為74萬(wàn)～75萬(wàn)，純合位點(diǎn)約占插入位點(diǎn)的92%（表4）。

表4 InDel數(shù)量統(tǒng)計(jì)

對(duì)InDel的長(zhǎng)度信息進(jìn)行統(tǒng)計(jì)，得到突變頻率較高的長(zhǎng)度為1～20 bp的InDel數(shù)量分布圖，從圖中可看出，4 個(gè)品種的插入和缺失長(zhǎng)度與數(shù)量趨勢(shì)一致，插入和缺失長(zhǎng)度主要集中在長(zhǎng)度為1～10 bp，隨著InDel長(zhǎng)度的增加，數(shù)量越少（圖3）。

圖3 InDel長(zhǎng)度分布圖

基因組注釋得到每個(gè)基因位置信息，使用Annovar注釋，統(tǒng)計(jì)基因的上下游區(qū)域、外顯子、內(nèi)含子、基因間區(qū)中變異位點(diǎn)的數(shù)量，同時(shí)判斷InDel 變異會(huì)否導(dǎo)致移碼突變。Indel分布趨勢(shì)與SNP基本一致，但未見(jiàn)位于上游/下游區(qū)域的Indel（表5）。

表5 InDel突變位點(diǎn)注釋信息統(tǒng)計(jì)

3.6 結(jié)構(gòu)性變異（SV）檢測(cè)分析

本研究中應(yīng)用BreakDancer 檢測(cè)結(jié)構(gòu)性變異（Structural variation，SV）數(shù)量和總長(zhǎng)度情況，包括長(zhǎng)度在50 bp 以上的長(zhǎng)片段序列的插入（INS）、缺失（DEL）、倒位（INV）、染色體內(nèi)易位（ITX）和染色體間易位（CTX）等，桂蠶8號(hào)的4個(gè)親本均有一定數(shù)量的結(jié)構(gòu)性變異（表6、表7）。4個(gè)親本均為缺失位點(diǎn)最多，插入位點(diǎn)最少，從相同類型的結(jié)構(gòu)性變異來(lái)看，日系親本“壯”（7MH）和“麗”（87H）插入位點(diǎn)比2個(gè)中系親本多，缺失位點(diǎn)比中系親本少，倒位數(shù)量相當(dāng)，染色體內(nèi)易位和染色體間易位略多于中系親本。

表6 各類型結(jié)構(gòu)性變異（SV）數(shù)量統(tǒng)計(jì)

表7 各種結(jié)構(gòu)性變異（SV）總長(zhǎng)度統(tǒng)計(jì)

從結(jié)構(gòu)性變異總長(zhǎng)度來(lái)看，插入總長(zhǎng)度亦為日系親本長(zhǎng)；缺失長(zhǎng)度相當(dāng)，但日系親本缺失數(shù)量較少，單個(gè)缺失位點(diǎn)長(zhǎng)度較長(zhǎng)；染色體內(nèi)易位（ITX）單位點(diǎn)平均長(zhǎng)度遠(yuǎn)大于染色體間易位（CTX）。

基于基因組注釋得到每個(gè)基因位置信息的gtf文件，使用Annovar注釋軟件將變異位點(diǎn)所在的基因信息注釋出來(lái)，并統(tǒng)計(jì)基因的上下游區(qū)域、外顯子、內(nèi)含子、基因間區(qū)中結(jié)構(gòu)性變異位點(diǎn)的數(shù)量，結(jié)構(gòu)性變異主要分布在基因間區(qū)、內(nèi)含子、外顯子、基因上游、基因下游，在可變剪切位點(diǎn)、UTR5、UTR3 及位于一個(gè)基因的上游同時(shí)也是另一個(gè)基因的下游區(qū)間亦有分布。

表8 結(jié)構(gòu)性變異（SV）位置統(tǒng)計(jì)

3.7 拷貝數(shù)變異（CNV）檢測(cè)分析

拷貝數(shù)變異（CNV）一般指長(zhǎng)度為1 kb 以上的基因組大片段的拷貝數(shù)增加或者減少，依靠特定位置的測(cè)序深度與平均測(cè)序深度來(lái)估算的。本研究使用Control-FREEC 軟件來(lái)檢測(cè)每個(gè)樣品的CNV變異，根據(jù)設(shè)定的家蠶染色體倍性為2，滑窗長(zhǎng)度為50 000 bp，自動(dòng)計(jì)算固定長(zhǎng)度區(qū)域的測(cè)序深度，并對(duì)其進(jìn)行標(biāo)準(zhǔn)化，從而得到CNV 區(qū)域的信息?；贑ontrol-FREEC的分析結(jié)果，分別對(duì)CNV偏高（gain）和偏低（loss）的區(qū)域數(shù)量和長(zhǎng)度進(jìn)行統(tǒng)計(jì)。

表9 拷貝數(shù)變異（CNV）數(shù)量統(tǒng)計(jì)

表10 拷貝數(shù)變異（CNV）類型統(tǒng)計(jì)

4 個(gè)親本拷貝數(shù)偏低的區(qū)域均超200 個(gè)，偏高的區(qū)域在79～100 個(gè)之間，所有品種在基因上游、外顯子、UTR5、基因間區(qū)均有拷貝數(shù)變異，其中以外顯子區(qū)域最多，在可變剪切位點(diǎn)、基因下游、上游/下游、UTR3區(qū)域未見(jiàn)有拷貝數(shù)變異。

3.8 所有變異類型的Circos圖

使用Circos 軟件將每個(gè)品種在染色體上的各種變異位點(diǎn)信息作到一張圖上，SNP 或者Indel 使用每1 000 000 bp頻率作Circos圖（圖4）。

圖4 所有變異類型的Circos圖

4 小結(jié)

通過(guò)對(duì)桂蠶8 號(hào)的4 個(gè)雜交親本“錦”（9MN）“繡”（芙H）“壯”（7MH）“麗”（87H）開(kāi)展全基因組重測(cè)序研究，獲得4 個(gè)全基因組單核苷酸多態(tài)性（SNP）、插入缺失（InDel）、結(jié)構(gòu)性變異（SV）和拷貝數(shù)變異（CNV），4個(gè)雜交親本SNP總數(shù)均超過(guò)600萬(wàn)個(gè)位點(diǎn)，插入缺失（Indel）總數(shù)約140 萬(wàn)；在結(jié)構(gòu)性變異（SV）方面，插入（INS）、缺失（DEL）、倒位（INV）、染色體內(nèi)易位（ITX）和染色體間易位（CTX）類型均有存在，主要分布在基因間區(qū)、內(nèi)含子、外顯子、基因上游、基因下游，在可變剪切位點(diǎn)、UTR5、UTR3及位于一個(gè)基因的上游同時(shí)也是另一個(gè)基因的下游區(qū)間亦有分布；拷貝數(shù)偏低的區(qū)域均超200 個(gè)，偏高的區(qū)域在79～100 個(gè)之間，所有品種在基因上游、外顯子、UTR5、基因間區(qū)均有拷貝數(shù)變異，其中以外顯子區(qū)域最多，在可變剪切位點(diǎn)、基因下游、上游/下游、UTR3區(qū)域未見(jiàn)拷貝數(shù)變異。