title: NormqPCR筆記對應文獻之一
date: 2019-10-01 12:00:00
type: "tags"
tags:
- miRNA
- qPCR
- 內(nèi)參
categories: - 生信工具
前言
最近在處理qPCR數(shù)據(jù)過程中，設計的qPCR陣列中采用了多個內(nèi)參，我對于多內(nèi)參的處理不太理解。后來找到了一篇文獻，就是講qPCR中的多內(nèi)參問題的。文獻信息如下：
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），雖然文獻比較古老，但是了解一下文獻中的思想還是比較有用的，以下是文獻的閱讀筆記（基本上就是翻譯原文外加自己檢索）。
文獻中涉及一定的線性代數(shù)知識，我不太懂，只能先暫時先記錄下來，等看懂的時候再補充一篇筆記。
研究背景
? 在生物學研究中，基因表達的分析有著重要的作用。對于基因表達模式的理解可以為復雜的調控網(wǎng)絡研究提供思路。最近開發(fā)的兩種檢測轉錄本豐度的方法得到了非常廣泛的應用（該文章是2002年寫的，這里的最近已經(jīng)是非常久遠了）。微陣列的方法（微陣列其實就是芯片法）可以同時檢測數(shù)以千計的基因表達情況，而RT-PCR的方法則是用于檢測許多不同的樣本中數(shù)量有限基因的表達情況（在細胞數(shù)量比較少的情況下用這種方法更好，說服力更強）。與傳統(tǒng)的方法例如northern-blot，核糖核酸酶保護，或者是競爭RT-PCR（competitive rt-PCR）相比，這兩種方法（微陣列與RT-qPCR）都有快速，高通量的優(yōu)點。但是這兩種方法都需要與傳統(tǒng)方法一樣，需要對目的mRNA進行均一化（normalization）。
以下是補充一些與核酸保護實驗相關的背景知識。
核酸酶保護實驗
? 糖核酸酶保護實驗（Ribonuclease protection assay，RPA）是一種mRNA定量分析方法。其基本原理是將標記的特異RNA探針（32P或生物素）與待測的RNA樣品液相雜交，標記的特異RNA探針按堿基互補的原則與目的基因特異性結合，形成雙鏈RNA；未結合的單鏈RNA經(jīng)RNA酶A或RNA酶T1消化形成寡核糖核酸，而待測目的基因與特異RNA探針結合后形成雙鏈RNA，免受RNA酶的消化，故該方法命名為RNA酶保護實驗。
? 對于32P標記的探針，雜交雙鏈進行變性聚丙酰胺凝膠電泳，用放射自顯影或磷屏成像系統(tǒng)檢測被保護的探針的信號; 對于生物素標記的探針，雜交雙鏈經(jīng)過變性聚丙酰胺凝膠電泳后電轉移至尼龍膜，采用鏈霉親和素－辣根過氧化物酶（Streptavidin-HRP）和化學發(fā)光底物與膜上biotin標記的探針結合，X射線膠片或化學發(fā)光圖像分析儀檢測雜交信號。
? 在利用RT-PCR對基因進行定量時，為了縮小誤差，我們需要對一些因素進行控制，例如原始樣品的數(shù)量（細胞或者是組織），酶的效率，組織或細胞的整個轉錄活動差異等因素。在控制條件的理想情況下，提取的高質量RNA的轉錄本數(shù)量的可重復結果與細胞數(shù)量是呈正比的，但是精確的細胞數(shù)目很難進行定量（尤其是當原始材料是組織時更是如此）。
? 另外一種常用的定量手段就是RNA的量（這里指的是數(shù)量，原文寫的是mass quantity），尤其是northern blot分析方法。但是這種方法也存在著一些爭議。因為這種方法沒有考慮到RNA的質量（這里側重的是RNA的完整性，原文寫的是quality）與相關酶的反應效率。此外，在某些情況下，這些因素也無法定量，例如從提取的微切割的組織（microdissected tissue）中得到的RNA數(shù)量極少，這個時候就無法對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亮就行了），如果這兩個條帶變?nèi)酰?S條件變亮，則說明RNA降解，如下所示：
? 不過電泳法是比較粗糙的方法，畢竟這篇文獻是2002年的，現(xiàn)在如果要進行測序建庫，要求就比較高，會采用Qubit檢測RNA的濃度，用安捷倫2100來檢測RNA的完整性，具體的檢測指標可以檢索相關資料。
但是在總的RNA中，mRNA的含量僅占7.5%。此外，還有報告指出，rRNA的轉錄還受到生物因素與藥物的影響。利用18S或28S作為RNA定量的標準缺陷還在于，在純的mRNA分子中，并不含有18S或28S，與靶mRNA轉錄本相比，18S或28S的含量太高。在利用RT-PCR進行分析時，很難排除基線(baseline)。
現(xiàn)在補充一些競爭性PCR的背景知識。
競爭性PCR
? 在競爭性PCR中，逆轉錄之前加入一個外源的RNA轉錄本（內(nèi)部標準RNA）作為樣品間差異的對照。內(nèi)部標準RNA被逆轉錄并同目的模板一起擴增，作為cDNA轉化和擴增效率差異的對照。通過將特異性的目的序列同已知濃度的內(nèi)部標準RNA一起擴增進行定量。通過比較由內(nèi)部標準獲得的信號和目的模板所獲得的信號可以確定目的模板的豐度。
? 到目前為止（也是2018年為止），內(nèi)參基因是最常用的對mRNA進行均一化的方法。內(nèi)參基因（internal control）通常指的是管家基因（housekeeping gene），這類基因并不會隨著組織或細胞的不同而變化， 或者是不會隨著實驗條件而變化。但是在許多的研究里，研究者使用的這些持續(xù)表達的內(nèi)參基因中，并沒有對這些基因的穩(wěn)定表達進行驗證。文獻中使用的這些管家基因偶爾會出現(xiàn)變化極大的現(xiàn)象。隨著RT-PCR靈敏度，可靠性以及使用范圍的增加，對于合適內(nèi)參基因的需求也日益迫切。些外，通過利用利用公共微陣列數(shù)量，這種均一化因素（normalization factor）也在廣泛的應用的微陣列定標系數(shù)（scaling factor）中得到了驗證。
需要了解的幾個公式
文獻末尾處提到了幾個公式，現(xiàn)在匯總如下所示：
公式1:單一內(nèi)參均一化錯誤值（Single control normalization error）
公式上面的英文注釋如下所示：
E值： 對于任何m個樣本，用real-time RT-PCR技術可以測量其n個內(nèi)參基因的表達水平，其表達水平為aij。對于每2個樣本（分別記為p和q）的組合來說，它們的兩個內(nèi)參基因j和k的每一種組合來說，可以計算出它們的單一內(nèi)參均一化錯誤值E（single control normalization error，公式1）。
當樣本p和q分別針對其內(nèi)參基因j或k進行均一化的時候，其樣本p和樣本q的表達總倍數(shù)差異就是錯誤值E。 n表示內(nèi)參的數(shù)目（注：既然是多內(nèi)參，這個n必然是大于等于2）。
aij表示的是內(nèi)參的表達水平（這是一組數(shù)據(jù)）。
j與k的范圍為1到n，p與q的值為1到m（都是閉區(qū)間），并且j不等于k，p不等于q。Rjkpq的倒數(shù)就是E值，公式1中，aqj是說在樣本q中，以內(nèi)參基因j為標準進行均一化后的數(shù)據(jù)，aqk則是說，以內(nèi)參基因k為標準進行均一化后的數(shù)據(jù)，aqj/aqk則是這二者之比，apk/apj為例，aqj/aqk與apk/apj的比值則是指p與q樣本中，a基因表達水平的誤差。
個人理解：這個E值其實就是每個樣本，每個內(nèi)參均一化后誤差的乘積。
公式2:內(nèi)參基因穩(wěn)定檢測值M
原文如下所示：
公式解釋：
M值：對于任意兩個內(nèi)參基因j與k，可以計算出由它們構成的一個Ajk數(shù)列（array），這個數(shù)列包含一組經(jīng)log2轉化來的比值（aij/aik），如果是m個樣本，那么Ajk中一共是m個值（公式2）。作者定義Vjk是內(nèi)參基因j和k的成對變異（pairwise variation）它是Ajk的標準差（公式3）。其中內(nèi)參基因j的穩(wěn)定表達值Mj就是所有成對基因變異Vjk的算術平均數(shù)（公式4）。
畫個圖，如下所示：
以上計算的只是內(nèi)參基因j的穩(wěn)定值。
研究思路
提取了不同組織或細胞的RNA進行檢測，然后檢測這些RNA中的幾個選定的內(nèi)參，再通過一定的算法進行比較。
實驗過程
（一）管家基因表達譜
? 本文研究了10個常用管家基因在人類13個不同組織中的表達，這10個管家基因如下所示：
管家基因列表。
注：選這個10個管家基因的思路是要盡量使這10個基因屬于不同的功能分類，這樣就能極大地了降低它們之間存在的共調控現(xiàn)象。
作者檢測了80個組織中的這10個管家基因的表達情況，這80個組織分別為：在34個神經(jīng)母細胞瘤（來源不同的實驗室以及不同的病人，NB1-34）；20個短期培養(yǎng)的常規(guī)成纖維細胞（FIB1-20）；13個常規(guī)的淋巴細胞樣本（LEU1-13）；9個骨髓樣本（BM1-9）；9個其他的人類正常組織（分別為heart，brain，fetal brain，lung，trachea，kidney，mammary gland，small intestine和uterus）。
（二）單一內(nèi)參均一化錯誤值（Single control normalization error）
Figure 1 -對于兩個理想的內(nèi)參基因來說，E值等于1，這個很好理解，我們還看一下前面的E公式，如下所示：
也就是說，理想的兩個內(nèi)參基因表達水平必然是相同的，因此E等于1。
事實上觀察到的E值通常大于1，這就構成了兩個樣本之間的E倍差異，具體的差異值取決用于均一化的特定管家基因。對所有的10個內(nèi)參基因兩兩配對（一共45種組合），以及所有樣本的兩兩配對（一共是865個組合）繪制E值，如figure 1所示。此外，通過分析重復運行的相同內(nèi)參基因繪制出了系統(tǒng)錯誤的分位數(shù)。E值的75百分位數(shù)和90百分位分別是3.0（范圍是2.1-3.9）和6.4（范圍是3.0-10.9）。
注：所有樣本的兩兩配對是指，同一個組織內(nèi)的兩兩配對，因此根據(jù)前面的內(nèi)容：
34個神經(jīng)母細胞瘤（來源不同的實驗室以及不同的病人，NB1-34）；20個短期培養(yǎng)的常規(guī)成纖維細胞（FIB1-20）；13個常規(guī)的淋巴細胞樣本（LEU1-13）；9個骨髓樣本（BM1-9）；9個其他的人類正常組織（分別為heart，brain，fetal brain，lung，trachea，kidney，mammary gland，small intestine和uterus）
排列組合數(shù)目就是865。
（三）基因表達穩(wěn)定性檢測以及所選基因的秩
? 在研究中，通常會認為基因的表達水平應該精選一個穩(wěn)定表達的基因進行均一化。但是為了驗證選定的這個內(nèi)參基因的穩(wěn)定表達，我們需要一種可靠的檢測手段手段來確定這個內(nèi)參基因確實是穩(wěn)定表達的，從而剔除掉非特異性的變異。為了解決這個循環(huán)問題（circular problem），作者以非均一化表達水平為基礎，提出了一種基因穩(wěn)定性檢測算法。
? 這種算法的原理是，無論實驗條件，無論細胞類型，兩種理想的內(nèi)參基因在所有樣本中的表達必然是相同的。而在實際條件下，兩種真實內(nèi)參基因的表達比值（ratio）的變異就會反映出一種（或者是兩種）內(nèi)參基因的表達并非是恒定不變的，這種比值的差異與基因表達的穩(wěn)定程度呈反比。對于每一個內(nèi)參基因來說，作者都研究了它與其他所有內(nèi)參基因兩兩的差異，這種差異用表達比值的對數(shù)轉換的標準差來表示，并且定義了一個內(nèi)參基因表達穩(wěn)定檢測值M作為一個特定基因與其他所有內(nèi)參基因的平均兩兩變異。
? 其中，最小M值的基因是表達最穩(wěn)定的內(nèi)參基因。假如選定的這些內(nèi)參基因不存在共調控，那么分步排除掉的最高M值的基因后，最終剩下的內(nèi)參基因就會產(chǎn)生一對組合，構成這種組合的兩個基因是穩(wěn)定表達的管家基因，即它們是樣品中最穩(wěn)定的內(nèi)參基因。
? 為了計算最佳的內(nèi)參基因組合，作者利用VBA開發(fā)了一個Excel插件，命名為geNorm，它能自動計算所有內(nèi)參基因的的M值（現(xiàn)在這個插件已經(jīng)整合到了qbase+中，網(wǎng)址為：https://www.qbaseplus.com/）。這個程序可以剔除那些不好的內(nèi)參基因（即最高M值的基因），并且重新計算剩余內(nèi)參基因的M值。利用geNorm，作者計算了5類組織中的10個內(nèi)參基因的M值，并按從大到小的順序排列（參見figure2與table 3）。此外，作者還計算了系統(tǒng)變異，其方法對相同的基因進行重復的RT-PCR實驗，這種變異反映了儀器自身，加樣，酶的固有變異。
Figure 2 橫坐標表示的是剩余的內(nèi)參數(shù)目。
注：qbase+這個軟件我并不熟悉，這是沒有重現(xiàn)數(shù)據(jù)的分析過程。但是在前一篇筆記中，我使用了NormqPCR這個R包來重現(xiàn)了數(shù)據(jù)分析過程。
（四）基于多內(nèi)參基因的幾何均數(shù)進行均一化因子計算
? 為了精確地評估基因的表達水平，作者推薦使用多個內(nèi)參基因進行均一化。將多個內(nèi)參取幾何均數(shù)。至于使用的內(nèi)參基因數(shù)目，需要在實際情況與精確性方面進行取舍。但是就實際情況來看，選擇太多的內(nèi)參基因浪費（例如10個太浪費了），如果太少，不為精確，作者推薦最少使用3個表達穩(wěn)定的內(nèi)參基因作為均一化因子（normalizaiton factor, NFn, n=3）。
? 作者還研究了一些情況，即那些利用了超過3種內(nèi)參基因用于均一化的情況，作者計算了相同組織類型內(nèi)的所有樣本，按照兩次均一化因子（NFn和NFn+1）的兩兩變異Vn/n+1（即aij=NFn,i和aik=NFn+1,n代表用于均一化的基因數(shù)目，范圍是3到9，i表示樣本，具體參照公式2和3）。我的理解是，在選定了幾個內(nèi)參，例如5個，對所有的樣本進行均一化，計算其變異，即V5，然后再排除一個內(nèi)參，剩了4個內(nèi)參，再對所有的樣本均一化，計算其變異，即V4，最后計算V4/V5）。這種算法研究的是，如果存在著一個很大的變異，就意味著，這個加入的內(nèi)參基因（我的理解是，排除的那個內(nèi)參）有著很強的效應，它應該優(yōu)先選擇用于均一化因子的計算（不太理解這種算法）。
? 對于所有的組織類型，用于計算的均一化因子是3個最穩(wěn)定的內(nèi)參基因，通過逐步添加剩余的7個內(nèi)參基因。逐步計算兩兩配對變異，對于連續(xù)的NFn和NFn+1均一化因子，這就反應了加入的第n+1基因的效應。
（如figure 3a所示）。從圖中可以明顯看出來，對于leukicytes, fibroblast和bone marrow來說，添加的第4個基因沒有特別明顯的影響（就是說低V3/4值），這也說明，在NF3和NF4值之間存在著很強的關聯(lián)，如figure 3b所示。以這些數(shù)據(jù)為基礎，作者采用了0.15作為閾值，低于0.15的內(nèi)參基因不需要添加進來。對于neuroblastoma和正常組織來說，分別需要加入1個和2個內(nèi)參基因（參見figure 3b）。最高的V8/9和V9/10值對于normal pool，neuroblastoma和leukocyte來說就確認了figure 2中逐步排除的最差內(nèi)參基因。這種分析表明了平均M值的在開始時的強烈下降，這就指出了，leukocyte的兩個內(nèi)參基因的表達異常，neuroblastoma和normal tissues一個內(nèi)參基因的不穩(wěn)定。此外，對于后兩種組織來說，需要包含另外的內(nèi)參基因，從而與內(nèi)參基因的高度變異保持一致。
驗證可能的RT-PCR均一化因子
? 為了驗證這種基因檢測算法的可靠性，即最低M值是最穩(wěn)定的表達的內(nèi)參基因，作者檢測了每個內(nèi)參基因的特異性變化，作為均一化之后表達的變異系數(shù)（variation coefficient）。對于合適的內(nèi)參基因來說，這個系數(shù)應該是最小的。作者以三個基因的幾何平均數(shù)為基礎，分別計算了三個不同的均一化因子，最低的M值（NF3（1-3）），最高的M值（NF3（8-10）和中等的M值（NF3（6-8））。接著，作者計算用最低變異系數(shù)，在每種組織類型內(nèi)，計算了這三個基因的平均基因特異性變異（figure 4a）。從計算結果可以看出，當數(shù)據(jù)用NF3（1-3）進行均一化時，這些基因特異性變異最小。
? 這就表明，基因穩(wěn)定性檢測能夠確定地確定表達最穩(wěn)定的內(nèi)參基因。為了確認高M值是一類表達不穩(wěn)定，或者是表達有差異的基因，作者分析了MYCN的表達水平，MYCN這個基因是一種在neuroblastoma中表達有差異的基因。在neuroblastoma中，MYN的M值是6.02，而B2M（表達最穩(wěn)定的內(nèi)參基因）的M值是2.17。進一步的觀察發(fā)現(xiàn)，用單一的內(nèi)參基因進行均一化會導致其它內(nèi)參基因更高的的基因特異性變異，這就進一步說明了利用多個管家基因進行均一化的優(yōu)勢。
? 為了說明最佳的個內(nèi)參基因不受細胞增殖的影響，作者分析了PCNA的表達水平（PCNA是增殖的標志），研究了四個最佳管家基因和標志基因PCNA之間的Spearman秩相關系數(shù)。從分析結果來看，管家基因之間存在著關聯(lián)（p值小于0.001，相關系數(shù)在0.6到0.76之間）。相比之下，PCNA和三個管家基因之間不存在聯(lián)系，而與HPRT1則存在著弱相關（p值為0.024，相關系數(shù)為0.43）。這些數(shù)據(jù)清楚地說明了，最穩(wěn)定的管家基因與細胞的增殖狀態(tài)無關。
? 為了進一步研究選定的內(nèi)參基因幾何均數(shù)的精確性，作者研究了公共數(shù)據(jù)庫中的微陣列數(shù)據(jù)，利用這些數(shù)據(jù)中的內(nèi)參基因的幾何均數(shù)進行了均一化，具體的選擇過程參見文獻，作者選定了5個內(nèi)參基因，其標準是M值小于0.7，其計算結果與微陣列的計算結果類似。
組織特異性管家基因的表達
? 為了比較13所測組織中的管家基因表達水平的異質性。作者剔除了差異最大的4個管家基因（這4個管家基因分別為B2M，RPL13A，ACTB與HMBS，參見Table 2）后，計算了6個內(nèi)參基因的幾何均數(shù)作為“偽內(nèi)參”，然后將這10個內(nèi)參基因再根據(jù)這個“偽內(nèi)參”進行均一化。Figure 5顯示了所有的10個基因的不同豐度類型，表達最高的ACTB基因是最低的HMBS的400倍。雖然內(nèi)參基因的豐度在不同的組織中類似，但還是觀察到了組織特異性表達的差異，例如，B2M的表達水平在leukocytes中比Fetal Brain中高出112倍，而Fibroiblast和Heart tisssue中，ACTB的差異有22倍。與這兩個基因相比，有2個基因的表達比較穩(wěn)定，例如UBC和HPRT1。
討論部分略。
參考資料：
