家蠶品種桂蠶8號親本全基因組重測序研究初報

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

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

1 分析方法

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

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

1.2 比對方法

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

1.3 SNP和InDel檢測方法

應(yīng)用GATK 軟件開展SNPs 和InDels 檢測，用Annovar軟件對SNPs和InDels進行注釋和統(tǒng)計［5］，檢測的過程如下。

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

（2）使用GATK的HaplotypeCaller工具進行檢測獲得所有位點的GVCF文件，使用GenotypeGVCFs工具進行變異檢測SNP和Indel變異位點，再進行過濾。

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

（4）基于基因組的gtf 文件分別對每個樣品的SNP和Indel位點進行Annovar［5］注釋。

1.4 SV檢測方法

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

1.5 CNV檢測方法

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

2 結(jié)果與分析

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

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

2.2 比對結(jié)果分析

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

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

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

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

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

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

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

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

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

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

2.4 插入缺失（InDel）分析

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

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

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

圖3 InDel長度分布圖

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

表5 InDel突變位點注釋信息統(tǒng)計

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

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

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

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

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

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

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

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

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

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

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

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

3.8 所有變異類型的Circos圖

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

圖4 所有變異類型的Circos圖

4 小結(jié)

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