使用BioC软件包分析qPCR数据

（本文已于2014.10.26更新）

BioC中有关qPCR分析的软件包有: HTqPCR, ddCt 和 qpcrNorm。其中HTqPCR软件包可用于Ct值分析，其功能包括：数据读取, 质量分析，归一化处理，数据可视化，feature（基因或miRNA）间Ct值的参数或非参数检验等，并引入limma的方法对多组试验之间的差异进行显著性检验，基本的统计分析功能都具备了。ddCt软件包能使用多内参基因计算基因表达量，和HTqPCR包配合使用基本能满足一般qPCR数据分析需要。

但是，HTqPCR和ddCt这两个软件包能处理的数据格式相对固定，它只能处理少数几种qPCR仪产生的数据。虽然软件包设计的函数具有读取非标准格式数据的功能，但实际应用过程中会有很多问题，对R语言不是很熟悉的用户很可能就卡在了数据读取步骤。因此，本文主要介绍qPCR非标准数据的读取和分析，软件包的其他功能请参考它们的说明文档。

library(BiocInstaller)
biocLite(c("HTqPCR", "ddCt"))

1 使用HTqPCR分析非标准化qPCR格式数据

library("HTqPCR")

1.1 readCtData函数

HTqPCR用于从文件中读取数据的函数为 readCtData ，它的使用语法为：

# 非运行代码
readCtData(files, path = NULL, n.features = 384, format = "plain", column.info,
    flag, feature, type, position, Ct, header = FALSE, SDS = FALSE, n.data = 1,
    samples, na.value = 40, sample.info, ...)

主要参数说明：

• files： 包含Ct值的一系列文本文件
• n.features：基因/探针/引物对的数量，如一个96孔板上用8对引物进行了扩增，n.features=8
• format 参数用于设置读取文件的格式，它预设了三种格式：
  • plain：用于读取ABI Taq-Man Low Density Arrays（TLDA）数据
  • OpenArray：也是ABI的格式
  • SDS：旧格式（没用过，不了解）
• column.info：文件数据列信息定义，用于指定探针所在位置、Ct值、基因名称、基因类型、样品编号、样品处理类型等所在的列
• flag, feature, type, position, Ct：直接对相应的信息进行赋值，它们的作用等同于column.info中的设置
• header（FALSE/TRUE）：文件是否包含表头
• n.data：向量，每个文件中包含的样品数
• na.value：NA值
• …：传递给read.table, read.csv的其他参数

1.2 自定义格式的整理和导入

1.2.1 qPCR数据整理

如果我们的数据格式完全符合HTqPCR预定义的格式，直接设置文件名、路径和格式进行读取就可以了。但除非你用的PCR仪跟他们指定的完全一样，而且Ct导出方法也一样，否则数据格式都会不那么“标准”。

其实“标准”是人定的，了解规则就很好办。readCtData函数无非是从文本文件中读入了下面的信息：

• PCR孔的位置
• Ct值
• 生物学处理，即treatment
• 样品/生物学重复的名称：sample
• 基因名称：gene
• 基因类型，即内参和检测基因，可以用ref和test表示
• Ct状态标记flag，可以用“good”或“bad”表示是否有问题，当然也可以用1/0表示

只要把数据处理成含有这些信息的文本表格就可以了，这很容易用表格处理软件产生，每一列保存一类变量。但Excel专家们要注意了：数据输入从第一个单元格开始，行和列之间别留空行或列，更不要使用单元格合并这些伎俩。

为了演示，我这做一组96孔板的qPCR数据吧：

set.seed(100)
well <- paste0(rep(LETTERS[1:8], each = 12), 1:12)
Ct <- rep(sample(20:32, 96/3, replace = TRUE), each = 3)
Ct <- round(rnorm(96, mean = Ct, sd = 0.5), 2)
treatment <- rep(rep(c("CK", "T1", "T2", "T3"), each = 3), 8)
samplex <- rep(paste0("S", 1:12), 8)
gene <- rep(paste0("Gene", 1:8), each = 12)
typex <- ifelse(grepl("A|D", well), "ref", "test")
## 对参考基因的Ct做一些处理
Ct[gene == "Gene1"] <- round(rnorm(12, mean = 18, sd = 0.1), 2)
Ct[gene == "Gene4"] <-