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

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

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

1 分析方法

1.1 全基因組重測序建庫測序和原始數據過濾方法

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

1.2 比對方法

應用BWA［3-4］的“mem-t 5-M”算法，將CleanData比對到西南大學研發的參考基因組SilkDB 3.0（https：//silkdb.bioinfotoolkits.net）上，生成比對結果文件sam 文件，用samtools 的sort 工具將sam 文件轉變為排序后bam 文件，用qualimap分析比對質量，再用samtools軟件過濾比對質量小于30的reads。

1.3 SNP和InDel檢測方法

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

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

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

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

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

1.4 SV檢測方法

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

1.5 CNV檢測方法

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

2 結果與分析

2.1 重測序數據產出與質量

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

2.2 比對結果分析

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

表1 比對質量情況

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

2.3 單核苷酸多態性（SNP）分析

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

表2 SNP數量及類型信息統計

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

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

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

表3 SNP突變注釋信息統計表

2.4 插入缺失（InDel）分析

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

表4 InDel數量統計

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

圖3 InDel長度分布圖

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

表5 InDel突變位點注釋信息統計

3.6 結構性變異（SV）檢測分析

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

表6 各類型結構性變異（SV）數量統計

表7 各種結構性變異（SV）總長度統計

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

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

表8 結構性變異（SV）位置統計

3.7 拷貝數變異（CNV）檢測分析

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

表9 拷貝數變異（CNV）數量統計

表10 拷貝數變異（CNV）類型統計

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

3.8 所有變異類型的Circos圖

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

圖4 所有變異類型的Circos圖

4 小結

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