发表于 | 分类于
| 字数统计: 4,526 | 阅读时长 ≈ 17

注：以前看综述的时候，我喜欢把原文逐字翻译为汉语，效率比较低，现在我就采用笔记的形式的来看综述，按照原文的结构，把关键点记录下来，再加上自己的理解即可。

文献信息

Hailing Shi, Jiangbo Wei, Chuan He,Where, When, and How: Context-Dependent Functions of RNA Methylation Writers, Readers, and Erasers,Molecular Cell,Volume 74, Issue 4,2019,Pages 640-650,
ISSN 1097-2765,https://doi.org/10.1016/j.molcel.2019.04.025.

摘要

细胞RNA会天然地经历各种化学修饰。这些经化学修饰的核苷酸又赋予了其自身结构的多样性，从而在各种有序的代谢与机体的功能方面发挥着各种调控功能，进而影响了基因的表达。随着对回贴(mapping)位点修饰的全转录组测序方法的发展，对精确修饰进行检测和定量的高灵敏质谱方法的进步，以及参与修饰的“效应器”(effectors，包括各种能改变各种修饰位点的酶，例如写入器(writers)和擦除器(erasers)和能识别化学修饰位点的读取器(readers)）理解的加深，最近几年有关RNA修饰的研究呈现暴发式地增长。然而由于RNA种类庞杂，表达水平有差异，因此这给RNA修饰的研究带来了很大的挑战，另外，RNA的修饰还与细胞类型，组织特异性，以及修饰效应器的定位有关，这又增加了研究难度。我们在这篇综述里，重点阐述了最近几年关于$N^6$-甲基化转换酶(m6A, $N^6-methyladenosine$)功能的研究进展，强调了环境(context)对RNA修饰调控和功能的重要性。

前言

早期研究发现，m6A主要发生在(G/A)(m6A)C序列上。最近的研究发现，m6A主偏向于出现在终止密码子和3’UTR上，因此我们可以推断，m6A可以影响基因的表达。
m6A途径的效应器(effectors)包括写入器(writers)，擦除器(erasers)和读取器(readers)，其中写入器往核苷酸上添加上甲基，擦除器反之，即去掉核苷酸上的甲基，读取器则是能够识别那些核苷酸上含有甲基的序列，它们的成员如下图的Figure 1所示：

机体的整个生物学环境会影响m6A通路，并且决定了m6A的功能。
一方面，不同的细胞，以及环境刺激能影响m6A这些效应器自身的表达水平，PTM以及细胞定位。例如，在多数细胞系中，m6A的去甲基化酶(demethylase)FTO主要位于细胞核，并且参与5%-10%的的mRNA的m6A去甲基化，但是在某些淋巴瘤细胞系中，这个比例达到40%。
另一方面，RNA分子自身的异质性又增加了m6A的研究难度，这是因为：
①RNA种类繁多，包括mRNA，tRNA，rRNA，以及其它大量的非编码RNA；
②不同类型的细胞中的RNA丰度不同；
③RNA存在着复杂的二级结构；
④RNA存在着不同的转录状态（例如非活跃时，RNA与组蛋白紧密结合，转录活跃时，则与组蛋白松开）；
⑤RNA中的顺式/反式作用元件（例如一些RNA结合蛋白，RBP, RNA binding proteins）会调控m6A效应子与RNA底物。因此说，m6A的活动与环境紧密相关。

m6A甲基化的写入器

m6A是通过一个复合物加到mRNA上的（Figure1），这个复合物有多个亚基，包括METTL3，METTL14，WTAP，VIRMA，ZC3H13，HAKAI，RBM15/15B。其中包括以下成员：

  • METTL3/14，全称是methyltransferase-like 3/14，即甲基转移酶3/14， 这两个甲基转移酶是这个复合物的核心成员，其中MELLT3起催化作用，MELLT14促进复合物与RNA结合。
  • WTAP，全称是Wilms tumor 1-associating protein，其功能是结合METTL3/14，诱导它们向底物招募以及定位。
  • VIRMA，全称是Vir like m6A methyltransferase associated，功能是将m6A与3’UTR结合。
  • ZC3H13，全称是zinc finger CCCH-type containing 13，功能是诱导写入器复合物向核内转移。
  • RBM15/15B全称是RNA binding motif protein 15/15B，功能是与尿嘧啶富集区域结合，促进某些RNAs的甲基化。

METTL3的功能

METTL3在各类细胞中的分布不同，例如在HeLa细胞中，METTL3/METTL14会形成二聚体结构与WTAP产生相互作用，然后进入细胞核形成斑点(speckle)。METTL3和WTAP中含有NLS(nuclear localization signals)，NLS的突变会导致复合物无法入核。METTL3还存在于一些细胞质中，例如在一些癌细胞系，包括HeLa，MDA-MB-231，MOLM13细胞系等。还有一些情况，例如在小鼠皮质神经元中，METTL14位于细胞质和细胞核。现在还不清楚为什么METTL3的定位为什么会出现如此差异。METTL3/METTL14的比例在不同的细胞系中也有差异，这可能会影响METTL3的定位。
写入器的招募导致了转录本上的m6A的特异性，这一招募过程可能是由TF或表观遗传标记介导的，如下下图的Figure 2所示：

METTL3在细胞核中的作用

当出现热休克时，METTL3就会定位到染色质的热休克相关的基因上，m6A就会被添加到热休克转录本上，从而诱导这些转录本在应激压力下被清除。UV诱导DNA破坏时，METTL3/14就会在2分钟内定位到UV诱导的损伤位点，同时伴随着m6A活动的增强。在人类多能干细胞的TGF-$\beta$信号转录通路转导方面，活化的SMAD2/3会与METTL3/14-WTAP相互作用，诱导m6A添加到特定的转录本上。在AML细胞中，一部分METTL3会通过METTL14非依赖方式与CCAAT增强子结合蛋白zeta(CEBPZ)的启动子区域结合。在HepG2细胞中，H3K36me3能通过与METTL14相互作用招募写入器，诱导mRNA的CDS和3’UTR上添加m6A。

METTL3在细胞质中的促进翻译作用

METTL3还能作为一个潜在的m6A读取器，在Figure 2B中，我们可以看到，在肺癌细胞中，METTL3此时并没有发挥催化作用，而是在细胞质中锚定到了3’UTR上，促进了一个报告mRNA的翻译。进一步的研究表明，METTL3是通过eIF3h相互作用来促进翻译的。
METTL3蛋白自身会通过PTM或与其它蛋白质相互作用来发挥调控功能，例如人类的METTL14能够诱导METTL3的丝氨酸399磷酸化。人类的METTL3会出现SUMOylation修饰，导致METTL3/14的活性降低。但这些现象的机制还不清楚。

序列和结构特异性诱导的甲基化

METTL3/14复合物偏爱RRACH模序(R=A或G；H=A，C或U)，METTL16则偏爱另外的序列与结构。METTL16被验证的两个特异性序列是U6(U6 small nuclear RNA, snRNA)和人类MAT2A(methionine adeno-syltransferase 2A)上的一个发夹结构(hp1)，MAT2A编码S-腺苷甲硫氨酸合成酶(S-adeno-sylmethionine synthetase, SAM synthetase)。有实验分析发现，METTL16偏爱一个形成茎环结构膨大处的UAC(m6A)GAGAA序列。EMTTL16介导的MAT2A甲基化是SAM稳态的负反馈信号。最近发现了ZCCH4是一个新的m6A读取器。

m6A读取器

m6A读取器通过作用于含有m6A的mRNA来发挥作用，Figure 1中可以把读取器分为三类。

第一类读取器含有YTH结构域（YTH的全是YT521-B homology），这些成员包括YTHDF1-3(YTH domain family 1-3)，YTHDC1-2(YTH domain containing 1-2)（参考Figure 1）。细胞质中的YTHDF2会通过招募CCR4-NOT腺嘌呤酶复合物来诱导靶转录本的部分降解。细胞质中的YTHDF1和YTHDF3能促进HeLa细胞中靶转录本的翻译。细胞核中的YTHDC1有多个作用，包括招募某些剪接因子调控mRNA的剪接，促进mRNA的输出，加速某些转录本的降解。YTHDC2调控mRNA的稳定，翻译以及精子形成。

第二类读取器是HNRNPs，全称是heterogeneous nuclear ribonucleoproteins，成员包括HNRNPC，HNRNPG与HNRNPA2B1，它们能调控靶转录本的剪接，加工，它们能与RNA的某些结构结合，这些结构是由m6A介导形成的，因为m6A会重构局部的RNA结构，调控RNA-蛋白质之间的相互作用，这一现象称为m6A开关(m6A-switch)。

第三类读取器含有相同的RNA结合结构域(RBD, RNA binding domains)，例如KH结构域(K homology)，RNA识别模序(RNA recognition motif)，以及RGG( arginine/glycine-rich)结构域，它们都能与含有m6A的mRNA结合，包括FMR1和IGF2BP1-3。FMR1(Fragile X mental retardation 1)含有3个KH结构域，1个RGG结构域，它有可能通过与YTHDF1和YTHDF2相互作用来影响RNA的转移，RNA的稳定。IGF2BP1-3(insulin-like growth factor 2 mRNA-binding proteins 1–3)能通过m6A的方式稳定靶转录本（YTHDF2与之相反，它是降解靶转录本）。IGF2BP蛋白功能的核心是KH3-4结构域，这是与m6A结合的关键区域。
Prrc2a是最近发现的一个读取器，它与髓鞘形成有关，具体机制不明。

为什么读取器会结合一些m6A位点

以下为作者的几个猜测。

第一，读取器有可能与其它的RBP相互作用，从而被招募到mRNA的不同区域。IGF2BP1-3和YTHDF2通过在不同位点来调控mRNA的稳定性，例如IGF2BP1-3结合到3’UTR，而YTHDF2结合到CDS区。而RBP对RNA的识别取决于多个因素，例如结合位点的序列，序列的侧翼以及RNA的二级结构。因此RBP与读取器蛋白质的相互作用有可能涉及了m6A位点的特异性。

第二，m6A位点的密度和序列也可能与之有关。

第三，读取器蛋白会在细胞区室的特定位点富集，因此会偏向与局部一些RNA类型结合。例如，YTHDF1-3，FMR1和HNRNPA2B1位于细胞的应激颗粒核心地区。在乳腺癌细胞中，YTHDF1-3，HNRNPK和IGF2BP2-3位于细胞的突起(protrusion)。
甲基化可以促进翻译或影响某一转录本的稳定性，这取决于揎细胞环境下哪个读取器占据主导地位，如Figure 3所示：

读取器对外界刺激的应答

在某些刺激下，例如热休克应激或病毒感染会诱导细胞质中的YTHDF进入细胞核。在热休克出现的几小时内，YTHDF2的转录本与蛋白水平都上升，并且多数YTHDF2转移到细胞核内。有人认为，YTHDF2核心的YTHDF2会与去甲基化酶FTO竞争，阻止那些热休克应答基因5‘UTR的去甲基化，从而增强这些基因在细胞质中非加帽方式的翻译(cap-independent translation)。体外实验发现，YTHDF2能抑制FTO的去甲基化活性。在Vero细胞感染实验中，细胞感染了enterovirus type 71后，YTHDF1和YTHDF2会上调，并且在感染12小时后分布在细胞质，感染24小时后分布在细胞核，这一现象机制不明。

在有丝分裂后期的细胞中，当需要功能蛋白时，YTHDF1的翻译促进作用就变得特别明显(Figure 3)。在损伤诱导的轴突再生的背根节(DRG)模型中，再生相关基因大量地被甲基化，在这个恢复过程中，YTHDF1是蛋白质强烈合成所必需的条件；相比之下，YTHDF1在损伤前的基础翻译增强方面，作用很小。这一现象与最初在HeLa细胞中研究YTHDF1时类似，YTHDF1会促进一些转录本的翻译。在关于记忆研究方面，有研究发现，YTHDF1的缺陷会导致海马依赖性神经元功能出现障碍，恢复YTHDF1的蛋白水平，则能修复这种问题。YTHDF1依赖的翻译增强会在一些癌细胞中持续出现，而在其它一些细胞中则是由刺激诱导的，这一过程是如何触发的，还不清楚。

读取器的不同功能

YTHDF1-3的结构类似，它们都含有N末端低复杂序列，和C末端的保守YTH结构域。但是它们的功能并不相同，在不同的细胞中，它们有的能促进蛋白翻译，有的能抑制蛋白翻译，具体的案例可以参考原文。并且在不同的m6A阅读器蛋白之间还存在着交叉作用或竞争作用，即使在FTO和阅读器之间，也有着类似的作用。因为这些蛋白本身就构成了一个有相互作用的网络结构，因此要理解它们的环境调控作用对于相关的研究有着重要意义。

m6A擦除器

m6A甲基化可以被FTO或ALKBH5逆转，这两个蛋白也就是m6A擦除器。FTO是第一个被发现的m6A擦除器。FTO能够去甲基化m6Am，以及部分mRNA上的cap-m6Am，如Figure 4所示：

CLIP-seq验证了很多FTO的结合位点，这些结合位点包括tRNAs，snRNA等。FTO的空间分布也能发挥调控作用，FTO的N末端有一个NLS，它能部分地分布在细胞核中，也能分布在细胞质中，FTO的分布在不同的细胞系中有所不同，例如在AML细胞中的分布和HEK，HeLa细胞中的分布就不同，这可能与肿瘤的发生有关，FTO在AML细胞中的这种分布也有可能是环境依赖的，因为2HG(2-hydroxy-glutarate)能抑制FTO。

CLIP-seq的测序数据分析表明，在某些细胞系中，FTO偏爱结合GAC-和/或GGAC模序。有研究发现，FTO的表达和定位可能被PTM调控（例如Lys-216）。
$cap-m^{6}A_{m}$位于mRNA的5’cap中的+1核糖上，它几乎是与内部m6A(internal m6A，这个内部指的是mRNA中间的部分)同时发现的。m6A可能占据了整个腺苷酸的0.1%-0.4%，它主要发生在第二个核苷酸2’-O-methylated($m^{7}GpppN_{m}$)的帽子结构附近。只有不到30%的$m^{7}GpppN_{m}$含有m6Am。m6Am的比例占到整个mRNA的m6A十分之一或更低（在一些细胞系中还要低）。m6A丰度要远高于$m6A_{m}$，因此，虽然FTO更偏爱结合$m6A_{m}$，但是FTO的主要靶点还是m6A。由于几乎所有的mRNA帽子结构都在胃镜核中被加帽蛋白结合，因此在细胞核中，cap-$m6A_{m}$不太可能被去甲基化，这也反映了FTO对$cap-m^{6}A_{m}$的局限，如Figure 4所示。现在研究发现的FTO还功能涉及：热休克应答，UV损伤应答，HCV感染等方面。

mRNA的$cap-m^{6}A_{m}$甲基化不影响转录本的稳定

最初研究认为FTO介导的$cap-m^6A_{m}$的去甲基化在功能上类似于DCP-2介导的去帽(decapping)化诱导mRNA的降解以及miRNA介导的mRNA降解。当敲低FTO后，FTO的靶mRNA出现了很多变化，这些效应被认为是由于$cap-m^6A_{m}$的去甲基化导致的，但是FTO敲低与那些只含有m6A或$m^6A_{m}$转录本水平之间的相关却不支持这个结论，并且发现只有内部m6A修饰的转录本的变化与FTO的敲低相关。

后来发现了PCIF1是mRNA的$cap-m^6A_{m}$甲基转移酶，它只加工$cap-m^6A_{m}$，却不加工内部的m6A(internal m6A)。还有研究发现，PCIF1添加的$m^6A_{m}$并不改变基因表达的水平或转录本的稳定性，这也不支持刚开始的观点，即FTO介导的$cap-m^6A_{m}$的去甲基化导致的mRNA降解。另外，有研究发现，当FTO敲低时，$cap-m^6A_{m}$被认为还能促进那些含有$cap-m^6A_{m}$mRNA的翻译效率，但也有一些研究发现了相反的现象，这一过程的调控也有可能是与环境有关。

ALKBH5的作用与底物

ALKBH5是第二个被发现的m6A去甲基化酶。ALKBH5在小鼠的精子形成中发挥了作用，ALKBH5和FTO的表达在不同组织中各异。例如，在小鼠中，Alkhb5主要表达在睾丸，而Fto主要表达在大脑，这可能与它们参与的不同生物学途径有关。在人源细胞系中敲低ALKBH5后，会诱导其靶RNA从细胞核转移到细胞质。几乎所有报道的ALKBH5功能研究都揭示了类似的分子途径，ALKBH5介导某些转录本的3’UTR m6A，其中就包括促进缺氧诱导的HIF依赖的乳腺癌干细胞表型，通过ALKBH5-FOXM1途径调节的胶质母细胞瘤增殖和肿瘤发生，调控雄性生殖细胞中长3’UTR mRNAs的剪接和稳定。

结论与展望

作者在结论部分中讨论了两个问题，第一个是m6A的调控，第二个是m6A未解决的一些问题。

其中m6A的调控功能可能受几个层面的影响，包括以下几个方面

第一，m6A的功能取决于细胞分化和细胞的发育状态。

第二，m6A的功能可能被环境因素或细胞信号触发。

第三，m6A效应子的定位也影响了m6A的功能。

第四，即使是同一个转录本，m6A精确的定位和其它的修饰也会影响其功能。

关于m6A的一些关键问题还没搞清楚，包括以下几个方面：

第一，m6A效应子针对转录本的选择性和m6A位点的选择性是如何实现的？

第二，m6A效应子（包括写入器，擦除器，读取器）是如何整合到不同的生物信号转导与调控过程的？

发表于 | 分类于
| 字数统计: 15,547 | 阅读时长 ≈ 56

前言

这篇文献是邵峰于2019年11月发表在《Nature Reviews Immunology》上的综述，对gasdermin的最近研究做了总结，对天然免疫有兴趣的同学可以读一读。

注：由于gasdermin没有相应的中文译名，我个人觉得可以译为“焦孔素”比较好。

文献信息

Broz, P., et al. (2019). “The gasdermins, a protein family executing cell death and inflammation.” Nature Reviews Immunology.

文献摘要

gasdermins是最近鉴定出来的一个具有成孔效应蛋白家族，这个家族的能造成膜通透性(permeabilization)与细胞焦亡(pyroptosis)，细胞焦亡是一种裂解性的，炎症性的细胞死亡形式。gasdermins含有一个细胞毒性的N末端结构域和一个C末端的抑制结构域，它们通过一 个柔性构件连接。这两个结构域之间的连接蛋白水解后，就会解除细胞毒性结构域的分子内抑制作用，使N末端结构域插入到细胞膜中，形成大寡聚孔，从而破坏细胞的离子稳态并诱导细胞死亡。gasdermins诱导的细胞焦亡在许多遗传病、自身炎症性疾病和癌症中扮演着重要的角色。篇综述讨论了gasdermin最新的研究进展，重点涉及gasdermin活化、孔形成和gasdermin诱导的膜通透化。

相关术语

  1. inflammasome（炎症小体）：一种多蛋白信号转导复合物，当检测到宿主来源或病原体来源的危险信号时，炎症小体就会在细胞质开始组装，促进细胞因子的释放，细胞的焦亡型死亡以及炎症。
  2. caspases（半胱天冬氨酸蛋白酶）：半胱氨酸依赖性天冬氨酸蛋白酶家族，它作用于底物上特定的天冬氨酸残基，这个家族的酶在细胞死亡和炎症中发挥着核心功能。
  3. necrosis（坏死）：细胞死亡方式之一，细胞经历这种死亡形式时，细胞膜丧失完整性，释放细胞成分。
  4. NLRP3：英文全称是Nucleotide- binding oligomerization domain, leucine- rich repeat and pyrin domain- containing 3，即核苷酸结合寡聚结构域，富含亮氨酸重复序列和含Pyrin结构域3，这是组装炎症小体的细胞质感受器。NLRP3是一个广谱感受器，它能识别病原体相关和损伤相关分子模式。
  5. Pannexin-1 channels（膜联蛋白-1通道）：一种细胞膜通道，可以通透离子和小分子代谢物，例如ATP。
  6. NETosis：中性粒细胞的一种特殊死亡方式，其特征是释放DNA，组蛋白和颗粒内容物到细胞外空间。
  7. Liposomes（脂质体）：由磷脂组成的，含有至少一个脂质双层的人造球形囊泡。
  8. Mitophagy（线粒体自噬）：在营养饥饿或线粒体应激的情况下，通过自噬选择性地去除线粒体。
  9. ASC foci（ASC焦点），英文全称为apoptosis-associated speck-like protein containing a caspase recruitment domain foci，即含有caspase募集结构域焦点的凋亡相关斑点样蛋白，这是一个多聚体蛋白，由炎症小体蛋白的同源寡聚引发的螺旋原纤维形成。
  10. Familial Mediterranean fever（家族性地中海热）：一种常染色体隐性遗传性自身炎症疾病，由MEFV基因突变引起，该蛋白编码蛋白Pyrin，疾病的特征是自限性发热和浆膜炎(Serositis)。
  11. Pyrin：由细菌诱导的RHO修饰和肌动蛋白细胞骨架动力学破坏的炎症小体形成的细胞质感受器。由MEFV基因编码。
  12. Peroxisomes（过氧化物酶体）：真核细胞的过氧化物酶体细胞器参与不同脂肪酸的分解代谢和活性氧的减少。
  13. NLRP1：英文全称为(Nucleotide- binding oligomerization domain, leucine- rich repeat and pyrin domain- containing 1，即核苷酸结合寡聚结构域，富含亮氨酸重复序列和含1个Pyrin结构域，这是一个细胞质感受器，组装成炎症小体。它通过N末端规蛋白酶体降解途径(N-end rule proteosomal degradation pathway)活化。
  14. NLRC4：英文全称为Nucleotide-binding oligomerization domain, leucine-rich repeat and caspase recruitment domain-containing 4，即核苷酸结合寡聚结构域，富含亮氨酸重复序列和含caspase募集结构域-4，一种炎症小体支架蛋白(scaffold protein)，与炎症小体感受器(NLR家族凋亡抑制蛋白(NAIP))结合时寡聚化，并形成活化caspase 1的平台。
  15. Toll样受体：位于内体和细胞膜上的一类模式识别受体，识别病原体相关和损伤相关的分子模式，并启动信号通路以诱导炎症。

前言

gasdermin基因家族最早于20世纪初报道，该家族的基因当时是作为小鼠皮肤突变的候选基因来进行研究的。gasdermin这个名称是基于gasdermin A(GSDMA) 蛋白只表达在小鼠的肠道(gastrointestinal tract )和皮肤上皮细胞(epithelium of the skin)中。早期的研究还发现，gasdermin与DFNA5（deafness autosomal dominant non-syndromic sensorineural 5）的N末端区域的序列非常相似，DFNA5这个蛋白于1998年发现，它与人类的非综合征型常染色体显性遗传性有关。基于这种同源特性，一些其它的gasdermin家族的成员以及gasdermin样蛋白也被陆续鉴定出来了，目前在人类中有6个同源基因( paralogous genes)，即GSDMA，GSDMB，GSDMA，GSDMD，GSDME（也被称为DFNA5）和PJVK（也被称为DFNB59）（Fig.1a)。

gasdermin家成员是基于序列同源进行鉴定的，它们表示出不同的组织表达特征（Box 1与Table 1）。但是在过去的15年里，这些蛋白的精确功能却并不清楚。尽管如此，在这些蛋白被鉴定后不久，人们就发现了这些蛋白与细胞活力和炎症之间存在着某些联系。例如在小鼠方面，Gsdama3（编码小鼠的GSDMA3蛋白）的9种突变体与小鼠的脱毛有关，并伴随着明显的干细胞耗竭，角化过度(hyperkeratosis)以及炎症。另一方面，敲除Gsdma3基因后，小鼠并不表现出肉眼可见的皮肤表型，这说明，Gsdma3基因的突变体增加了某些功能。表明gasdermins潜在细胞毒性的最直接证据的研究指出，使用某些突变的，C末端截短式的人源GSDME能够导致酵母细胞的细胞周期停滞，以及人类细胞（凋亡型）的死亡。不过，gasdermins是否是细胞死亡程序体(cell death programme)的组成部分，以及这些蛋白控制的细胞死亡类型仍然不清楚。

2015年，两项独立的研究揭示了gasdermin功能的机制，随后不久，第三项研究证实了gasdermin的作用机制。这些研究表明，GSDMD是焦亡的唯一执行者(BOX 2)。自从这一发现以来，大量的文献开始描述gasdermins在炎症生物学、细胞死亡等方面的作用。在这篇综述中，我们讨论了对gasdermin激活和调节、gasdermin孔隙的组装和gasdermin蛋白家族相关的生物学功能的最新进展。

Fig.1：gasdermin家族

Fig1:gasdermin蛋白家族。

A. 系统发育树显示了人，小鼠和大鼠gasdermin(GSDM)家族蛋白的差异。人gasdermin蛋白(GSDM)家族由6个基因编码，其总序列相似性范围为23.9~49.4%。这些家族可以被进一步细分，因为GSDME和pejvakin(PJVK)也属于耳聋相关基因(DFN)，这两类蛋白质序列聚在一起，远离其他人类gasdermins(GSDMA-GSDMD)。从进化角度来讲，GSDME和PJVK是最古老的gasdermin成员，在低等脊椎动物和一些无脊椎动物中也发现了相似的序列，但在蠕虫和苍蝇中没有发现。GSDMA基因序列存在于哺乳动物以及鸟类和爬行动物中，而GSDMB、GSDMC和GSDMD基因仅存在于哺乳动物基因组中，并且与GSDMA密切相关，这表明它们是通过基因复制产生的。小鼠和大鼠缺乏GSDMB，但除了Gsdmd，Gsdme和Pjvk外，小鼠还具有三个GSDMA同源物(Gsdma1-Gsdma3)和四个GSDMC同源物(Gsdmc1-Gsdmc4)。从结构上讲，gasdermins由两个不同的结构域组成，通过一个柔性连接区连接，除了PJVK呈现出一个较小的C-末端结构域(CT)。N末端gasdermin(GSDM NT)结构域在所有家族成员中显示出最高的序列相似性(相似性范围为28.8%至50.5%)。相比之下，C-末端gasdermin(GSDM CT)结构域呈现可变长度和较低的相似性(范围从1.3%到46.3%)。GSDM NT具有固有的成孔/诱导焦亡活性，而GSDM CT与GSDM NT相互作用，从而在没有活化信号的情况下抑制其活性。标度表示序列中每个氨基酸的替换次数。上述系统发育树由欧洲生物信息学研究所(EMBL-EBI)Clustal Omega tool 从UniProt序列生成的，并由Figtree软件版本1.4.3绘制。**

B. 人GSDMD和GSDME的结构域结构，其特征在于连接区中的caspase裂解序列(顶部)和膜插入的N-末端结构域(NT；底部)。

Box1：gasdermins的鉴定和表达模式

  • gasdermin A(GSDMA；也称为GSDM，GSDM1或FKSG9)首先从小鼠皮肤中克隆发现，其表达主要局限于人的食管、膀胱和皮肤的上皮细胞。T淋巴细胞也能表达出可检测到的GSDMA蛋白。小鼠Gsdma1、Gsdma2和Gsdma3的表达也仅限于上皮和皮肤，包括表皮、毛囊和胃。
  • GSDMB(也称为GSDML，PP4052或PRO2521)是通过使用GSDMA序列作为诱饵进行数据库同源性搜索发现的。GSDMB是gasdermin家族中最具分散特征的成员(另见图1a)，它不存在于小鼠和大鼠基因组中，尽管一些啮齿动物物种具有GSDMB同源的序列。GSDMB的表达主要在气道上皮，食道，胃，肝脏，小肠和结肠等组织中也检测到。已经在人类中检测到不同的GSDMB剪接可变体，其中一个转录本编码为GSDMB，在域间连接区(interdomainlinker)(由外显子6编码)中具有caspase 1裂解位点。过表达实验表明，caspase 1可以切割这种异构体并诱导裂解性细胞死亡。
  • GSDMC(也称为MLZE)首先被鉴定为在转移性小鼠黑色素瘤细胞中表达上调的基因，后来被鉴定为gasdermin家族的成员。小鼠基因组包含四个Gsdmc同源基因。GSDMC的表达仅限于食管、皮肤、脾脏和阴道。人工截短的N端GSDMC(GSDMC NT)能够诱导焦亡，但是什么信号可以激活GSDMC以及如何激活GSDMC仍然不清楚。
  • GSDMD(也称为GSDMDC1，DFNA5L或FKSG10)首先通过在人类基因组数据库中搜索GSDMA的同源物而鉴定出来的。GSDMD在不同的人类组织以及白细胞的不同亚群中广泛表达。GSDMD直系同源基因仅存在于哺乳动物基因组中，并且都包含一个大的中央结构域，具有caspase 1和小鼠caspase 11或人caspase 4和5的切割位点。然而，值得注意的是，在较低等的脊椎动物，例如斑马鱼中，caspase a(caspy)或caspase b(Caspy2)是人类caspase 1和人类caspase 4或5的同源物，据报道它们能诱导细胞死亡并参与免疫。这表明GSDMD的功能同源物可能存在于这些低等脊椎动物中，尽管不能排除存在其他替代的细胞死亡途径。在哺乳动物中，caspase 1切割前体pro-IL-1β，产生成熟的有生物活性的IL-1β细胞因子。同时，由caspase1、小鼠caspase11或人caspase4或5切割GSDMD导致形成高度裂解的GSDMD NT蛋白片段，其允许释放成熟的IL-1β(Fig. 2)。IL-1β的低级脊椎动物序列缺少保守的caspase 1切割位点，但是，尽管如此，鱼类caspase a，caspase b或IL-1β的抑制在感染期间对宿主是有害的。Gsdmd和pro-IL-1β在哺乳动物中都具有caspase 1切割位点，因此可以赋予炎症小体对IL-1β信号传导的两个关键步骤的调控，即其加工和释放。相比之下，这两个步骤可能由低等脊椎动物中的同源蛋白酶(如caspase a或caspase b)和GSDMD同源物控制。
  • GSDME(也称为ICERE-1或DFNA5)最初被克隆为常染色体显性非综合征性听力损失的候选基因，后来被发现与gasdermins具有序列和结构上的相似性。GSDME在不同的人类细胞和组织中有不同的表达，包括脑、子宫内膜、胎盘和肠等。在小鼠和人类中，GSDME由caspase 3加工，它直接诱导或在细胞出现凋亡形态后间接诱导焦亡。
    GSDME在不同种类的低等脊椎动物中也有表达，例如，GSDME的两个同源基因(GsdmEa和GsdmEb)可以在硬骨鱼(bony fish)中找到。Caspase 3切割位点在斑马鱼GsdmEa中存在，但在GsdmEb中不存在，这表明GsdmEa可以被认为是GSDME的功能同源物。目前尚不清楚GsdmEb是否由鱼类Caspy或caspy2加工，但如果是这样，GsdmEb可能作为哺乳动物GSDMD的功能同源物(见上文)。有趣的是，斑马鱼中GsdmEb的缺失会导致耳朵的半规管畸形，这表明GsdmEb可能导致与人类GSDME相关的听力损失。
  • Pejvakin(PJVK；也称为DFNB59或GSDMF)是另一种与耳聋相关的突变蛋白，但最初是从人类睾丸中克隆得到的。PJVK与GSDME高度相似，PJVK同源基因存在于早期脊椎动物和无脊椎动物中，这表明gasdermin蛋白家族可能是从这些祖先进化而来的。PJVK在睾丸中表达较高，但在其他组织中也广泛表达，包括内耳的毛细胞和听觉系统的其他细胞(TABLE 1)。到目前为止，还不清楚PJVK是否能够被蛋白酶处理，以及它的N末端或全长PJVK是否形成膜孔。

Fig.2：gasdermin D的经典和非经典炎性小体活化

Fig. 2。gasdermin D在经典和非经典的炎症小体活化方面的作用。
A. 模式识别受体例如pyrin，AIM2，NAIP-NLRC4，NLRP3和NLRP1产生信号后，炎症小体复合物就会组装起来。这识别受体能识别PAMP，内源性危险信号或细胞死亡，损伤或感染引起的细胞稳态改变。
这些受体通过同型或异型PYD/CARD(caspase活化和招募结构域)相互作用来招募含有caspase募集结构域(ASC)的衔接蛋白凋亡相关speck样蛋白的，从而进一步招募pro-caspase 1或直接招募pro-caspase 1。caspase 1在炎症小体中被活化，并且caspase 1加工gasdermin D(GSDMD)以及pro-IL-1β和pro-IL-18等细胞因子。当细胞膜被GSDMD孔通透化时，细胞经历溶裂解、促炎性细胞死亡（焦亡），促进成熟的IL-1β和IL-18的释放。在没有细胞裂解的情况下，GSDMD孔还可以直接释放细胞因子，如IL-1α，危险分子，如高迁移率族蛋白-1(HMGB1)，甚至包括caspase 1在内的整个炎症复合物。

B. 非经典的炎症小体通路会导致小鼠caspase 11的活化（人类则是caspase 4和caspase 5）。来自革兰氏阴性菌的脂多糖(LPS)结合并诱导这些caspase的寡聚化和活化，使它们能够切割GSDMD。在第一步中，GSDMD孔允许钾释放，导致NLRP3炎症小体的活化和IL-1β/IL-18的成熟。在第二个步骤中，GSDMD孔导致焦亡，从而驱动成熟细胞因子的释放。由GSDMD孔引起的膜损伤通过运输(ESCRT)机制所需的内体分选复合体进行修复。

NAIP, NLR family apoptosis inhibitory protein; NLR , nucleotide- binding oligomerization domain, leucine- rich repeat- containing protein; NLRC4, NLR with CARD domain- containing 4; NLRP1/3, NLR with pyrin domain- containing 1/3.

TABLE 1:gasdermin表达谱

gasdermin活化

在20世纪80年代和90年代，焦亡被认为是由毒素刺激或病原体感染诱导，并导致巨噬细胞由caspase-1介导的细胞死亡形式。后来的研究表明，促炎caspases（这是在炎症复合体中活化的一组proteases）的主要效应机制就是焦亡。炎症小体的活化有两种不同的途径，即经典炎症途径和非经典炎症途径，它们能识别病原体来源或宿主来源的危险信号，并开启人源caspase 1和caspase 4的活化或鼠源caspase 11的活化。这些caspases能够在GSDMD的中央连接区（人类中的是FLTD，小鼠中的是LLSD）切割GSDMD(Fig. 1b)，产生两个片段，一个是31kDa 的N末端GSDMD NT片段，该片段有着内在的成孔活性，另外一个是22kDa的C末端GSDMD CT片段，它能与GSDMD NT结合，作为抑制后者的抑制剂。GSDMD NT本身的表达就能诱导细胞焦亡，而GSDMD CT的过表达则会阻止细胞死亡。

这些结果就构建了一个caspase介导的GSDMD切割的模型，切割后生成的GSDMD NT从GSDMD CT的抑制中释放出来（也就是说，GSDMD NT与膜磷脂的结合后，分子内的相互作用就会很容易破坏细胞膜），因此诱导细胞的焦亡。其它的一些发现也支持了这个模型，gasdermin家族的成员多数都含有类似的结构（除了PJVK，pejvakin含有一个截短的C末端结构域），也就是说含有一个N末端细胞毒性结构域和一个C末端抑制结构域，这两个结构域由一个构件连接，异位表达GSDMA，GSDMB，GSDMC或GSDME的N末端结构诱导会诱导坏死(necrosis)，这种形式的细胞死亡在形态类似于GSDMD诱导的焦亡。因此，gasdermin家族作为一组新的细胞死亡效应因子就被发现了，这些效应因子由其能诱导N末端焦亡的形式来发挥作用。

GSDMD中含有caspase 1切割位点，而在其它的gasdermins成员中，除了GSDMB的一个较小的剪接突变体外，其余的成员并不含有这个位点，不过它们可能含有其它的protease切割的位点。事实上，人源和鼠源GSDME的连接区(Fig. 1b)含有一个caspase 3剪切基序(motif)，另外有研究指出，当受到凋亡刺激后，活化的caspase 3会在这个位点切割GSDME(Fig. 3)。根据这一发现，有人指出，当细胞已经以进入凋亡程序后，GSDME会导致细胞出现焦亡，丧失膜完整性。后来的研究则不同意该假设，并表明，尽管有些细胞在自然条件下含有高水平的GSDME，并表现出一些凋亡的特征，但是GSDME也还能直接诱导焦亡；另外还发现，Gsdme缺陷的巨噬细胞仍能通过某些未知程序，出现二级坏死(secondary necrosis)，进入凋亡程序。有意思的是，GSDME不是唯一在凋亡诱导后被活化的gasdermin成员。使用药物诱导或病原体诱导的激酶TAK1或抑制凋亡(IAP)，能够诱导GSDMD的切割，并且这一过程并不依赖于caspase 1和caspase 11。有人指出，在这种条件下，GSDMD直接由caspase 8切割(Fig. 3)，这与caspase 8可以在体外加工GSDMD的事实是相符的，尽管这一过程比较慢。也有人认为，由caspase 8驱动的GSDMD孔隙形成过程是造成凋亡细胞钾外流和NLRP3炎性小体活化的原因，但其它的发现指出，这一过程是由 pannexin-1通道活化的caspase驱动介导的。有意思的是，凋亡细胞中的GSDMD的活性受到caspase 3的负调控，caspase 3通过切割GSDMS的N末端结构域来发挥作用，进而产生失活片段。caspase 8活化GSDMD的生理功能和caspase 3的负调控作用还有待进一步确定。然而现在已经明确，GSDMD可以在“调亡”刺激后造成裂解性细胞死亡和炎症。

现在人们认识到，导致gasdermin活化的通路已经扩展到了caspase家族之外。两项研究表明，在活化的中性粒细胞中，细胞中的弹性蛋白酶(elastase)能切割GSDMD，这是一种丝氨酸蛋白酶，它对中性粒细胞的成熟和抗菌功能非常重要。虽然弹性蛋白酶在炎性caspases作用基序的上游几个位点处切割GSDMD，但是这种切割仍然能产生功能性孔形成片段。GSDMD在中性粒细胞中的精确机制现在仍然有争议。已经有研究发现，Gsdmd缺陷的小鼠对大肠杆菌更有抵抗性，这可能是由于中性粒细胞的寿命延长所致。另一方面，有人认为，GSDMD与中性粒细胞的死亡（注：原文中使用的术语是NeTosis，这是一种中性粒细胞的死亡方式）的有关。因此，GSDMD在中性细胞中同时发挥着好坏双重作用，这要取决于感染的类型以及是否需要中性粒细胞存活或死亡来限制病原体。

虽然诱导gasdermins活化的唯一已知生理机制就是蛋白水解，但是有几条证据表明，移除C末端gasdermin(GSDM CT)结构域并不是gasdermin活化的绝对条件。某些突变会导致自身抑制性结构域的相互作用被破坏，从而触发gasdermin的活化，这表明，C末端结构域本身的存在并不会干预孔的形成。GSDMA3的晶体结构揭示了结构域间相互作用的关键特征是GSDMA3 NT的α1螺旋和β1-β2发夹环，它们深入GSDMA3 CT的疏水核的凹槽里。GSDMA，GSDMA3，GSDMC，GSDMD和GSDME的C-末端结构域的疏水核心突变都会导致焦亡，并且导致脱毛的Gsdma3突变体的几个映射到C-末端结构域和结构域间相互作用界面。因此，生理信号通路有可能通过解除自身抑制，例如通过磷酸化或其他翻译后修饰，以类似的方式诱导gasdermin的活化。

Fig. 3：gasdermin被凋亡型caspase诱导的活化

Fig. 3 “凋亡型”caspases诱导gasdermins的活化。

某些外部刺激，例如肿瘤坏死因子(TNF)或Toll样受体3和4的配体(TLR3和TLR4)，都能抑制凋亡抑制因子(IAP)蛋白，第二线粒体来源的胱氨酸酶活化剂(second mitochondria-derived activator of caspase,SMAC)模拟化合物或病原体诱导的TAK1激酶抑制可以促进激酶受体相互作用蛋白1(RIP1)依赖的胞质caspase 8激活复合物(复合物IIb/坏死性凋亡小体(Ripoptosome))的组装。活性caspase 8通过活化效应性caspase，如caspase 3来驱动细胞凋亡，但也切割gasdermin D(GSDMD)以产生活性N-末端片段。活性GSDMD可以诱导孔形成，但仅局限在天冬氨酸D87(小鼠D88)处受到caspase 3依赖的切割，从而产生GSDMD的非活性p20/p10片段。然而，caspase 3也可以切割和激活GSDME，从而在GSDME高表达的细胞中将细胞凋亡转化为焦化。GSDME诱导的焦亡不同于继发性坏死，继发性坏死以GSDMD/GSDME非依赖性方式进行。

FADD, Fas associated via death domain; SMAC, second mitochondrial- derived activator of caspases; TNFR1, tumour necrosis factor receptor 1; TRAM, TRIF- related adaptor molecule; TRIF, TIR domain-containing adapter- inducing IFNβ.

gasdermin的孔形成机制

焦亡的特征在于由gasdermin孔形成，进而引起的细胞膜破裂。体外结合试验表明，gasdermin N末端可以直接与膜脂相互作用。有研究报道，GSDMD NT优先靶向作用于酸性磷脂，例如磷脂酰肌醇和心磷脂，但也可以弱结合到磷脂酸和磷脂酰丝氨酸上。 其它的gasdermins例如GSDME，GSDMA和鼠源GSDMA3的N末端结构域表现出类似的脂质结合特性，这表明整个gasdermin家族具有共同的膜靶向机制。磷脂酰肌醇仅存在于细胞膜的细胞质小叶中。与此一致的是，N末端gasdermin(GSDM NT)结构域只能诱导细胞内的焦亡，在细胞外添加活化的gasdermin并不会造成细胞的膜裂解。心磷脂具有类似于磷酸肌醇的带负电荷的头部结构，它们存在于真核生物和细菌的线粒体内膜中。在大肠杆菌中表达GSDM NT结构域会表现出严重的毒性，重组GSDM NT蛋白可以裂解巨大芽孢杆菌(Bacillus megaterium)的原生质体。已经有研究表明，GSDMA3 NT和GSDMD NT单独或当不受分子内抑制时可以破坏线粒体，并且当细菌暴露于重组的GSDMD NT时，细菌的生成会受到抑制，不过还有待进一步的研究才能确定gasdermin的成孔结构域是如何进入内膜的。除了磷脂能够与GSDM NT结构域发生强烈的相互作用外，至于其它的膜脂结构，即使它们没有表现出与GSDM NT的特异性和强烈的结合，也有GSDM NT结构域通过影响膜的物理特性来发挥作用。例如，鞘磷脂的存在可以极大地促进GSDMD NT与脂质体的结合，而将胆固醇包在脂膜中则显著降低了与GSDMD NT的结合。虽然对于大多数gasdermins来说，只能用游离的GSDM NT结构域观察到脂质结合，但全长GSDMB则表现出与单独的GSDMB NT结构域类似的脂质结合能力，这表明GSDMB CT结构域不妨碍GSDMB的脂质结合功能。 除了磷酸肌醇结合外，GSDMB还与硫脂能特异结合，但这种结合的生理相关性仍有待确定。

考虑到所有GSDM NT结构域之间序列的高度相似性，因此可以推测出大多数（有可能不是全部）gasdermins使用类似的机制在膜中形成孔(Fig. 4)。全长GSDMA3的高分辨率晶体结构和最近鉴定的从重组脂质体(liposomes)中提取的GSDMA3 NT孔洞的低温电镜结构提供了一个很好的模板，让人们详细地来理解gasdermin形成的孔洞。在GSDMA3晶体结构中，GSDMA3 NT结构呈现出延伸的扭曲β片状结构，两侧有几个螺旋，这代表了一种在已知的成孔蛋白中未发现的新的球状折叠。

从β片的一端延伸出一个长环，连接到紧密并列在GSDM NT结构域一侧的螺旋GSDM CT结构域。在GSDM NT结构域内，螺旋α1和位于β片层结构凹侧的短β发夹与GSDM CT结构域具有密切的相互作用。在β片层的另一端，由两个柔性环持有的短螺旋从球状褶皱中突出，它们与GSDM CT结构域的另一部分相互作用(Fig. 4)。两个结构域间的相互作用将全长GSDMA3锁定到自抑制状态。当自抑制被破坏时，GSDM CT结构域从凹面释放，解除对GSDM NT的抑制，而GSDM NT结构域则以形成膜孔。

与自抑制状态相比，孔内的GSDMA4 NT结构显示出构象发生了强烈的变化，这种变化主要存在于两个结构元件中。在自抑制结构中与GSDM CT结构域接触的短螺旋连同其侧翼环一起重新折叠成两条β链，形成一个β发夹。在邻近区域，另一条β链及其侧翼环也重新折叠成β链，并形成一个额外的β发夹。四个新形成的β链各自与核心β片层结构中的现有β链合并。两个反向平等的β发夹形成了一个长的四链双亲性β片层结构，这个结构从gasdermin N末端结构域的核心球状折叠延伸中伸出去。

这些构象变化会产生三个寡聚界面，驱动GSDM NT结构域形成环状孔结构，并伴随着两亲性β片层结构绑定在一起，组装成一个插入膜中的β桶结构。β桶或寡聚界面上的突变，例如GSDMD中的E15K和L192D突变会严重操作GSDM NT结构域的成孔活性。每个GSDMA3孔包含大约26-28个GSDM NT原聚体(promoters)，不过多数都是27倍原聚体(Fig. 4)。GSDMA3孔的内径约为18 nm，外径约为28nm。GSDMD孔的并不均一，它们的内径在10到20 nm之间，这表明由不同GSDM NT结构域形成的孔中存在着环境依赖的化学计量。GSDMD孔的尺寸足够大，这样在经典炎症小体中活化过程中生成的成熟的IL-1β就能流出来。

导致孔形成的构象变化是由GSDM NT结构域与膜磷脂的结合触发的。在GSDMA3孔的低温电子显微镜结构中观察到了可能的脂质结合位点。在螺旋α1和插入膜的β片层之间的根部和带正电荷的口袋中电子密集，心磷脂的头部基团中也有类似的电子密度，也被用于重建GSDMA3孔。该口袋的表面被自抑制GSDMA3结构中的GSDM CT域完全掩蔽。由于结构域间切割不能解锁自抑制相互作用(GSDMD和GSDME)，因此仍有待研究带电的磷脂头如何能够进入完全掩埋的口袋从而触发随后的构象变化。或者说，在GSDM NT结构域中可能存在另一个脂质结合位点，该位点具有与膜中磷脂完全结合的启动作用。在GSDMA3孔中观察到的结构变化使人想起膜攻击复合体穿孔素样/胆固醇依赖性细胞溶解素(MACPF/CDC)家族中所观察到的结构变化，尽管两种类型的孔形成蛋白的总体结构非常不同。对于MACPF/CDC家族来说，单体成孔结构域在随后的构象变化之前寡聚成可溶性的前孔结构(pre-pore)，从而介导膜插入和成熟孔的形成。利用高分辨率原子力显微镜，最近的一项研究分析了GSDMD孔隙形成的动态过程，研究发现GSDMD NT单体嵌入脂细胞膜并组装成弧形或狭缝状中间寡聚物，然后生长为环状跨膜孔(Fig. 4)。该过程是连续的，并且不涉及前孔结构的过渡阶段，这不同于MACPF/CDC家族的特征。

Fig.4：gasdermin膜的插入和孔形成机制

Fig. 4 gasdermin膜插入和孔形成的机制。N-末端gasdermin(GSDM NT)和C-末端gasdermin(GSDM CT)结构域之间的相互作用使蛋白质保持在自抑制状态。在GSDMNT结构域内，螺旋α1和位于β片层结构凹侧的短β-发夹与GSDM CT结构域产生紧密的相互作用。此外，从β片层的一端延伸出一个长环路与GSDM CT连接。一旦自身抑制被破坏，通过哺乳动物caspase 1，小鼠caspase 11或人caspase 4/5和GSDME切割gasdermin D(GSDMD)的机制，GSDM CT从凹面释放，从而释放GSDM NT，后者用于膜插入和孔形成。膜靶向需要具有带负电荷头部基团的磷脂，就如在细胞膜的内叶上发现的那样。与自抑制状态相比，GSDM NT的孔构象显示出剧烈的构象变化，这涉及到新的β链与扭曲的β片层结构合并的重新折叠。这些变化还会产生新的寡聚界面，驱动跨膜的β桶，即GSDM NT孔的组装。从结构角度来看，gasdermin可以在没有结构域间切割的情况下通过破坏分子内抑制的机制来活化。

gasdermin孔形成后的结果

在细胞膜上形成的gasdermin孔通常会导致由于孔诱导的膜裂解而导致的细胞坏死。考虑到不同的信号和细胞环境，不同细胞或不同细胞群体的gasdermin孔的数量和大小可能不同，这就导致了不同的细胞死亡依赖或细胞死亡无关的细胞效应。

亚溶解孔形成

虽然在实验体系中，gasdermin孔最终会导致焦亡，但是细胞裂解并非总是这些孔的主要功能。到目前为止，这种“非依赖裂解”的功能仅在GSDMD孔方面提到，不过其它gasdermin家族的成员也有可能有类似的机制。例如，研究表明，GSDMD孔的形成可能由于细胞类型，GSDMD表达水平，活化和时间以及对抗机制的效率而有所不同。例如在小鼠巨噬细胞中，用细菌肽聚糖的N-乙酰氨基葡萄糖片段或OxPAPC能诱导NLRP3炎性小体的活化，即在没有出现细胞裂解的情况下，还是会出现GSDMD依赖的活细胞中IL-1β的释放。类似的，在炎性小体活化时，LPS能诱导人单核细胞释放IL-1，以及诱导中性粒细胞非裂解形式(lysis-independent)的IL-1释放。这些研究表明不仅表明了会出现亚裂解形式的GSDMD孔的形成，并且亚裂解形式的GSDMD孔也许是一种直接的，非常规形式的释放的IL-1β和IL-18的释放途径，从而可以使在细胞不发生裂解时能释放这两种细胞因子。

gasdermin孔的这种非溶解依赖的功能或许也有可能释放其它蛋白质或参与信号通路的调控。GSDMD相对较大尺寸的孔可以直接释放IL-1β和IL-18（这些细胞因子的分子直径约为5nm）以及其它小的胞质蛋白，例如小GTP酶，半乳糖凝集素或半胱氨酸内肽酶抑制剂cystatins。此外，由于gasdermin孔是大的、非选择性的膜通道，由gasdermin孔引起的离子流出甚至在细胞死亡之前就可以对细胞信号通路产生深远的影响。例如，由GSDMD孔引起的钾流出在LPS诱导的非典型炎症途径活化后触发NLRP3的活化。由于钾驱动的NLRP3活化是一种细胞内在机制，因此钾外流和NLRP3活化必定在细胞经历GSDMD驱动的焦亡之前发生。类似地，有人提出，在嗜肺军团菌以亚裂解形式活化AIM2炎性小体后，GSDMD依赖的钾流出活化NLRP3，并且钾流出损害IFN-I应答，而与最终细胞的死亡无关。

那么，细胞是如何调节gasdermin的活化或孔形成的水平，以及膜上孔形成的过程是可逆的么？干扰素调节转录因子IRF2能强烈上调GSDMD的表达，这可能使处于休眠状态的细胞保持低水平的GSDMD，从而避免细胞死亡。此外，不同细胞类型和活化触发器之间的caspase的活性有很大差异，这也可能是导致亚裂解型GSDMD活化的条件。此外，对成孔毒素或机械或激光诱导的膜损伤的研究表明，细胞膜损伤不是一个最终事件，并确定了几种类型的膜修复机制，它们可以在几秒钟或几分钟内恢复膜的完整性。参与修复机制的一些蛋白被招募到细胞膜上，从而对Ca2+通过GSDMD孔的内流产生应答，并促进含有受损膜的囊泡的出芽以及释放。在中性粒细胞或亚裂解型炎性小体活化过程中，ESCRT或其它的膜修复系统是否活跃还需要进一步的证明。

需要注意的是，ESCRT或其它膜出芽机制释放的囊泡可能代表了非常规的蛋白质分泌的替代途径。在早期关于炎性小体活化的研究发现，在炎性小体活化后不久，外泌体脱落增加，并且可以在炎性小体活化的细胞外泌体中发现成熟的IL-1β。最近的研究还发现，类似的外泌体形成与外泌体介导的IL-1β的释放依赖于GSDMD。这些囊泡是否能以非特异性的试释放胞质蛋白或IL-1β和其它蛋白质是否优先包装到这些囊泡中还有待进一步的研究。总体而言，根据GSDMD的活化水平，gasdermin孔作为非常规蛋白分泌的主要调节器，通过直接膜转位(孔功能)、囊泡释放(诱导膜修复)或通过膜裂解的被动释放(焦亡)来促进无导(leaderless)蛋白的释放。

焦亡型细胞死亡

在大多数情况下，GSDMD加工水平和GSDMD孔形成的增加最终将克服调节机制并诱导细胞死亡。这种由炎症caspase的活化控制的特定类型的导致细胞坏死的开启类型最初被称为焦亡；然而，由于所有的GSDM NT结构域都可以在没有caspase活化的情况下诱导焦亡，因此术语“焦亡”已经被重新定义为一种以gasdermin依赖型的细胞死亡形式(Box 2)。在细胞培养中，焦亡细胞的特征是广泛的胞膜空泡化，随后膜膨胀，最终膜完整性丧失，这可能是由于渗透溶解导致的。导致细胞裂解的精确事件尚不完全清楚，可能gasdermin在细胞器膜，例如线粒体或核膜上形成了孔，这一事件有助于诱导细胞死亡以及导致焦亡细胞发生形态学上的改变。例如，GSDMD NT与心脏脂质结合，并被发现它线粒体膜为靶点，并且能增强活性氧的生成。GSDMA3 NT也能够损伤线粒体并诱导线粒体自噬。在中性粒细胞中，GSDMD已被证明能与核膜结合并破坏它，从而在NETosis期间促进DNA的分泌。因此，即使在高效的细胞膜修复的条件下，影响细胞内细胞器的gasdermin孔仍可能导致细胞死亡。

体内的焦亡型细胞死亡证据源于对于自身炎症冷卟啉相关周期性综合征(CAPS)患者的研究，研究发现，在患者的炎症爆发时，能在患者的血液中系统地检测到炎症寡聚体，这表明潜在的细胞焦亡伴随着活化的炎症小体。从冷卟啉相关周期性综合征患者体内分离单核细胞，使用炎症小体激动剂来刺激该细胞，就能发现一小部分细胞表现出炎症小体标记物的阳性结果，这说明，并非所有的单核细胞都会出现焦亡。由GSDMD驱动的焦亡也被证明能够影响冷卟啉相关周期性综合征的小鼠模型的疾病严重程度。至于焦亡在体内参与疾病以及宿主防御的功能方面，还需要更多的研究。

焦亡通常指的是一种炎症形式的细胞死亡，发生焦亡的细胞释放出大量的分子，这些分子能发生“找我”的信号，并且呈递“吃我“的信号，例如磷脂酰丝氨酸。因此，可以推断，通过胞吐(efferocytosis)作用可以有效地除去焦亡细胞的残留物。然而，焦亡形式的细胞死亡与坏死形式的细胞死亡不同，在多数情况下，焦亡是由炎症caspases引发的，因此焦亡的过程伴随着IL-1β和其它IL-1家族成员的释放。这些细胞因子的释放有可能是导致GSDMD依赖性的焦亡表型，并使其与坏死区分开来。但是，现在不能排除焦亡的细胞也能释放其它独特的危险信号分子，这些信号分子是在细胞器破裂或caspase活性的过程产生的。一些研究已经鉴定了与焦亡相关的分泌组(secretome)，鉴定了在caspase 1活化和细胞膜通透人时时释放的900多个蛋白质。然而，目前还没有系统的研究将焦亡的分泌组和免疫学结果与其他类型的细胞死亡进行比较。总之，现在可以假设焦亡能导致不同的免疫学结果，导致低水平或高水平的炎症反应(Fig. 5)，并且这种炎症应答将取决于活化的gasdermin的类型，活化的机制(蛋白酶驱动或其他驱动)，执行焦亡的细胞环境和细胞类型，以及gasdermin损伤介导的膜修复效率。

Fig.5：gasdermin孔形成以及焦亡的免疫学结果

Fig. 5 gasdermin孔隙形成和焦亡的免疫学结果。焦亡细胞释放大量细胞内分子，这些分子可以通过充当警报器和“找到我”信号来活化免疫系统。
A. 如果gasdermin孔被修改并终止了gasdermin活化信号，那么细胞内容物的释放可能是短暂的，并且仅限于能够穿过gsdermin孔的小分子(损伤相关的分子模式(DAMP))。
B. 在病原体或激活NF-κB的损伤相关信号存在的情况下，核苷酸结合寡聚结构域、富含亮氨酸重复和含Pyrin结构域的3(NLRP3)炎症小体就会活化焦亡，并与caspase 1的活化和促炎细胞因子(IL-1β，IL-18)的释放以及穿过gasdermin毛孔的小细胞内蛋白(DAMP)有关。在这种情况下，如果没有修复细胞膜上的gasdermin孔，则随着促炎细胞因子释放的激增，以及大量细胞内成分(如炎症小体寡聚体)的释放，最终会以导致促炎性焦亡。因此，焦亡最有可能导致不同的免疫学结果，导致低水平或高水平的炎症反应。

GSDMD, gasdermin D; PAMP, pathogen- associated molecular pattern.

Box 2：焦亡：一种由gasdermin诱导的坏死型细胞死亡

焦亡这个术语最初是根据其形态学特征和对caspase 1的严格而特异性的要求而定义的，它被定义为“caspase 1依赖性坏死”。焦亡这个词的英文是pyroptosis，它最初来源于希腊语，其中pyro的意思是“火”或者“发热”，ptosis的意思是“落”（注：这里的“落”可以理解为树叶落下这个动作，跟凋亡的英语apoptosis中的后缀ptosis一样），这两个部分构成了pyroptosis，意思是强调这种细胞死亡形式的特征是伴随着促炎发生以及成熟的IL-1β和IL-18的释放。随后在2002年鉴定了炎症小体复合物，以及在2011年发现了非典型炎症小体通路，焦亡又被重新定义为炎症小体依赖性细胞死亡以及炎性小体的效应机制。然而，现在的研究也明确发现，焦亡并不一定伴随着成熟的IL-1β和IL-18的释放，因为即使在没有caspase 1的情况下，小鼠caspase 11和人caspase 4诱导的具有所有特征的形态学特征也会诱发焦亡。这一点通过鉴定GSDMD作为焦亡的唯一执行者而进一步得确认。并且实验表明，GSDMD或其他gasdermins的N-末端结构域的表达足以诱导焦亡，而不需要caspase的活化。

更多的报道已经开始从炎症caspase和炎性小体中研究焦亡和gasdermin活化。例如，研究表明，中性粒细胞弹性蛋白酶和caspase 8可以切割和激活GSDMD以导致细胞死亡，并且caspase 3加工GSDME也可以导致焦亡型细胞死亡。由于在所有情况下，细胞死亡在形态上类似于焦亡，并且依赖于gasdermin家族成员的活化，因此似乎很明显，焦亡这个术语需要重新定义。因此，我们建议将“焦亡”定义为“gasdermin诱导的坏死型细胞死亡”，并将此术语应用于所有可通过诱导膜通透化而导致细胞死亡的gasdermin家族成员。我们建议可以将这个术语独立地使用于涉及了gasdermin激活的实际机制和它发生的细胞类型中(We also suggest to use this term independently of the actual mechanism of gasdermin activation and the cell type in which it occurs.)。我们承认上游信号事件或受影响的细胞类型可以改变焦亡刺激炎症或免疫反应的能力（由caspase 1引起的焦亡涉及成熟的IL-1β的释放，而GSDME依赖的焦亡可能不会），这与活化机制有关，与作为细胞死亡执行者的gasdermins的功能无关。

gasdermins在疾病方面的作用

在不同的病理生理条件下，gasdermin家族发挥着重要的，以及各种功能的焦亡作用，因此gasdermin家族可能在机体的健康和疾病中发挥了重要作用。事实上，事实上，第一个gasdermin基因就是因为其在小鼠皮肤病中的作用才被克隆出来的。随着2015年首次发生gasdermin在焦亡和炎症小体活化中的作用，现在已经涌现了很多研究gasdermin活化与各种疾病关系的文献。

gasdermin A

GSDMA的在胃肠道和皮肤中高表达，但在原发性胃癌和胃癌细胞系中不表达。早期的研究表明，恢复人胃癌细胞系中的GSDMA表达可以增加其对转化生长因子-β(TGF-β)诱导的凋亡的敏感性。然而，GSDMA在癌细胞中的促死亡活性是否需要切割以及GSDMA在健康胃上皮细胞中的生理功能尚不清楚。如上所述，GSDMA3 NT在培养的细胞中的表达后，它会通过细胞膜通透化(permeabilization)来诱导细胞死亡。有趣的是，在这一过程中，还存在LC3-II的上调，这是自噬途径的标记，这表明GSDMA3诱导的自噬可能与自噬成分的活化同时发生，但其潜在机制以及GSDMA3 NT诱导自噬的生理相关性仍有待确定。GSDMA3 NT诱导的细胞死亡可以通过共表达GSDMA3 CT而被抑制，这一现象类似于该结构域在其他gasdermins中发现的抑制功能。Gsdma3的不同突变体与皮肤相关表型相关，其中就包括角质形成和脱发。然而，Gsdma3-/-小鼠没有表现出明显的皮肤发育异常，这表明这些突变赋予了突变体蛋白新的功能。皮肤中Gsdma3的生理功能似乎与毛发的发育有关。Gsdma3中的功能获得突变(gain-of-function)被发现能破坏GSDMA3的自身抑制，从而使GSDMA3 NT结构域诱导焦亡。在过表达GSDMA3 NT的细胞中能观察到线粒体活性的下降。然而，目前尚不清楚GSDMA3 NT和/或GSDMA3的功能获得突变是否可以直接靶向作用于线粒体膜并形成孔，从而导致线粒体衰竭并促进线粒体自噬和焦亡。

gasdermin B

现已发现GSDMB与肿瘤进展有关，其在胃癌、宫颈癌和乳腺癌以及肝癌中的表达上升。在人类表皮生长因子受体2(HER2；也称为ERBB2)阳性的乳腺癌患者中，肿瘤细胞中GSDMB基因表达增加与预后不良、生存率降低和转移增加有关，也与HER2靶向治疗的不良治疗反应有关；同时，GSDMB被发现与ERBB2有共表达。GSDMB在人类中有几种剪接体，GSDMB异构体2(isoform)与乳腺癌细胞中的促瘤性(pro-tumorigenic)和促转移性表型(pro-metastatic)密切相关。目前尚不清楚GSDMB如何促进癌细胞存活的，因为GSDMB NT在培养的细胞中过度表达时能够诱导细胞焦亡。

不同的全基因组关联研究也揭示了某些GSDMB单核苷酸多态性与疾病易感性增加之间存在着相关性，如哮喘、克罗恩病和溃疡性结肠炎(ulcerative colitis)。全长GSDMB和GSDMB NT结构域都可以与磷脂酰肌醇和巯基糖脂硫脂结合。事实上，GSDMB也被认为在硫脂的细胞运输方面发挥作用，研究表明GSDMB CT内的单核苷酸多态性可以导致GSDMB结构的改变，从而影响细胞内的硫脂水平。他研究发现，GSDMB的剪接异构体(splice variant)与哮喘风险较低有关。这个异构体缺少编码caspase 1裂解位点的外显子，因此影响GSDMB诱导气道上皮细胞的焦亡能力。要充分了解GSDMB在癌症、感染和自身免疫性疾病中的作用还需要进一步的深入研究。

gasdermin C

GSDMC最初被发现在转移性黑色素瘤细胞中高表达，因此最初被命名为黑素瘤衍生亮氨酸拉链额外核因子(MLZE)。最近有报道称，GSDMC敲除后能降低结直肠癌细胞系的增殖，这为GSDMC的促肿瘤作用提供了进一步的证据。然而，另一项研究表明，GSDMC在许多食管鳞状细胞癌中的表达受到抑制，这表明它可能是肿瘤抑制基因。因此，到目前为止还没有关于GSDMC是否促进或抑制癌症发展，以及这一功能是否需要其N-末端成孔结构域活化的研究报道。

gasdermin D

体外研究表明，GSDMD是非典型炎症途径下游细胞死亡的唯一执行者。然而，值得注意的是，由于存在ASC焦点刺激的备用细胞死亡程序或更混杂有caspase 1的蛋白水解活性，因此GSDMD的缺失并不能完全消除经典炎症小体活化后的细胞死亡。在活体实验方面，研究表明Gsdmd-/-小鼠与Casp11-/-小鼠一样，对LPS诱导的致死性具有类似的保护作用。此外，Casp11-/-和Gsdmd-/-小鼠肠沙门氏菌亚种(Salmonella enterica subsp)或肠血清鼠伤寒沙门氏菌(TyphimuriumΔSifa)和流产布鲁氏菌(Brucella Abortus)的感染均表现出易感性增强，这两类细菌都能活化非经典炎症小体通路。相比之下，在感染新凶手弗朗西丝氏菌(Francisella novicida)方面，Gsdmd-/-小鼠比CASP1-/-或AIM2-/-小鼠抵抗力更强。同样，据报道，感染鼠伤寒沙门氏菌(S. Typhimurium)的gsdmd-/-小鼠的腹腔IL-1β水平高于casp1-/-casp11-/-小鼠。然而，Gsdmd缺乏症完全保护了携带有Pyrin家族性地中海热(familial Mediterranean fever)相关的Mefv V276A等基因的小鼠免受自身炎症疾病的发展，并且在表达有NLRP3的D301N功能获得突变的小鼠中完全消除了所有新生儿期发病多系统炎性疾病(neonatal onse multisystem inflammatory disease,NOMID)相关的炎症症状。因此，GSDMD可能在体内发挥环境依赖性作用，并且在某些情况下，Gsdmd-/-动物可能参与替代的细胞死亡途径，从而允许保护动物免受典型炎症小体活化后的微生物攻击。

gasdermin E

GSDME最初报道为与人类非综合征性听力损失相关的基因。所有已知的GSDME突变跳过了外显子8，从而产生具有细胞毒性的截短蛋白。由于GSDME NT具有重组后活性和跳过外显子8导致GSDME CT阻遏结构域的破坏，因此可以假设GSDME相关的听力损失是由细胞死亡引起的。然而，由于GSDME在多个组织中表达，目前尚不清楚为什么内耳细胞在GSDME自动活化时优先经历细胞死亡。一种合理的解释是，由于移码突变，从而导致了编码与听力有关的某个蛋白的尾部出现了41个残基的非自然序列，使得正常的蛋白被截短，而这个截短的蛋白对蛋白酶体的降低非常敏感；因此，该蛋白的剩余低水平仅足以损害听觉系统。或者说，GSDME的不同表达水平和不同细胞类型中截短的突变体的降解途径存在差异，这都有可能导致一些细胞对耳聋突变更具抵抗力。

GSDME的生理活化机制可能涉及caspase 3的切割。激活caspase 3的化疗药物，如拓扑替康(topotecan)、依托泊苷(etoposide)、顺铂(cisplatin)和CPT-11，随后被证明可诱导那些高表达GSDME的细胞系的焦亡型死亡，而这些药物则诱导GSDME阴性细胞的凋亡。因此，表达低水平GSDME的小鼠骨髓来源的巨噬细胞在诱导凋亡时不会经历GSDME依赖的焦亡，即使这一过程存在着GSDME的加工。然而，GSDME可能是化疗治疗中一些副作用出现的原因，因为与野生型小鼠相比，Gsdme-/-小鼠更能抵抗顺铂注射引起的组织损伤和体重减轻。因此，GSDME有可能以细胞类型特异性的方式将细胞凋亡转化为细胞焦亡型死亡。

Pejvakin

PJVK突变也与人类和小鼠的听力损伤有关。与导致常染色体显性听力损失的GSDME突变不同，PJVK的所有已知突变都会造成常染色体隐性听力损伤，这种受损是由于功能障碍的外毛细胞(outer hair cells)引发的。有趣的是，研究表明Pjvk-/-小鼠表现出早发性进行性听力损失的听觉表型，类似于携带PJVK突变的患者。这与所有其他gasdermin家族成员不同，在gasdermin家族中，敲除小鼠没有显示出肉眼可见的表型，这表明PJVK的突变实际上导致了功能丧失型表型(loss-of-function)。虽然尚未证明PJVK NT结构域在孔形成方面的功能，但其生理功能可能需要PJVK即使在没有刺激的情况下也持续活跃，即形成孔。事实上，PJVK CT结构域要短得多，并且与其他gasdermin家族成员没有明显的同源性(参考Supplementary Figure 1)，因此PJVK CT结构域可能无法充当抑制性结构域。有趣的是，有人提出PJVK定位于内毛细胞中的过氧化物酶体的膜上，并且它能直接招募自噬蛋白LC3B来驱动在噪声过度暴露时引起的氧化应激之后，由自噬介导的受损过氧化物酶体(Pexophagy)的清除。PJVK驱动的pexophagy随后就是过氧化物酶体的增殖，从而保护听觉毛细胞免受氧化损伤。事实上，Pjvk-/-小鼠表现出过氧化物酶体功能障碍和抗氧化防御受损的迹象，并且Pjvk-/-毛细胞中的过氧化物酶体在听觉诱发后就显示出结构异常。有人认为，这有助于Pjvk-/-小鼠和携带有PJVK突变的人对声音的异常易感性，从而导致进行性听力损失。

抑制gasdermins与治疗

鉴于GSDMD在炎症小体诱导的细胞死亡和细胞因子释放中的关键作用，针对gasdermin孔的抑制正在成为抗炎治疗的新靶点。抑制焦亡作用的思路来自于对这种类型的细胞死亡的初次研究，研究表明细胞的裂解可以被渗透保护剂(osmoprotectant)或高浓度的甘氨酸(在毫摩尔范围内)阻断。然而，这一策略并不阻止分子直接通过GSDMD孔隙(如摄取染料或IL-1β释放)，因此不适合于靶向GSDMD的相关炎症。第一种被发现能够阻止焦亡和IL-1β释放的化合物是石榴多酚(Penicalagin)，这是一种在石榴中发现的复杂的抗氧化剂多酚。石榴多酚能可逆地抑制细胞膜通透性和caspase 1活化后成熟的IL-1β和IL-18的释放，在低微摩尔范围内具有半最大抑制浓度(IC50)，而不影响GSDMD NT的产生。石榴多酚不影响NLRP3或AIM2炎性小体活化，但能阻断细胞膜的流动性，并可能干扰GSDMD NT正确地插入细胞膜，阻止其寡聚化和/或孔的形成。最近，一份报告表明，石榴多酚影响NlrP1和NLRC4炎性小体的活化；然而，由于石榴多酚以干扰细胞膜流动性和外源颗粒和蛋白质的摄取，因此在刺激前应用时，它也可能干扰细菌性炎症小体活化剂(如鞭毛蛋白或炭疽致死因子)向细胞的传递。因此，需要进一步的研究来了解石榴多酚的作用机制及其对GSDMD的潜在特异性。体外实验表明，经石榴多酚清洗后，可以观察到快速的细胞溶解，这种溶解可以被镧系元素(La3+和Gd3+)所阻断。这些金属化学元素也被发现能抑制巨噬细胞中在焦亡前细胞膜孔的形成。然而，镧系元素是否影响细胞膜上的GSDMD NT寡聚化尚不清楚，因为它们不能阻断IL-1β释放。

最近报道了几种化合物能直接靶向GSDMD。例如，necrosulfonamide(NSA)是一种先前已知的半胱氨酸反应性药物，此化合物能够抑制坏死性凋亡(necroptosis)中人混合系列蛋白激酶样结构域(MLKL,mixed-lineage kinase domain-like pseudo- kinase)，它也能抑制人源和鼠源细胞的焦亡。NSA与GSDMD结合并抑制细胞膜上的GSDMD NT的寡聚化，但不影响Toll样受体信号传导、炎症小体的活化、GSDMD的裂解或细胞因子的成熟。在低微摩尔范围内，NSA能阻断炎症小体活化时的细胞膜通透性和IL-1β释放，并保护小鼠免受内毒素诱导的小鼠感染。从机制上讲，NSA能共价修饰GSDMD的Cys191，这是孔形成所必需的残基。有趣的是，NSA不影响GSDME NT诱导的细胞死亡，这与GSDME中同样的位点缺少一个半胱氨酸的效应是一致的。另一种化合物是LDC7559，据报道该化合物可抑制在小鼠或人类细胞中发生的弹性蛋白酶依赖性的NETosis或焦亡，但其作用机制至今尚不清楚。

结论与展望

自从GSDMD被鉴定为焦亡的执行者以来，我们对这个新兴的细胞死亡执行者家族的了解取得了快速的进展。对全长gasdermins的结构和诱变以及GSDMA3 NT孔隙的结构研究已经揭示了自抑制和膜插入的机制。此外，大量研究表明，在炎症小体诱导的焦亡之外，gasdermins还参与细胞死亡，例如在出现典型的凋亡形态后发生后的NETosis和坏死死亡过程中。另一方面，对gasdermin孔隙的非裂解功能的鉴定扩展了gasdermins的已知功能，提示了gasdermin在非常规蛋白质分泌中的关键作用。

很明显，这些研究只是冰山的一角，在未来几年中，gasdermin家族将成为免疫、癌症治疗和其他领域的核心参与者。然而，许多问题仍未得到解答。例如，了解不同gasdermins的活化机制，哪些细胞类型能产生这些活性gasdermins，以及它们将引发什么生物学效应，这都是非常重要的问题。目前尚不清楚不同的GSDMB和GSDMC是否具有不同的作用(GSDMB和GSDMC的功能在很大程度上是未知的)，以及是否所有的GSDMB都具有炎症和宿主防御微生物感染方面的功能。

随着知识的进步，gasdermin活性的药理调节有可能成为治疗具有炎症成分的不同疾病的关键。特异性gasdermin阻滞剂的表征以及验证焦亡是一个可行的药理学靶点才刚刚开始。在许多体内模型和临床试验中，已经用化学抑制剂或阻断抗体证实了阻断炎性小体-IL-1途径在临床方面的意义，针对IL-1的阻断抗体已经被批准用于治疗自身炎症和慢性炎症性疾病。新型gasdermin阻滞剂的开发不仅对理解焦亡孔在不同疾病中的作用具有重要意义，而且还将为开发新的炎症性疾病治疗方案铺平道路。

发表于 | 分类于
| 字数统计: 1,739 | 阅读时长 ≈ 6

前言

在学习统计学的时候，我遇到过这么一个问题，也就是说，样本的标准差公式，如下所示：

很多统计学书上都提到，在样本标准差的计算公式中，平方根中的分子是 $n-1$，而总体标准差则是 $n$ 。其理由是为了校正样本变异性而做出的调整，这是对总体标准差的无偏估计。

但是，为什么说这是一种无偏估计，很多书中并没有提及，或者说是只用了很粗略的语言简单地说了一下，其实也没必过于纠结这个问题，记住就行。但是，如果实在是想弄明白这个问题，网上也有人给出了证明过程，但是证明过程对于没有数学基础的人来讲，还是有点难的，这个完整的证明过程的可以参考知乎上的这个帖子《为什么样本方差(sample variance)的分子是n-1》。

最近我看到了一本统计学的书《行为科学统计》（第七版）作者：[美]FrederickJ Gravetter，这本书中对这个问题的描述很清楚，通过用举例子的方式说明了一下（并非严格证明），为什么在样本标准差中，使用 $n-1$ 是对总体方差的无偏估计。

另外说明一下，《行为科学统计》这本书原本就是给社会学的学生学习统计学准备的，里面的语言浅显易懂，没有复杂的公式，对于数学功底差的学生来说，非常友好，最新一版已经到了第9版。

背景知识

  1. 离差：数据到平均数的距离，例如对于一个 $\mu = 50$ 的分布来说，如果你的一个数据是 $X=53$ ，那么离差就是 $X-\mu = 53 - 50 = 3$。如果数据是45，那么离差就是 $45-50=-5$。
  2. 离均差平方和（SS，sum of squares of deviation from mean）：由于离差有正有负，最终所有离均差的和即 $(X-\mu)$ 为0，因此离均差的和无法描述一组数据的变异大小。因此将离均差平方后相加得到平方和$Var(X)=E(X-\mu)^2$，这就是离均差平方和（sum of squares of deviations from mean, SS）。
  3. 方差：方差定义为离均差平方和的平均数，如下所示：

  1. 标准差：方差的平方根。

计算过程

先来看一组数据，即1, 9, 5, 8, 7，我们把这个数据当作是总体，现在我们计算它的离差，离差的平方等，如下所示：

1
2
3
4
5
6
a <- c(1,9, 5, 8, 7) # 原始数据
a - mean(a) # 离差
(a- mean(a))^2 # 离差的平方
sum((a- mean(a))^2) # 离均差平方和
sum((a- mean(a))^2)/length(a) # 方差
sqrt(sum((a- mean(a))^2)/length(a)) #标准差

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
> a <- c(1,9, 5, 8, 7) # 原始数据
> a - mean(a) # 离差
[1] -5  3 -1  2  1
> (a- mean(a))^2 # 离差的平方
[1] 25  9  1  4  1
> sum((a- mean(a))^2) # 离差的平方和
[1] 40
> sum((a- mean(a))^2)/length(a) # 方差
[1] 8
> sqrt(sum((a- mean(a))^2)/length(a)) #标准差
[1] 2.828427

总体方差与样本方差

总体方差的计算公式就是方差的定义：

总体标准差就是总体方差的平方根，如下所示：

样本方差与总体方差有所不同，为了校正样本变异性，我们需要对方差和标准有效期的公式做出调整，此时样本方差公式中的分母是 $n-1$， 如下所示：

样本标准差的公式如下所示：

这里要注意的是，公式使用了 $n-1$ 来代替 $n$ ，这是为了校正样本变异性的偏误做出的调整，调整的结果使所得的结果变大，从而使样本方差成为对总体方差精确的无偏估计（如果是n，则是有偏估计）。

下面我们用不太严谨的一个案例来说明一下为什么样本方差中的分母是 $n-1$ 。

举例说明为什么是n-1

现在我们设计一个N=6的总体，它的元素为0, 0, 3, 3, 9, 9，现在我们计算可知它的总体均数为 $\mu=4$， 方差 $\sigma^2 = 14$。

现在我们从这个总体中选择一个 $n=2$ 的样本，我们选出所有可能的组合，并计算出其平均数，有偏方差，无偏方差，如下所示：

样本编号 第1个数 第2个数 平均数 有偏的方差(n) 无偏的方差(n-1)
1 0 0 0 0 0
2 0 3 1.5 2.25 4.5
3 0 9 4.5 20.25 40.5
4 3 0 1.5 2.25 4.5
5 3 3 3 . 0
6 3 9 6 9 18
7 9 0 4.5 20.25 40.5
8 9 3 6 9 18
9 9 9 9 0 0
总和 36 63 126

现在我们观察平均数这一列，原始的总体均数为 $\mu = 4$。虽然没有一个样本的均数恰好为4，但是如果考虑整组样本，将会发现，9个样本的平均数总和为36，因此样本均数数的平均数为 36/9=4，此时样本平均数恰好等于总体平均数。根据定义，这是一个无偏的统计量，也就是说，样本精确地代表了总体。

现在我们考虑用除以n得到的存在偏误的样本方差这一列。原始的总体方差是 $\sigma^2 = 14$。 然而，9个样本方差的总和为63， 这使得63/9=7。注意，这些样本方差的平均值不等于总体方差，也就是说，如果用除以n得到的样本方差，得出的结果不能精确估计总体方差， 也就是说，这些样本方差低估了总体方差，因此是存在偏误的统计量。

现在我们再考虑除了n-1得到的样本详这一列，虽然总体方差为$\sigma^2=14$，然而没有一个样本的方差恰好等于14。但是，如果考虑整组样本方差，将会发现这9个值总和为126，因此方差的平均值为126/9=14。因此，样本方差的平均值恰好等于总体方差。也就是说，样本方差（此时是使用了n-1来代替n）是对总体方差的一个精确的、无偏的估计。

结论就是，样本平均数和样本方差（使用n-1）都是无偏估计的例子。这个事实使样本平均数和样本方差在推论统计方面变得非常重要。虽然没有单个样本恰好具有与总体一样的平均数和方差，但是，样本平均数和样本方差的平均值确实提供了对相应总体参数的精确估计。

参考资料

  1. 行为科学统计 作者: [美] F. J. Gravetter / [美] L. B. Wallnau，出版社: 中国轻工业出版社，原作名: Statistics for the behavioral sciences，译者: 王爱民 / 李悦，出版年: 2008-7

发表于 | 分类于
| 字数统计: 410 | 阅读时长 ≈ 2

前言

这是一棵树（网名）同学在群里分享的教程，有关R语言多线程的内容。

关于多线程是什么，可以搜一下，网上很多资料。

R的一个比较大的缺点就是它只能单线程操作，例如如果你的电脑有200个处理器，R一次只能使用一个，这个时候工作效率会降低。为了解决这个问题，可以调用R的多线程。

R中与多线程有关的包是parallel。

parallel包

paralle包是R中的基础包，不用额外安装，直接调用即可，它的使用如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Load parallel package
library(parallel)
# Check how many threads and cores in your computer
detectCores()
detectCores(logical = F)
# Open multi-threads
# You will see 4 R programs run in your taks managers simultaneously
cl <- makeCluster(3)
# Close cluster
# Run following command will close 4 R program
stopCluster(cl)

detectCores()这个命令用来查看你的电脑有几个核心，如果添加上logical=F参数，则表示查看电脑有几个物理核心，如果是logical=T则表示有几个逻辑核心。

makecluster()函数用于创建几个并行运算。

我们先来看一下运行这个函数之前的任务管理器，如下所示：

只有一个R，现在我们使用并行计算，即运行这个cl <- makeCluster(3)命令，再看一下任务管理器，如下所示：

现在再关闭这个并行计算，即stopCluster(cl)，如下所示：

使用并行计算的思路就是：打开并行计算，计算，关闭并行计算。

并行计算的分配

当我们使用makeCluster(4)打开了4个R后，如何分配计算任务呢？

发表于 | 分类于
| 字数统计: 3,838 | 阅读时长 ≈ 18

前言

这个教程是一棵树（网名）演示的爬虫笔记，他的Github为：https://github.com/yikeshu0611

爬取科学网，网址为：http://fund.sciencenet.cn/search/smallSubject?subject=H0101&yearStart=2018&yearEnd=2018&submit=list

内容如下所示：

爬虫的思路就是：

  1. 读取网页；
  2. 提取数据。

R包

使用rvest包中的read_html()函数提取网页中的内容。

读取国自然操作

1. 读取网页

  1. 安装并加载rvest包；
  2. 将网址赋值给url；
  3. 使用read_html()函数读取，如下所示：
1
2
3
4
install.packages("rvest")
library(rvest)
url='http://fund.sciencenet.cn/search/smallSubject?subject=H0101&yearStart=2018&yearEnd=2018&submit=list'
read_html(url)

结果如下所示：

1
2
3
4
5
> read_html(url)
{html_document}
<html xmlns="http://www.w3.org/1999/xhtml">
[1] <head>\n<meta http-equiv="Content-Type" content="text/html; charset=UTF-8" ...
[2] <body>\n        <div id="hd">\n    <div class="wp">\n                <ifra ...

上面的结果就是网页结构，其中<head>是网页的头部，<body>是网页的主体。

2. 读取数据

读取数据则是要定位从哪里开始读取，还是先回到网页部分，如下所示：

把鼠标放到上面的题目上，然后单击右键，选择审查元素（chrome浏览器），如下所示：

在上面网址那一行单击右键，复制->Xpath，如下所示：

复制的内容是//*[@id="resultLst"]/div[1]/p/a。现在复制另外一个题目的xpath，内容为//*[@id="resultLst"]/div[2]/p/a。

从这两个内容上我们可以大概知道，id="resultLst"对应了<div id="resultLst"--------</div>，如下所示：

在上面的html代码中，我们哦可以发现，这一部分有2个div，第一个题目是div[1]，第2个题目是div[2]，依次类推，最后一个就是div[last()]（不要过于纠结这种写法）。

现在我们看第1个div，div下面是p节点，如下所示：

p节点下面又有2个节点，b和a，b节点那里是1，就是项目前面的标号，如下所示：

a节点下面是href="..."，如下所示：

我们可以看到，在a节点现在有2个内容，第1个是链接，第2个是文本，也就是标题，我们的目标就是这个项目标题，现在我们从div那个节点开始，来写这个标题的地址，这个网址的结果如下所示：

在rvest包中，网页的定位是使用html_nodes()函数，现在我们定位第1个标题的位置，现在将读取的网页赋值给content，来定位网页中的某个东西，例如标题1，如下所示：

1
2
content <- read_html(url)
html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[1]/p/a')

结果如下所示：

1
2
3
4
5
> content <- read_html(url)
> 
> html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[1]/p/a')
{xml_nodeset (1)}
[1] <a href="http://fund.sciencenet.cn/project/50 ...

标题的xpath地址赋值给xpath，上面的结果就是相应的内容，里面就是一个文本，我们使用html_text()函数来提取这些内容，并将定位的内容赋值给location，然后再提取，如下所示：

1
2
location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[1]/p/a')
html_text(location)

结果如下所示：

1
2
3
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[1]/p/a')
> html_text(location)
[1] "CFTR/EGFR反馈环路调控肺液清除功能在支气管肺发育不良发病中的作用和分子机制"

现在我们提取第2个标题,如下所示：

1
2
location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[2]/p/a')
html_text(location)

运行结果如下所示：

1
2
3
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[2]/p/a')
> html_text(location)
[1] "II型肺泡上皮细胞（AT2）在重症流感肺泡损伤修复过程中的参与作用及调控机制"

其中改变的就是div[2]这个参数。如果是要提取最后1个题目，就是div[last()]，如下所示：

1
2
3
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]/div[last()]/p/a')
> html_text(location)
[1] "II型肺泡上皮细胞（AT2）在重症流感肺泡损伤修复过程中的参与作用及调控机制"

如果是100个题目，不能这么干，我们需要把div[1]这类字符串删掉，如下所示：

1
2
3
4
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]//p/a')
> html_text(location)
[1] "CFTR/EGFR反馈环路调控肺液清除功能在支气管肺发育不良发病中的作用和分子机制"
[2] "II型肺泡上皮细胞（AT2）在重症流感肺泡损伤修复过程中的参与作用及调控机制"

现在就提取了所有的题目。

第二个任务：提取姓名。

现在我们再来提取作者的姓名，例如赵冬莹，前面的xpath操作和前面的相同，即//*[@id="resultLst"]/div[1]/div/p[1]/span[1]/i，现在再来复制第2个名字，即//*[@id="resultLst"]/div[2]/div/p[1]/span[1]/i，复制2个名字主要是为了找到规律。现在我们把div[1]删掉，如下所示：

1
2
location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]//div/p[1]/span[1]/i')
html_text(location)

运行结果如下所示：

1
2
3
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]//div/p[1]/span[1]/i')
> html_text(location)
[1] "赵冬莹" "李辉"

第三个任务：提取标题部分的网址，这个网址，就是标题后面链接的网址，有时候，我们需要爬取二级页面，就地需要获得二级页面的网址，我们看到这个网址不是文本，它虽然和标题在同一个位置，都是a节点下面，但是我们使用html_text()命令并没有将其提取出现，因为这个函数认为它不是文本，而是链接，对应的是herf="----------------"这种格式，如下所示：

现在我们要提取某一个具体的网页(html)属性(attribute)内容，此时我们使用html_attr()命令，例如我们要提取超链接，就写成html_attr("href")，所以，如果我们要提取标题处的链接，就需要先定位到标题那里，然后使用html_attr()函数，如下所示：

1
2
location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]//p/a')
html_attr(location, name = 'href')

结果如下所示：

1
2
3
4
> location <- html_nodes(x = content, xpath = '//*[@id="resultLst"]//p/a')
> html_attr(location, name = 'href')
[1] "http://fund.sciencenet.cn/project/509194"
[2] "http://fund.sciencenet.cn/project/509195"

读取Pubmed

现在来讲一下大致思路：第一，找到网址；第二，定位，也就是说从哪个地方开始抓取数据；第三步，抓取数据。

  1. 打开pubmed，https://www.ncbi.nlm.nih.gov/pubmed，输入circulation，点击搜索如下所示：

  1. 加载rvest包，输入目标网址，如下所示：
1
2
3
## Crawl pubmed
## Setup 1
page_content <- read_html(x = 'https://www.ncbi.nlm.nih.gov/pubmed/?term=circulation')
  1. 像前面一样，右键xpath，如下所示：

其中，一个rprt对应的就是左侧的蓝色阴影部分，一共有20个这样的结构（其实就是一页中的20个结果），如下所示：

我们再回到第1个题目，右键，copy xpath，如下所示：

1
2
3
4
5
6
7
8
## Crawl pubmed
## Setup 1
page_content <- read_html(x = 'https://www.ncbi.nlm.nih.gov/pubmed/?term=circulation')
## Step 2: Crawl content
#xpath
node = '//*[@id="maincontent"]/div/div[5]/div[1]/div[2]/p/a'
html_nodes(x = page_content, xpath = node)

结果如下所示：

1
2
3
> html_nodes(x = page_content, xpath = node)
{xml_nodeset (1)}
[1] <a href="/pubmed/31603578" ref="ordinalpos=1&amp;ncbi_uid=31603578&amp;link_ ...

读取上面内容中的文本部分，如下所示：

1
2
3
4
5
6
7
8
## Step 2: Crawl content
#xpath
### 2.1 Location
node = '//*[@id="maincontent"]/div/div[5]/div[1]/div[2]/p/a'
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_text(x = nodes_content, trim = T)

现在我们提取第2个题目的xpath，再将它与第1个比较，如下所示：

1
2
//*[@id="maincontent"]/div/div[5]/div[1]/div[2]/p/a
//*[@id="maincontent"]/div/div[5]/div[2]/div[2]/p/a

我们可以找到它们的共同规律，即//*[@id="maincontent"]/div/div[5]//div[2]/p/a，如下所示：

1
2
3
4
5
6
7
8
### 2.1 Location
node = '//*[@id="maincontent"]/div/div[5]//div[2]/p/a'
# 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_text(x = nodes_content, trim = T)

运行结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
> ### 2.1 Location
> node = '//*[@id="maincontent"]/div/div[5]//div[2]/p/a'
> nodes_content <- html_nodes(x = page_content, xpath = node)
> 
> ### 2.2 Crawl data
> html_text(x = nodes_content, trim = T)
 [1] "Variation of organ position in snakes."                                                                                                                                   
 [2] "Similar articles"                                                                                                                                                         
 [3] "Pathogenesis of Cluster Headache: From Episodic to Chronic Form, the Role of Neurotransmitters and Neuromodulators."                                                      
 [4] "Similar articles"                                                                                                                                                         
 [5] "Haemodynamic evaluation of the new pulsatile-flow generation method in vitro."                                                                                            
 [6] "Similar articles"                                                                                                                                                         
 [7] "Stacking as a Key Property for Creating Nanoparticles with Tunable Shape: The Case of Squalenoyl-Doxorubicin."                                                            
 [8] "Similar articles"                                                                                                                                                         
 [9] "[Study on absorbed components of Aconitum kusnezoffii under Yunnan Baiyao compatibility in effect of activating blood circulation and removing blood stasis]."            
[10] "Similar articles"                                                                                                                                                         
[11] "[Application of quality constant method in grade evaluation of Aurantii Fructus Immaturus pieces]."                                                                       
[12] "Similar articles"                                                                                                                                                         
[13] "[Study on regularity of Tibetan medicine in treatment of gZav-Grib disease (apoplexy sequelae) based on HIS clinical medical records]."                                   
[14] "Similar articles"                                                                                                                                                         
[15] "Exercise responses in children and adults with a Fontan circulation at simulated altitude."                                                                               
[16] "Similar articles"                                                                                                                                                         
[17] "The Course of Immune Stimulation by Photodynamic Therapy: Bridging fundamentals of photochemically-induced Immunogenic Cell Death to the Enrichment of T Cell Repertoire."
[18] "Similar articles"                                                                                                                                                         
[19] "Design of Dense Brush Conformation Bearing Gold Nanoparticles as Theranostic Agent for Cancer."                                                                           
[20] "Similar articles"                                                                                                                                                         
[21] "Liquid biopsy: one cell at a time."                                                                                                                                       
[22] "Similar articles"                                                                                                                                                         
[23] "The Evolving Role of Neutrophils in Liver Transplant Ischemia-Reperfusion Injury."                                                                                        
[24] "Similar articles"                                                                                                                                                         
[25] "Predictors of Balloon Guide Catheter Assistance Success in Stent-retrieval Thrombectomy for an Anterior Circulation Acute Ischemic Stroke."                               
[26] "Similar articles"                                                                                                                                                         
[27] "High mobility group box protein 1 neutralization therapy in ovine bacteremia: Lessons learned from an ovine septic shock model incorporating intensive care support."     
[28] "Similar articles"                                                                                                                                                         
[29] "Low human and murine Mcl-1 expression leads to a pro-apoptotic plaque phenotype enriched in giant-cells."                                                                 
[30] "Similar articles"                                                                                                                                                         
[31] "Upper-tropospheric bridging of wintertime surface climate variability in the Euro-Atlantic region and northern Asia."                                                     
[32] "Similar articles"                                                                                                                                                         
[33] "Discordances between predicted and actual risk in obese patients with suspected cardiac ischaemia."                                                                       
[34] "Similar articles"                                                                                                                                                         
[35] "Tumor-Associated Release of Prostatic Cells into the Blood after Transrectal Ultrasound-Guided Biopsy in Patients with Histologically Confirmed Prostate Cancer."         
[36] "Similar articles"                                                                                                                                                         
[37] "Circulating tumor cells exhibit metastatic tropism and reveal brain metastasis drivers."                                                                                  
[38] "Similar articles"                                                                                                                                                         
[39] "Permanent pacing post-Fontan is not associated with reduced long-term survival."                                                                                          
[40] "Similar articles"

一共有40个结果，里面有相同的字符串，即"Similar articles"，这个不是我们需要的，这因为定位出现了问题，也就是说node = '//*[@id="maincontent"]/div/div[5]//div[2]/p/a'这段代码有问题，现在我们查看原题目与Similar articles的元素，如下所示：

其中，红框是我们要爬取的题目，而蓝框则similar articles的内容，因此我们需要把蓝框的内容给剔掉，只爬取到class="title"这个字段就行，也就是说添加上p[@class="title"]，如下所示：

1
2
3
4
5
6
7
8
### 2.1 Location
node = '//*[@id="maincontent"]/div/div[5]//div[2]/p[@class="title"]/a'
# 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_text(x = nodes_content, trim = T)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
> ### 2.1 Location
> node = '//*[@id="maincontent"]/div/div[5]//div[2]/p[@class="title"]/a'
> # 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
> 
> nodes_content <- html_nodes(x = page_content, xpath = node)
> 
> ### 2.2 Crawl data
> html_text(x = nodes_content, trim = T)
 [1] "Variation of organ position in snakes."                                                                                                                                   
 [2] "Pathogenesis of Cluster Headache: From Episodic to Chronic Form, the Role of Neurotransmitters and Neuromodulators."                                                      
 [3] "Haemodynamic evaluation of the new pulsatile-flow generation method in vitro."                                                                                            
 [4] "Stacking as a Key Property for Creating Nanoparticles with Tunable Shape: The Case of Squalenoyl-Doxorubicin."                                                            
 [5] "[Study on absorbed components of Aconitum kusnezoffii under Yunnan Baiyao compatibility in effect of activating blood circulation and removing blood stasis]."            
 [6] "[Application of quality constant method in grade evaluation of Aurantii Fructus Immaturus pieces]."                                                                       
 [7] "[Study on regularity of Tibetan medicine in treatment of gZav-Grib disease (apoplexy sequelae) based on HIS clinical medical records]."                                   
 [8] "Exercise responses in children and adults with a Fontan circulation at simulated altitude."                                                                               
 [9] "The Course of Immune Stimulation by Photodynamic Therapy: Bridging fundamentals of photochemically-induced Immunogenic Cell Death to the Enrichment of T Cell Repertoire."
[10] "Design of Dense Brush Conformation Bearing Gold Nanoparticles as Theranostic Agent for Cancer."                                                                           
[11] "Liquid biopsy: one cell at a time."                                                                                                                                       
[12] "The Evolving Role of Neutrophils in Liver Transplant Ischemia-Reperfusion Injury."                                                                                        
[13] "Predictors of Balloon Guide Catheter Assistance Success in Stent-retrieval Thrombectomy for an Anterior Circulation Acute Ischemic Stroke."                               
[14] "High mobility group box protein 1 neutralization therapy in ovine bacteremia: Lessons learned from an ovine septic shock model incorporating intensive care support."     
[15] "Low human and murine Mcl-1 expression leads to a pro-apoptotic plaque phenotype enriched in giant-cells."                                                                 
[16] "Upper-tropospheric bridging of wintertime surface climate variability in the Euro-Atlantic region and northern Asia."                                                     
[17] "Discordances between predicted and actual risk in obese patients with suspected cardiac ischaemia."                                                                       
[18] "Tumor-Associated Release of Prostatic Cells into the Blood after Transrectal Ultrasound-Guided Biopsy in Patients with Histologically Confirmed Prostate Cancer."         
[19] "Circulating tumor cells exhibit metastatic tropism and reveal brain metastasis drivers."                                                                                  
[20] "Permanent pacing post-Fontan is not associated with reduced long-term survival."

此时，如果想读取链接，则如下所示：

1
2
3
4
5
6
7
8
### 2.1 Location
node = '//*[@id="maincontent"]/div/div[5]//div[2]/p/a'
# 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_attr(x = nodes_content, name = "href")

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
> html_attr(x = nodes_content, name = "href")
 [1] "/pubmed/31603578"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31603578"
 [3] "/pubmed/31603552"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31603552"
 [5] "/pubmed/31603372"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31603372"
 [7] "/pubmed/31603305"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31603305"
 [9] "/pubmed/31602894"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602894"
[11] "/pubmed/31602882"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602882"
[13] "/pubmed/31602864"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602864"
[15] "/pubmed/31602790"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602790"
[17] "/pubmed/31602649"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602649"
[19] "/pubmed/31602558"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602558"
[21] "/pubmed/31602399"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602399"
[23] "/pubmed/31602356"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602356"
[25] "/pubmed/31602354"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602354"
[27] "/pubmed/31602200"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31602200"
[29] "/pubmed/31601924"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601924"
[31] "/pubmed/31601851"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601851"
[33] "/pubmed/31601728"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601728"
[35] "/pubmed/31601564"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601564"
[37] "/pubmed/31601552"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601552"
[39] "/pubmed/31601284"                                 "/pubmed?linkname=pubmed_pubmed&from_uid=31601284"

此时读取了40项内容，又多了一倍，还是跟前面的问题一样，此时需要过滤一下（核心内容就是正则表达式），如下所示：

1
2
3
4
5
6
7
8
9
10
11
### 2.1 Location
node = '//div[@class="rprt"]/div[@class="rslt"]/p[@class="title"]/a'
# 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_attr(x = nodes_content, name = "href")

结果如下所示：

1
2
3
4
5
> html_attr(x = nodes_content, name = "href")
 [1] "/pubmed/31603578" "/pubmed/31603552" "/pubmed/31603372" "/pubmed/31603305" "/pubmed/31602894" "/pubmed/31602882"
 [7] "/pubmed/31602864" "/pubmed/31602790" "/pubmed/31602649" "/pubmed/31602558" "/pubmed/31602399" "/pubmed/31602356"
[13] "/pubmed/31602354" "/pubmed/31602200" "/pubmed/31601924" "/pubmed/31601851" "/pubmed/31601728" "/pubmed/31601564"
[19] "/pubmed/31601552" "/pubmed/31601284"

其实我们可以发现，node = '//div[@class="rprt"]/div[@class="rslt"]/p[@class="title"]/a'这一句中的最后一部分node = '//p[@class="title"]/a'其实是唯一标记的，也就是说在这个路径中没有重复，因此我们还可以改一下代码，把这个字符前面的都删掉，如下所示：

1
2
3
4
5
6
### 2.1 Location
node = '//p[@class="title"]/a'
# 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
nodes_content <- html_nodes(x = page_content, xpath = node)
### 2.2 Crawl data
html_attr(x = nodes_content, name = "href")

运行一下，如下所示：

1
2
3
4
5
6
7
8
9
10
> ### 2.1 Location
> node = '//p[@class="title"]/a'
> # 这里要注意//与/的区别，/表示绝对路径，//表示相对路径
> nodes_content <- html_nodes(x = page_content, xpath = node)
> ### 2.2 Crawl data
> html_attr(x = nodes_content, name = "href")
 [1] "/pubmed/31603578" "/pubmed/31603552" "/pubmed/31603372" "/pubmed/31603305" "/pubmed/31602894" "/pubmed/31602882"
 [7] "/pubmed/31602864" "/pubmed/31602790" "/pubmed/31602649" "/pubmed/31602558" "/pubmed/31602399" "/pubmed/31602356"
[13] "/pubmed/31602354" "/pubmed/31602200" "/pubmed/31601924" "/pubmed/31601851" "/pubmed/31601728" "/pubmed/31601564"
[19] "/pubmed/31601552" "/pubmed/31601284"

结果是一样的，这里要学习的就是唯一标记符，使用这种方法非常高效（核心就是找到唯一的节点）。

简化操作之管道

上面介绍的这个爬虫过程，都需要找到网址，输入节点，比较麻烦，因此可以采取管道（%>%）来简化操作，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
## pipeline operation
page_content <- read_html(x = 'https://www.ncbi.nlm.nih.gov/pubmed/?term=circulation')
# title
node <-  '//p[@class="title"]/a'
page_content %>%
  html_nodes(xpath = node) %>%
  html_text(trim = T)
# Link
page_content %>%
  html_nodes(xpath = node) %>%
  html_attr(name = 'href')

简化操作之函数

上面的操作还能继续简化，也就是写成一个函数，R中构建函数的内容可以找本书看看，这里直接放代码，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
## Contruction function for web site
read_html.my <- function(url){
  library(rvest)
  page_concent <<- read_html(x = url)
  # note symbol "<<-" is transform local variable into global variable
}
read_html.my('https://www.ncbi.nlm.nih.gov/pubmed/?term=circulation')
# Read text from web site
html_text.my <- function(xpath){
  page_concent %>%
    html_nodes(xpath = xpath) %>%
    html_text(trim = T)
}
# Read link from web site
html_href.my <- function(xpath){
  page_concent %>%
    html_nodes(xpath = xpath) %>%
    html_attr(name = 'href')
}
html_text.my('//p[@class="title"]/a')
html_href.my('//p[@class="title"]/a')

总结

涉及到的知识点大概如下所示：

  1. 网页的构成（xpath，html，css，绝对路径(/)与相对路径(//，节点，内容）；
  2. 正则表达式；
  3. R中函数的构建（局部变量，变局变量，局部变量切换为全局变量<<-）；
  4. 管道操作(%>%)。

发表于 | 分类于
| 字数统计: 5,618 | 阅读时长 ≈ 20

前言

最近在处理qPCR数据过程中，设计的qPCR阵列中采用了多个内参，我对于多内参的处理不太理解。后来找到了一篇文献，就是讲qPCR中的多内参问题的。文献信息如下：

1
Vandesompele, J., et al. (2002). "Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes." Genome Biology 3(7): research0034.0031.

此篇文献引用量很高（到2019年10月为止，引用量约为16000），虽然文献比较古老，但是了解一下文献中的思想还是比较有用的，以下是文献的阅读笔记（基本上就是翻译原文外加自己检索）。

文献中涉及一定的线性代数知识，我不太懂，只能先暂时先记录下来，等看懂的时候再补充一篇笔记。

研究背景

​ 在生物学研究中，基因表达的分析有着重要的作用。对于基因表达模式的理解可以为复杂的调控网络研究提供思路。最近开发的两种检测转录本丰度的方法得到了非常广泛的应用（该文章是2002年写的，这里的最近已经是非常久远了）。微阵列的方法（微阵列其实就是芯片法）可以同时检测数以千计的基因表达情况，而RT-PCR的方法则是用于检测许多不同的样本中数量有限基因的表达情况（在细胞数量比较少的情况下用这种方法更好，说服力更强）。与传统的方法例如northern-blot，核糖核酸酶保护，或者是竞争RT-PCR（competitive rt-PCR）相比，这两种方法（微阵列与RT-qPCR）都有快速，高通量的优点。但是这两种方法都需要与传统方法一样，需要对目的mRNA进行均一化（normalization）。

以下是补充一些与核酸保护实验相关的背景知识。

核酸酶保护实验

​ 糖核酸酶保护实验（Ribonuclease protection assay，RPA）是一种mRNA定量分析方法。其基本原理是将标记的特异RNA探针（32P或生物素）与待测的RNA样品液相杂交，标记的特异RNA探针按碱基互补的原则与目的基因特异性结合，形成双链RNA；未结合的单链RNA经RNA酶A或RNA酶T1消化形成寡核糖核酸，而待测目的基因与特异RNA探针结合后形成双链RNA，免受RNA酶的消化，故该方法命名为RNA酶保护实验。

​ 对于32P标记的探针，杂交双链进行变性聚丙酰胺凝胶电泳，用放射自显影或磷屏成像系统检测被保护的探针的信号; 对于生物素标记的探针，杂交双链经过变性聚丙酰胺凝胶电泳后电转移至尼龙膜，采用链霉亲和素－辣根过氧化物酶（Streptavidin-HRP）和化学发光底物与膜上biotin标记的探针结合，X射线胶片或化学发光图像分析仪检测杂交信号。

​ 在利用RT-PCR对基因进行定量时，为了缩小误差，我们需要对一些因素进行控制，例如原始样品的数量（细胞或者是组织），酶的效率，组织或细胞的整个转录活动差异等因素。在控制条件的理想情况下，提取的高质量RNA的转录本数量的可重复结果与细胞数量是呈正比的，但是精确的细胞数目很难进行定量（尤其是当原始材料是组织时更是如此）。

​ 另外一种常用的定量手段就是RNA的量（这里指的是数量，原文写的是mass quantity），尤其是northern blot分析方法。但是这种方法也存在着一些争议。因为这种方法没有考虑到RNA的质量（这里侧重的是RNA的完整性，原文写的是quality）与相关酶的反应效率。此外，在某些情况下，这些因素也无法定量，例如从提取的微切割的组织（microdissected tissue）中得到的RNA数量极少，这个时候就无法对RNA进行定量。

​ 其中最大的争议还在于，用于均一化时所考虑的总RNA量，因为总的RNA主要是由rRNA组成的，它并不代表mRNA。这里大概介绍一下使用总RNA量来作为RNA质控的方法。

​ 当我们提取了RNA后，在进行质控时，通常是使用1%的琼脂糖凝胶电泳来跑200ng到400ng的总RNA，然后观察RNA的条带。良好的RNA（这里指的是动物组织或细胞的样本）中可以看到3个条带，从上到下依次是28S、18S、5S（这是核糖体RNA的主要成分），其中28S条带与18S条带非常清楚，5S条件非常弱，28S和18S比值可以判断RNA样本的完整性，通常这个比值是1.5（肉眼不太好判断，通常来说就是28S比18S亮就行了），如果这两个条带变弱，5S条件变亮，则说明RNA降解，如下所示：

​ 不过电泳法是比较粗糙的方法，毕竟这篇文献是2002年的，现在如果要进行测序建库，要求就比较高，会采用Qubit检测RNA的浓度，用安捷伦2100来检测RNA的完整性，具体的检测指标可以检索相关资料。

但是在总的RNA中，mRNA的含量仅占7.5%。此外，还有报告指出，rRNA的转录还受到生物因素与药物的影响。利用18S或28S作为RNA定量的标准缺陷还在于，在纯的mRNA分子中，并不含有18S或28S，与靶mRNA转录本相比，18S或28S的含量太高。在利用RT-PCR进行分析时，很难排除基线(baseline)。

现在补充一些竞争性PCR的背景知识。

竞争性PCR

​ 在竞争性PCR中，逆转录之前加入一个外源的RNA转录本（内部标准RNA）作为样品间差异的对照。内部标准RNA被逆转录并同目的模板一起扩增，作为cDNA转化和扩增效率差异的对照。通过将特异性的目的序列同已知浓度的内部标准RNA一起扩增进行定量。通过比较由内部标准获得的信号和目的模板所获得的信号可以确定目的模板的丰度。

​ 到目前为止（也是2018年为止），内参基因是最常用的对mRNA进行均一化的方法。内参基因（internal control）通常指的是管家基因（housekeeping gene），这类基因并不会随着组织或细胞的不同而变化， 或者是不会随着实验条件而变化。但是在许多的研究里，研究者使用的这些持续表达的内参基因中，并没有对这些基因的稳定表达进行验证。文献中使用的这些管家基因偶尔会出现变化极大的现象。随着RT-PCR灵敏度，可靠性以及使用范围的增加，对于合适内参基因的需求也日益迫切。些外，通过利用利用公共微阵列数量，这种均一化因素（normalization factor）也在广泛的应用的微阵列定标系数（scaling factor）中得到了验证。

需要了解的几个公式

文献末尾处提到了几个公式，现在汇总如下所示：

公式1:单一内参均一化错误值（Single control normalization error）

公式上面的英文注释如下所示：

E值： 对于任何m个样本，用real-time RT-PCR技术可以测量其n个内参基因的表达水平，其表达水平为aij。对于每2个样本（分别记为p和q）的组合来说，它们的两个内参基因j和k的每一种组合来说，可以计算出它们的单一内参均一化错误值E（single control normalization error，公式1）。

当样本p和q分别针对其内参基因j或k进行均一化的时候，其样本p和样本q的表达总倍数差异就是错误值E。 n表示内参的数目（注：既然是多内参，这个n必然是大于等于2）。

aij表示的是内参的表达水平（这是一组数据）。

j与k的范围为1到n，p与q的值为1到m（都是闭区间），并且j不等于k，p不等于q。Rjkpq的倒数就是E值，公式1中，aqj是说在样本q中，以内参基因j为标准进行均一化后的数据，aqk则是说，以内参基因k为标准进行均一化后的数据，aqj/aqk则是这二者之比，apk/apj为例，aqj/aqk与apk/apj的比值则是指p与q样本中，a基因表达水平的误差。

个人理解：这个E值其实就是每个样本，每个内参均一化后误差的乘积。

公式2:内参基因稳定检测值M

原文如下所示：

公式解释：

M值：对于任意两个内参基因j与k，可以计算出由它们构成的一个Ajk数列（array），这个数列包含一组经log2转化来的比值（aij/aik），如果是m个样本，那么Ajk中一共是m个值（公式2）。作者定义Vjk是内参基因j和k的成对变异（pairwise variation）它是Ajk的标准差（公式3）。其中内参基因j的稳定表达值Mj就是所有成对基因变异Vjk的算术平均数（公式4）。

画个图，如下所示：

以上计算的只是内参基因j的稳定值。

研究思路

提取了不同组织或细胞的RNA进行检测，然后检测这些RNA中的几个选定的内参，再通过一定的算法进行比较。

实验过程

（一）管家基因表达谱

​ 本文研究了10个常用管家基因在人类13个不同组织中的表达，这10个管家基因如下所示：

管家基因列表。

注：选这个10个管家基因的思路是要尽量使这10个基因属于不同的功能分类，这样就能极大地了降低它们之间存在的共调控现象。

作者检测了80个组织中的这10个管家基因的表达情况，这80个组织分别为：在34个神经母细胞瘤（来源不同的实验室以及不同的病人，NB1-34）；20个短期培养的常规成纤维细胞（FIB1-20）；13个常规的淋巴细胞样本（LEU1-13）；9个骨髓样本（BM1-9）；9个其他的人类正常组织（分别为heart，brain，fetal brain，lung，trachea，kidney，mammary gland，small intestine和uterus）。

（二）单一内参均一化错误值（Single control normalization error）

Figure 1 -对于两个理想的内参基因来说，E值等于1，这个很好理解，我们还看一下前面的E公式，如下所示：

也就是说，理想的两个内参基因表达水平必然是相同的，因此E等于1。

事实上观察到的E值通常大于1，这就构成了两个样本之间的E倍差异，具体的差异值取决用于均一化的特定管家基因。对所有的10个内参基因两两配对（一共45种组合），以及所有样本的两两配对（一共是865个组合）绘制E值，如figure 1所示。此外，通过分析重复运行的相同内参基因绘制出了系统错误的分位数。E值的75百分位数和90百分位分别是3.0（范围是2.1-3.9）和6.4（范围是3.0-10.9）。

注：所有样本的两两配对是指，同一个组织内的两两配对，因此根据前面的内容：

34个神经母细胞瘤（来源不同的实验室以及不同的病人，NB1-34）；20个短期培养的常规成纤维细胞（FIB1-20）；13个常规的淋巴细胞样本（LEU1-13）；9个骨髓样本（BM1-9）；9个其他的人类正常组织（分别为heart，brain，fetal brain，lung，trachea，kidney，mammary gland，small intestine和uterus）

排列组合数目就是865。

（三）基因表达稳定性检测以及所选基因的秩

​ 在研究中，通常会认为基因的表达水平应该精选一个稳定表达的基因进行均一化。但是为了验证选定的这个内参基因的稳定表达，我们需要一种可靠的检测手段手段来确定这个内参基因确实是稳定表达的，从而剔除掉非特异性的变异。为了解决这个循环问题（circular problem，这里的循环问题就是说，我要选择稳定的内参，就需要“内参”进行均一化，但是选择均一化的这个“内参”又得是稳定的，因此是一个不断循环的问题），作者以非均一化表达水平为基础，提出了一种基因稳定性检测算法。

​ 这种算法的原理是，无论实验条件，无论细胞类型，两种理想的内参基因在所有样本中的表达必然是相同的。而在实际条件下，两种真实内参基因的表达比值（ratio）的变异就会反映出一种（或者是两种）内参基因的表达并非是恒定不变的，这种比值的差异与基因表达的稳定程度呈反比。对于每一个内参基因来说，作者都研究了它与其他所有内参基因两两的差异，这种差异用表达比值的对数转换的标准差来表示，并且定义了一个内参基因表达稳定检测值M作为一个特定基因与其他所有内参基因的平均两两变异。

​ 其中，最小M值的基因是表达最稳定的内参基因。假如选定的这些内参基因不存在共调控，那么分步排除掉的最高M值的基因后，最终剩下的内参基因就会产生一对组合，构成这种组合的两个基因是稳定表达的管家基因，即它们是样品中最稳定的内参基因。

​ 为了计算最佳的内参基因组合，作者利用VBA开发了一个Excel插件，命名为geNorm，它能自动计算所有内参基因的的M值（现在这个插件已经整合到了qbase+中，网址为：https://www.qbaseplus.com/）。这个程序可以剔除那些不好的内参基因（即最高M值的基因），并且重新计算剩余内参基因的M值。利用geNorm，作者计算了5类组织中的10个内参基因的M值，并按从大到小的顺序排列（参见figure2与table 3）。此外，作者还计算了系统变异，其方法对相同的基因进行重复的RT-PCR实验，这种变异反映了仪器自身，加样，酶的固有变异。

Figure 2 横坐标表示的是剩余的内参数目。

注：qbase+这个软件我并不熟悉，这是没有重现数据的分析过程。但是在前一篇笔记中，我使用了NormqPCR这个R包来重现了数据分析过程。

（四）基于多内参基因的几何均数进行均一化因子计算

这一部分就是说明要选择几个内参来进行均一化。

​ 为了精确地评估基因的表达水平，作者推荐使用多个内参基因进行均一化。将多个内参取几何均数。至于使用的内参基因数目，需要在实际情况与精确性方面进行取舍。但是就实际情况来看，选择太多的内参基因浪费（例如10个太浪费了），如果太少，不为精确，作者推荐最少使用3个表达稳定的内参基因作为均一化因子（normalizaiton factor, NFn, n=3）。

​ 作者还研究了一些情况，即那些利用了超过3种内参基因用于均一化的情况，作者计算了相同组织类型内的所有样本，按照两次均一化因子（NFn和NFn+1）的两两变异Vn/n+1（即aij=NFn,i和aik=NFn+1,n代表用于均一化的基因数目，范围是3到9，i表示样本，具体参照公式2和3）。我的理解是，在选定了几个内参，例如5个，对所有的样本进行均一化，计算其变异，即V5，然后再排除一个内参，剩了4个内参，再对所有的样本均一化，计算其变异，即V4，最后计算V4/V5）。这种算法研究的是，如果存在着一个很大的变异，就意味着，这个加入的内参基因（我的理解是，排除的那个内参）有着很强的效应，它应该优先选择用于均一化因子的计算（不太理解这种算法）。

​ 对于所有的组织类型，用于计算的均一化因子是3个最稳定的内参基因，通过逐步添加剩余的7个内参基因。逐步计算两两配对变异，对于连续的NFn和NFn+1均一化因子，这就反应了加入的第n+1基因的效应。
（如figure 3a所示）。从图中可以明显看出来，对于leukicytes, fibroblast和bone marrow来说，添加的第4个基因没有特别明显的影响（就是说低V3/4值），这也说明，在NF3和NF4值之间存在着很强的关联，如figure 3b所示。以这些数据为基础，作者采用了0.15作为阈值，低于0.15的内参基因不需要添加进来。对于neuroblastoma和正常组织来说，分别需要加入1个和2个内参基因（参见figure 3b）。最高的V8/9和V9/10值对于normal pool，neuroblastoma和leukocyte来说就确认了figure 2中逐步排除的最差内参基因。这种分析表明了平均M值的在开始时的强烈下降，这就指出了，leukocyte的两个内参基因的表达异常，neuroblastoma和normal tissues一个内参基因的不稳定。此外，对于后两种组织来说，需要包含另外的内参基因，从而与内参基因的高度变异保持一致。

验证可能的RT-PCR均一化因子

​ 为了验证这种基因检测算法的可靠性，即最低M值是最稳定的表达的内参基因，作者检测了每个内参基因的特异性变化，作为均一化之后表达的变异系数（variation coefficient）。对于合适的内参基因来说，这个系数应该是最小的。作者以三个基因的几何平均数为基础，分别计算了三个不同的均一化因子，最低的M值（NF3（1-3）），最高的M值（NF3（8-10）和中等的M值（NF3（6-8））。接着，作者计算用最低变异系数，在每种组织类型内，计算了这三个基因的平均基因特异性变异（figure 4a）。从计算结果可以看出，当数据用NF3（1-3）进行均一化时，这些基因特异性变异最小。

​ 这就表明，基因稳定性检测能够确定地确定表达最稳定的内参基因。为了确认高M值是一类表达不稳定，或者是表达有差异的基因，作者分析了MYCN的表达水平，MYCN这个基因是一种在neuroblastoma中表达有差异的基因。在neuroblastoma中，MYN的M值是6.02，而B2M（表达最稳定的内参基因）的M值是2.17。进一步的观察发现，用单一的内参基因进行均一化会导致其它内参基因更高的的基因特异性变异，这就进一步说明了利用多个管家基因进行均一化的优势。

​ 为了说明最佳的个内参基因不受细胞增殖的影响，作者分析了PCNA的表达水平（PCNA是增殖的标志），研究了四个最佳管家基因和标志基因PCNA之间的Spearman秩相关系数。从分析结果来看，管家基因之间存在着关联（p值小于0.001，相关系数在0.6到0.76之间）。相比之下，PCNA和三个管家基因之间不存在联系，而与HPRT1则存在着弱相关（p值为0.024，相关系数为0.43）。这些数据清楚地说明了，最稳定的管家基因与细胞的增殖状态无关。

​ 为了进一步研究选定的内参基因几何均数的精确性，作者研究了公共数据库中的微阵列数据，利用这些数据中的内参基因的几何均数进行了均一化，具体的选择过程参见文献，作者选定了5个内参基因，其标准是M值小于0.7，其计算结果与微阵列的计算结果类似。

组织特异性管家基因的表达

​ 为了比较13所测组织中的管家基因表达水平的异质性。作者剔除了差异最大的4个管家基因（这4个管家基因分别为B2M，RPL13A，ACTB与HMBS，参见Table 2）后，计算了6个内参基因的几何均数作为“伪内参”，然后将这10个内参基因再根据这个“伪内参”进行均一化。Figure 5显示了所有的10个基因的不同丰度类型，表达最高的ACTB基因是最低的HMBS的400倍。虽然内参基因的丰度在不同的组织中类似，但还是观察到了组织特异性表达的差异，例如，B2M的表达水平在leukocytes中比Fetal Brain中高出112倍，而Fibroiblast和Heart tisssue中，ACTB的差异有22倍。与这两个基因相比，有2个基因的表达比较稳定，例如UBC和HPRT1。

讨论部分略。

参考资料：

  1. PCR和RT-PCR基础（强烈推荐）
  2. 核糖核酸酶保护实验的基本原理
  3. 大家好，给大家介绍一下，这是RNA质检结果解读指南

发表于 | 分类于
| 字数统计: 8,176 | 阅读时长 ≈ 43

前言

NormqPcR这个包用于对qPCR数据进行均一化(normalization，这里需要注意的是，这里的normalization只是使用内参来校正目标基因的表达的意思，与统计学中的归一化(Normalization)，中心化(center)，标准化(scales)之类的术语没有关系)。这篇笔记主要是翻译了这个包的英文文档，本文档以一个小型案例(vignette)来演示了这个包的主要功能，笔记主要内容如下所示：

  • 合并技术重复(technical replicates)（了解即可）；
  • 处理低Cq值(undetermined values)（了解即可）；
  • 计算用于均一化的最佳管家基因(housekeeping genes)（这个最重要，其余的很多函数其实了解即可，这一部分要了解geNorm和NormFinder的算法原理）；
  • 使用管家基因来对数据进行均一化(normalization)（了解）。

注：在这篇笔记中，管家基因(housekeeping genes)与内参(reference)指的是同一个意思。

合并技术重复

我们先使用ReadqPCR包来导入原始数据，这个原始数据中含有技术重复，其后缀都是_TechRep.n，其中n表示总的重复数。读取数据过程如下所示：

1
2
3
4
5
6
7
8
9
# BiocManager::install("ReadqPCR")
# BiocManager::install("NormqPCR")
# input raw data for qPCR
library(NormqPCR)
library(ReadqPCR)
path <- system.file("exData", package = "ReadqPCR")
qPCR.example.techReps <- file.path(path, "qPCR.techReps.txt")
qPCRBatch.qPCR.techReps <- read.qPCR(qPCR.example.techReps)
rownames(exprs(qPCRBatch.qPCR.techReps))[1:8]

计算结果如下所示：

1
2
3
4
> rownames(exprs(qPCRBatch.qPCR.techReps))[1:8]
[1] "gene_aj_TechReps.1" "gene_aj_TechReps.2" "gene_al_TechReps.1"
[4] "gene_al_TechReps.2" "gene_ax_TechReps.1" "gene_ax_TechReps.2"
[7] "gene_bo_TechReps.1" "gene_bo_TechReps.2"

上面是读取数据的一些信息，从中我们可以发现，里面有技术重复，例如gene_aj_TechReps.1与gene_aj_TechReps.2就是做了2个技术重复。

在我们计算之前，我们还要计算技术重复Cq的算术平均值(arithmetic mean)，这一步我们使用combineTechReps()函数来计算一下，这个函数会生成一个qPCRBatch对象，这样就计算出了技术重复的算术平均值，如下所示：

1
2
3
4
5
6
7
# Cal arithmetic mean of the raw qPCR Cq value
# combineTechReps function will return a new qPCRBatch object
combinedTechReps <- combineTechReps(qPCRBatch.qPCR.techReps)
combinedTechReps
exprs(qPCRBatch.qPCR.techReps)
exprs(combinedTechReps)
apply(exprs(qPCRBatch.qPCR.techReps)[1:2,],2,mean)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
> combinedTechReps <- combineTechReps(qPCRBatch.qPCR.techReps)
> combinedTechReps
qPCRBatch (storageMode: lockedEnvironment)
assayData: 8 features, 4 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: four one three two
  varLabels: sample
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  
> exprs(qPCRBatch.qPCR.techReps)
                       four      one    three      two
gene_aj_TechReps.1 18.02512 20.72960 21.87109 17.51878
gene_aj_TechReps.2 17.59243 19.58196 21.71061 18.40029
gene_al_TechReps.1 20.30233 22.23137 19.30596 21.32223
gene_al_TechReps.2 20.56011 21.66608 20.07169 19.53765
gene_ax_TechReps.1 20.67044 19.09129 21.90034 21.15897
gene_ax_TechReps.2 20.44186 18.22294 21.28660 21.31048
gene_bo_TechReps.1 20.62939 18.77342 19.70848 21.08782
gene_bo_TechReps.2 20.70826 17.95030 19.53260 21.53267
gene_cp_TechReps.1 23.50768 20.39530 22.06632 22.99578
gene_cp_TechReps.2 22.73511 19.64446 20.74793 22.34119
gene_dx_TechReps.1 21.16744 19.57131 21.62546 21.65467
gene_dx_TechReps.2 21.36217 20.21180 21.41619 20.73930
gene_ei_TechReps.1 19.90058 19.78343 20.26527 19.83862
gene_ei_TechReps.2 19.67431 20.97861 19.82184 19.16033
gene_jy_TechReps.1 21.60261 19.23546 18.22053 22.21470
gene_jy_TechReps.2 21.04445 18.75998 19.31650 20.94050
> exprs(combinedTechReps)
            four      one    three      two
gene_aj 17.80877 20.15578 21.79085 17.95954
gene_al 20.43122 21.94873 19.68883 20.42994
gene_ax 20.55615 18.65711 21.59347 21.23473
gene_bo 20.66882 18.36186 19.62054 21.31024
gene_cp 23.12139 20.01988 21.40713 22.66848
gene_dx 21.26481 19.89156 21.52082 21.19699
gene_ei 19.78745 20.38102 20.04355 19.49948
gene_jy 21.32353 18.99772 18.76851 21.57760
> apply(exprs(qPCRBatch.qPCR.techReps)[1:2,],2,mean)
    four      one    three      two 
17.80877 20.15578 21.79085 17.95954

从上面我们可以看到：

  1. 做了2个技术重复；
  2. 使用combineTechReps()函数计算了技术重复的算术平均值；
  3. 又用apply()函数核对了一下，没问题。

处理未检测值

在qPCR实验中，当Cq值超过一定阈值后，qPCR仪有可能会将这个读数标记为Undetermined，也就是数值太小（Cq值越大，表示信号越弱），在qPCRBatch对象中，这个值就会被标记为NA。不同的人有不同的处理策略，常见的处理策略有两种。

第一种就是将那些大于38的Cq值认为是无法检测值，并将其标记为40，这是Taqman软件中的常见做法。而在我们的这个案例中，我们不会将其替换为40，而是替换为NA值，现在我们读入一个Taqman的qPCR数据，这批数据是用96孔板测的，其中，实验组(mia)4个重复，对照组(no-mia)4个重复，如下所示：

1
2
3
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "/example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)

现在我们只看某个样本的检测值，即Ccl20.Rn00570287_m1这个数据，如下所示：

1
exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]

结果如下所示：

1
2
3
4
5
> exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]
    fp1.day3.v     fp2.day3.v   fp5.day3.mia   fp6.day3.mia   fp.3.day.3.v 
            NA             NA       35.74190       34.05922       35.02052 
  fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia 
            NA       35.93689       36.57921

现在我们使用replaceAboveCutOff()函数将大于35的数值替换为NA，如下所示：

1
2
3
# Use replaceAboveCutOff method in order to replace anything above 35 with NA
qPCRBatch.taqman.replaced <- replaceAboveCutOff(qPCRBatch.taqman,newVal = NA, cutOff = 35)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]

结果如下所示：

1
2
3
4
5
> exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
    fp1.day3.v     fp2.day3.v   fp5.day3.mia   fp6.day3.mia   fp.3.day.3.v 
            NA             NA             NA       34.05922             NA 
  fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia 
            NA             NA             NA

但是对于有些同学来说，我不想把大于35的数值替换为NA，而是想把所有的NA都替换为40，那么可以使用replaceNAs()函数，如下所示：

1
2
3
# Use replaesNAs methods to replace all NAs with 40
qPCRBatch.taqman.replaced <- replaceNAs(qPCRBatch.taqman, newNA = 40)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]

结果如下所示：

1
2
3
4
5
> exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
    fp1.day3.v     fp2.day3.v   fp5.day3.mia   fp6.day3.mia   fp.3.day.3.v 
      40.00000       40.00000       35.74190       34.05922       35.02052 
  fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia 
      40.00000       35.93689       36.57921

有的时候我们会遇到这样的一样种情况，在一次实验中，例如我们做了96孔板的qPCR实验（做了好几块96孔板），某个孔的检测器中检测的数值大于设定的阈值，而有的孔的检测器检测到的数值则没有大于设定的阈值。有的同学可能就会认为，即然一些孔的检测器检测到的值设为NAs，那么这个孔的检测器在读取其它数值时，也应该设为NAs。因为否则一个检测器的计数出现了异学，那么会导致某个基因表达水平的计算出现错误估计（这一段不太理解，原话如下）：

In addition, the situation sometimes arises where some readings for a given detector are above a given cycle threshold, but some others are not. The user may decide for example that if a given number of readings are NAs, then all of the readings for this detector should be NAs. This is important because otherwise an unusual reading for one detector might lead to an inaccurate estimate for the expression of a given gene.

对于不同的样本类型这个过程应该分开处理，因为你也许想比较不同样本之间的某个特定的基因。因此有必要使用比较矩阵来指定每个样本类型的重复数。因此，在上面的案例中，我们一共有两个类型，分别是实验组与对照组，每组4个重复，我们可以使用contrastMatrix()函数和sampleMaxMatrix()函数来创建比较矩阵，如下所示：

1
2
3
4
5
6
7
8
9
10
11
# Use contrastMatrix() and sampleMaxMatrix() to generate contrast matrix
sampleNames(qPCRBatch.taqman) # Check group information
a <- c(0,0,1,1,0,0,1,1) # one for each sample type, with 1 representing
b <- c(1,1,0,0,1,1,0,0) # position of sample type in samplenames vector
contM <- cbind(a,b)
colnames(contM) <- c("case","control") # set the names of each sample type
rownames(contM) <- sampleNames(qPCRBatch.taqman) # set row names
contM
sMaxM <- t(as.matrix(c(3,3))) # now make the contrast matrix
colnames(sMaxM) <- c("case","control") # make sure these line up with samples
sMaxM

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> # Use contrastMatrix() and sampleMaxMatrix() to generate contrast matrix
> sampleNames(qPCRBatch.taqman) # Check group information
[1] "fp1.day3.v"     "fp2.day3.v"     "fp5.day3.mia"   "fp6.day3.mia"  
[5] "fp.3.day.3.v"   "fp.4.day.3.v"   "fp.7.day.3.mia" "fp.8.day.3.mia"
> a <- c(0,0,1,1,0,0,1,1) # one for each sample type, with 1 representing
> b <- c(1,1,0,0,1,1,0,0) # position of sample type in samplenames vector
> contM <- cbind(a,b)
> colnames(contM) <- c("case","control") # set the names of each sample type
> rownames(contM) <- sampleNames(qPCRBatch.taqman) # set row names
> contM
               case control
fp1.day3.v        0       1
fp2.day3.v        0       1
fp5.day3.mia      1       0
fp6.day3.mia      1       0
fp.3.day.3.v      0       1
fp.4.day.3.v      0       1
fp.7.day.3.mia    1       0
fp.8.day.3.mia    1       0
> sMaxM <- t(as.matrix(c(3,3))) # now make the contrast matrix
> colnames(sMaxM) <- c("case","control") # make sure these line up with samples
> sMaxM
     case control
[1,]    3       3

有关更多的算法细节可以参考limmma包中的内容，那个包中创建差异基因比较的矩阵也是类似方法。

例如，如果用户决定将那些4个重复中，有3个是NA的值都设为NA的话，可以使用makeAllNewVal()函数，如下所示：

1
qPCRBatch.taqman.replaced <- makeAllNewVal(qPCRBatch.taqman, contM,sMaxM, newVal=NA)

现在我们来看一下Ccl20.Rn00570287 m1这个样本的检测值，对照组的值都为NA，因为原始的4个技术重复中，有3个是NA，有一个是35，因此全设为NA。不过实验组的数据全部进行了保留，如下所示：

1
2
exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]

结果如下所示：

1
2
3
4
5
6
7
8
9
10
> exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]
    fp1.day3.v     fp2.day3.v   fp5.day3.mia   fp6.day3.mia   fp.3.day.3.v 
            NA             NA       35.74190       34.05922       35.02052 
  fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia 
            NA       35.93689       36.57921 
> exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
    fp1.day3.v     fp2.day3.v   fp5.day3.mia   fp6.day3.mia   fp.3.day.3.v 
            NA             NA       35.74190       34.05922             NA 
  fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia 
            NA       35.93689       36.57921

前面这两部分内容（合并技术重复与处理未检测值）了解一下就行，平时也用不到，最重要的是下面的内容。

选择最稳定内参基因

这一部分涉及选择最稳定的内参基因，其中用的是geNorm算法，这个算法的相关文献信息如下所示：

Jo Vandesompele, Katleen De Preter, Filip Pattyn et al. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biology 2002. 3(7):research0034.1-0034.11. http://genomebiology.com/2002/3/7/research/0034/

geNorm

先加载数据，如下所示：

1
2
3
options(width = 68)
data(geNorm)
str(exprs(geNorm.qPCRBatch))

数据结构如下所示：

1
2
3
4
5
> str(exprs(geNorm.qPCRBatch))
 num [1:10, 1:85] 0.0425 0.0576 0.1547 0.1096 0.118 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:10] "ACTB" "B2M" "GAPD" "HMBS" ...
  ..$ : chr [1:85] "BM1" "BM2" "BM3" "BM4" ...

现在我们对选择的内参基因进行排序。我们在这里使用selectHKs()函数来实现geNorm算法中的逐步处理过程，这里的数据都是文献中材料与方法中提到的数据。

第一步，我们计算所有内参基因稳定检测值M后，排除那些M值最高的内参基因；

第二步，再次计算M值，再排除掉最高M值的内参基因，直到剩下两个内参基因（也就是参数minNrHK=2，如下所示：

1
2
tissue <- as.factor(c(rep("BM", 9), rep("FIB", 20), rep("LEU", 13),rep("NB", 34), rep("POOL", 9)))
res.BM <- selectHKs(geNorm.qPCRBatch[,tissue == "BM"], method = "geNorm",Symbols = featureNames(geNorm.qPCRBatch),minNrHK = 2, log = FALSE)

计算过程中，这个函数会给出不断计算的M值，结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
res.BM <- selectHKs(geNorm.qPCRBatch[,tissue == "BM"], method = "geNorm",Symbols = featureNames(geNorm.qPCRBatch),minNrHK = 2, log = FALSE)
###############################################################
Step 1:
stability values M:
    HPRT1     YWHAZ    RPL13A       UBC      GAPD      SDHA       TBP 
0.5160313 0.5314564 0.5335963 0.5700961 0.6064919 0.6201470 0.6397969 
     HMBS       B2M      ACTB 
0.7206013 0.7747634 0.8498739 
average stability M:	0.63628545246682
variable with lowest stability (largest M value):	ACTB
Pairwise variation, (9/10):	0.0764690052563778
###############################################################
Step 2:
stability values M:
    HPRT1    RPL13A     YWHAZ       UBC      GAPD      SDHA       TBP 
0.4705664 0.5141375 0.5271169 0.5554718 0.5575295 0.5738460 0.6042110 
     HMBS       B2M 
0.6759176 0.7671985 
average stability M:	0.582888329316757
variable with lowest stability (largest M value):	B2M
Pairwise variation, (8/9):	0.0776534266912183
###############################################################
Step 3:
stability values M:
    HPRT1    RPL13A      SDHA     YWHAZ       UBC      GAPD       TBP 
0.4391222 0.4733732 0.5243665 0.5253471 0.5403137 0.5560120 0.5622094 
     HMBS 
0.6210820 
average stability M:	0.530228279613623
variable with lowest stability (largest M value):	HMBS
Pairwise variation, (7/8):	0.0671119963410967
###############################################################
Step 4:
stability values M:
    HPRT1    RPL13A     YWHAZ       UBC      SDHA      GAPD       TBP 
0.4389069 0.4696398 0.4879728 0.5043292 0.5178634 0.5245346 0.5563591 
average stability M:	0.499943693933222
variable with lowest stability (largest M value):	TBP
Pairwise variation, (6/7):	0.0681320232188603
###############################################################
Step 5:
stability values M:
    HPRT1    RPL13A       UBC     YWHAZ      GAPD      SDHA 
0.4292808 0.4447874 0.4594181 0.4728920 0.5012107 0.5566762 
average stability M:	0.477377523800525
variable with lowest stability (largest M value):	SDHA
Pairwise variation, (5/6):	0.0806194432580746
###############################################################
Step 6:
stability values M:
      UBC    RPL13A     HPRT1     YWHAZ      GAPD 
0.4195958 0.4204997 0.4219179 0.4424631 0.4841646 
average stability M:	0.437728198765878
variable with lowest stability (largest M value):	GAPD
Pairwise variation, (4/5):	0.0841653121631615
###############################################################
Step 7:
stability values M:
   RPL13A       UBC     YWHAZ     HPRT1 
0.3699163 0.3978736 0.4173706 0.4419220 
average stability M:	0.406770625156432
variable with lowest stability (largest M value):	HPRT1
Pairwise variation, (3/4):	0.097678269387021
###############################################################
Step 8:
stability values M:
      UBC    RPL13A     YWHAZ 
0.3559286 0.3761358 0.3827933 
average stability M:	0.371619241507029
variable with lowest stability (largest M value):	YWHAZ
Pairwise variation, (2/3):	0.113745049966055
###############################################################
Step 9:
stability values M:
   RPL13A       UBC 
0.3492712 0.3492712 
average stability M:	0.349271187472188
> res.BM
$ranking
       1        1        3        4        5        6        7 
"RPL13A"    "UBC"  "YWHAZ"  "HPRT1"   "GAPD"   "SDHA"    "TBP" 
       8        9       10 
  "HMBS"    "B2M"   "ACTB" 
$variation
      9/10        8/9        7/8        6/7        5/6        4/5 
0.07646901 0.07765343 0.06711200 0.06813202 0.08061944 0.08416531 
       3/4        2/3 
0.09767827 0.11374505 
$meanM
       10         9         8         7         6         5 
0.6362855 0.5828883 0.5302283 0.4999437 0.4773775 0.4377282 
        4         3         2 
0.4067706 0.3716192 0.3492712

以上只是BM这类组织的计算结果。现在我们把其它的组织都计算出来，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
res.POOL <- selectHKs(geNorm.qPCRBatch[,tissue == "POOL"],
	method = "geNorm",
	Symbols = featureNames(geNorm.qPCRBatch),
	minNrHK = 2, trace = FALSE, log = FALSE)
res.FIB <- selectHKs(geNorm.qPCRBatch[,tissue == "FIB"],
	method = "geNorm",
	Symbols = featureNames(geNorm.qPCRBatch),
	minNrHK = 2, trace = FALSE, log = FALSE)
res.LEU <- selectHKs(geNorm.qPCRBatch[,tissue == "LEU"],
	method = "geNorm",
	Symbols = featureNames(geNorm.qPCRBatch),
	minNrHK = 2, trace = FALSE, log = FALSE)
res.NB <- selectHKs(geNorm.qPCRBatch[,tissue == "NB"],
	method = "geNorm",
	Symbols = featureNames(geNorm.qPCRBatch),
	minNrHK = 2, trace = FALSE, log = FALSE)

我们就获得了基因的下面排序信息，这个就是文献中的Table 3，文献中的信息如下所示：

我们自己的计算结果：

1
2
3
4
5
6
ranks <- data.frame(c(1, 1:9),
	res.BM$ranking, res.POOL$ranking,
	res.FIB$ranking, res.LEU$ranking,
	res.NB$ranking)
names(ranks) <- c("rank", "BM", "POOL", "FIB", "LEU", "NB")
ranks

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
ranks
   rank     BM   POOL    FIB    LEU     NB
1     1 RPL13A   GAPD   GAPD    UBC   GAPD
2     1    UBC   SDHA  HPRT1  YWHAZ  HPRT1
3     2  YWHAZ   HMBS  YWHAZ    B2M   SDHA
4     3  HPRT1  HPRT1    UBC   GAPD    UBC
5     4   GAPD    TBP   ACTB RPL13A   HMBS
6     5   SDHA    UBC    TBP    TBP  YWHAZ
7     6    TBP RPL13A   SDHA   SDHA    TBP
8     7   HMBS  YWHAZ RPL13A  HPRT1   ACTB
9     8    B2M   ACTB    B2M   HMBS RPL13A
10    9   ACTB    B2M   HMBS   ACTB    B2M

现在绘制出每个细胞类型的平均基因表达稳定值M的折线图，这个对应于文献中的Figure 2，如下所示：

代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(RColorBrewer)
mypalette <- brewer.pal(5, "Set1")
matplot(cbind(res.BM$meanM, res.POOL$meanM, res.FIB$meanM,
res.LEU$meanM, res.NB$meanM), type = "b",
ylab = "Average expression stability M",
xlab = "Number of remaining control genes",
axes = FALSE, pch = 19, col = mypalette,
ylim = c(0.2, 1.22), lty = 1, lwd = 2,
main = "Figure 2 in Vandesompele et al. (2002)")
axis(1, at = 1:9, labels = as.character(10:2))
axis(2, at = seq(0.2, 1.2, by = 0.2), labels = seq(0.2, 1.2, by = 0.2))
box()
abline(h = seq(0.2, 1.2, by = 0.2), lty = 2, lwd = 1, col = "grey")
legend("topright", legend = c("BM", "POOL", "FIB", "LEU", "NB"),
fill = mypalette)

绘制每个细胞类型的成对方差(pairwise variation)图，这个对应文献中的Figure 3，如下所示：

代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
mypalette <- brewer.pal(8, "YlGnBu")
barplot(cbind(res.POOL$variation, res.LEU$variation, res.NB$variation,
res.FIB$variation, res.BM$variation), beside = TRUE,
col = mypalette, space = c(0, 2),
names.arg = c("POOL", "LEU", "NB", "FIB", "BM"),
ylab = "Pairwise variation V",
main = "Figure 3(a) in Vandesompele et al. (2002)")
legend("topright", legend = c("V9/10", "V8/9", "V7/8", "V6/7",
"V5/6", "V4/5", "V3/4", "V2/3"),
fill = mypalette, ncol = 2)
abline(h = seq(0.05, 0.25, by = 0.05), lty = 2, col = "grey")
abline(h = 0.15, lty = 1, col = "black")

注：文献中推荐的成对方差值的截止值为0.15，因此低于此值后，就不需要使用其它的内参基因了。

NormFinder

第二种选择内参的基因的方法是NormFinder，这种方法来源于以下文献：

Andersen, C. L., et al. (2004). “Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets.” Cancer Res 64(15): 5245-5250.

现在我们计算Table 3，如下所示：

第一步，导入数据：：

1
2
3
## NormFinder
data(Colon)
Colon

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
> data(Colon)
> Colon
qPCRBatch (storageMode: lockedEnvironment)
assayData: 13 features, 40 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: I459N 90 ... I-C1056T (40 total)
  varLabels: Sample.no. Classification
  varMetadata: labelDescription
featureData
  featureNames: UBC UBB ... TUBA6 (13 total)
  fvarLabels: Symbol Gene.name
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation: Table 1 in Andersen et al. (2004)

继续计算，如下所示：

1
2
3
Class <- pData(Colon)[,"Classification"]
res.Colon <- stabMeasureRho(Colon, group = Class, log = FALSE)
sort(res.Colon) # see Table 3 in Andersen et al (2004)

结果如下所示：

1
2
3
4
5
> sort(res.Colon) # see Table 3 in Andersen et al (2004)
      UBC      GAPD      TPT1       UBB     TUBA6     RPS13      NACA 
0.1821707 0.2146061 0.2202956 0.2471573 0.2700641 0.2813039 0.2862397 
     CFL1      SUI1      ACTB      CLTC     RPS23  FLJ20030 
0.2870467 0.3139404 0.3235918 0.3692880 0.3784909 0.3935173

再看另外的一个数据集Bladder，如下所示：

1
2
3
4
5
6
7
data(Bladder)
Bladder
grade <- pData(Bladder)[,"Grade"]
res.Bladder <- stabMeasureRho(Bladder,
                              group = grade,
                              log = FALSE)
sort(res.Bladder)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
> data(Bladder)
> Bladder
qPCRBatch (storageMode: lockedEnvironment)
assayData: 14 features, 28 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: 335-6 1131-1 ... 1356-1 (28 total)
  varLabels: Sample.no. Grade
  varMetadata: labelDescription
featureData
  featureNames: ATP5B HSPCB ... FLJ20030 (14 total)
  fvarLabels: Symbol Gene.name
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation: Table 1 in Andersen et al. (2004) 
> grade <- pData(Bladder)[,"Grade"]
> res.Bladder <- stabMeasureRho(Bladder,
+                               group = grade,
+                               log = FALSE)
Warning message:
In .local(x, ...) : Argument 'group' is transformed to a factor vector.
> sort(res.Bladder)
    HSPCB      TEGT     ATP5B       UBC     RPS23     RPS13      CFL1 
0.1539598 0.1966556 0.1987227 0.2033477 0.2139626 0.2147852 0.2666129 
 FLJ20030      TPT1       UBB     FLOT2      GAPD    S100A6      ACTB 
0.2672918 0.2691553 0.2826051 0.2960429 0.3408742 0.3453435 0.3497295

我们也可以使用geNorm算法来计算，如下所示：

1
2
3
4
selectHKs(Colon, 
          log = FALSE, 
          trace = FALSE,
          Symbols = featureNames(Colon))$ranking

计算结果如下所示：

1
2
3
4
5
6
7
8
> selectHKs(Colon, 
+           log = FALSE, 
+           trace = FALSE,
+           Symbols = featureNames(Colon))$ranking
         1          1          3          4          5          6          7 
   "RPS23"     "TPT1"    "RPS13"     "SUI1"      "UBC"     "GAPD"    "TUBA6" 
         8          9         10         11         12         13 
     "UBB"     "NACA"     "CFL1"     "CLTC"     "ACTB" "FLJ20030"

再来看Baldder的数据，如下所示：

1
2
3
4
selectHKs(Bladder,
	log = FALSE,
	trace = FALSE,
	Symbols = featureNames(Bladder))$ranking

结果如下所示：

1
2
3
4
5
6
7
8
> selectHKs(Bladder,
+           log = FALSE,
+           trace = FALSE,
+           Symbols = featureNames(Bladder))$ranking
         1          1          3          4          5          6          7 
    "CFL1"      "UBC"    "ATP5B"    "HSPCB"     "GAPD"     "TEGT"    "RPS23" 
         8          9         10         11         12         13         14 
   "RPS13"     "TPT1" "FLJ20030"    "FLOT2"      "UBB"     "ACTB"   "S100A6"

我们也可以使用NormFinder算法的分步计算方式工来看多个内参基因的情况，也就是文献中附录中的Average control gene部分，selectHKs函数中含有这种算法的参数，如下所示：

1
2
3
4
5
6
7
8
Class <- pData(Colon)[,"Classification"]
selectHKs(Colon, 
	group = Class, 
	log = FALSE, 
	trace = TRUE,
	Symbols = featureNames(Colon), 
	minNrHKs = 12,
	method = "NormFinder")$ranking

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
> Class <- pData(Colon)[,"Classification"]
> selectHKs(Colon, 
+           group = Class, 
+           log = FALSE, 
+           trace = TRUE,
+           Symbols = featureNames(Colon), 
+           minNrHKs = 12,
+           method = "NormFinder")$ranking
###############################################################
Step 1:
stability values rho:
      UBC      GAPD      TPT1       UBB     TUBA6     RPS13      NACA 
0.1821707 0.2146061 0.2202956 0.2471573 0.2700641 0.2813039 0.2862397 
     CFL1      SUI1      ACTB      CLTC     RPS23  FLJ20030 
0.2870467 0.3139404 0.3235918 0.3692880 0.3784909 0.3935173 
variable with highest stability (smallest rho value):	UBC
###############################################################
Step 2:
stability values rho:
     GAPD      TPT1       UBB      NACA      CFL1     RPS13     TUBA6 
0.1375298 0.1424519 0.1578360 0.1657364 0.1729069 0.1837057 0.1849021 
     SUI1      ACTB     RPS23  FLJ20030      CLTC 
0.2065531 0.2131651 0.2188277 0.2359623 0.2447588 
variable with highest stability (smallest rho value):	GAPD
###############################################################
Step 3:
stability values rho:
     TPT1      NACA       UBB     RPS13      CFL1     TUBA6  FLJ20030 
0.1108474 0.1299802 0.1356690 0.1411173 0.1474242 0.1532953 0.1583031 
     SUI1      ACTB     RPS23      CLTC 
0.1586250 0.1682972 0.1686139 0.1926907 
variable with highest stability (smallest rho value):	TPT1
###############################################################
Step 4:
stability values rho:
       UBB      TUBA6       ACTB       CFL1      RPS13       SUI1       CLTC 
0.09656546 0.09674897 0.10753445 0.10830099 0.11801680 0.12612399 0.12773131 
      NACA   FLJ20030      RPS23 
0.13422958 0.14609897 0.16530522 
variable with highest stability (smallest rho value):	UBB
###############################################################
Step 5:
stability values rho:
     RPS13       SUI1      TUBA6       NACA   FLJ20030       CFL1       ACTB 
0.09085973 0.09647829 0.09943424 0.10288912 0.11097074 0.11428399 0.11495336 
     RPS23       CLTC 
0.12635109 0.13286210 
variable with highest stability (smallest rho value):	RPS13
###############################################################
Step 6:
stability values rho:
      ACTB      TUBA6       CFL1   FLJ20030       NACA       CLTC       SUI1 
0.09215478 0.09499893 0.09674032 0.10528784 0.10718604 0.10879846 0.11368091 
     RPS23 
0.13134766 
variable with highest stability (smallest rho value):	ACTB
###############################################################
Step 7:
stability values rho:
      SUI1       NACA   FLJ20030      RPS23      TUBA6       CFL1       CLTC 
0.08281504 0.08444905 0.08922236 0.09072667 0.10559279 0.10993755 0.13142181 
variable with highest stability (smallest rho value):	SUI1
###############################################################
Step 8:
stability values rho:
      NACA       CFL1      TUBA6   FLJ20030       CLTC      RPS23 
0.08336046 0.08410148 0.09315528 0.09775742 0.10499056 0.10554332 
variable with highest stability (smallest rho value):	NACA
###############################################################
Step 9:
stability values rho:
      CFL1      TUBA6       CLTC   FLJ20030      RPS23 
0.07222968 0.07722737 0.08440691 0.09831958 0.12735605 
variable with highest stability (smallest rho value):	CFL1
###############################################################
Step 10:
stability values rho:
  FLJ20030      TUBA6       CLTC      RPS23 
0.08162006 0.08189011 0.10705192 0.11430674 
variable with highest stability (smallest rho value):	FLJ20030
###############################################################
Step 11:
stability values rho:
     TUBA6       CLTC      RPS23 
0.06978897 0.08069582 0.13702726 
variable with highest stability (smallest rho value):	TUBA6
###############################################################
Step 12:
stability values rho:
     CLTC     RPS23 
0.1199009 0.1245241 
variable with highest stability (smallest rho value):	CLTC
         1          2          3          4          5          6          7 
     "UBC"     "GAPD"     "TPT1"      "UBB"    "RPS13"     "ACTB"     "SUI1" 
         8          9         10         11         12 
    "NACA"     "CFL1" "FLJ20030"    "TUBA6"     "CLTC" 
Warning message:
In .local(x, ...) : Argument 'group' is transformed to a factor vector.

在Bladder数据集中，我们发现，排列前面的2个内参基因是HSPCB与RPS13，也就是Figure 1中的内容，如下所示：

1
2
3
4
5
6
7
8
grade <- pData(Bladder)[,"Grade"]
selectHKs(Bladder, 
	group = grade, 
	log = FALSE, 
	trace = FALSE,
	Symbols = featureNames(Bladder), 
	minNrHKs = 13,
	method = "NormFinder")$ranking

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
> grade <- pData(Bladder)[,"Grade"]
> selectHKs(Bladder, 
+           group = grade, 
+           log = FALSE, 
+           trace = FALSE,
+           Symbols = featureNames(Bladder), 
+           minNrHKs = 13,
+           method = "NormFinder")$ranking
         1          2          3          4          5          6          7 
   "HSPCB"    "RPS13"      "UBC"    "RPS23"    "ATP5B"     "TEGT"      "UBB" 
         8          9         10         11         12         13 
"FLJ20030"     "CFL1"   "S100A6"    "FLOT2"     "ACTB"     "TPT1"

使用多内参进行均一化

单内参的ΔCq方法

单内参均一化的方法就是，使用某个基因的Cq值减去内参的Cq值。NormqPCR中就有相关的计算函数，如下所示：

1
2
3
4
5
6
7
8
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "example.txt")
qPCR.example <- file.path(path, "qPCR.example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)
hkgs<-"Actb-Rn00667869_m1"
qPCRBatch.norm <- deltaCq(qPCRBatch = qPCRBatch.taqman, hkgs = hkgs, calc="arith")
head(exprs(qPCRBatch.taqman)) # raw data
head(exprs(qPCRBatch.norm)) # Normlazation data

这里我们使用了Actb作为内参，计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
> head(exprs(qPCRBatch.taqman)) # raw data
                     fp1.day3.v fp2.day3.v fp5.day3.mia fp6.day3.mia
Actb.Rn00667869_m1     20.92301   21.00262     20.88438     20.38190
Adipoq.Rn00595250_m1   20.93906   20.88610     23.81790     22.92289
Adrbk1.Rn00562822_m1         NA         NA     27.45101     27.02446
Agtrl1.Rn00580252_s1   25.82239   26.03846     27.28175     26.06274
Alpl.Rn00564931_m1     33.45495   32.81128     33.91955     32.62145
B2m.Rn00560865_m1      21.66457   21.89334     22.92485     22.61651
                     fp.3.day.3.v fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia
Actb.Rn00667869_m1       20.94854     21.48838       20.97863       21.10281
Adipoq.Rn00595250_m1     20.76957     20.92512       23.43714       23.83928
Adrbk1.Rn00562822_m1           NA           NA       24.71573       27.97638
Agtrl1.Rn00580252_s1     26.16934     25.91374       25.77341       26.44801
Alpl.Rn00564931_m1       33.34335     33.26128       33.08863       33.35799
B2m.Rn00560865_m1        21.45406     22.36598       22.90620       23.00608
> head(exprs(qPCRBatch.norm)) # Normlazation data
                     fp1.day3.v fp2.day3.v fp5.day3.mia fp6.day3.mia
Actb.Rn00667869_m1     0.000000   0.000000     0.000000     0.000000
Adipoq.Rn00595250_m1   0.016052  -0.116520     2.933523     2.540987
Adrbk1.Rn00562822_m1         NA         NA     6.566628     6.642561
Agtrl1.Rn00580252_s1   4.899380   5.035841     6.397364     5.680837
Alpl.Rn00564931_m1    12.531942  11.808657    13.035166    12.239549
B2m.Rn00560865_m1      0.741558   0.890717     2.040470     2.234605
                     fp.3.day.3.v fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia
Actb.Rn00667869_m1       0.000000     0.000000       0.000000       0.000000
Adipoq.Rn00595250_m1    -0.178971    -0.563263       2.458509       2.736475
Adrbk1.Rn00562822_m1           NA           NA       3.737100       6.873568
Agtrl1.Rn00580252_s1     5.220796     4.425364       4.794776       5.345202
Alpl.Rn00564931_m1      12.394802    11.772896      12.110000      12.255186
B2m.Rn00560865_m1        0.505516     0.877598       1.927563       1.903269

生成的结果是一个qPCRBatch对象，很多R包例如heatmap就兼容这类对象。

多内参的ΔCq方法

如果使用前面提到的geNorm或NormFinder算法找到了多个内参，此时想要用这多个内参进行均一化，那么就需要计算这几个内参的均值，构成一个“伪内参“(pseudo-housekeeper)进行均一化。现在我们使用三个内参，分别为GAPDH，beta-2-microglobulin和Beta-actin来进行均一化，如下所示：

1
2
3
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
qPCRBatch.norm <- deltaCq(qPCRBatch = qPCRBatch.taqman, hkgs = hkgs, calc="arith")
head(exprs(qPCRBatch.norm))

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
> head(exprs(qPCRBatch.norm))
                     fp1.day3.v fp2.day3.v fp5.day3.mia fp6.day3.mia
Actb.Rn00667869_m1   -1.2998917 -1.2816963    -1.380296   -1.5106197
Adipoq.Rn00595250_m1 -1.2838397 -1.3982163     1.553227    1.0303673
Adrbk1.Rn00562822_m1         NA         NA     5.186332    5.1319413
Agtrl1.Rn00580252_s1  3.5994883  3.7541447     5.017068    4.1702173
Alpl.Rn00564931_m1   11.2320503 10.5269607    11.654870   10.7289293
B2m.Rn00560865_m1    -0.5583337 -0.3909793     0.660174    0.7239853
                     fp.3.day.3.v fp.4.day.3.v fp.7.day.3.mia fp.8.day.3.mia
Actb.Rn00667869_m1     -1.1644617   -1.1714227      -1.323712      -1.286277
Adipoq.Rn00595250_m1   -1.3434327   -1.7346857       1.134797       1.450198
Adrbk1.Rn00562822_m1           NA           NA       2.413388       5.587291
Agtrl1.Rn00580252_s1    4.0563343    3.2539413       3.471064       4.058925
Alpl.Rn00564931_m1     11.2303403   10.6014733      10.786288      10.968909
B2m.Rn00560865_m1      -0.6589457   -0.2938247       0.603851       0.616992

单内参2^-ΔΔCq方法

使用2^-ΔΔCq方法哦可以计算两个样本类型之间的相对表达值。默认情况下，同一个样本的目标基因与其内参基因配对进行计算，此时的标准差就体现在不同基因与不同内参上。但是，如果参数中设定了paired=FALSE，那么标准差的计算公式就是s=sqrt(s1^2+s2^2)，这里的s1是目标基因的标准差，s2是内参基因的标准差，但是，如果当你要比较的基因与内参基因是来源于同一个样本时，不推荐这种方法，例如使用Taqman板的情况就是如此，但是，如果要考虑到完整性(completeness)以及要使用不同的生物学重复中获得的内参基因时，可以考虑使用这种计算方法，例如我们最初设计的内参数据并不好时，当采用事后检验(post hoc)时，就采用这种方法，或者是当NormqPCR用于处理那些单独板孔中的扩增数据时，也是采用这种方法。

现在我们读取原始数据，如下所示：

1
2
3
4
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "example.txt")
qPCR.example <- file.path(path, "qPCR.example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)

deltaDeltaCq()函数在进行计算时需要一个比较矩阵。在这个矩阵中，列是样本名（例如case或control），这与limma包中的数据类似。这个比较矩阵的列中只包括1或0这两个数据，用于指定分类变量，如下所示：

1
2
3
4
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
contM

结果如下所示：

1
2
3
4
5
6
7
8
9
10
> contM
               interestingPhenotype wildTypePhenotype
fp1.day3.v                        0                 1
fp2.day3.v                        0                 1
fp5.day3.mia                      1                 0
fp6.day3.mia                      1                 0
fp.3.day.3.v                      0                 1
fp.4.day.3.v                      0                 1
fp.7.day.3.mia                    1                 0
fp.8.day.3.mia                    1                 0

此时，我们可以指定一个内参基因来进行均一化，然后再看一个某个基因在实验组(case)和对照组(control)的比值，计算结果包括以下信息（按列显示）：

  1. 每个基因的表达值；
  2. 实验组的ΔCq值：
  3. Cq值的sd；
  4. 对照组的ΔCq值；
  5. 对照组的ΔCq值的sd；
  6. 2^-ΔΔCq值，这是对照组与实验组的差异；
  7. 和8. 分别是1和0时的sd值，如下所示：
1
2
3
4
5
6
7
8
9
10
hkg <- "Actb-Rn00667869_m1"
ddCq.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, 
                            maxNACase=1, maxNAControl=1,
                            hkg=hkg, 
                            contrastM=contM, 
                            case="interestingPhenotype",
                            control="wildTypePhenotype", 
                            statCalc="geom",
                            hkgCalc="arith")
head(ddCq.taqman)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
> head(ddCq.taqman)
                    ID 2^-dCt.interestingPhenotype interestingPhenotype.sd
1   Actb.Rn00667869_m1                   1.000e+00               0.000e+00
2 Adipoq.Rn00595250_m1                   1.587e-01               2.280e-02
3 Adrbk1.Rn00562822_m1                   2.602e-02               3.266e-02
4 Agtrl1.Rn00580252_s1                   2.300e-02               1.014e-02
5   Alpl.Rn00564931_m1                   1.892e-04               4.770e-05
6    B2m.Rn00560865_m1                   2.464e-01               2.498e-02
  2^-dCt.wildTypePhenotype wildTypePhenotype.sd              2^-ddCt 2^-ddCt.min
1                1.000e+00            0.000e+00 1                             NA
2                1.171e+00            2.131e-01 0.135541545192243             NA
3                       NA                   NA +                             NA
4                3.434e-02            8.584e-03 0.669721905042939             NA
5                2.298e-04            6.107e-05 0.823327272466571             NA
6                5.965e-01            7.668e-02 0.413128242070071             NA
  2^-ddCt.max
1          NA
2          NA
3          NA
4          NA
5          NA
6          NA

也可以使用每个单独的样本/孔(samples/wells)方法来计算平均值。这里我们单独计算sd，然后将它们汇总起来。因此，同一个样本的内参与检测值的本领信息氷是人干活的，这同时会增加变异，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
hkg <- "Actb-Rn00667869_m1"
ddCqAvg.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, 
                               maxNACase=1, 
                               maxNAControl=1,
                               hkg=hkg, 
                               contrastM=contM, 
                               case="interestingPhenotype",
                               control="wildTypePhenotype", 
                               paired=FALSE, statCalc="geom",
                               hkgCalc="arith")
head(ddCqAvg.taqman)

运行结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
> hkg <- "Actb-Rn00667869_m1"
> ddCqAvg.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, 
+                                maxNACase=1, 
+                                maxNAControl=1,
+                                hkg=hkg, 
+                                contrastM=contM, 
+                                case="interestingPhenotype",
+                                control="wildTypePhenotype", 
+                                paired=FALSE, statCalc="geom",
+                                hkgCalc="arith")
There were 12 warnings (use warnings() to see them)
> 
> head(ddCqAvg.taqman)
                    ID 2^-dCt.interestingPhenotype interestingPhenotype.sd
1   Actb.Rn00667869_m1                   1.000e+00               0.000e+00
2 Adipoq.Rn00595250_m1                   1.587e-01               2.280e-02
3 Adrbk1.Rn00562822_m1                   2.602e-02               3.266e-02
4 Agtrl1.Rn00580252_s1                   2.300e-02               1.014e-02
5   Alpl.Rn00564931_m1                   1.892e-04               4.770e-05
6    B2m.Rn00560865_m1                   2.464e-01               2.498e-02
  2^-dCt.wildTypePhenotype wildTypePhenotype.sd              2^-ddCt 2^-ddCt.min
1                1.000e+00            0.000e+00 1                             NA
2                1.171e+00            2.131e-01 0.135541545192243             NA
3                       NA                   NA +                             NA
4                3.434e-02            8.584e-03 0.669721905042939             NA
5                2.298e-04            6.107e-05 0.823327272466571             NA
6                5.965e-01            7.668e-02 0.413128242070071             NA
  2^-ddCt.max
1          NA
2          NA
3          NA
4          NA
5          NA
6          NA

多内参2^-ΔΔCq方法

如果想使用大于1个内参进行均一化，例如，我们想使用NormFinder/geNorm算法来进行均一化。我们就可以使用几何平均数(geometric mean)来构建一个伪内参(pseudo-housekeeper)来实现这一目的。现在我们使用GAPDH，beta-2-microglobulin和Beta-actin这三个内参来进行均一化，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
qPCRBatch.taqman <- read.taqman(taqman.example)
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
ddCq.gM.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman,
                               maxNACase=1,
                               maxNAControl=1,
                               hkgs=hkgs, 
                               contrastM=contM,
                               case="interestingPhenotype",
                               control="wildTypePhenotype", 
                               statCalc="arith",
                               hkgCalc="arith")
head(ddCq.gM.taqman)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
> qPCRBatch.taqman <- read.taqman(taqman.example)
Warning message:
In read.taqman(taqman.example) :
  Incompatible phenoData object. Created a new one using sample name data derived from raw data.
> contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
> colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
> rownames(contM) <- sampleNames(qPCRBatch.taqman)
> hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
> ddCq.gM.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman,
+                                maxNACase=1,
+                                maxNAControl=1,
+                                hkgs=hkgs, 
+                                contrastM=contM,
+                                case="interestingPhenotype",
+                                control="wildTypePhenotype", 
+                                statCalc="arith",
+                                hkgCalc="arith")
There were 12 warnings (use warnings() to see them)
> head(ddCq.gM.taqman)
                    ID 2^-dCt.interestingPhenotype interestingPhenotype.sd
1   Actb.Rn00667869_m1                   2.594e+00                 0.09819
2 Adipoq.Rn00595250_m1                   4.083e-01                 0.24929
3 Adrbk1.Rn00562822_m1                   4.182e-02                 1.45844
4 Agtrl1.Rn00580252_s1                   5.520e-02                 0.63719
5   Alpl.Rn00564931_m1                   4.767e-04                 0.42589
6    B2m.Rn00560865_m1                   6.367e-01                 0.05413
  2^-dCt.wildTypePhenotype wildTypePhenotype.sd              2^-ddCt 2^-ddCt.min
1                2.345e+00             0.071373 1.10638851325547       1.034e+00
2                2.713e+00             0.201905 0.150497255530234      1.266e-01
3                       NA                   NA +                             NA
4                7.878e-02             0.333840 0.700597907024805      4.505e-01
5                5.242e-04             0.386280 0.909381199520663      6.769e-01
6                1.390e+00             0.163975 0.457939394245865      4.411e-01
  2^-ddCt.max
1    1.184310
2    0.178884
3          NA
4    1.089636
5    1.221662
6    0.475448

我们也可以使用那些在样本间方差相同的内参来进行均一化，这个类似于前面提到的第二个deltaDeltaCq方法，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
qPCRBatch.taqman <- read.taqman(taqman.example)
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
ddAvgCq.gM.taqman <-deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, 
                                 maxNACase=1, maxNAControl=1,
                                 hkgs=hkgs, 
                                 contrastM=contM, 
                                 case="interestingPhenotype",
                                 control="wildTypePhenotype",
                                 paired=FALSE, 
                                 statCalc="arith",
                                 hkgCalc="arith")
head(ddAvgCq.gM.taqman)

计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
> qPCRBatch.taqman <- read.taqman(taqman.example)
Warning message:
In read.taqman(taqman.example) :
  Incompatible phenoData object. Created a new one using sample name data derived from raw data.
> contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
> colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
> rownames(contM) <- sampleNames(qPCRBatch.taqman)
> hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
> ddAvgCq.gM.taqman <-deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, 
+                                  maxNACase=1, maxNAControl=1,
+                                  hkgs=hkgs, 
+                                  contrastM=contM, 
+                                  case="interestingPhenotype",
+                                  control="wildTypePhenotype",
+                                  paired=FALSE, 
+                                  statCalc="arith",
+                                  hkgCalc="arith")
There were 12 warnings (use warnings() to see them)
> head(ddAvgCq.gM.taqman)
                    ID 2^-dCt.interestingPhenotype interestingPhenotype.sd
1   Actb.Rn00667869_m1                   2.594e+00                  0.3849
2 Adipoq.Rn00595250_m1                   4.083e-01                  0.4822
3 Adrbk1.Rn00562822_m1                   4.182e-02                  1.4545
4 Agtrl1.Rn00580252_s1                   5.520e-02                  0.6905
5   Alpl.Rn00564931_m1                   4.767e-04                  0.5846
6    B2m.Rn00560865_m1                   6.367e-01                  0.2777
  2^-dCt.wildTypePhenotype wildTypePhenotype.sd              2^-ddCt 2^-ddCt.min
1                2.345e+00               0.3574 1.10638851325547       8.473e-01
2                2.713e+00               0.2495 0.150497255530234      1.077e-01
3                       NA                   NA +                             NA
4                7.878e-02               0.2813 0.700597907024805      4.341e-01
5                5.242e-04               0.3689 0.909381199520663      6.064e-01
6                1.390e+00               0.4576 0.457939394245865      3.778e-01
  2^-ddCt.max
1    1.444684
2    0.210221
3          NA
4    1.130625
5    1.363762
6    0.555126

NRQ值计算

导入数据，如下所示：

1
2
3
path <- system.file("exData", package = "ReadqPCR")
qPCR.example <- file.path(path, "qPCR.example.txt")
Cq.data <- read.qPCR(qPCR.example)

合并技术重复，并计算SD，如下所示：

1
Cq.data1 <- combineTechRepsWithSD(Cq.data)

加载数据集的效应(efficiencies)，并将它们添加到数据集中，如下所示：

1
2
3
4
5
Effs <- file.path(path, "Efficiencies.txt")
Cq.effs <- read.table(file = Effs, row.names = 1, header = TRUE)
rownames(Cq.effs) <- featureNames(Cq.data1)
effs(Cq.data1) <- as.matrix(Cq.effs[,"efficiency",drop = FALSE])
se.effs(Cq.data1) <- as.matrix(Cq.effs[,"SD.efficiency",drop = FALSE])

计算数据集的NRQ(normalized relative quantities)，我们只考虑两个特征值作为内参基因，如下所示：

1
2
3
4
5
6
res <- ComputeNRQs(Cq.data1, hkgs = c("gene_az", "gene_gx"))
## NRQs
exprs(res)
## SD of NRQs
se.exprs(res)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
> res <- ComputeNRQs(Cq.data1, hkgs = c("gene_az", "gene_gx"))
Warning message:
In .local(qPCRBatch, ...) : This function is still experimental!
> ## NRQs
> exprs(res)
            caseA     caseB  controlA  controlB
gene_ai 1.9253072 1.3586729 0.6479659 0.8749479
gene_az 1.0567118 1.1438982 1.0331980 0.9134997
gene_bc 1.1024935 0.7193500 0.7030487 1.2140836
gene_by 1.5102316 0.9573047 0.7527082 1.6008850
gene_dh 1.2982037 1.0722522 0.9623335 0.9392871
gene_dm 0.6590246 1.1690720 1.2475372 0.9366210
gene_dq 0.7541955 0.7036408 0.8327917 1.6165326
gene_dr 2.2192305 1.0581211 0.7026411 0.6900584
gene_eg 0.9366671 0.5800339 0.8313720 1.1848856
gene_er 0.5269062 0.9375427 0.6953326 2.3195978
gene_ev 1.4622280 2.3457021 0.9038912 1.1454535
gene_fr 1.4954763 1.6200792 0.9641192 0.7295680
gene_fw 0.6944248 0.8051075 1.5698382 0.7978611
gene_gx 0.9463318 0.8742037 0.9678687 1.0946911
gene_hl 1.0009372 1.4015267 0.7683665 0.7713712
gene_il 1.4632019 1.2595559 0.7216891 0.9318860
gene_iv 1.7263335 1.2275001 1.5464212 0.8881605
gene_jr 0.8984351 0.9834026 0.8754813 0.6637941
gene_jw 1.4655948 0.9340184 1.0505200 1.5504136
gene_qs 0.6730225 0.7610418 1.0665938 3.5329891
gene_qy 0.5287127 1.5722670 1.0615326 3.3252907
gene_rz 0.8690600 1.5588299 0.7287288 1.4812753
gene_sw 0.5975288 1.2406438 0.6982954 1.6007333
gene_vx 0.6942254 0.7168408 2.0253177 1.3190943
gene_xz 0.7668030 1.0218209 0.6136038 1.6729352
> ## SD of NRQs
> se.exprs(res)
            caseA     caseB  controlA  controlB
gene_ai 1.3996554 0.8787290 0.4855882 1.0034912
gene_az 0.6832730 0.7601966 0.8971054 0.6031927
gene_bc 0.7225348 0.4746146 0.9626570 1.0478385
gene_by 1.1522746 0.6116269 0.6088836 2.0409211
gene_dh 1.2483072 0.7889984 0.6165041 0.9767947
gene_dm 0.4711409 0.7780238 0.8476053 0.7294405
gene_dq 0.7023561 0.4849899 0.5813310 1.4067670
gene_dr 1.4407662 1.0804211 0.4543153 0.5149367
gene_eg 0.7355269 0.5497433 0.5588801 1.0938601
gene_er 0.4301195 0.6119514 0.4471454 1.5115897
gene_ev 1.0094209 2.4267114 0.6337126 0.7782519
gene_fr 1.6760391 1.1119157 0.6226081 0.5040967
gene_fw 0.5041070 0.9131565 1.1153268 0.7234551
gene_gx 0.6046042 0.9027816 0.6713914 1.7394961
gene_hl 0.7633174 0.9123997 1.0000329 0.5005813
gene_il 1.4621406 0.9540445 0.5678634 0.6067147
gene_iv 1.2668346 0.8039841 0.9995225 0.8171996
gene_jr 0.5749672 0.6786989 0.5595295 0.4405919
gene_jw 0.9626606 0.7890401 0.7194378 1.1512512
gene_qs 0.5830335 0.5309990 0.6828952 2.8253799
gene_qy 0.5947918 1.1294199 0.6794829 2.1713340
gene_rz 0.5846751 1.7926435 0.4911506 2.1424580
gene_sw 0.6284440 0.8062083 0.9638307 1.8398593
gene_vx 0.5285361 1.0126959 1.3861226 0.8683886
gene_xz 0.5231477 0.9270275 0.3972901 1.3643840

参考资料

  1. NormqPCR: Functions for normalisation of RT-qPCR data

发表于 | 分类于
| 字数统计: 4,416 | 阅读时长 ≈ 16

试剂配置

电泳液(10X)

称量
Tris 30.3 g
Glycine 144 g
SDS 10g

转膜液(10X)

成分 称量
Tris 58g
Glycine 29g
SDS（电泳级） 3.8g

TBS(10X)

成分 称量
NaCl 80g
1M Tris-HCl((PH8.0) 24g

使用时，取100mL TBS(10X)，加入900mL去离子水，再加入1mL 的Twen20（Sigma，货号：T2700-1L）。

其它试剂

封闭液

5%脱脂奶（BD:232100），或索莱宝（货号D8340）

2.5g脱脂奶粉加入20ml的TBST Wash Buffer中，搅拌溶解，4℃或-20℃保存。

封闭液也可以直接使用商业公司的产品，例如碧云天的QuickBlock封闭液，货号：P0228

一抗二抗稀释液

一抗稀释液：碧云天Western一抗稀释液（货号：P0023A）

二抗稀释液：碧云天Western二抗稀释液（货号：P0023D）

或用5%脱脂奶粉

常用内参一抗与二抗

一抗：鼠源小鼠actin抗体（碧云天，货号：AA128）

一抗：兔源小鼠Gapdh抗体（碧云天，货号：AF1186）

兔抗鼠Vinculin(CST，货号：13901S)

兔抗鼠beta-tubulin(CST，货号：15115S)

羊抗兔二抗（鼎国，货号：SA00001-2）

羊抗小鼠二抗（碧云天：货号：A0216）

羊抗兔二抗（CST，货号：7074S）。

(4) 蛋白Marker

ThermoFisher protein Marker（货号：26619）

其它试剂

BCA蛋白浓度测定试剂盒（碧云天，货号：P0010）
Twen20（Sigma，货号：T2700-1L）。

实验步骤

（一）样本裂解

  1. 将2×106的细胞悬液离心后，除去上清，加入RIPA溶液40µL（根据细胞的数目而定），反复吹打，随后放冰上静置10分钟；（注：RIPA是一种传统的细胞组织快速裂解液。成分为 ①50 mMTris-HCl(pH 7.4)② 150 mM NaCl③1% NP-40）；
  2. 12000g，离心10min，4℃，取上清，检测蛋白含量。

（二）蛋白浓度测定(BCA法)

试剂盒：碧云天（BCA蛋白浓度测定试剂盒，货号：P0010）或BCA蛋白浓度测定试剂盒(增强型，货号：P0012)
①BCA试剂A②BCA试剂B③蛋白标准品（5mg/mL BSA，BSA为牛血清白蛋白，bovine serum albumin）。BCA试剂A和B室温保存，蛋白标准于-20℃冻存（提前可以配置为梯度BSA溶液，于负20度保存）；④5XSDS变性剂（碧云天，货号：P0015）

检测原理 碱性条件下，蛋白将Cu2+还原为Cu+， Cu+与BCA试剂形成紫颜色的络合物，测定其在562nm处的吸收值，并与标准曲线对比，即可计算待测蛋白的浓度。

检测流程：

  1. 根据样品数量，按试剂A：试剂B=50：1的比例进行稀释（检测的时候，此工作的体积为180µL，样品以及标准品的体积为20µL，因此要按照这个体积以及相当样品数来进行稀释），BCA工作液室温24小时内稳定（标准品的梯度为8个，样品是N个，则BCA工作液的体积一般是200x(8+N+1)µL）。注：根据试剂盒的说明进行操作，例如我使用的P0010试剂盒，里面是要求加200µL的反应液。
  2. 取提前配置好的蛋白标准品梯度溶液，即0µg/µL（编号1），0.025µg/µL（编号2），0.05µg/µL（编号3），0.1µg/µL（编号4），0.2µg/µL（编号5），0.3µg/µL（编号6），0.4µg/µL（编号7），0.5µg/µL（编号8）。
  3. 在96孔板相应的孔中加入200µL的BCA工作液，然后将编号为1到8的标准品加入到这些孔中，体积为20µL；
  4. 将样品稀释10或20倍，然后将20µL的样品加入到样品孔中（注：在这一步的操作中，通常是先在96孔的孔中提前加入9μL或19μL的双蒸水，然后再加入1μL的样本，最后加入200μL的反应液）；
  5. 在37℃条件下放置30分钟；
  6. 测定OD562nm处的吸光度，根据标准曲线计算出蛋白浓度，这一步注意，通常我们是使用y = ax+b来计算的，也就是x轴是吸光度，y轴是蛋白浓度，计算后再乘以稀释倍数即可。

（三）蛋白变性

  1. 将准备上样的所有蛋白样品调至等浓度（一般上样量最多是20µg，体积是10µL，因此蛋白的终浓度控制在2µg/µL，但根据Biorad公司的说明，1.5mm的胶，大梳子每个孔最多上样66uL，但推荐最大上样是50μL，小梳子最多一个孔上样40μL），蛋白与5×SDS－PAGE Loading Buffer以4:1 混匀（体积比），1×loading buffer补充至等体积，SDS要在不超过37度的水浴中进行融化，用完放到负20度即可。

建议：蛋白样品在加入SDS之前最好统一稀释至浓度为2.5µg/uL（因为最终上样体积为10uL，上样蛋白为20ug，其中在变性前，1次实验的蛋白样本为8µL，需要加入含有5XSDS的Loading Buffer为2uL，因此蛋白浓度最终就是2µg/µL）。每次制备的样品需要按3次实验进行分装，即每个试管中装24uL，每支的浓度为2.5µg/µL，使用的时候，直接加6uL的Loading Buffer，一次上样为10uL;如果细胞数量较多，蛋白浓度较大，可以将5XSDS稀释为1XSDS,并对蛋白样本进行稀释后再变性。

  1. 100℃干热或水浴（封口）10min充分变性，12000g离心30s后，置于冰上或4℃保存。

（四）SDS-PAGE胶配置

  1. 清洗玻璃板：根据实验需要选择一定厚度的玻璃板 （1.0mm，1.5mm），注意一定要将玻璃板洗净（洗洁精），最后用无水乙醇冲洗，与胶接触的一面向上置于干净的纸巾晾干或用电吹风吹干。
  2. 装置玻璃板：将两块玻璃板装好，低的玻璃板朝内，注意玻璃底部平整，防止漏胶，四个螺丝钉一起拧紧，防止玻璃破裂，再固定于制胶板上，可同时做一块胶（对面有塑料挡板）或两块胶，在加入分离胶之前，可以加入无水乙醇，观察一下是否漏液；如使用Bio-Rad, 将两块玻璃板装好，低的玻璃板朝外，夹好玻璃。
  3. 配分离胶：分离胶的其它成分可事先配好（事先配置好的意思就是，先把除了10% AP（过硫酸铵）及TEMED外的其它成分加好，在实验之前再加入AP与TEMED），避光存放于4℃，可至少存放1个月，临用前取出室温平衡（否则凝胶过程产生的热量会使低温时溶解于储存液中的气体析出而导致气泡，有条件者可真空抽吸3min），加入10%AP及TEMED即可，也可临时按比例配制，注意胶用枪吹打混匀（鉴定的Rab32蛋白是按照10%SDS-PAGE的胶进行配置，体积为5mL，不同的蛋白分子量需要不同浓度的胶，通常是分子量越大的蛋白，需要浓度越低，同时胶也越脆弱）注：试剂盒的A盒（即带有TEMED的那个试剂盒）放4℃保存；B盒常温下保存。
  4. 灌胶：将凝胶溶液平缓地注入两层玻璃极中，灌入分离胶后（可预先插入梳子，梳子底下约0.5cm处标记预备配制的分离胶高度，拔出梳子）应立即封胶，可用在液面上小心注入一层无水乙醇，以阻止氧气进入凝胶溶液中，然后静置，封胶后切记勿动，分离胶的上层距离低的玻璃板的最上层距离约为1.5~2cm，配置分离胶的体积约为5mL，用掉约4.5mL。
  5. 待分离胶凝后（30min左右出现一条明显的分离线）将无水乙醇倒掉，用滤纸或纸巾吸干。
  6. 配浓缩胶： 同前按比例配制浓缩胶，一个板子约2mL，以连续平稳的液流注入凝胶溶液，然后小心插入梳子（不同厚度的玻璃板匹配不同厚度的梳子），注意不得在齿尖留有气泡，静置1h以保证完全聚合。

(5) 电泳

电泳：上样前将胶板下的气泡赶走，用注射器吸取电泳缓冲液冲孔，可根据实验需要设定阳性对照（即肯定可以表达你的目的蛋白的组织或者细胞，同时可以提示你一抗的质量）或阴性对照（即肯定不会表达你的目的蛋白的组织或者细胞），样品两侧的泳道用等体积的1×loading buffer上样，Marker（Ferments）也用1×loading buffer调整至与样品等体积，上样完毕后宜尽快电泳（注意电极正负）防止扩散，以初始电压为80V进行稳压电泳，当溴芬兰泳动至分离胶时，这个过程约为20min~30min，改电压为120 V继续稳压电泳，最终的蓝色缓冲液条件不要超过底部支架的绿色部分，电泳的过程约为60min~90min，整个过程要注意电泳的条件位置，不同分子量的蛋白电泳的时间与位置不同。

(6) 转膜

  1. 取胶：在目的条带充分分离后（参考Marker的分离效果），或电泳至溴酚染料前沿下至凝胶末端处可结束电泳，将胶卸下，保留30-100KD或分子量范围更广些的胶（以便以后杂其他感兴趣的蛋白），在转移液（转移液即1X的transfer buffer，要现用现配）中浸泡。
  2. 转膜：转膜液是800mL的1X转膜液加上200mL的甲醇，根据胶的大小剪裁稍大的PVDF膜，在左上角剪角，甲醇中充分浸泡5min（活化PVDF膜）后转入转移液中，去掉PVDF膜的上下两层塑料封纸，将厚的白色滤纸及海绵也浸入转移液中。15min-30min后按如下顺序安装转移装置（如图1所示）：

a.平放底部电极（负极），放一张海绵垫片。
b.在海绵垫片上放置1张用转移缓冲液浸泡过的滤纸，然后用一玻璃棒作滚筒以挤出所有气泡，必要时可滴加转膜液润湿。
c.取出浸在转膜液中的凝胶平放于滤纸上，排除所有气泡。
d.把PVDF膜放在聚丙烯酰胺凝胶上，在硝酸纤维素滤膜与聚丙烯酰胺凝胶之间不留有气泡。
e. 把最后1张滤纸放在PVDF膜上方，同样须确保不留气泡。
f.放上另一张海绵垫片，盖上阳极板，夹紧。保证对凝胶有一定的压力。加入预冷的转膜液及内置的冰盒，盖好盖子，连接电源，（注意正负极）电转移40min, 300mA（根据实验需要摸索最佳转膜条件）。电转过程中，尽量保持低温环境，可放在冰水里转膜。电转结束后，断开电源，拆卸转移装置，逐一掀去各层。将凝胶转移至盛有考马斯亮蓝染液的托盘中，进行染色，可检查蛋白质转移是否完全。要观察膜上蛋白的情况可将膜在丽春红中染色，然后蒸馏水冲洗。（电泳上槽液及转膜液可回收）。

  1. 如果不染色，可以跳过后面几步，直接进行封闭。染色和脱色 ：电泳完毕可通过染色和脱色评定其结果的优劣.对凝胶中分离成不同条带的蛋白进行检测.此外,根据不同的研究目的,也可将凝胶电印迹后的胶或膜染色。
  2. 考马斯亮蓝染色：将凝胶取出放入塑料盒中倒入少量的染色液,以浸没凝胶为宜,在摇床上震荡,直到凝胶上出现明显的条带为止,一般5-6小时或者4℃过夜;然后倾去染色液，用脱色液洗涤震荡,其间要不断换脱色液,直至脱色液颜色稍清亮.,凝胶背景清楚为止。（染色液可回收）
  3. 丽春红染色：详细方法见说明书。

(7) 膜的封闭

在进行抗体杂交之前，需要先对转印膜进行封闭，防止免疫试剂的非特异性吸附。封闭一般采用异源性蛋白质或去污剂，常用的有0.2%Tween-20，20%BSA，10%马血清以及5% Non-fat milk等，至于选择哪一类封闭液，首先应考虑与检测试剂相适应，如采用葡萄球蛋白A（SPA）作用检测试剂，就不能用全血清封闭，其次是尽可能使非特异着色背静浅，本实验室常用封闭液5%Non-fat milk。封闭过程如下所示（与胶接触的PVDF膜的一面是正面，要在正面的左上角剪下一小块，作为区分正反面的依据，前面剪过的则忽略此步）：

  1. 洗转印膜：用1xTBS-T室温漂洗3次×10min（漂洗1次也行），以尽量洗去转印膜上的SDS，防止影响后面的抗体结合。
  2. 取漂洗的转印膜，放入5% Non-fat milk的封闭液内，摇床震动，室温封闭2h，也可在4℃过夜， 或使用商用封闭液。
  3. 封闭后，使用1×TBS-T室温漂洗3次×10min。

(8) 抗体杂交

抗体表面主要采用间接法。即先加入未标记特异性抗体（Ab1，即一抗）与膜上抗原结合，再加入标记的抗抗体（Ab2，即二抗）进行杂交检测，标记Ab2 物质有放射性核素酶及生物素等。 （抗体稀释倍数需要摸索通常一抗浓度1：2000到1:1000，具体的稀释浓度可以参考一抗的说明书，二抗浓度1:1000）。杂交过程如下所示：

  1. 封闭后的杂交膜放入塑料盒或杂交袋中，加入抗体稀释液稀释的Ab1（宜设立内参照，beta-actin，beta-tubulin等），4℃孵育过夜或室温（22-25℃）摇动孵育2h。（一抗可回收3-5次）；
  2. 1×TBS-T洗膜3次×10min；
  3. 抗体稀释液稀释的Ab2（抗兔-HRP或抗鼠-HRP，根据一抗种属来源决定，注意避光）室温2h，随后洗膜3×10min。
  4. ECL显色：将显色剂A、B按1：1 混匀，将杂交后PVDF膜放入显色盒中，加上混合好的显色液，约1-5min用纸巾吸去印迹膜边缘或者边角部分多余的显色液，在暗室中使胶片曝光1-5min。注：不同的显色仪参数不同，需要阅读一下说明书，找到最优曝光条件。

WB常见问题

  1. 如果我想观察2个蛋白，但这两个蛋白的分子量离得很近，怕剪膜剪坏，如何操作？

    答：通常这种情况，可以先孵A蛋白，然后将A蛋白的抗体洗掉，再进行封闭操作，随后孵B蛋白的抗体。

    洗抗体通常会使用抗体去除剂来进行操作，这样通常可以重复利用膜3-5次。抗体去除剂通常有中性，酸性，碱性，要根据自己的目的来进行选择，官网通常就有，例如碧云天的官网《碧云天生产的各种Western一抗二抗去除液的主要特点、差异和选择》。

  2. 如果我想用5%的脱脂奶粉来配制一抗和二抗，但由于担心不好保存怎么办？

    答：如果要用奶粉来配制抗体，通常会加入叠氮化钠。实验室配制的叠氮钠是10%的储备液，常温放置备用，配制一抗时每1mL加入叠氮化钠10μL,能使一抗保存时间一个月左右，使用叠氮化钠时一定要小心，此试剂有毒。

  3. 转膜时间与蛋白分子量有没有关系？

    答：通常如果使用Bio-Rad的标准湿式转膜装置，可以设定转膜电流为300-400mA，转膜时间为30-60分钟。也可以在15-20mA转膜过夜。具体的转膜时间要根据目的蛋白的大小而定，目的蛋白的分子量越大， 需要的转膜时间越长，目的蛋白的分子量越小，需要的转膜时间越短。

  4. 转膜的效率与分子量大小有关系没？

    答：如果是大蛋白（大于100 kDa），需要注意以下几点：

    1. 对于大蛋白而言，其在凝胶电泳分离迁移较慢，而从凝胶转出也非常慢，因此对于这种大分子量蛋白应该用低浓度的凝胶，8% 或更低，但因低浓度的胶非常易碎，所以操作时需十分小心。
    2. 大蛋白易在凝胶里形成沉淀，从而影响转膜。转膜时在电转移缓冲液加入SDS 至终浓度0.1%，避免出现这种情况。甲醇易使SDS 从蛋白上脱失，因此相应降低甲醇的浓度至10% 或更低，以防止蛋白沉淀。
    3. 降低电转移缓冲液中甲醇的比例，这样可以促进凝胶的膨胀，更易于大蛋白的转出。
    4. 只有使用硝酸纤维素膜时，甲醇才是必需的。如果是PVDF 膜，甲醇可以不必加入电转移缓冲液中，但转膜前PVDF 需用甲醇活化。
    5. 选择湿式，4℃ 转膜过夜，以取代半干式转膜。

参考资料

  1. Western Blot实验——蛋白转膜和染色篇

发表于 | 分类于
| 字数统计: 1,519 | 阅读时长 ≈ 8

性能增强概览

软件系统和数据分析的可扩展性(scalability)的概念是复杂的，精确的度量(precise metrics)和标准现在还不好处理。可以参考一下来自计算机界的 Weinstock and Goodenough 的一个概述。基于Bioconductor的分析工作流的可扩展性(scalability)可以沿着不同的路径进行，Bioconductor的两个核Mike Lawrence和Martin Morgan在工作就提到了许多有意的想法。

本文档仅解决Bioconductor直接支持的并行计算，主要通过两种方法实现：

  • 共享内存并行性。使用内存映像的完整副本对R进程进行任意次数的分叉，并且对每个映像独立地进行计算。选定的结果将返回到主进程。可以有效地在多环境中实现。
  • 分布式并行性。可以运行不同操作系统的独立处理器可以运行R的兼容实例。R或集群调度程序可以执行作业控制(job control)。

简要说明

现在我们来说明一下共享内存方法：

1
2
3
4
5
6
7
8
9
10
system.time( lapply(1:8, function(x)Sys.sleep(1) ) )
##    user  system elapsed 
##   0.004   0.001   8.020
library(parallel)
detectCores()
## [1] 4
options(mc.cores=4)
system.time( mclapply(1:8, function(x)Sys.sleep(1) ) )
##    user  system elapsed 
##   0.007   0.022   2.031

上面只是一个无意义的计算，我们就是为了说明并行计算的原理，通过并行计算，我们加快了计算速度，计算过程的时间缩短为原来的四分之一。

RNA-seq实验中的BAM文件

我们将使用 Zarnack et al. 2013 等的人研究中的一些BAM文件。2013年的这项假设是，有一类蛋白质，即异质核糖蛋白C1和C2，也就是指的hnRNP C(heterogeneous nuclear ribonucleoproteins C1 and C2)，这个蛋白会阻止Alu因子整合到成熟到转录本上。这样的整合会导致蛋白功能障碍，进而导致疾病；如果想了解更深入的信息，可以去查看文献原文。

我们下面所使用的包含有8个BAM文件，它们已经与chr14染色体进行了比对。其中4个文件是正常的HeLa细胞的读长文件(reads)（阴性对照），另外4个文件是敲低了hnRNP C后的HeLa的读长文件，如下所示：

1
2
3
4
5
6
7
8
9
10
11
library(RNAseqData.HNRNPC.bam.chr14)
dir(system.file("extdata", package="RNAseqData.HNRNPC.bam.chr14"))
##  [1] "ERR127302_chr14.bam"     "ERR127302_chr14.bam.bai"
##  [3] "ERR127303_chr14.bam"     "ERR127303_chr14.bam.bai"
##  [5] "ERR127304_chr14.bam"     "ERR127304_chr14.bam.bai"
##  [7] "ERR127305_chr14.bam"     "ERR127305_chr14.bam.bai"
##  [9] "ERR127306_chr14.bam"     "ERR127306_chr14.bam.bai"
## [11] "ERR127307_chr14.bam"     "ERR127307_chr14.bam.bai"
## [13] "ERR127308_chr14.bam"     "ERR127308_chr14.bam.bai"
## [15] "ERR127309_chr14.bam"     "ERR127309_chr14.bam.bai"
## [17] "tophat2_out"

这些都是由Tabix建立索引的BAM文件。其中*.bam 是原始的BAM文件，而*.bam.bai则是添加了索引后的BAM文件，从上面我们可以看到有8个bam.bai文件，前面四个是敲低组，后面四个是对照组。

使用BiocParallel实现隐式并行

为了促进可靠的执行程序的和谐发展，Bioconductor开发了BiocParallel包。在下面的案例中，我们将读取HNRNP C案例数据集中的reads读取到由外显子地址定义的bins中（这里我不清楚bin是什么，列出了原文）：

Consider the following example: we will count reads into bins defined by exon addresses in the HNRNPC example dataset.

RNAseqData.HNRNPC.bam.chr14 这个包含有BAM文件的绝对路径的一个向量，我将其分配给fns，如下所示：

1
2
3
4
library(RNAseqData.HNRNPC.bam.chr14)
fns = RNAseqData.HNRNPC.bam.chr14_BAMFILES
length(fns)
## [1] 8

Here we establish the exon bins into which we will count reads.

1
2
3
4
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
seqlevels(txdb) = "chr14"
ebg = exonsBy(txdb, by="gene")

Now we use the summarizeOverlaps function from the GenomicAlignments package to count reads into the exon bins. We’ll time the counting from a single file, and then time the counting from all files at once.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
library(GenomicAlignments)
# summarizeOverlaps uses bplapply to iterate over files
s1 = system.time(i1 <- summarizeOverlaps( ebg, fns[3] ))
s1
##    user  system elapsed 
##   2.975   0.178   3.178
# show implicit config
BiocParallel::bpparam()
## class: MulticoreParam
##   bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB
##   bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
##   bptimeout: 2592000; bpprogressbar: FALSE
##   bpRNGseed: 
##   bplogdir: NA
##   bpresultdir: NA
##   cluster type: FORK
s2 = system.time(i2 <- summarizeOverlaps( ebg, fns ))
s2
##    user  system elapsed 
##   0.271   0.162  10.480

This is not a thorough way of measuring speedup but it shows reasonable enhancement. In the second computation, we did approximately 8 times as much computation, but the clock time elapsed increased only by a factor of (3.3). The default behavior of BiocParallel on the MacBook Air used to produce this document is to pick up the value of the option mc.cores and use this as the number of workers in a MulticoreParam configuration; if options()$mc.cores is NULL, 2 workers are specified.

When BiocParallel is attached, the summarizeOverlaps function will iterate over the files using bplapply from the BiocParallel package. That function will check the R session for specific parallelization configuration information. If it finds none, it will check for multiple cores and make arrangements to use them if present. The “check” occurs via the function bpparam.

The default situation on a MacBook Air running MacOSX 10.9.5 with an Intel core i7 processor, (two physical cores with two logical cores each, allowing for four concurrent threads) as follows.

1
2
3
4
5
6
7
8
9
10
11
12
options(mc.cores=NULL)
library(BiocParallel)
register(MulticoreParam())
bpparam()
## class: MulticoreParam
##   bpisup: FALSE; bpnworkers: 2; bptasks: 0; bpjobname: BPJOB
##   bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
##   bptimeout: 2592000; bpprogressbar: FALSE
##   bpRNGseed: 
##   bplogdir: NA
##   bpresultdir: NA
##   cluster type: FORK

This identifies an object called MulticoreParam which is used to configure the behavior of bplapply and other utilities of the BiocParallel package. There are various configuration object classes that can be used.

1
2
3
4
5
6
7
8
9
‘SnowParam’: distributed memory computing
‘MulticoreParam’: shared memory computing
‘BatchJobsParam’: scheduled cluster computing
‘DoparParam’: foreach computing
‘SerialParam’: non-parallel execution

We need to use register to determine the type of concurrent computation that will be performed.

If process size is large, we may want to leave some cores idle. We can accomplish that by using register. Here we’ll check for any advantage of using all logical cores.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(BiocParallel)
register(MulticoreParam(workers=1))
system.time(i3 <- summarizeOverlaps( ebg, fns ))
##    user  system elapsed 
##  15.653   1.286  17.070
register(MulticoreParam(workers=2))
system.time(i3 <- summarizeOverlaps( ebg, fns ))
##    user  system elapsed 
##   0.103   0.029   8.818
register(MulticoreParam(workers=3))
system.time(i3 <- summarizeOverlaps( ebg, fns ))
##    user  system elapsed 
##   0.051   0.023   9.543
register(MulticoreParam(workers=4))
system.time(i3 <- summarizeOverlaps( ebg, fns ))
##    user  system elapsed 
##   0.056   0.026   8.774
all.equal(i3,i2)  # check that the results do not change
## [1] TRUE

Note that in this environment, despite increasing the number of CPUs linearly do not appreciably reduce run time after moving to two cores. This due mainly to communication costs that typically increase with the number of CPUs, and this phenomenon will be environment-specific.

In summary, it is very easy to perform embarrassingly parallel tasks with R, and this carries over to genomic data analysis thanks to BiocParallel. There are some strategic considerations concerning control of memory consumption and communication costs, and full mastery of the topic involves attention to profiling and benchmarking methods that can be addressed in an advanced software development course.

参考资料

  1. Architecture: considerations on high performance computing with Bioconductor