全基因组siRNA高通量筛选调控Wnt/β-catenin信号转导的关键基因
High-throughput Genome-wide siRNA Screen for Key Genes Regulating Wnt/β-CateninSignaling   

材料与试剂

1. 细胞培养384微孔板 (Perkin Elmer, catalog number: 6007460)
2. L细胞 (ATCC, catalog number: CRL-2648)
3. 胎牛血清 (Invitrogen, catalog number: 10099141) 
4. 多聚甲醛 (Alfa Aesar, catalog number: 30525894)
5. DMEM培养液 (Invitrogen, catalog number: C11995500CP) 
6. siRNA文库 (Dharmacon,100 nM) 
7. Opti-MEM (Invitrogen, catalog number: 31985070) 
8. 转染试剂RNAiMAX (Invitrogen, catalog number: 13778150) 
9. PBST (PBS添加 0.1% triton-X100) 
10. DAPI (Sigma, catalog number: D9542) 
11. Ctrl siRNA (Dharmacon, catalog number: D-001220) 
12. 针对小鼠LRP6的siRNA (Dharmacon, catalog number: M-040651)
13. 针对小鼠APC的siRNA (Dharmacon, catalog number: M-043292) 
14. 靶向β-catenin 的siRNA 5'-ACCAUGCAGAAUACAAAUGAU-3'(GenePharma) 
15. β-catenin 抗体 (BD Transduction Laboratories, catalog number: 610154) 
16. β-tubulin抗体 (Cell Signaling Technology, catalog number: 2146S) 
17. Cy3荧光标记的二抗 (Jackson, catalog number: 115-165-062)
18. PBST (见溶液配方)
19. Wnt3a纯化 (见溶液配方)
    

仪器设备

1. Biomek FX自动化液体处理平台体(Beckman Coulter 公司)
2. Multidrop微孔板分液器 (Thermo Fisher公司) 
3. ELX405微孔板洗板机 (BioTek公司) 
4. 高内涵细胞分析仪Opera LX (PerkinElmer公司) 
5. 细胞培养箱 (Thermo Fisher公司)
   

实验步骤

1. 检测L细胞中β-catenin对Wnt信号的响应。
   L细胞在含10%胎牛血清的DMEM，37 °C，5% CO2的条件下培养，细胞隔天传代，保持传代时的密度为85%至90%之间。0.1 μl转染试剂RNAiMAX稀释于10 μl的Opti-MEM中，加入到10 μl 100 nM的Ctrl siRNA中，室温放置20 min，然后，加入L细胞悬液30 μl，放入培养箱培养72h后取出，此时，加入纯化的10 μl Wnt3a (稀释于Opti-MEM。Wnt3a 是一种分泌性糖蛋白，是经典Wnt信号通路的主要配体，它和其细胞膜上受体Frizzled以及共受体LRP5/6结合后，可抑制细胞内由Axin、APC 、GSK3以及β-TrCP等蛋白质形成的降解复合物对β-catenin的降解，从而促使β-catenin在细胞质中大量累积，部分累积的β-catenin进入细胞核，与核内的TCF4/LEF1家族相互作用，最终启动下游靶基因的表达。加入Wnt3a 2h后取出细胞培养板，准备进行免疫荧光染色。首先加入20 μl的16%的多聚甲醛固定15 min (多聚甲醛终浓度为4%，Alfa Aesar, 30525894)，将固定好的细胞用PBST (PBS 加 0.1% triton-X100) 于微孔板洗板机 (BioTeck ELX405) 中洗涤3次，并通透20min，然后，加入β-catenin的抗体 (BD,610154,1:1000)，4°C孵育过夜。微孔板洗板机洗涤3次后加入Cy3荧光标记的二抗 (Jackson, 115-165-062) 和DAPI (Sigma, D9542)，室温孵育1h后，用微孔洗板机洗涤后，利用Opera LX的 20X–Air-LUCPLFLN 物镜 (NA=0.45, Perkin Elmer) 对β-catenin免疫染色共聚焦成像，并用Acapella软件分析图像并进行定量。
2. 对细胞核和细胞质中β-catenin的定量分析。
   捕捉β-catenin和DAPI的染色结果后，细胞核区域根据DAPI染色由机器进行识别，定义为核的边界到核中心，胞浆区域定义为细胞核周围依照核的形状圈定的核外第三个像素至第八个像素之间的带状区域，如图1所示。根据相应定义区域中平均荧光强度，我们就能得到每个细胞的细胞核与细胞质中的β-catenin的定量值。每个孔随机取4个视野内所有细胞做分析，这样我们就得到了该孔样品中所有细胞核和细胞质的平均荧光强度。我们把β-catenin的水平定义为细胞核和细胞质的平均荧光强度。细胞核和细胞质的比值Ratio (Nuc/Cyto) (Nuc: Nucleus,Cyto: Cytoplasm) 来指示β-catenin的核质定位。
   
   
   图1. 细胞核和细胞质中β-catenin的定量分析. 细胞核区域根据DAPI染色由机器Acapella软件进行识别，核区域为0至infinity (即核的边界到核中心)；胞浆区域定义为细胞核周围依照核的形状圈定的带状区域，胞浆区域为-3至-8 (即核外第三个像素至第八个像素)
   
   
3. β-catenin蛋白质荧光信号特异性验证。
   为了证明我们的实验系统中得到的荧光信号特异地来源于β-catenin蛋白质，我们在L细胞中，用β-catenin的siRNA敲低β-catenin后，检测β-catenin的免疫荧光信号强度。同时，为了确认实验系统中β-catenin的累积是特异性地由Wnt3a诱导引起的，我们用siRNA对Wnt信号通路调控因子APC和LRP6进行沉默，检测其对β-catenin荧光强度的影响。β-catenin或者APC和LRP6的敲低实验和随后的荧光检测具体实验过程和步骤1相同 (将步骤1中的siRNA control替换成β-catenin或者APC和LRP6的siRNA，其他步骤相同)。
4. siRNA文库筛选前的整板均一性测试。
   我们用siRNA缓冲液代替siRNA，按照筛选的标准流程进行实验，检测384孔板相同处理的各孔间的荧光信号值的均一性，根据Zhang et al. (1999) 的报道(1)，Z=1-3* (实验组SD+对照组SD) / (实验组平均值-对照组平均值) 计算整块板的Z值。Z值介于0-1之间的实验系统适合应用于高通量筛选。
5. 全基因组单基因siRNA文库筛选
   在准备工作完成后，我们开始对小鼠全基因组siRNA文库进行高通量功能缺失筛选。我们对小鼠全基因组中19,059个基因进行筛选，独立进行三次生物学重复筛选，筛选具体步骤叙述如下。
   1)

   用液体工作站Biomek FX，取10 μl小鼠全基因组siRNA文库 (100 nM) 预先加入到384孔板。
   2)

   用Multidrop微孔板分液器加入10 μl稀释于Opti-MEM的转染试剂RNAiMax (1:100稀释比)，室温放置20min，然后加入L细胞悬液30 μl，培养72h。
   3)

   加入纯化的稀释于Opti-MEM的10 μlWnt3a蛋白刺激细胞2h，取出细胞培养板加入20 μl的16%的多聚甲醛固定15min (多聚甲醛终浓度为4%)。
   4)

   将固定好的细胞用PBST (PBS + 0.1% triton-X100) 于微孔板洗板机 (BioTeck ELX405) 中洗涤3次，并通透20min，然后加入β-catenin的抗体 (1:1000)，4°C孵育过夜。
   5)

   用微孔板洗板机洗涤3次后加入Cy3荧光标记的二抗和DAPI，室温孵育1h后，用微孔洗板机洗涤后，利用Opera LX的 20X-Air-LUCPLFLN 物镜 (NA=0.45, Perkin Elmer) 对β-catenin免疫染色共聚焦成像，并用Acapella软件分析图像并进行定量。
6. 生物信息学分析筛选数据
   我们将获得的高内涵筛选数据用Bioconductor的R软件包OperaMate(2)完成标准的处理流程，包含B-分数正规化 (B Score Normalization) (3)和质量控制两部分。其中，质量控制采用自助抽样法 (Bootstrapping Method) 构建β-catenin水平在重复实验中的平均值和标准差比值的经验分布，并以此剔除标准差相对均值显著偏大的实验 (p-value<0.05)。而后，基于多重t检验和t分数法筛选显著变化的候选因子。多重t检验主要对照β-catenin水平在处理组和对照组的区别，其产生的p值通过错误发现率 (False Discovery Rate) 方法(4)进行校正。T分数法则仅考虑处理组β-catenin水平，用偏态t分布 (skew t distribution) 模拟该水平特征，从而建立t分数指标 (t-score)。选用偏态t分布的主要原因是蛋白质数据多为不对称并具有较重尾部的数据，而这正符合该分布的特征(5)。同时，我们的数据可以通过偏态t分布进行很好的拟合(6)。t分数指标定义为2*(1-cdf(|ts-1|+1))，其中ts是检验统计量，cdf是拟合的偏态t分布的累积分布函数。在实验中，我们选取t-score<0.1作为显著变化候选因子的筛选标准。由于siRNA具有脱靶效应，而目前没有恰当的阴性对照siRNA，所以我们利用中值法对筛选数据进行了标准化处理，并通过标准化消除各板间信号值的差异，使得所有筛选数据能够进行平行比较。
7. 基因家族siRNA组合文库的筛选
   1)

   在自然选择和生物进化的过程，伴随着多种多样的变异，形成了复杂和多样化的生物体。在漫长的进化过程中，为了应对变异产生的有害缺失突变（deleterious mutation），生物体逐渐进化出表达多种相似功能蛋白质的能力来保证生物体功能的正常运行，即功能冗余（functional redundancy）。因而，功能冗余是生物体为了克服突变等功能缺失而衍生出的一种保护机制，由于进化过程中产生的这种保护机制，生物体大量存在功能相似的基因。而事实上，在我们的功能缺失筛选中，功能冗余却是一个很大的干扰。基因功能冗余使得基因间存在着代偿效应，单个基因功能的沉默可能会被其它功能相似的基因补偿，这种代偿效应常常导致沉默单个基因无法观察到表型的变化。而这些基因往往属于同一家族，具有类似的功能，我们将这些基因称为一个功能基因簇。针对功能基因簇的基因沉默，有助于消除基因之间代偿效应，更可能观察到表型的变化，从而发现基因（或基因家族）的新功能，因此，我们决定建立一个靶向基因家族的siRNA组合文库。对基因进行家族分类
   我们所购买的小鼠全基因组siRNA文库来自于Dharmacon公司。为了降低脱靶效应并提高沉默效率，针对每个基因的siRNA是4条不同序列的siRNA组合，分别针对目的基因不同靶序列。如果我们要沉默一个基因家族，则需要混合针对这个基因家族所有成员的多条siRNA并进行转染。由于Wnt信号通路以及其他一些典型信号通路的关键组分大部分拥有2个或3个基因家族成员，综合考虑沉默效率问题，我们设定基因家族siRNA组合文库中一个家族最多包含3个基因成员。在我们的这个研究以前并没有我们所定义的这种基因家族数据库，于是，我们根据蛋白质功能和序列的相似性开发了一个计算分析流程对基因进行家族分类。这个流程的大概思路如图2A示。首先，利用Pfam(7)对来自GenBank(8)的蛋白质序列进行解析，根据功能相关性等因素对基因进行家族分类，把基因归类为超家族，每个超家族都拥有共同的结构域。然后，再分别对每一个超家族根据蛋白质序列的相似性进行序列联配和进化树构建(9)。最后，利用进化树反映出的序列相似性关系，将大的超家族划分为小的家族，并使每个家族的成员数至多为3个。图2B将一个激酶超家族拆分为若干成员数至多为3的家族的实例。这个流程的步骤详述如下。
   
   
   图2.基因家族siRNA组合文库的建立 (A) 基因家族分类流程示意图，首先，根据蛋白质的注释将基因分成超家族，红、绿、黄表示不同的结构域，再分别对每一个超家族进行序列联配，进化树构建并进行进一步的细化分类。(B) Gsk3α/β所在的激酶家族分类举例。在这个分类树中，两个序列最相似的蛋白质被分为一个家族，如果有第三个蛋白质的序列和他们最相近，则被分于同一个基因家族。
   
   第一步，获取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://w ww.bioperl.org/wiki/Main_Page) 中的相关模块构建系统发育树。首先，用Bio::AlignIO模块读入clustal格式的序列比对。其次，用Bio::Align::ProteinStatistics模块计算两两序列间的距离，采用Kimura方法来近似 PAM距离。在一些情况下，两个分歧较远的蛋白质可能没有公共区域能比对上，我们将他们之间的距离设置为1。最后，基于距离矩阵用Bio::Tree::DistanceFactory模块构建系统发育树，构树方法选择为非加权组平均法 (UPGMA)。UPGMA是一种较常用的聚类分析方法，可以得到有根树。得到的树用Bio::TreeIO模块保存为newick格式。
   第五步，基于系统发育树反应出的序列远近关系，将整个系统树拆分为较小的家族，每个家族的成员数不超过3。系统树拆分为较小家族的步骤详见下一节。
   2)

   系统树拆分为较小家族
   我们构建了一种标号算法来将整个系统树拆分为较小的家族，该算法的输入为结点数大于3的有根树，输出为每一个家族的基因成员列表。算法调用了Bio::Tree::Node模块的相关函数对树进行操作，步骤祥述如下。
   步骤1，初始化标号：对每个节点，用group标号标记该节点已归类入的家族号。由于初始所有节点都未归类，group均设为0。
   步骤2，叶节点标号：遍历每个叶节点，如果该叶节点已经归类 (grouped !=0)，则跳过。如果该叶节点还未归类 (grouped = 0)，则取得该叶节点的直接祖先节点。根据直接祖先节点的另一个子节点是否是叶节点，又分为两种情况。 (1) 如果该直接祖先节点的另一个子节点也是一个叶节点，则将这两个叶节点归到一个家族 (首先归类的两个成员)，将他们和直接祖先的group设置为同一个家族号 (例如为a)；同时，如果该直接祖先不是根节点则取得它的更上一级祖先节点，如果该祖先节点还有一个直接叶节点，则将该叶节点作为第三个成员归入，将它和二级祖先节点的group也设置为a。(2) 如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点还未归类，跳过。如果该直接祖先节点的另一个子节点是一个中间节点且该中间节点已经归类，则只将这一个叶节点作为第一个成员，把它和直接祖先节点的group设置为同一个家族号 (例如为b)；同时如果该祖先节点不是根节点，则再取得该祖先节点的更上一级祖先节点；如果更上一级结点还有一个直接叶节点，就将这个叶节点作为第二个成员，把它和二级祖先节点的group设为b；以同样的方法可以找更高一级祖先节点的直接叶节点作为第三个成员，把它和三级祖先节点的group也设为b。
   步骤3，中间节点标号：检查经过步骤2) 后group还未标号的中间节点，如果该节点存在一个group为0的子节点，跳过。如果该节点的两个子节点的group均不为0，该节点的 group设置为任意子节点的 group号，表示该节点的所有下级节点均已归类。
   循环步骤2和3，直至所有叶节点均已归类。
   3)

   基因家族siRNA文库筛选
   完成了基因家族siRNA组合文库的构建后，我们针对基因家族siRNA组合文库进行了高通量功能缺失筛选。我们采用和单基因筛选相同的流程及数据分析方法，同样进行了独立的三次生物学重复筛选。
   

结果与分析

1. L细胞中Wnt诱导β-catenin水平上升和入核。
   L细胞是一种来源于小鼠皮下结缔组织的成纤维细胞，在这个细胞中，β-catenin蛋白质表达水平非常低，因而几乎观察不到细胞膜上β-catenin的这一特征使得我们能够不受膜定位β-catenin的干扰，从而能较准确地定量细胞质中的β-catenin免疫染色信号。在我们的实验系统中， L细胞在没有Wnt3a刺激下，β-catenin水平很低，并且细胞膜上的定位也很少，而在Wnt3a刺激下，β-catenin在细胞质中大量累积，并进入细胞核中，如图3所示。
   
   
   图3. 免疫染色检测β-catenin蛋白质水平 L细胞在Wnt3a刺激2h或不刺激条件下，β-catenin和DAPI染色结果图，图像由Opera LX system捕捉。
   
   为了确认实验系统中β-catenin的累积是特异性地由Wnt3a诱导引起的，我们对Wnt信号通路的调控因子APC和LRP6进行了siRNA敲减。APC和LRP6是已知的Wnt信号通路的重要调控因子，其中，APC是负向调控因子，而LRP6则是介导Wnt信号从细胞膜向细胞质传递的关键蛋白质。在我们的实验系统中，沉默APC可以显著增加β-catenin的稳定性，而沉默LRP6则导致Wnt3a诱导的β-catenin累积减少，这一结果也很好地证明了我们实验系统中β-catenin的累积是特异性地由Wnt3a诱导引起的，如图4所示。
   
   
   图4. 沉默APC、LRP6和β-catenin，在Wnt3a刺激或不刺激条件下，检测β-catenin蛋白质水平的相应变化
   
   
2. 系统稳定性检测结果。
   为了检测筛选系统的稳定性，我们进行了siRNA文库筛选前的整板均一性检测。我们利用siRNA 缓冲液代替siRNA，按照筛选的标准实验流程进行实验，检测384孔板相同处理的各孔间的荧光信号值的均一性，计算整块板的Z值，实验结果如图5所示。Z值介于0-1之间的实验系统适合应用于高通量筛选，据计算，在我们的实验体系中，针对细胞核和细胞质β-catenin定量系统中，Z值分别为0.61和0.63，均大于0.5，达到了高通量筛选的要求。
   
   
   图5. 整板筛选均一性验证. 利用siRNA buffer代替siRNA按照标准筛选流程进行实验，用每个孔的编号 (横坐标) 和细胞核、细胞质的β-catenin平均荧光强度或核质比 (纵坐标) 进行作图，计算Z值检测整板的均一性。
   
   
3. 筛选数据标准化处理。
   我们对小鼠全基因组siRNA文库进行单基因高通量功能缺失筛选后，利用BioConductor的R软件包OperaMate(2)，对筛选数据进行了均一化处理，具体步骤如前所述。由于siRNA可能具有脱靶效应，而目前没有恰当的阴性对照siRNA，所以我们利用中值法对筛选数据进行了标准化处理，并通过标准化消除各板间信号值的差异，使得所有筛选数据能够进行平行比较，处理前后结果如图6A，B所示。
   
   
   图6. 单基因siRNA文库筛选数据标准化. 筛选数据标准化处理前 (A) 后 (B) 热图。热图以渐变色标记各板孔相应信号强度。
   
   筛选数据显示β-catenin水平的分布呈现类似正态分布的T分布，如图7所示。根据T分布的规律，以及处理组与对照组两组数据间的p-value，我们设定变化显著候选因子的阈值为t-score<0.1 和p-value<0.05，详细数据处理过程参照实验方法。
   
   
   图7. 单基因筛选数据T分布图. 分别以细胞质、细胞核的β-catenin平均荧光强度或核质比 (横坐标) 和数据密度 (纵坐标) 作图。数据整体呈现类似正态分布的T分布，大部分数据都集中于1附近，表征没有变化，大于1信号强度代表升高，小于1代表信号强度降低。
   
   
4. 全基因组单基因siRNA文库筛选结果。
   我们对小鼠全基因组siRNA文库进行高通量功能缺失筛选，对其结果进行分析后 (步骤详见实验方法)，发现此次筛选得到的候选因子中包含了很多已知的Wnt信号通路的调控子，如表1所示。
   
   表1. 筛选结果中包含部分已知的Wnt信号通路调控子
   

   Barcode	Well	GeneSymbol	Gene ID	Intensity
   DSIMGA22	B06	Prmt1	15469	1.69
   DSIMGA13	B19	APC	11789	3.46
   DSIMGA40	C07	Axin1	12005	1.63
   DSIMGA55	H10	LECT2	16841	2.04
   DSIMGA44	H18	2610109H07RIK	70433	1.48
   DSIMGA35	I03	UBE2R2	67615	1.59
   DSIMGA47	J04	NKD1	93960	1.61
   DSIMGA41	J21	Podnl1	244550	1.43
   DSIMGA35	K05	Huwe1	59026	1.49
   DSIMGA05	O13	PRDM5	70779	1.63
   DSIMGA01	K18	LRP6	16974	0.47
   DSIMGA03	G16	MSX2	17702	0.51
   DSIMGA05	J07	DDX3X	13205	0.44
   DSIMGA09	D09	WDR5	140858	0.65
   DSIMGA10	I21	PIAS4	59004	0.37
   DSIMGA12	C12	CATNB	12387	0.28
   DSIMGA14	M18	RUVBL1	56505	0.45
   DSIMGA18	H19	Calcoco1	67488	0.41
   DSIMGA20	B06	STAU1	20853	0.52
   DSIMGA27	G08	FZD7	14369	0.57
   DSIMGA29	D07	Trrap	100683	0.55
   DSIMGA29	F19	Cdc2l1	12537	0.55
   DSIMGA29	L04	CSNK2B	13001	0.4

   
   
   但是，出乎意料的是，在我们的筛选系统中沉默Wnt信号通路关键成员比如Dvl、β-TrCP和GSK3等，并没有检测到相应的β-catenin变化。这个意外的现象提示我们，是否存在同家族基因成员间的代偿效应，从而导致沉默单个基因观察不到表型？事实上，Wnt信号通路中存在的基因功能冗余现象已有报道((Holder and Lewis, 2003; Huang da et al., 2009; Mao et al., 2016)，据此我们推测，针对这部分已知Wnt信号调控子的阴性结果很有可能是由于基因功能冗余造成的。为了验证我们的猜想，我们进行了基因家族组合沉默实验。通过比较我们发现，同时沉默Dvl家族成员Dvl1/2/3时，Wnt3a诱导的β-catenin水平显著下降，而单独沉默任何一个成员，都不能观察到明显的变化，如图8A所示。类似的结果也在β-TrCP1/2家族中呈现，如图8B所示。这些结果说明，针对具有类似功能的基因家族 (功能基因簇) 进行整体沉默有助于消除基因之间代偿效应，更可能帮助我们观察到表型的变化，从而发现基因 (或基因家族) 的新功能。
   
   
   图8. 基因家族组合沉默测试. (A) 沉默Dvl家族能抑制Wnt3a诱导的β-catenin水平，但单独沉默任何一个Dvl1/2/3成员都观察不到明显的变化。(B) 沉默β-TrCP1/2家族促进β-catenin累积，但单基因沉默敲低观察不到明显的变化。
   
   
5. 功能基因簇图谱的建立
   我们把小鼠全基因组19,059个基因分成了5850 个基因家族和4779个单基因。其中，三个基因和两个基因组成的基因家族的基因数量分别占全基因组基因数的41%和34%。图9统计了所有超家族包含的家族基因数。在完成了基因家族分类之后，我们运用cherry-picking和pooling等程序，使用Beckman液体工作站在无菌的条件下完成了基因家族siRNA组合文库的构建。
   
   
   图9. 基因家族分类总览. 纵轴指示含有不同基因数目的超家族，柱子左边的数字代表超家族基因成员数，柱子右边的数字代表相应成员数超家族在基因组中的数量。红、绿、黄分别表示3个，2个，1个基因成员的基因家族。
   
   
6. 基因家族siRNA文库筛选结果
   完成了基因家族siRNA组合文库的构建后，我们针对基因家族siRNA组合文库进行了高通量功能缺失筛选。我们采用和单基因筛选相同的流程及数据分析方法，同样进行了独立的三次生物学重复筛选，筛选数据中值标准化如图10所示，数据整体分布如图11所示。
   
   
   图10. 基因家族筛选数据标准化(A，B) 筛选数据标准化处理前 (A) 后 (B) 热图。热图是以渐变色标记各板孔相应信号强度。
   
   
   图11.基因家族筛选数据T分布图 分别以细胞质，细胞核的β-catenin平均荧光强度或核质比 (横坐标) 和数据密度 (纵坐标) 作图。数据整体呈现类似正态分布的T分布，大部分数据都集中于1附近，表征没有变化，大于1信号强度代表升高，小于1代表信号强度降低。
   
   根据前文所述的数据处理方法，我们对处理后的单基因和基因家族筛选结果进行了初步比较分析。单基因和基因家族筛选总览结果如图12散点火山图所示，其中红色标记的变化候选因子同时符合β-catenin水平变化强度t-score<0.1和统计学上显著因子p-value<0.05这两个标准。从图上可以看出，沉默Wnt/β-catenin信号通路关键成员家族，如β-TrCP1/2、Dvl1/2/3和Gsk3α/β等，都能显著影响Wnt3a诱导的β-catenin水平，相反，在单基因筛选中他们都未出现在候选因子名单中。这个结果说明，基于基因功能簇的筛选策略能够克服基因功能冗余导致的假阴性现象。
   
   
   图12. 单基因和基因家族筛选结果散点火山图. 以log2 (β-catenin intensity)为横坐标，-log10 (p-value)为纵坐标， p-value根据处理组和对照组的三次独立实验计算得出，同时符合变化强度 (t-score<0.1) 和统计学显著差异 (p-value<0.05) 标准的候选因子标记为红色。部分Wnt信号通路关键成员标记为黄色背景。
   
   与此同时，我们对单基因筛选和基因家族筛选数据进行了进一步整合分析，结果如图13所示。具体来说，我们把基因家族筛选结果中具有统计学意义变化的候选家族与该家族各个成员的单基因筛选结果进行比较，并把这些家族归为三大类：至少有一个单基因具有显著变化；至少有一个单基因变化，但变化程度弱；单基因成员不变。从图13中可以看出，大部分基因家族的变化都是由家族中至少一个单基因成员的变化所引起的，这个结果说明两次筛选数据具有很高的一致性。事实上我们对后两类家族候选因子更感兴趣，很有可能的一个情况是，在单基因筛选中由于功能冗余导致它们真实的功能被掩盖而无法被发现。
   
   
   图13. 单基因和基因家族筛选整合分析饼图. 第一组：至少有一个单基因具有显著变化 (inhibition t-score<0.1，promotion t-score<0.2)；第二组：至少有一个单基因变化，但变化程度弱 (inhibition t-score<0.2 ，promotion t-score<0.3)；第三组：单基因成员不变。
   
   为了验证筛选结果的可靠性，我们对第三组分类中促进或抑制Wnt3a诱导β-catenin累积最显著的前十位候选因子进行了Western Blot实验验证，结果如图14所示，其中60%的Western Blot验证结果与筛选结果一致，这说明我们的筛选系统具有较高的可靠性。
   
   
   图14. Western Blot实验验证Opera筛选数据. L细胞转染针对相应基因或基因家族的siRNA，Wnt3a刺激或不刺激条件下，Western Blot检测-catenin水平。Western Blot验证结果和筛选结果一致的基因家族标记为斜体。
   
   最后，我们分别把单基因筛选和基因家族筛选得到的影响Wnt3a诱导β-catenin水平变化显著的候选因子进行了信号通路分析，结果显示基因家族筛选策略更能帮助富集Wnt信号通路以及肿瘤相关信号通路的因子，如图15所示。
   
   
   图15. 信号通路富集分析. 利用DAVID功能注释软件(13)对单基因筛选及基因家族筛选得到候选因子进行KEGG和Biocarta信号通路分析。
   
   以上一系列针对两种筛选结果的分析比较，充分体现了基因家族筛选策略的优越性。总的来说，我们以检测Wnt3a诱导的β-catenin蛋白质的稳定性和核质分布为出发点，证明了基因家族筛选能有效地减少基因冗余导致的假阴性现象，这种全新的针对基因家族的功能缺失筛选的概念和策略同样适用于其他的以功能缺失为基础的筛选系统。
   

溶液配方

1. PBST (PBS and 0.1% triton-X100) 
2. Wnt3a纯化流程参照http://web.stanford.edu/group/nusselab/cgi-bin/wnt/purific ation
   

致谢

本章内容改编自毛丽博士的毕业论文以及发表的文章 (Mao et al., 2016). 本课题获得国家自然科学基金 (31230044，31530094)，以及中科院先导项目 (XDB19000000) 的支持。

参考文献

1. Brideau, C., Gunter, B., Pikounis, B. and Liaw, A. (2003). Improved statistical methods for hit selection in high-throughput screening. J Biomol Screen 8: 634-647.
2. Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate - a practical and powerful approach to multiple testing.J Roy Stat Soc B Met 57: 289-300.
3. Coordinators, N. R. (2014) .Database resources of the national center for biotechnology information. Nucleic Acids Res 42: D7-17.
4. Doble, B. W., Patel, S., Wood, G. A., Kockeritz, L. K. and Woodgett, J. R. (2007). Functional redundancy of GSK-3alpha and GSK-3beta in Wnt/beta-catenin signaling shown by using an allelic series of embryonic stem cell lines.Dev Cell 12: 957-971.
5. Etheridge, S. L., Ray, S., Li, S., Hamblet, N. S., Lijam, N., Tsang, M., Greer, J., Kardos, N., Wang, J., Sussman, D. J., Chen, P. and Wynshaw-Boris, A. (2008).Murine dishevelled 3 functions in redundant pathways with dishevelled 1 and 2 in normal cardiac outflow tract, cochlea, and neural tube development. PLoS Genet 4: e1000259.
6. Finn, R. D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R. Y., Eddy, S. R., Heger, A., Hetherington, K., Holm, L., Mistry, J., Sonnhammer, E. L., Tate, J. and Punta, M. (2014).Pfam: the protein families database. Nucleic Acids Res 42: D222-230.
7. Guardavaccaro, D., Kudo, Y., Boulaire, J., Barchi, M., Busino, L., Donzelli, M., Margottin-Goguet, F., Jackson, P. K., Yamasaki, L. and Pagano, M. (2003).Control of meiotic and mitotic progression by the F box proteinβ-Trcp1 in vivo.Dev Cell 4: 799-812.
8. Gentleman, R. C., Carey, V. J., Bates, D. M., Bolstad, B., Dettling, M., Dudoit, S., Ellis, B., Gautier, L., Ge, Y., Gentry, J., Hornik, K., Hothorn, T., Huber, W., Iacus, S., Irizarry, R., Leisch, F., Li, C., Maechler, M., Rossini, A. J., Sawitzki, G., Smith, C., Smyth, G., Tierney, L., Yang, J. Y. and Zhang, J. (2004). Bioconductor: open software development for computational biology and bioinformatics.Genome Biol 5: R80.
9. Hansen, B. E. (1994).Autoregressive conditional density estimation.IntEcon Rev 35: 705-730.
10. Holder, M. and Lewis, P. O. (2003).Phylogeny estimation: traditional and Bayesian approaches. Nat Rev Genet 4: 275-284.
11. Huang da, W., Sherman, B. T. and Lempicki, R. A. (2009).Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 37: 1-13.
12. Mao, L., Liu, C., Wang, Z., Niu, X., Xue, L., Zhou, Z., Cai, Z., Yu, M., Li, Y., Wu, D. and Li, L. (2016). A genome-wide loss-of-function screening method for minimizing false-negatives caused by functional redundancy. Cell Res 26:1067-1070.
13. Shahrezaei, V. and Swain, P. S. (2008).Analytical distributions for stochastic gene expression. Proc Natl Acad Sci U S A 105: 17256-17261.
14. Zhang, J. H., Chung, T. D. and Oldenburg, K. R. (1999). A simple statistical parameter for use in evaluation and validation of high throughput screening assays. J Biomol Screen 4: 67-73.