(本文已于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"] <-