建立克服基因功能冗余的全基因组功能缺失筛选方法
技术领域
本发明属于生物技术领域，具体地说，本发明涉及建立克服基因功能冗余的全基因组功能缺失筛选方法及其应用。
背景技术
探索基因及其表达蛋白质在各种生理和病理过程中的作用是生命科学领域中永恒的研究主题，而高通量的功能缺失筛选成为了最强有力的研究工具之一。随着科学技术的发展，功能缺失筛选从早期的基于诱导DNA突变的筛选进入到最近10几年来兴起的大规模全基因组RNA干扰筛选时代。RNA干扰技术原理是利用小RNA片段靶向目的基因序列而沉默基因表达的[1]。目前RNA干扰(RNA interference,RNAi)文库被广泛应用于功能缺失型基因筛选[2-4]。这种文库主要分为siRNA和shRNA两类,二者的原理都是通过RNA干扰降低目的基因的表达来引起表型改变。siRNA为人工合成的短RNA片段,能在微孔培养板中进行彼此独立的表型筛选,配合高通量液体工作站等自动化机器，能稳定快速简便地获取目标信号或者图像等表型信息,因而被广泛应用。shRNA则是通过构建到慢病毒等表达载体上，通过病毒侵染方式进行混合型文库筛选,再通过微阵列芯片或者深度测序技术对筛选后富集的shRNA进行分析。混合型shRNA文库筛选具有简便、经济等特点,但是观察表型仅局限于细胞生长。最近，规律间隔成簇短回文重复序列(clustered regularlyinterspaced palindromic repeats，CRISPR-Cas9)这一强大的基因编辑工具，如风暴一般席卷了整个基因组工程领域，它是利用特异靶向目的基因的gRNA指导Cas9蛋白切割特定的目的DNA序列，从而起到基因编辑的作用[5,6]。并且，和shRNA相似的混合型gRNA文库筛选也已经应用于生物研究[7-10]。相信随着科学技术的进步，在不久的将来，shRNA和gRNA文库也能结合高通量自动化仪器进行多样化的表型筛选。
综上所述的功能缺失筛选都面临一个共同问题，即基因功能冗余所导致的假阴性表型。而事实上，基因功能冗余在基因组中又是一个非常普遍的现象。因为在进化过程中，功能冗余是生物体为了克服突变等功能缺失而衍生出的一种保护机制[11-13]。但是在功能缺失筛选中，功能冗余却会带来很大的干扰。因为基因间存在的代偿效应，往往使得沉默单个基因观察不到表型的变化。因此本领域中致力于开发能够克服基因功能冗余所导致的假阴性问题，高效、准确的鉴定基因功能的方法。
发明内容
本发明的目的在于提供建立克服基因功能冗余的全基因组功能缺失筛选方法。
本发明的第一方面，提供了一种构建靶向基因家族的组合siRNA文库的方法，所述方法包括步骤：
(1)提供一蛋白组群；
(2)根据所述蛋白组群中的各蛋白序列信息，进行基于结构域的多序列比对，把具有相同结构域的蛋白归为一类，形成蛋白超家族；
(3)对所述蛋白超家族中蛋白种类＞n的蛋白超家族进行拆分获得蛋白家族；对所述蛋白超家族中蛋白种类≤n的蛋白超家族不进行拆分，直接归类为蛋白家族；从而实现各蛋白家族中家族成员数均≤n；
(4)提供针对各蛋白家族中各成员的siRNA，将针对同一蛋白家族中各成员的siRNA组成siRNA集，针对不同蛋白家族的siRNA集构成了所述靶向基因家族的组合siRNA文库；
其中n为2、3、4、或5。
在另一优选例中，n为2、或3。
在另一优选例中，所述蛋白组群包括≥200种蛋白，优选地≥500种蛋白，优选地≥1000种蛋白，优选地≥2000种蛋白，优选地≥5000种蛋白。
在另一优选例中，所述蛋白组群包括同一物种的70％～100％的蛋白种类。
在另一优选例中，所述物种为哺乳动物，优选地为鼠、或人。
在另一优选例中，所述步骤(1)中，所述蛋白组群中的各蛋白均具有对应的天然或非天然siRNA(优选地，所述天然或非天然siRNA为已经报道的)。
在另一优选例中，所述步骤(2)中，对于包含多个结构域的蛋白，经基于结构域的多序列比对后，根据比对结果的统计显著性(e-value)，留取统计显著性最小的结构域，把具有相同的所述统计显著性最小的结构域的蛋白归类到一起，形成所述蛋白超家族。
在另一优选例中，所述步骤(3)中，对所述蛋白超家族中蛋白种类＞n的蛋白超家族进行拆分获得蛋白家族的具体步骤包括：
(a)对蛋白质超家族中的各个蛋白进行多重序列比对；
(b)根据步骤(a)的比对结果构建系统发育树；
(c)基于所述系统发育树反应出的序列远近关系，将系统树拆分为蛋白家族，每个蛋白家族的成员数≤n。
在另一优选例中，所述步骤(c)中，使用标号算法将系统树拆分为较小的蛋白家族，具体步骤包括：
1)初始化标号：对系统树的每个节点，用group标号标记该节点已归类入的蛋白家族号；将初始所有节点的group标号均设为0；
2)叶节点标号：遍历每个叶节点，如果该叶节点已经归类，则跳过；如果该叶节点还未归类，则取得该叶节点的直接祖先节点；
根据直接祖先节点的另一个子节点是否是叶节点，又分为两种情况:
2.1)如果该直接祖先节点的另一个子节点也是一个叶节点，则将这两个叶节点归到一个蛋白家族(首先归类的两个成员)，将他们和直接祖先的group设置为同一个家族号，例如为a；
同时，如果该直接祖先不是根节点则取得它的更上一级祖先节点，如果该祖先节点还有一个直接叶节点，则将该叶节点作为第三个成员归入，将它和二级祖先节点的group也设置为a。
2.2)如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点还未归类，则跳过；如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点已经归类，则只将这一个叶节点作为第一个成员，把它和直接祖先节点的group设置为同一个家族号，例如为b；
同时如果该祖先节点不是根节点，则再取得该祖先节点的更上一级祖先节点；如果更上一级结点还有一个直接叶节点，就将这个叶节点作为第二个成员，把它和二级祖先节点的group设为b；以同样的方法找更高一级祖先节点的直接叶节点作为第三个成员，把它和三级祖先节点的group也设为b；
3)中间节点标号：检查经过步骤2)后group还未标号的中间节点，如果该节点存在一个group为0的子节点，则跳过；
如果该节点的两个子节点的group均不为0，该节点的group设置为任意子节点的group号，表示该节点的所有下级节点均已归类。
4)循环2)和3)步骤直至所有叶节点均已归类。
在另一优选例中，所述方法还包括步骤(5)，将所述siRNA文库作用于生物样本(包括：微生物、植物或动物细胞、植物或动物组织、植物体或动物体等)，然后检测所述生物样本的表型变化。
本发明的第二方面，提供了一种靶向基因家族的组合siRNA文库，其中，所述组合siRNA文库包括t个siRNA集，所述各siRNA集分别包括针对一蛋白家族的siRNA，所述蛋白家族包括m个蛋白成员，其中m为≤n的正整数，且n为2-5的正整数；
其中，所述siRNA文库中siRNA集的数量t≥20；优选地≥50；更优选地≥100；最优选地≥500；如≥1000、≥2000、≥5000；
并且所述siRNA文库中至少30％(优选地至少40％，更优选地至少50％)的siRNA集的各自m值为2或3或4。
在另一优选例中，所述siRNA文库中60-100％，较佳地70-99％，更佳地80-90％的siRNA集的各自m值为2或3。
在另一优选例中，所述蛋白家族通过以下方法归类：
(1)提供一蛋白组群；
(2)根据所述蛋白组群中的各蛋白序列信息，进行基于结构域的多序列比对，把具有相同结构域的蛋白归为一类，形成蛋白超家族；
(3)对所述蛋白超家族中蛋白种类＞n的蛋白超家族进行拆分获得蛋白家族；对所述蛋白超家族中蛋白种类≤n的蛋白超家族不进行拆分，直接归类为蛋白家族；从而实现各蛋白家族中家族成员数均≤n。
在另一优选例中，所述的siRNA文库包括选自下组的一个或多个siRNA集：
(1)针对VPS4A基因的siRNA，针对VPS4B基因的siRNA，和针对SPG4基因的siRNA；
(2)针对BBS4基因的siRNA，和针对ST13基因的siRNA；
(3)针对DVL3基因的siRNA，针对DVL1基因的siRNA，和针对DVL2基因的siRNA；
(4)针对Gsk3a基因的siRNA，和针对GSK3b基因的siRNA；和
(5)针对Fbxw11基因的siRNA，和针对BTRC基因的siRNA。
应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。
附图说明
图1显示了本发明筛选系统的建立过程。
图2显示了全基因组siRNA文库的筛选。
图3显示了根据本发明的对基因进行家族分类的流程。
图4显示了单基因和基因家族筛选结果。
图5显示了沉默BBS4/ST13家族，Vps4a/Vps4b/Spg4家族影响Wnt3a诱导的β-catenin累积。
具体实施方式
本发明人通过广泛而深入的研究，获得一种建立克服基因功能冗余的全基因组功能缺失筛选方法，实验结果表明，所述方法能够克服基因功能冗余所导致的假阴性问题，高效、准确的鉴定出功能基因簇。
在描述本发明之前，应当理解本发明不限于所述的具体方法和实验条件，因为这类方法和条件可以变动。还应当理解本文所用的术语其目的仅在于描述具体实施方案，并且不意图是限制性的，本发明的范围将仅由所附的权利要求书限制。
除非另外定义，否则本文中所用的全部技术与科学术语均具有如本发明所属领域的普通技术人员通常理解的相同含义。如本文所用，在提到具体列举的数值中使用时，术语“约”意指该值可以从列举的值变动不多于1％。例如，如本文所用，表述“约100”包括99和101和之间的全部值(例如，99.1、99.2、99.3、99.4等)。
虽然在本发明的实施或测试中可以使用与本发明中所述相似或等价的任何方法和材料，本文在此处例举优选的方法和材料。
具体地，本发明人根据基因序列和蛋白质功能的相似性建立了一个全基因组基因家族的siRNA组合文库，借助Opera高内涵分析系统(该系统购自PerkinElmer公司)直观地对细胞核和细胞质中β-catenin进行定量分析，对小鼠全基因组siRNA文库和基于功能基因簇siRNA文库进行了高通量筛选。对比两个筛选，本发明人发现基于基因家族的siRNA组合文库能排除基因功能冗余所导致假阴性实验结果这一干扰，找到了一些影响β-catenin稳定性的候选因子。这种新的功能缺失筛选策略同样适用于其他的功能缺失筛选。
功能基因簇
沉默基因的功能可能会被另外功能相似的基因补偿，而这些基因往往属于同一家族，具有类似的功能，本文中将这些基因称为一个功能基因簇(基因家族)。针对功能基因簇的基因沉默，有助于消除基因之间代偿效应，更可能观察到表型的变化，从而发现基因(或基因家族)的新功能。本发明中，将功能基因簇中各个基因所编码的蛋白的集合称为蛋白(超)家族。
RNA干扰
如本文所用，术语“RNAi”(RNA interference,RNA干扰)是指在进化过程中高度保守的、由双链RNA(dsRNA)诱发的、高效特异性降解具有互补配对序列的RNA的现象。由于使用RNAi技术可以特异性关闭特定基因的表达，所以该技术已被广泛用于探索基因功能和传染性疾病及肿瘤的基因治疗等领域。dsRNA介导的RNAi现象在真菌、果蝇、拟南芥、锥虫、水螅、涡虫、斑马鱼等多种真核生物中均有发现，而且在植物中的转录后基因沉默(posttranscriptional gene silencing，PTGS)、共抑制(cosuppression)及RNA介导的病毒抗性、真菌的抑制(quelling)现象也均属于RNAi在不同物种的表现形式。
如本文所用，术语“siRNA”(Small interfering RNA,siRNA)是指一种小RNA分子(约21-25个核苷酸)，可由Dicer(RNA酶Ⅲ家族中对双链RNA具有特异性的酶)从其前体(比如dsRNA、shRNA等)加工而成，也可由化学方法合成或由其它蛋白加工产生。siRNA是siRISC的主要成员，激发与之序列互补的目标RNA被迅速切割降解，导致目标基因的沉默，因此成为RNAi中的关键功能分子。
如本文所用，术语“siRNA前体”是指可以在哺乳动物细胞中被加工产生siRNA的RNA分子，具体地说，是由Dicer或其它类似蛋白选择性加工从而产生成熟的siRNA，进而实施RNAi。类似地，如本文所用，术语“表达盒”是指包含本发明核酶增强型shRNA的编码序列以及与所述编码序列操作性相连的启动子和终止信号的表达盒，所述表达盒在转录后产生本发明的核酶增强型shRNA；而如本文所用，术语“构建物”是包含所述表达盒的构建物。
如本文所用，术语“shRNA”是short hairpin RNA的缩写，即，“短发夹RNA”。shRNA包括两个短反向互补序列，中间由一顶端环(loop)序列分隔的，组成发夹结构，通常由细胞内源的RNA聚合酶III(RNA polymerase III)启动子控制转录，shRNA序列的末端连接5-6个T作为RNA聚合酶Ⅲ的转录终止子。shRNA也可以由其它RNA聚合酶的启动子转录产生。
在活体中产生“小干扰RNA”(siRNA)的一种办法是，将siRNA序列作为“短发夹”的一部分克隆进质粒载体中。当送入动物体内时，该发夹序列被表达出来，形成一个带有顶端环结构的“双链RNA”(shRNA)，被细胞内的Dicer蛋白所识别和加工，产生有功能的siRNA。
RNAi筛选策略和Wnt/β-catenin信号通路
RNAi筛选已经被应用于经典的Wnt/β-catenin信号通路的研究中[14-17]。Wnt/β-catenin信号通路调控了许多生命过程，包括生物体的生长、发育、疾病、衰老与死亡等；也包括细胞形态与功能的分化与维持、免疫、应激、细胞癌变与细胞凋亡等[18,19]。在没有Wnt信号刺激时，以Axin和APC为骨架，GSK3，CK1，以及β-TrCP等蛋白质形成降解复合物，识别细胞质中游离的β-catenin。随后，在CK1和GSK3作用下，对β-catenin进行磷酸化修饰，进而在β-TrCP介导下进行泛素化修饰、降解。在发育的特定时间段，某些组织或细胞群体分泌的Wnt蛋白质，结合受体Frizzled家族成员，和共受体低密度脂蛋白质LRP5/6，将信号传递至细胞内，降解复合物的功能被抑制，促使β-catenin在细胞质中大量累积。部分累积的β-catenin进入细胞核，与核内的TCF4/LEF1家族相互作用，启动下游靶基因的表达。像其他的生物过程一样，功能冗余现象也出现于Wnt/β-catenin信号通路中[20-23]。比如说，人类基因组中受体Fz家族有10个成员(Fz1-10)，受体LRP有2个成员(LRP5/6)，3个Dvl(Dvl1-3)，2个GSK3(GSK3α/β),2个Axin(AXIN1/2)和2个β-TrCP(β-TRCP1/2)。
材料和方法
1.高内涵筛选
L细胞(ATCC:CRL-2648)使用含有10％胎牛血清(Gibco)的DMEM(Invitrogen)培养，37℃，CO2浓度为5％。小鼠全基因组siRNA文库(Dharmacon)用液体工作站(Beckman Coulter Biomek FX)取10μl(浓度为100nM)预先加入到384孔板。实验时用Multidrop微孔板分液器(Thermo Fisher)加入10μl 1:100稀释于Opti-MEM的转染试剂RNAiMax，室温放置20分钟，加入L细胞悬液30μl，在细胞培养箱培养48-72小时后取出，加入纯化的Wnt3a(纯化流程参照
http://web.stanford.edu/group/nusselab/cgi-bin/wnt/purification)刺激2小时后加入16％的多聚甲醛固定15分钟(多聚甲醛终浓度为4％，AlfaAesar,30525894))。固定好的细胞用微孔板洗板机(BioTeck ELX405)，PBST(PBS and 0.1％triton-X100)清洗，并通透20分钟，再用β-catenin的抗体(BD，610154)1:10004℃过夜孵育。微孔板洗板机(BioTeck ELX405)清洗后加入Cy3荧光标记的二抗(Jackson,115-165-062)和DAPI(Sigma,D9542)，1小时后用微孔板洗板机清洗后，利用Opera LX的20X–Air-LUCPLFLN物镜(NA＝0.45,Perkin Elmer)对β-catenin免疫染色共聚焦成像并用Acapella软件分析图像并进行定量。Ctrl siRNA和LRP6、APC、BBS4,、ST13、Vps4a、Vps4b和Spg4siRNA都购买于Dharmacon(D-001220、M-040651、M-043292、M-054691、M-056945、M-046156、M-044487、M-058588)，靶向5’-ACCAUGCAGAAUACAAAUGAU-3’(SEQ ID NO.15)的β-catenin siRNA合成于GenePharma。
2.生物信息学分析筛选数据
实验产生的高内涵筛选数据基于Bioconductor的R软件包OperaMate[24]，完成标准的处理流程，包含B-分数正规化(B Score Normalization)[25]和质量控制两部分。其中质量控制采用自助抽样法(Bootstrapping Method)构建β-catenin水平在重复实验中的平均值和标准差比值的经验分布，并以此剔除标准差相对均值显著偏大的实验(p-value<0.05)。而后，基于多重t检验和t分数法筛选显著变化的候选因子。多重t检验主要对照β-catenin水平在处理组和对照组的区别，其产生的p值通过错误发现率(False Discovery Rate)方法[26]进行校正。T分数法则仅考虑处理组β-catenin水平，用偏态t分布(skew t distribution)模拟该水平特征，从而建立t分数指标(t-score)。选用偏态t分布的主要原因是蛋白质数据多为不对称并具有较重尾部的数据，而这正符合该分布的特征[27]。同时，QQ图显示，本发明人的数据的确可以通过偏态t分布进行很好的拟合[28]。t分数指标定义为2*(1-cdf(|ts-1|+1))，其中ts是检验统计量，cdf是拟合的偏态t分布的累积分布函数。在实验中，本发明人选取t-core<0.1作为显著变化候选因子的筛选标准。
3.Western Blot,RT-PCR和定量实时PCR(quantitative real-time PCR)
对于Western Blot实验,2×SDS loading buffer裂解L细胞并100℃煮10分钟。β-catenin和β-tubulin抗体分别来源于BD TransductionLaboratories,610154，和Cell Signal ing Technology,2146。对于反转录实验，L细胞直接用TRIzol(Invi trogen)裂解，总RNA用酚氯仿抽提后以ol igo(dT)为引物，用superscriptTM III first strand sythesi s system(Invitrogen)试剂盒反转录制备cDNA文库。cDNA适当稀释后，使用Quantitative SYBR green PCR ki t(Takara SYBR premix Ex Taq)试剂盒配制实时定量PCR反应体系。实时定量PCR反应使用的仪器是ABI 7500Fastreal-time PCR system(Appl ied Biosystems)。实验中使用的定量PCR引物序列如下：GAPDH引物为5'-GCCTGCTTCACCACCTTC-3'(SEQ ID NO.1)和5'-CAAGGTCATCCATGACAACT-3'(SEQ ID NO.2)；CTNNB1引物为5'-TGCAGTTCGCCTTCACTATG-3'(SEQ ID NO.3)和5'-ACTAGTCGTGGAATGGCACC-3'(SEQ ID NO.4)；BBS4引物为5'-TGAAAACTCAGGTTCCTGCATC-3'(SEQ ID NO.5)和5'-CCTTCCAGGCGAAAAATCAGTG-3'(SEQ ID NO.6)；ST13引物为5'-TCGGGCCTTCGTGAAGATG-3'(SEQ ID NO.7)和5'-GTAGCAGGTGGTACTTTCCCC-3'(SEQ ID NO.8)；Vps4a引物为5'-ACGGTGGAATGATGTAGCTGG-3'(SEQ ID NO.9)和5'-CCAAAGAGGAGTATGCCTCGC-3'(SEQ ID NO.10)；Vps4b引物为5'-CACAAGGTGATAAAGCCAAGCA-3'(SEQ ID NO.11)和5'-GGTCGCTCTATAACAATGGCAC-3'(SEQ ID NO.12)；Spast引物为5'-AACCTGACATGCCGCAATG-3'(SEQ ID NO.13)和5'-GGACAGTTTTTGATCGAGGCAAT-3'(SEQ ID NO.14)。
本发明的主要优点在于：
(1)本发明提供的全基因组功能缺失筛选方法，能够克服基因功能冗余所导致的假阴性问题；
(2)本发明提供的全基因组功能缺失筛选方法，能够高效、准确的鉴定出功能基因簇；
(3)本发明提供的全基因组功能缺失筛选方法，能够减少全基因组基因沉默功能筛选的工作量。
下面结合具体实施例，进一步详陈本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明详细条件的实验方法，通常按照常规条件如美国Sambrook.J等著《分子克隆实验室指南》(黄培堂等译，北京：科学出版社，2002年)中所述的条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数按重量计算。以下实施例中所用的实验材料和试剂如无特别说明均可从市售渠道获得。本文中所涉及的计算机程序，如无特别说明，均为开源程序或者本领域技术人员可以从公开的渠道获得该程序。
实施例1.筛选系统的建立
Opera高内涵筛选系统能精确直观地检测免疫染色的样品，因而本发明人利用这个系统来检测细胞核和细胞质中β-catenin的分布并进行定量。根据细胞的贴壁性，形态，对Wnt3a信号的响应情况，本发明人在测试的多种细胞中选了小鼠的L细胞用于筛选。
在本发明的实验系统中，小鼠L细胞在Wnt3a刺激下，β-catenin在细胞质中累积，并进入细胞核中，如(图1，a)所示。细胞核区域根据DAPI染色由机器进行识别，胞浆区域定义为细胞核周围依照核的形状圈定的带状区域，如(图1，b)所示。这样，细胞核与细胞质中的β-catenin水平就能根据相应的区域定量出的平均荧光强度得到，每个孔核质的定量是随机取的4个视野中所有细胞核质的平均荧光强度。
本发明中把β-catenin水平定义为细胞核和细胞质的平均荧光值。本发明的实验结果显示，敲低β-catenin观察到了染色强度极大地减少，这验证了荧光信号是特异性地来源于β-catenin蛋白质。同时，敲低APC可以极大地稳定β-catenin，而当敲低LRP6时，Wnt3a稳定的β-catenin则减少，以此也很好地验证了本发明系统中稳定的β-catenin是由Wnt3a刺激引起的，如(图1，c)所示。
在进行siRNA文库筛选前，进行了384孔板整板均一性验证，利用siRNAbuffer代替siRNA，按照筛库的标准实验流程进行实验，观察相同处理的各孔间是否有很好的均一性，同时计算整块板的Z值。实验结果如(图1，d)所示。根据Zhang等[29]的报道，Z值介于0-1之间符合筛库的条件，而本发明人实验的细胞核和细胞质的Z值分别为0.61和0.63，均大于0.5，很好地达到了筛库的要求。
图1显示了本发明筛选系统的建立过程。
a：免疫染色检测β-catenin水平。L细胞在Wnt3a刺激2小时或不刺激条件下，β-catenin和DAPI染色结果图，图像由Opera LX system捕捉。
b：细胞核和细胞质定量区域的界定。细胞核区域根据DAPI染色由机器Acapella软件进行识别，核区域为0至infnity(即核的边界到核中心)；胞浆区域定义为细胞核周围依照核的形状圈定的带状区域，胞浆区域为-3至-8(即核外第三个像素至第八个像素)。
c：敲低APC，LRP6和β-catenin，在Wnt3a刺激或不刺激条件下，观察到β-catenin水平相应的变化。
d：整板筛选均一性验证，利用siRNA buffer代替siRNA按照标准筛选流程进行实验，用每个孔的编号(横坐标)和细胞核或细胞质的β-catenin平均荧光强度(纵坐标)进行作图，观察整板均一性，同时计算整块板的Z值。
实施例2.全基因组单基因siRNA文库筛选
接下来，本发明人利用建立好的Opera高内涵筛选系统进行筛选。分别敲低小鼠全基因组19059个基因，免疫染色检测Wnt3a诱导细胞核和细胞质中β-catenin的变化，独立进行了三次生物学重复筛选。并对筛选后得到的数据利用BioConductor bioinformatics package，OperaMate[24]进行均一化处理及分析，筛选流程及数据处理示意图如(图2，a)所示。筛选大数据β-catenin水平的分布呈T分布，本发明人根据T分布的规律，以及处理组与对照组两组数据间p-value设定了具有变化显著的阈值：t-score<0.1和p-value<0.05，详细见实验方法。本发明人在这批筛选数据中，确实发现了很多已经报道的Wnt信号通路的调控子，部分列于表1，然而，与此同时本发明人也发现，敲低Wnt信号通路关键的成员DVL,β-TrCP,和GSK3等并没有看到相应的β-catenin变化。根据已有的报道，发现Wnt信号通路中存在很多基因功能冗余的现象，因而，本发明人推断本发明人筛选数据中出现的这些阴性结果是由于基因功能冗余导致的。事实上，当本发明人同时敲低DVL1/2/3这个家族时，才能观察到Wnt3a诱导的细胞核和胞浆中β-catenin水平的下降，而单独敲低任何一个，都看不到明显的变化，如(图2，b)所示，类似的现象同时也在β-TrCP1/2这个家族中呈现，如(图2，c)所示。所以，对功能基因簇的基因进行沉默，有助于消除基因之间代偿效应，更可能观察到表型的变化。
图2显示了全基因组siRNA文库的筛选。
a：筛选流程及数据处理示意图。
b，c：功能基因簇siRNA组合敲低效应验证。b，Dvl家族的敲低能抑制Wnt3a诱导的β-catenin水平，但是单独敲低任何一个Dvl1/2/3观察不到明显的变化。c，β-TrCP1/2家族敲低能累积β-catenin，但是单基因的敲低观察不到明显的变化。
表1 全基因组单基因siRNA文库筛选
实施例3.基因家族组合siRNA文库的建立
为了解决功能缺失筛选中基因功能冗余所造成的假阴性这个问题，本发明人决定建立一个靶向基因家族的siRNA文库。本发明人所购买的来自于Dharmacon的小鼠全基因组siRNA文库中，每个基因都是针对目的基因不同靶序列的4条siRNA的组合(购买的文库中4条siRNA都是针对一个基因的，只是针对这个基因上不同靶向序列，以确保敲低效率，而本发明的文库是针对不同靶基因的)，因而当本发明人要进行功能基因簇siRNA组合敲低一个基因家族表达水平的时候，是多条混合的siRNA。比如像Dvl这样的基因家族由Dvl1、Dvl2和Dvl3组成，β-TrCP家族由β-TrCP1和β-TrCP2组成。在本发明人的系统中，这样的siRNA组合并不会影响单个基因的敲低效率，于是本发明人设定一个基因家族最多含有3个基因成员(在本发明的系统中，通过Dvl基因家族验证3个基因的siRNA组合并不会影响单个基因的敲低效率，所以设定一个基因家族最多含有3个基因成员)。由于当前并没有本发明人所定义的这种基因家族数据库，于是，本发明人根据蛋白质功能和基因序列的相似性开发了一个计算分析流程对基因进行家族分类。这个流程的大概思路如(图3，a)所示。首先，利用Pfam[30]对来自GenBank[31]的蛋白质序列进行解析，根据功能相关性等因素对基因进行家族分类，把基因归类为超家族，每个超家族都有共同的结构域。然后，再分别对每一个超家族根据蛋白质序列的相似性进行序列联配和构建进化树[32]。最后，利用进化树反映出的序列相似性关系，将大的超家族划分为小的家族，每个家族的成员数至多为3个。(图3，b)所示就是一个基于激酶超家族的进化树，将其拆分为成员数至多为3的的家族的例子。这个流程的步骤可详述如下：
第一步，获取Dharmacon Mouse Genome siRNA Library中每条siRNA对应的靶基因的蛋白序列。方法是，对于给定的靶基因GI号，通过NCBI的E-utilities接口(http://www.ncbi.nlm.nih.gov/books/NBK25500)下载相应的基因信息(xml格式)。E-utilities允许以http协议批量查询和下载数据。然后，解析得到的xml文件，从其中的<IUPACaa>字段中提取相应的蛋白质序列。
第二步，将蛋白质序列批量提交到Pfam网站(http://pfam.xfam.org/search)上进行结构域的注释。Pfam数据库是一个蛋白质家族大集合，基于结构域的多序列比对和隐马尔可夫模型来构建。Pfam网站返回的结果包括蛋白质上比对到的所有结构域的列表和其统计显著性(e-value)。采用阈值(1e-4)对列表进行过滤，只保留具有统计显著性的结构域。对于包含多个结构域的蛋白质，本发明人只取e-value最小的的结构域，因为该结构域的特征最为显著。最后，把具有相同结构域的蛋白质序列归类到一起，形成蛋白质超家族(fasta格式)。
第三步，一个蛋白质超家族的成员数往往大于3，所以需要对这些大的家族进行拆分。为此，本发明人先调用ClustalW程序(http://www.clustal.org/)对每一个蛋白质超家族进行多重序列比对(参数默认)。ClustalW是一种渐进的多序列比对方法，先将多个序列两两比对构建距离矩阵；然后根据距离矩阵计算产生系统进化指导树，对关系密切的序列进行加权；然后从最紧密的两条序列开始，逐步引入临近的序列并不断重新构建比对，直到所有序列都被加入为止。比对好的蛋白质序列以clustal格式保存。
第四步，对成员数大于3的蛋白质超家族，利用Bioperl软件包(http://www.bioperl.org/wiki/Main_Page)中的相关模块构建系统发育树。首先，用Bio::AlignIO模块读入clustal格式的序列比对。其次，用Bio::Align::ProteinStatistics模块计算两两序列间的距离，采用Kimura方法来近似PAM距离。在一些情况下，两个分歧较远的蛋白质可能没有公共区域能比对上，本发明人将他们之间的距离设置为1。最后，基于距离矩阵用Bio::Tree::DistanceFactory模块构建系统发育树，构树方法选择为非加权组平均法(UPGMA)。UPGMA是一种较常用的聚类分析方法，可以得到有根树。得到的树用Bio::TreeIO模块保存为newick格式。
第五步，基于系统发育树反应出的序列远近关系，将整个系统树拆分为较小的家族，每个家族的成员数不超过3。本发明人构建了一种标号算法来实现这一目的，该算法的输入为结点数大于3的有根树，输出为每一个家族的基因成员列表。算法调用了Bio::Tree::Node模块的相关函数对树进行操作，步骤祥述如下：
1)初始化标号：对每个节点，用group标号标记该节点已归类入的家族号。由于初始所有节点都未归类，group均设为0。
2)叶节点标号：遍历每个叶节点，如果该叶节点已经归类(grouped！＝0)，则跳过。如果该叶节点还未归类(grouped＝0)，则取得该叶节点的直接祖先节点。根据直接祖先节点的另一个子节点是否是叶节点，又分为两种情况:
2.1)如果该直接祖先节点的另一个子节点也是一个叶节点，则将这两个叶节点归到一个家族(首先归类的两个成员)，将他们和直接祖先的group设置为同一个家族号(例如为a)；同时，如果该直接祖先不是根节点则取得它的更上一级祖先节点，如果该祖先节点还有一个直接叶节点，则将该叶节点作为第三个成员归入，将它和二级祖先节点的group也设置为a。
2.2)如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点还未归类，跳过。如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点已经归类，则只将这一个叶节点作为第一个成员，把它和直接祖先节点的group设置为同一个家族号(例如为b)；同时如果该祖先节点不是根节点，则再取得该祖先节点的更上一级祖先节点；如果更上一级结点还有一个直接叶节点，就将这个叶节点作为第二个成员，把它和二级祖先节点的group设为b；以同样的方法可以找更高一级祖先节点的直接叶节点作为第三个成员，把它和三级祖先节点的group也设为b。
3)中间节点标号：检查经过步骤2)后group还未标号的中间节点，如果该节点存在一个group为0的子节点，跳过。如果该节点的两个子节点的group均不为0，该节点的group设置为任意子节点的group号，表示该节点的所有下级节点均已归类。
4)循环2)和3)步骤直至所有叶节点均已归类。
本发明人通过这些分析计算，建立了一个较全面的功能基因簇图谱，如(图3，c)统计了超家族包含的家族的基因数。总体上，本发明人把小鼠全基因组19059个基因分成了5850个基因家族和4779个单基因。其中，三个基因和两个基因组成的基因家族分别占41％和34％。本发明人的基因家族siRNA文库是由Beckman液体工作站在无菌的条件下，利用cherry-picking和pooling等程序对Dharmacon Mouse Genome siRNA Library进行组合，耗时三个星期完成的。
图3，基因家族siRNA组合文库的建立
a：基因家族分类流程示意图，首先，根据蛋白质的注释将基因分成超家族，红、绿、黄表示不同的结构域，再分别对每一个超家族进行序列联配和构建进化树进行进一步的细化分类。
b：Gsk3α/β所在的激酶家族分类举例。
c：小鼠全基因组基因家族分类统计图。
实施例4.基因家族siRNA文库筛选
本发明人利用基因家族siRNA文库进行了筛选，同样进行了独立的三次生物学重复，筛选结果处理和单基因筛选一样。并且，本发明人使用和单基因筛选一样的阈值t-score<0.1和p-value<0.05，即符合β-catenin变化强度和统计学上显著这两个标准找到了一些阳性候选因子，如图4，a单基因和基因家族筛选散点火山图红色标记点所示。从图上可以看出，敲低Wnt/β-catenin信号通路关键成员的家族，β-TrCP1/2、Dvl1/2/3、Gsk3α/β，在基因家族筛选中都能显著影响β-catenin水平，而在单基因筛选中并未呈现出明显的变化。本发明人对单基因和基因家族筛选数据进行了整合分析，如(图4，b)所示，具体来说，是把基因家族筛选数据中具有统计学意义变化的候选家族与单基因筛选的结果进行比较，并把这些家族归为三大类：至少有一个单基因具有显著变化；至少有一个单基因变化，但变化程度弱；单基因成员不变。从(图4，b)中可以看出，大部分基因家族的变化都是由家族中至少有一个单基因成员的变化所引起的，这也暗示两次筛选数据的一致性很高。其中，本发明人对后两类数据更感兴趣，因为，这类候选因子很有可能是单基因筛选中的假阴性因子。本发明人分别把单基因筛选和基因家族筛选中β-catenin水平变化显著的基因进行功能分析，发现基因家族筛选更能富集Wnt信号通路以及肿瘤相关的信号通路因子，如(图4，c)所示。以上对这些数据的分析都很好地显现了本发明人基因家族筛选策略的优势。
图4，单基因和基因家族siRNA文库筛选对比分析
a：单基因和基因家族筛选结果散点火山图。以log2(β-cateninintensity)为横坐标，-log10(p-value)为纵坐标，其中p-value是根据处理组和对照组的三次实验计算出，红色标记是符合变化强度(t-score<0.1)和统计学上显著(p-value<0.05)这两个标准的候选因子。一些Wnt信号通路中关键的成员用黄色背景标记出。
b：单基因和基因家族筛选整合分析饼图。第一组：至少有一个单基因具有显著变化(inhibition t-score<0.1，promotion t-score<0.2)；第二组：至少有一个单基因变化，但变化程度弱(inhibition t-score<0.2，promotiont-score<0.3)；第三组：单基因成员不变。
c:信号通路富集分析。利用DAVID功能注释软件[33]对单基因筛选及基因家族筛选候选因子进行KEGG和Biocarta信号通路分析。
下表2中列出了部分通过上述的基因家族siRNA文库筛选获得的基因家族siRNA文库。
表2 基因家族siRNA文库筛选
经过进一步地验证，针对上述各蛋白家族(基因家族)中每个成员的单一siRNA抑制均无法观察到功能的变化，而在同时抑制该蛋白家族则可以观察到明显的功能变化。
实施例5.进一步验证基因家族筛选中的两个候选家族
本发明人挑选了两个家族，BBS4/ST13家族和Vps4a/Vps4b/Spg4家族进一步验证，这两个家族均属于单基因成员不变而基因家族变化(即，敲低家族中的单基因并不影响该基因家族的该功能，而只有同时敲除整个基因家族成员，才能观察到功能上的变化)的这一类候选家族因子。本发明人利用Opera重现了两次筛选的现象，如图(图5，ab)所示，并用Western blot实验验证了Opera实验系统如图(图5，cd)所示，确实这两个家族必须要同时敲低家族成员，才能看到相应的β-catenin水平的变化，单个无显著变化。另外，本发明人对这两个家族的siRNA进行了敲低效率验证，同时发现它们不影响β-catenin的mRNA水平，如(图5，e)所示。
图5，沉默BBS4/ST13家族，Vps4a/Vps4b/Spg4家族影响Wnt3a诱导的β-catenin累积
L细胞分别转染图示中各单基因siRNA或者基因家族siRNA，利用免疫染色Opera实验系统(a，b)或者Western Blot实验(c，d)检测在Wnt3a刺激2小时或不刺激的条件下，β-catenin的水平。各处理组siRNA的敲低效率及β-catenin的mRNA水平如图e所示。
总结与讨论
(1)在本发明中，本发明人建立了一个克服基因功能冗余的全基因组功能缺失筛选方法，以检测Wnt3a诱导的β-catenin蛋白质的稳定性为出发点，通过对比分析小鼠全基因组单基因siRNA文库和基因家族siRNA文库的筛选结果，证明了相比于常规的单基因siRNA文库，基于基因家族siRNA组合文库的筛选能排除同一家族基因间的代偿效应所导致的假阴性实验结果这一干扰。这种新的功能缺失筛选策略同样适用于其他的功能缺失筛选。
(2)对单基因和基因家族筛选数据进行整合分析，本发明人对常规单基因siRNA文库筛选可能忽略的，单基因成员都不变而基因家族变化的这一类候选基因较感兴趣。本发明人对这类数据中的BBS4/ST13和Vps4a/Vps4b/Spg4这两个家族进行了进一步验证。BBS4主要功能是参与纤毛的形成，和它具有类似序列和功能的BBS6/10/12[34-36]被报道可以发挥分子伴侣的功能影响蛋白质的折叠，而BBS4目前为止，无相关功能报道。有趣的是，ST13是一个已报道的分子伴侣[37]，所以本发明人的筛选结果就提示BBS4也可能具有分子伴侣功能，但是他的分子伴侣功能和ST13功能是冗余的，所以单独敲低BBS4，并不能显现它的分子伴侣功能。对于Vps4这个家族来说，它们是一个ATPase，是ESCRT复合物的关键组成成分[38]，并且ESCRT这个复合物影响细胞自体吞噬(autophagy)过程[39]，而Gao等工作发现自体吞噬会通过促进Dvl蛋白质的降解抑制Wnt信号通路[40]。所以，这些工作解释了本发明人敲低Vps4a/Vps4b/Spg4这个家族，观察到β-catenin水平上升这个现象。
(3)对于全基因组范围内筛选Wnt信号通路调控因子的相关工作已经有很多报道了，但是这类工作都是检测Wnt信号通路下游报告基因的转录水平，而本发明人检测的信号是Wnt3a诱导的内源β-catenin水平，相对来说可以排除一些非特异转录的干扰。
同时本发明人对细胞核和细胞质中的β-catenin定量，可以进一步分析调控β-catenin核质分布的因子。当然，本发明人的筛选阳性结果中还包含了很多间接的影响因子以及脱靶的假阳性因子，确实，这里面出现了很多house-keeping功能的因子，包括，转录、翻译、蛋白酶体以及代谢相关的分子这些都是大规模筛选不可避免的。虽然本发明人的基因家族siRNA文库的筛选能消除基因功能冗余的假阴性现象，但是同时本发明人也发现了一个不足之处。比如说，对于Fzd家族来说，在L细胞中高表达的Fzd7(数据未展示)，在单基因筛选结果中是阳性结果，但是包含Fzd7这一基因家族在基因家族筛选中却没有变化。本发明人推测这可能是由于脱靶效应或者是家族中其他成员相反效应的中和(Fzd1单基因筛选结果和Fzd7相反)所造成的。所以综合分析单基因，基因家族以及其他类型，比如说过表达等筛选结果，再进行第二次筛选得到阳性的概率会大大增加。
在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
参考文献
1.Fire,A.,et al.,Potent and specific genetic interference by double-stranded RNA in Caenorhabditiselegans.Nature,1998.391(6669):p.806-11.
2.Diehl,P.,D.Tedesco,and A.Chenchik,Use of RNAi screens to uncover resistance mechanisms incancer cells and identify synthetic lethal interactions.Drug Discov Today Technol,2014.11:p.11-8.
3.Gao,S.,et al.,Applications of RNA interference high-throughput screening technology in cancerbiology and virology.Protein Cell,2014.5(11):p.805-15.
4.Karlsson,C.,J.Rak,and J.Larsson,RNA interference screening to detect targe Table moleculesin hematopoietic stem cells.Curr Opin Hematol,2014.21(4):p.283-8.
5.Cong,L.,et al.,Multiplex genome engineering using CRISPR/Cas systems.Science,2013.339(6121):p.819-23.
6.Mali,P.,et al.,RNA-guided human genome engineering via Cas9.Science,2013.339(6121):p.823-6.
7.Koike-Yusa,H.,et al.,Genome-wide recessive genetic screening in mammalian cells with alentiviral CRISPR-guide RNA library.Nat Biotechnol,2014.32(3):p.267-73.
8.Zhou,Y.,et al.,High-throughput screening of a CRISPR/Cas9 library for functional genomics inhuman cells.Nature,2014.509(7501):p.487-91.
9.Konermann,S.,et al.,Genome-scale transcriptional activation by an engineered CRISPR-Cas9complex.Nature,2015.517(7536):p.583-8.
10.Parnas,O.,et al.,A Genome-wide CRISPR Screen in Primary Immune Cells to Dissect RegulatoryNetworks.Cell,2015.
11.Brookfield,J.F.,Genetic redundancy.Adv Genet,1997.36:p.137-55.
12.Nowak,M.A.,et al.,Evolution of genetic redundancy.Nature,1997.388(6638):p.167-71.
13.Wagner,A.,Selection and gene duplication:a view from the genome.Genome Biol,2002.3(5):p.reviews1012.
14.Major,M.B.,et al.,New regulators of Wnt/beta-catenin signaling revealed by integrativemolecular screening.Sci Signal,2008.1(45):p.ra12.
15.Tang,W.,et al.,A genome-wide RNAi screen for Wnt/beta-catenin pathway components identifiesunexpected roles for TCF transcription factors in cancer.Proc Natl Acad Sci U S A,2008.105(28):p.9697-702.
16.Simons,M.,et al.,Electrochemical cues regulate assembly of the Frizzled/Dishevelled complex atthe plasma membrane during planar epithelial polarization.Nat Cell Biol,2009.11(3):p.286-94.
17.Conrad,W.,et al.,FAM129B is a novel regulator of Wnt/beta-catenin signal transduction inmelanoma cells.F1000Res,2013.2:p.134.
18.Wang,J.,T.Sinha,and A.Wynshaw-Boris,Wnt signaling in Mammalian development:lessonsfrom mouse genetics.Cold Spring Harb Perspect Biol,2012.4(5).
19.Clevers,H.and R.Nusse,Wnt/beta-catenin signaling and disease.Cell,2012.149(6):p.1192-205.
20.Doble,B.W.,et al.,Functional redundancy of GSK-3alpha and GSK-3beta in Wnt/beta-cateninsignaling shown by using an allelic series of embryonic stem cell lines.Dev Cell,2007.12(6):p.957-71.
21.Schwab,K.R.,et al.,Pygo1 and Pygo2 roles in Wnt signaling in mammalian kidney development.BMC Biol,2007.5:p.15.
22.Etheridge,S.L.,et al.,Murine dishevelled 3 functions in redundant pathways with dishevelled 1and 2 in normal cardiac outflow tract,cochlea,and neural tube development.PLoS Genet,2008.4(11):p.e1000259.
23.Satoh,W.,et al.,Sfrp1,Sfrp2,and Sfrp5 regulate the Wnt/beta-catenin and the planar cellpolarity pathways during early trunk formation in mouse.Genesis,2008.46(2):p.92-103.
24.Gentleman,R.C.,et al.,Bioconductor:open software development for computational biology andbioinformatics.Genome Biol,2004.5(10):p.R80.
25.Brideau,C.,et al.,Improved statistical methods for hit selection in high-throughput screening.JBiomol Screen,2003.8(6):p.634-47.
26.Benjamini,Y.and Y.Hochberg,Controlling the False Discovery Rate-a Practical and PowerfulApproach to Multiple Testing.Journal of the Royal Statistical Society Series B-Methodological,1995.57(1):p.289-300.
27.Shahrezaei,V.and P.S.Swain,Analytical distributions for stochastic gene expression.Proc NatlAcad Sci U S A,2008.105(45):p.17256-61.
28.Hansen,B.E.,Autoregressive Conditional Density Estimation.International Economic Review,1994.35(3):p.705-730.
29.Zhang,J.H.,T.D.Chung,and K.R.Oldenburg,A Simple Statistical Parameter for Use inEvaluation and Validation of High Throughput Screening Assays.J Biomol Screen,1999.4(2):p.67-73.
30.Finn,R.D.,et al.,Pfam:the protein families database.Nucleic Acids Res,2014.42(Databaseissue):p.D222-30.
31.Coordinators,N.R.,Database resources of the National Center for Biotechnology Information.Nucleic Acids Res,2014.42(Database issue):p.D7-17.
32.Holder,M.and P.O.Lewis,Phylogeny estimation:traditional and Bayesian approaches.Nat RevGenet,2003.4(4):p.275-84.
33.Huang da,W.,B.T.Sherman,and R.A.Lempicki,Bioinformatics enrichment tools:paths towardthe comprehensive functional analysis of large gene lists.Nucleic Acids Res,2009.37(1):p.1-13.
34.Kim,J.C.,et al.,MKKS/BBS6,a divergent chaperonin-like protein linked to the obesity disorderBardet-Biedl syndrome,is a novel centrosomal component required for cytokinesis.J Cell Sci,2005.118(Pt 5):p.1007-20.
35.Stoetzel,C.,et al.,BBS10 encodes a vertebrate-specific chaperonin-like protein and is a majorBBS locus.Nat Genet,2006.38(5):p.521-4.
36.Stoetzel,C.,et al.,Identification of a novel BBS gene(BBS12)highlights the major role of avertebrate-specific branch of chaperonin-related proteins in Bardet-Biedl syndrome.Am J HumGenet,2007.80(1):p.1-11.
37.Johnson,B.D.,et al.,Hop modulates Hsp70/Hsp90 interactions in protein folding.J Biol Chem,1998.273(6):p.3679-86.
38.Wollert,T.,et al.,The ESCRT machinery at a glance.J Cell Sci,2009.122(Pt 13):p.2163-6.
39.Rusten,T.E.and H.Stenmark,How do ESCRT proteins control autophagy？J Cell Sci,2009.122(Pt 13):p.2179-83.
40.Gao,C.,et al.,Autophagy negatively regulates Wnt signalling by promoting Dishevelleddegradation.Nat Cell Biol,2010.12(8):p.781-90.
