使用NormqPCR包来选择多内参—科研工具箱

请关注公众号【叨客学习资料】 在使用网站的过程中有疑问,请来公众号进行反馈哦

前言

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

1. 合并技术重复(technical replicates)(了解即可);

2. 处理低Cq值(undetermined values)(了解即可);

3. 计算用于均一化的最佳管家基因(housekeeping genes)(这个最重要,其余的很多函数其实了解即可,这一部分要了解geNorm和NormFinder的算法原理);

4. 使用管家基因来对数据进行均一化(normalization)(了解)。

注:在这篇笔记中,管家基因(housekeeping genes)与内参(reference)指的是同一个意思。
由于原文档有点长,有一些东西了解一下就行,如果想要看完整的笔记,可以单击阅读原文直接看我博客的内容。公众号里只放第三部分内容(即如何计算最稳定的内参)。

选择最稳定内参基因

这一部分涉及选择最稳定的内参基因,其中用的是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方法

先加载数据,这里所用的数据都是前文提到的那篇文献里的数据,如下所示:

# BiocManager::install(\"ReadqPCR\")
# BiocManager::install(\"NormqPCR\")
# Following two packages is used to input raw data for qPCR
library(NormqPCR)
library(ReadqPCR)

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

数据结构如下所示:

> str(exprs(geNorm.qPCRBatch))
 num [1:101:850.0425 0.0576 0.1547 0.1096 0.118 ...
 - attr(*, \"dimnames\")=List of 2
  ..

nbsp;: chr [1:10\"ACTB\" \"B2M\" \"GAPD\" \"HMBS\" ...nbsp;: chr [1:85\"BM1\" \"BM2\" \"BM3\" \"BM4\" ...

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

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

第二步,再次计算M值,再排除掉最高M值的内参基因,直到剩下两个内参基因(也就是参数minNrHK=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值,结果如下所示:

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这类组织的计算结果。现在我们把其它的组织都计算出来,如下所示:

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]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

我们自己的计算结果:

ranks <- data.frame(c(11: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

结果如下所示:

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,如下所示:

图片[2]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

代码如下所示:

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.21.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.21.2, by = 0.2), labels = seq(0.21.2, by = 0.2))
box()
abline(h = seq(0.21.2, by = 0.2), lty = 2, lwd = 1, col = \"grey\")
legend(\"topright\", legend = c(\"BM\"\"POOL\"\"FIB\"\"LEU\"\"NB\"),
fill = mypalette)
图片[3]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

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

图片[4]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

代码如下所示:

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(02),
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.050.25, by = 0.05), lty = 2, col = \"grey\")
abline(h = 0.15, lty = 1, col = \"black\")

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

图片[5]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

如何将自己的数据导入分析?


输入数据格式如下:Data

qPCRBatch <- read.qPCR(“Data路径”)
rownames(exprs(qPCRBatch))
res <- selectHKs(qPCRBatch, method = \"geNorm\",Symbols = featureNames(qPCRBatch),minNrHK = 2, log = FALSE)
rank<-res$ranking

获得排名


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,如下所示:

图片[6]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网

第一步,导入数据::

## NormFinder
data(Colon)
Colon

结果如下所示:

> 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

继续计算,如下所示:

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

结果如下所示:

> 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,如下所示:

data(Bladder)
Bladder
grade <- pData(Bladder)[,\"Grade\"]
res.Bladder <- stabMeasureRho(Bladder,
                              group = grade,
                              log = FALSE)
sort(res.Bladder)

计算结果如下所示:

> 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算法来计算,如下所示:

selectHKs(Colon, 
          log = FALSE, 
          trace = FALSE,
          Symbols = featureNames(Colon))$ranking

计算结果如下所示:

> 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的数据,如下所示:

selectHKs(Bladder,
    log = FALSE,
    trace = FALSE,
    Symbols = featureNames(Bladder))$ranking

结果如下所示:

> 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函数中含有这种算法的参数,如下所示:

Class <- pData(Colon)[,\"Classification\"]
selectHKs(Colon, 
    group = Class, 
    log = FALSE, 
    trace = TRUE,
    Symbols = featureNames(Colon), 
    minNrHKs = 12,
    method = \"NormFinder\")$ranking

结果如下所示:

> 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中的内容,如下所示:

图片[7]-使用NormqPCR包来选择多内参—科研工具箱-叨客学习资料网
grade <- pData(Bladder)[,\"Grade\"]
selectHKs(Bladder, 
    group = grade, 
    log = FALSE, 
    trace = FALSE,
    Symbols = featureNames(Bladder), 
    minNrHKs = 13,
    method = \"NormFinder\")$ranking

结果如下所示:

> 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\" 

参考资料

  1. NormqPCR: Functions for normalisation of RT-qPCR data
  2. 2.
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
请输入有效评论哦,肆意灌水或者乱打评论是不会通过的,会影响您评论后获得资源哦~~
提交
头像

昵称

取消
昵称表情

    暂无评论内容