中國(guó)野生小鼠來(lái)源1號(hào)染色體替換系缺失突變的發(fā)掘及功能注釋

徐福意，晁天柱，胡世賢，李凱，周宇荀，周梁良，肖君華

(東華大學(xué)生物科學(xué)與技術(shù)研究所，上海 201620)

中國(guó)野生小鼠來(lái)源1號(hào)染色體替換系缺失突變的發(fā)掘及功能注釋

徐福意，晁天柱，胡世賢，李凱，周宇荀，周梁良，肖君華*

(東華大學(xué)生物科學(xué)與技術(shù)研究所，上海 201620)

目的 基于野生小家鼠來(lái)源1號(hào)染色體替換系群體(population of specific chromosome 1 substitution strains, PCSSs)中18個(gè)品系的全基因組重測(cè)序結(jié)果，鑒定1號(hào)染色體上的缺失突變并對(duì)其進(jìn)行功能注釋。方法 采用Illumina二代測(cè)序平臺(tái)獲取18個(gè)品系的全基因組序列信息，通過(guò)SpeedSeq軟件鑒定缺失突變，進(jìn)一步利用SnpEff軟件完成功能注釋。結(jié)果 在18個(gè)品系的1號(hào)染色體上共檢測(cè)到13 803個(gè)缺失突變。缺失長(zhǎng)度從51 bp到70 kb不等，其中長(zhǎng)度<500 bp的缺失突變約占總數(shù)的50%。多數(shù)缺失突變位于內(nèi)含子區(qū)(50.361%)和基因間隔區(qū)(28.745%)。發(fā)現(xiàn)31個(gè)蛋白編碼基因含有功能性缺失，其中有3個(gè)基因和人類(lèi)疾病相關(guān)，7個(gè)基因參與了11條KEGG通路。結(jié)論 PCSSs的1號(hào)染色體上含有豐富的缺失突變，是在研究復(fù)雜性狀的重要遺傳標(biāo)記。

缺失突變鑒定；二代測(cè)序；功能注釋?zhuān)籔CSSs

現(xiàn)有近交系小鼠品系匱乏的遺傳多樣性是其在多基因復(fù)雜性狀研究中進(jìn)展緩慢的主要原因[1]。以具有豐富遺傳多樣性的野生小家鼠為供體建立的PCSSs，克服了近交系小鼠在數(shù)量性狀基因座(quantitative trait locus,QTL)研究中的局限性，是QTL精確定位、基因鑒定的理想資源[2]。缺失突變作為結(jié)構(gòu)變異的一種重要類(lèi)型，其大小從1 bp到幾十kb不等，廣泛分布于人和小鼠基因組中[3,4]。研究表明，缺失突變與許多人類(lèi)疾病相關(guān)[3,5]。在小鼠基因組中，缺失突變已被證明和血清尿素濃度等[6]性狀相關(guān)。但迄今，PCSSs的野生小家鼠來(lái)源1號(hào)染色體中缺失突變的遺傳多態(tài)性未被解析。

隨著二代測(cè)序技術(shù)的不斷發(fā)展，低成本、高通量的全基因組測(cè)序得以實(shí)現(xiàn)。通過(guò)全基因組測(cè)序，不僅可以檢測(cè)單個(gè)位點(diǎn)的遺傳突變，如單核苷酸多態(tài)性，也可以實(shí)現(xiàn)全基因組范圍內(nèi)復(fù)雜的結(jié)構(gòu)變異檢測(cè)，包括插入、缺失、倒位、重復(fù)等。目前，針對(duì)二代測(cè)序數(shù)據(jù)檢測(cè)結(jié)構(gòu)變異主要有5種策略[7]：(1)Read pair(也稱(chēng)為Pair-end Mapping，簡(jiǎn)稱(chēng)PEM)；(2)Split read(簡(jiǎn)稱(chēng)SR)；(3)Read depth(簡(jiǎn)稱(chēng)RD)；(4)基于de novo組裝的方法；(5)結(jié)合上述方法的綜合策略(combinatorial approach，CA)。而SpeedSeq軟件內(nèi)置的lumpy-sv結(jié)合了結(jié)構(gòu)變異檢測(cè)的前三種策略，具有較高的靈敏度和準(zhǔn)確度，顯著優(yōu)于其他軟件，廣泛應(yīng)用于結(jié)構(gòu)變異的檢測(cè)[8]。

本研究以18個(gè)“野生小家鼠來(lái)源1號(hào)染色體替換系”品系為樣本，采用Illumina二代測(cè)序平臺(tái)獲取其全基因組序列信息。利用SpeedSeq 和SnpEff 軟件完成1號(hào)染色體上缺失突變的挖掘及其功能注釋?zhuān)媪私夂驮u(píng)估缺失突變?cè)谶@一群體中的分布特征，進(jìn)一步完善該群體的遺傳變異信息，為該群體用于相關(guān)遺傳學(xué)研究奠定了基礎(chǔ)。

1 材料與方法

1.1 實(shí)驗(yàn)動(dòng)物

實(shí)驗(yàn)動(dòng)物來(lái)源于本實(shí)驗(yàn)室自建的野生小家鼠來(lái)源1號(hào)染色體替換系。動(dòng)物實(shí)驗(yàn)遵守1988年動(dòng)物管理?xiàng)l例，并得到東華大學(xué)實(shí)驗(yàn)動(dòng)物委員會(huì)批準(zhǔn)。實(shí)驗(yàn)在東華大學(xué)生物科學(xué)與技術(shù)研究所屏障動(dòng)物實(shí)驗(yàn)設(shè)施進(jìn)行[ SYXK (滬) 2014- 0022]。收集18個(gè)成年1號(hào)染色體替換系的鼠尾組織，-20℃保存?zhèn)溆谩?/p>

1.2 DNA提取和測(cè)序

DNA提取采用Axygen(愛(ài)思進(jìn)生物技術(shù)有限公司)基因組DNA抽提試劑盒。以0.8%瓊脂糖凝膠電泳、全自動(dòng)紫外與可見(jiàn)光分析儀FR-200A(上海復(fù)日科技實(shí)驗(yàn)技術(shù)研究所)和NanoDrop 2000c超微量分光光度計(jì)(Thermo Fisher Scientific, 美國(guó))確定DNA質(zhì)量和濃度。采用隨機(jī)打斷方式，選取大小為400～500 bp的DNA片段建庫(kù)。以讀長(zhǎng)(reads)為150 bp的Illumina 高通量雙末端(pair-end)測(cè)序平臺(tái)(藥明康德新藥開(kāi)發(fā)有限公司)完成全基因組重測(cè)序。

1.3 缺失突變檢測(cè)及功能注釋

首先，采用NGS QC軟件[9]對(duì)原始數(shù)據(jù)(fastq文件)進(jìn)行質(zhì)控，保留質(zhì)量值Q>20堿基數(shù)占整條reads堿基數(shù)>70%的reads；其次利用SpeedSeq軟件[10]中的align模塊將高質(zhì)量的數(shù)據(jù)比對(duì)到C57BL/6J小鼠參考基因組(GRCm38/mm10)，并生成bam文件；之后通過(guò)SpeedSeq軟件中的sv模塊鑒定缺失序列，獲得初步的缺失數(shù)據(jù)集；最后，過(guò)濾掉質(zhì)量值低于100、SU(support reads)<4、長(zhǎng)度>50 bp以及鑒定為缺失一個(gè)拷貝的數(shù)據(jù)，最終生成高質(zhì)量的缺失數(shù)據(jù)集。

缺失突變的功能注釋利用SnpEff軟件[11]完成，基因注釋信息版本為Ensembl 82。此外，利用MGI中的Human-Mouse: Disease Connection 數(shù)據(jù)庫(kù)[12]和DAVID: Functional Annotation Tools[13]分別完成人類(lèi)疾病基因注釋和KEGG 通路注釋。

2 結(jié)果

2.1 測(cè)序

測(cè)序結(jié)果表明，18個(gè)1號(hào)染色體替換系均獲得了大量、有效、高質(zhì)量的測(cè)序數(shù)據(jù)，適用于后續(xù)生物信息學(xué)研究。如表1所示，18個(gè)樣本測(cè)序數(shù)據(jù)量均超過(guò)90 Gb；質(zhì)控后的高質(zhì)量數(shù)據(jù)比對(duì)到參考基因組的數(shù)據(jù)百分比均在99%以上；各品系的平均測(cè)序深度均達(dá)到或高于30倍。18個(gè)樣本的原始數(shù)據(jù)(fastq 文件)均已上傳到NCBI (National Center for Biotechnology Information) 的SRA(Sequence Reads Archive)數(shù)據(jù)庫(kù)(編號(hào)：SRP066591)。

表1 18個(gè)品系DNA測(cè)序及1號(hào)染色體缺失檢測(cè)統(tǒng)計(jì)

2.2 缺失檢測(cè)

從長(zhǎng)度分布來(lái)看，有超過(guò)半數(shù)的缺失小于500 bp(51.9%，7164個(gè))；長(zhǎng)度在0.5～1 kb和1～2 kb的缺失分別占15.5%和11.6%；而大于10 kb的只有122個(gè)，說(shuō)明缺失突變數(shù)量與缺失長(zhǎng)度呈負(fù)相關(guān)。但長(zhǎng)度在6～7 kb的有972個(gè)，占整體數(shù)量的7%，高于其他相近長(zhǎng)度(圖1)。此外，對(duì)缺失突變?cè)?8個(gè)小鼠品系中的頻率分布統(tǒng)計(jì)發(fā)現(xiàn)(圖2)，約40%的缺失頻率小于0.1，即這些缺失僅存在于單個(gè)品系中；另近1300個(gè)缺失頻率在0.9～1.0之間，即同時(shí)存在于17或18個(gè)品系中；頻率在0.5～0.6之間的最少，僅占總數(shù)的2%(278個(gè))。

圖1 缺失長(zhǎng)度頻率分布Fig.1 Size distribution of the identified deletions

圖2 缺失在18個(gè)品系中的頻率分布Fig.2 Frequency distribution of identified deletions among 18 lines

2.3 功能注釋

功能注釋結(jié)果如表2所示，多數(shù)缺失位于內(nèi)含子區(qū)(49.655%)和基因間隔區(qū)(27.955%)，其次為基因上下游，其余區(qū)域所占比例均接近或小于1%。共發(fā)現(xiàn)605個(gè)缺失對(duì)基因功能可能產(chǎn)生功能性的影響(loss of function deletion)，其中包括99個(gè)外顯子缺失、15個(gè)移碼突變、188個(gè)剪切受體或供體突變、4個(gè)起始或終止密碼子突變以及199個(gè)轉(zhuǎn)錄本ablation。這605個(gè)缺失突變共計(jì)影響205個(gè)基因，其中包括31個(gè)蛋白編碼基因(其中21個(gè)基因含有外顯子缺失突變)、9個(gè)LincRNA、10個(gè)miRNA、3個(gè)snoRNA和7個(gè)snRNA。此外，還有部分缺失位于3’ UTR和5’ UTR區(qū)(0.36%)，這些缺失突變可能通過(guò)影響基因表達(dá)調(diào)控序列而發(fā)揮作用。

2.4 蛋白編碼基因疾病與通路注釋

31個(gè)含缺失突變的蛋白編碼基因的注釋結(jié)果顯示，其中3個(gè)基因與人類(lèi)疾病相關(guān)，分別為Gigyf2，Ptpn14和Cfh。如圖3中所示，在B6-Chr1QP和B6-Chr1YX小鼠中，與人類(lèi)帕金森病相關(guān)(OMIM:607688)的Gigyf2基因，其第28個(gè)外顯子含有一634 bp的缺失。與后鼻孔閉鎖和淋巴水腫相關(guān)(OMIM:613611)，屬于蛋白酪氨酸磷酸酶家族的Ptpn14基因，其第3個(gè)外顯子完全丟失(B6-Chr1BLD小鼠)。而與補(bǔ)體因子H缺乏(OMIM:126700)、溶血性尿毒綜合征(OMIM:235400)以及年齡相關(guān)性

黃斑變性(OMIM:610698)有關(guān)的Cfh基因的剪切受體區(qū)域存在序列缺失。31個(gè)基因的KEGG 通路注釋結(jié)果如表3所示，共發(fā)現(xiàn)7個(gè)基因參與了11條代謝通路，包括補(bǔ)體和凝血級(jí)聯(lián)、藥物代謝、嗅覺(jué)轉(zhuǎn)導(dǎo)以及細(xì)胞因子-細(xì)胞因子受體相互作用等。

表2 缺失突變功能注釋結(jié)果

注：圖中上部分為Gigyf2基因示意圖，下部分為B6-Chr1QP和B6-Chr1YX品系在Gigyf2基因上的缺失突變。圖3 Gigyf2基因在B6-Chr1QP和B6-Chr1YX小鼠中的缺失突變Note. Show in upper panel is part of Gigyf2 gene model based on Ensembl,lower panel are deletions in B6-Chr1QPand B6-Chr1YXlines.Fig.3 An exon loss of Gigyf2 gene in B6-Chr1QPand B6-Chr1YXlines

3 討論

18個(gè)1號(hào)染色體替換系小鼠的基因組中，除1號(hào)染色體來(lái)源于不同供體小鼠外，其他基因組背景均來(lái)源于近交系品系C57BL/6J。因此，本研究只針對(duì)1號(hào)染色體進(jìn)行缺失的鑒定與分析。此外，因18個(gè)小鼠品系已近交化，因此，不可能出現(xiàn)單拷貝缺失，利用這一特征，我們對(duì)鑒定的缺失做了進(jìn)一步的過(guò)濾，最終獲得了高質(zhì)量的缺失數(shù)據(jù)集。

作為遺傳變異的一種類(lèi)型，缺失可通過(guò)諸如外顯子缺失、移碼框改變等方式改變蛋白結(jié)構(gòu)，使其喪失原有功能，也可通過(guò)改變基因拷貝數(shù)或基因轉(zhuǎn)錄調(diào)控序列而影響其表達(dá)量。本次研究，在18個(gè)品系的1號(hào)染色體上注釋到605個(gè)功能性缺失，共影響205個(gè)基因，其中大部分基因?yàn)榧倩?，只?1個(gè)為蛋白編碼基因。通過(guò)對(duì)31個(gè)基因做進(jìn)一步的疾病與通路注釋?zhuān)l(fā)現(xiàn)3個(gè)基因與人類(lèi)疾病相關(guān)，7個(gè)基因參與了11條已知通路，提示這些基因的相關(guān)表型可能在該群體中存在差異。此外，我們還發(fā)現(xiàn)了一些lincRNA、miRNA等存在功能性的缺失。作為非編碼RNA，二者在基因轉(zhuǎn)錄調(diào)控方面發(fā)揮著重要作用[14,15]。

表3 KEGG 通路注釋結(jié)果

中國(guó)野生小家鼠包含兩個(gè)亞種M.m.musculus和M.m.castaneus[16]，而KM鼠是上世紀(jì)50年代引入的Swiss小鼠，與現(xiàn)有近交系小鼠來(lái)源于同一亞種M.m.domesticus[17]，本研究結(jié)果揭示，B6-Chr1KM品系1號(hào)染色體上缺失突變的數(shù)量(1510個(gè))與其他17個(gè)品系相差較遠(yuǎn)，說(shuō)明野生小家鼠缺失突變的多樣性高于實(shí)驗(yàn)室小鼠。

總之，本次研究，利用二代測(cè)序技術(shù)，結(jié)合生物信息學(xué)分析手段準(zhǔn)確解析了18個(gè)1號(hào)染色體替換系缺失突變的遺傳多樣性，進(jìn)一步完善了該群體的遺傳變異信息，為該群體用于相關(guān)遺傳學(xué)研究奠定了基礎(chǔ)。

[1] Flint J, Valdar W, Shifman S, et al. Strategies for mapping and cloning quantitative trait genes in rodents[J]. Nat Rev Genet,2005, 6(4):271-286.

[2] Xiao J, Liang Y, Li K, et al. A novel strategy for genetic dissection of complex traits: the population of specific chromosome substitution strains from laboratory and wild mice[J]. Mammalian Genome, 2010, 21(7-8):370-376.

[3] Stankiewicz P, Lupski JR.Structural variation in the human genome and its role in disease[J]. Ann Rev Med,2010, 61(1):437-455.

[4] Yalcin B, Wong K, Bhomra A, et al. The fine-scale architecture of structural variants in 17 mouse genomes[J], Genome Biol, 2012, 13(3):405-413.

[5] Weischenfeldt J, Symmons O, Spitz F, et al. Phenotypic impact of genomic structural variation: insights from and for human disease[J]. Nat Rev Genet, 2013, 14(2):125-138.

[6] Keane TM, Wong K, Adams DJ, et al. Identification of structural variation in mouse genomes.[J]. Front Genet，2014, 5:192.

[7] Zhao M, Wang Q, Wang Q, et al. Computational tools for copy number variation (CNV) detection using next-generation sequencing data: features and perspectives[J]. BMC Bioinformatics, 2013, 14(S11):1-16.

[8] Layer RM, Chiang C, Quinlan AR, et al. LUMPY: a probabilistic framework for structural variant discovery[J]. Genome Biol,2012, 15(6):1-19.

[ 9 ] Patel RK, Jain M. NGS QC Toolkit: a toolkit for quality control of next generation sequencing data[J]. PLoS One, 2012, 7(2): e30619

[10] Chiang C, Layer RM, Faust GG, et al. SpeedSeq: ultra-fast personal genome analysis and interpretation[J]. Nat Methods, 2015, 12(10):966-968.

[11] Cingolani P, Platts A, Wang LL, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3[J]. Fly, 2012, 6(2):80-92.

[12] Eppig JT, Blake JA, Bult CJ, et al. The Mouse Genome Database (MGD): facilitating mouse as a model for human biology and disease[J]. Nucleic Acids Res,2014, 43(D1):726-736.

[13] Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources[J]. Nat Protoc, 2009, 4(1):44-57.

[14] Ulitsky I, Bartel DP. lincRNAs: Genomics, evolution, and mechanisms[J]. Cell, 2013, 154(1):26-46.

[15] Jonas S, Izaurralde E. Towards a molecular understanding of microRNA-mediated gene silencing[J]. Nat Rev Genet,2015, 16(7):421-433.

[16] Boursot P, Auffray JC, Brittondavidian J, et al. The evolution of house mice[J]. Ann Rev Ecol Syst,1993, 24(4):119-152.

[17] 岳秉飛, 劉雙環(huán), 劉殿峰,等. 五個(gè)昆明小鼠封閉群遺傳生化位點(diǎn)比較研究[J]. 實(shí)驗(yàn)動(dòng)物科學(xué)與管理, 2003, 20(z1):58-62.

Deletion detection and functional annotation of chromosome 1 substitution strains from Chinese wild mice

XU Fu-yi, CHAO Tian-zhu, HU Shi-xian, LI Kai, ZHOU Yu-xun, ZHOU Liang-liang, XIAO Jun-hua*

(Institute of Biological Sciences and Biotechnology, Donghua University, Shanghai 201620, China)

Objective Deletion detection and annotation of 18 lines from the population of specific chromosome 1 substitution strains (PCSSs)derived from Chinese wild mice based on whole genome re-sequencing data. Methods Whole genome re-sequencing of the 18 lines were performed on the Illumina Hiseq platform. SpeedSeq software was used to detect the deletion after read alignment. Further annotation was obtained using SnpEff software. Results 13 803 deletions were identified among the 18 lines, the length of deletion was ranged from 51bp to 70 kb, among them nearly 50% were less than 500 bp. Through functional annotation,we found most of the variants were located in intronic (50.361%) and intergenic (28.745%) regions. However, we also identified 31 protein coding genes harboring loss-of-function deletions. Among them, 3 genes were associated with human diseases, 7 genes were participated in 11 KEGG pathways. Conclusion The chromosome 1 of PCSSs harbors abundant deletion mutations which can be used as genetic markers in genetic studies.

Deletion detection; Next-generation sequencing;Function annotation; PCSSs

XIAO Jun-hua，E-mail: xiaojunhua@dhu.edu.cn

國(guó)家自然科學(xué)基金面上項(xiàng)目(編號(hào)：31371257)；上海市科技支撐計(jì)劃(編號(hào)：15140900500)。

徐福意(1987-)，男，博士研究生，研究方向：醫(yī)學(xué)分子遺傳學(xué)。E-mail： xufuyiphd@gmail.com

肖君華(1968-)，男，教授，研究方向：醫(yī)學(xué)分子遺傳學(xué)。E-mail： xiaojunhua@dhu.edu.cn

Q95-33

A

1005-4847(2017) 01-0020-05

10.3969/j.issn.1005-4847.2017.01.004

2016-06-28