R 中的数据挖掘算法/序列挖掘/DEGSeq
使用 DEGseq
首先，我们必须选择组织数据的方式。有两种方法可以做到这一点：使用基于 MA 图的随机抽样模型方法或使用基于 MA 图的技术重复方法。
基于 MA 图的随机抽样模型方法
MA 图是一种统计分析工具，已用于检测和可视化微阵列数据的强度依赖比率。设 C1 和 C2 表示从两个样本获得的映射到特定基因的读取次数，其中 Ci ~ 二项式(ni, pi)，I = 1,2，其中 ni 表示映射读取的总数，pi 表示读取来自该基因的概率。我们定义 M = log2C1 - log2C2，以及 A = (log2C1 + log2C2)/2。可以证明，在随机抽样假设下，给定 A = a（a 是 A 的观测值）的 M 的条件分布遵循近似正态分布。对于 MA 图上的每个基因，我们对 H0: p1 = p2 与 H1: p1 ≠ p2 进行假设检验。然后，可以基于条件正态分布分配一个 P 值。
基于 MA 图的技术重复方法
尽管据报道测序平台具有低背景噪声，但技术重复仍然可以提供信息，用于质量控制和估计由于不同机器或平台造成的变异。基于 MA 图的另一种方法通过比较数据中的技术重复（如果可用）来估计噪声水平。在这种方法中，首先在两个技术重复的 MA 图上沿 A 轴应用滑动窗口，以估计对应于不同表达水平的随机变异。通过局部回归对强度依赖噪声水平进行平滑估计，并在正态分布假设下转换为以 A 为条件的 M 的局部标准差。然后，局部标准差用于识别两个样本之间基因表达的差异。
多重检验校正
对于上述方法，针对每个基因计算的 P 值会调整为 Q 值，以进行多重检验校正。用户可以设置 P 值或错误发现率 (FDR) 阈值来识别差异表达基因。
选择数据后，您可以应用 R 包 DEGseq。有五种不同的方法来运行该程序。它们是 DEGexp、DEGseq、readGeneExp、samWrapper 和 getGeneExp。
描述
当用户已经拥有基因表达值（例如映射到每个基因的读取次数）时，此函数用于识别差异表达基因。
用法
DEGexp( geneExpFile1, geneCol1=1, expCol1=2, depth1=rep(0, length(expCol1) ), groupLabel1=" geneExpFile2, geneCol2=1, expCol2=2, depth2=rep(0, length(expCol2) ), groupLabel2=" header=TRUE, sep="", method=c("LRT", "CTR", "FET", "MARS", "MATR", "FC"), pValue=1e-3, zScore=4, qValue=1e-3, foldChange=4, thresholdKind=1, outputDir="none", normalMethod=c("none", "loess", "median"), replicate1="none", geneColR1=1, expColR1=2, depthR1=rep(0, length(expColR1) ), replicate2="none", geneColR2=1, expColR2=2, depthR2=rep(0, length(expColR2) ), rawcount=TRUE )
参数
| 参数 | 描述 |
|---|---|
| geneExpFile1 | 包含样本 1（或 method="CTR" 时的重复 1）重复的基因表达值的文件。 |
| geneCol1 | geneExpFile1 中的基因 ID 列。 |
| expCol1 | geneExpFile1 中样本 1（数值向量）重复的表达值列。注意：每列对应于样本 1 的一个重复。 |
| depth1 | 样本 1（数值向量）每个重复的唯一映射到基因组的读取总数，默认值：将映射到所有注释基因的读取总数作为每个重复的深度。 |
| groupLabel1 | 图中组 1 的标签。 |
| GeneExpFile2 | 包含样本 2（或 method="CTR" 时的重复 2）重复的基因表达值的文件。 |
| geneCol2 | geneExpFile2 中的基因 ID 列。 |
| expCol2 | geneExpFile2 中样本 2（数值向量）重复的表达值列。注意：每列对应于样本 2 的一个重复。 |
| depth2 | 样本 2（数值向量）每个重复的唯一映射到基因组的读取总数，默认值：将映射到所有注释基因的读取总数作为每个重复的深度。 |
| groupLabel2 | 图中组 2 的标签。 |
| header | 一个逻辑值，指示 geneExpFile1 和 geneExpFile2 是否包含变量名称作为第一行。 |
| sep | 字段分隔符字符。如果 sep = ""（read.table 的默认值），则分隔符为空格，即一个或多个空格、制表符、换行符或回车符。 |
| method | 识别差异表达基因的方法。可能的方法包括：• "LRT"：似然比检验 (Marioni 等人，2008)，• "CTR"：检查技术重复之间的变异是否可以用随机抽样模型 (Wang 等人，2009) 解释，• "FET"：Fisher 精确检验 (Joshua 等人，2009)，• "MARS"：基于 MA 图的随机抽样模型方法 (Wang 等人，2009)，• "MATR"：基于 MA 图的技术重复方法 (Wang 等人，2009)，• "FC"：MA 图上的倍数变化阈值。 |
| pValue | p 值阈值（对于方法：LRT、FET、MARS、MATR）。仅在 thresholdKind=1 时使用。 |
| zScore | z 分数阈值（对于方法：MARS、MATR）。仅在 thresholdKind=2 时使用。 |
| qValue | q 值阈值（对于方法：LRT、FET、MARS、MATR）。仅在 thresholdKind=3 或 thresholdKind=4 时使用。 |
| thresholdKind | 阈值类型。可能的类型包括：• 1：p 值阈值，• 2：z 分数阈值，• 3：q 值阈值 (Benjamini 等人，1995)，• 4：q 值阈值 (Storey 等人，2003)。 |
| foldChange | MA 图上的倍数变化阈值（对于方法：FC）。 |
| outputDir | 输出目录。 |
| normalMethod | 归一化方法："none"、"loess"、"median"。建议："none"。 |
| replicate1 | 包含重复批次 1（仅在 method="MATR" 时使用）的基因表达值的文件。注意：重复 1 和重复 2 是样本的两个（组）技术重复。 |
| GeneColR1 | 重复批次 1（仅在 method="MATR" 时使用）的表达文件中的基因 ID 列。 |
| expColR1 | 重复批次 1（数值向量）（仅在 method="MATR" 时使用）的表达文件中的表达值列。 |
| depthR1 | 重复批次 1（数值向量）中每个重复的唯一映射到基因组的读取总数，默认值：将映射到所有注释基因的读取总数作为每个重复的深度（仅在 method="MATR" 时使用）。 |
| ReplicateLabel1 | 图中重复批次 1 的标签（仅在 method="MATR" 时使用）。 |
| replicate2 | 包含重复批次 2（仅在 method="MATR" 时使用）的基因表达值的文件。注意：重复 1 和重复 2 是样本的两个（组）技术重复。 |
| geneColR2 | 重复批次 2（仅在 method="MATR" 时使用）的表达文件中的基因 ID 列。 |
| expColR2 | 重复批次 2（数值向量）（仅在 method="MATR" 时使用）的表达文件中的表达值列。 |
| depthR2 | 重复批次 2（数值向量）中每个重复的唯一映射到基因组的读取总数，默认值：将映射到所有注释基因的读取总数作为每个重复的深度（仅在 method="MATR" 时使用）。 |
| ReplicateLabel2 | 图中重复批次 2 的标签（仅在 method="MATR" 时使用）。 |
| rawCount | 一个逻辑值，指示基因表达值是基于原始读取计数还是归一化值。 |
示例
> library(DEGseq)
> geneExpFile <- system.file("data", "GeneExpExample5000.txt",
+ package = "DEGseq")
> if (geneExpFile == "") {
+ zipFile <- system.file("data", "Rdata.zip", package = "DEGseq")
+ if (zipFile != "") {
+ unzip(zipFile, "GeneExpExample5000.txt", exdir = tempdir())
+ geneExpFile <- file.path(tempdir(), "GeneExpExample5000.txt")
+ }
+ }
> layout(matrix(c(1, 2, 3, 4, 5, 6), 3, 2, byrow = TRUE))
> par(mar = c(2, 2, 2, 2))
> DEGexp(geneExpFile1 = geneExpFile, expCol1 = c(7, 9, 12, 15, 18), groupLabel1 = "kidney", geneExpFile2 = geneExpFile, expCol2 = c(8, 10, 11, 13, 16), groupLabel2 = "liver", method = "MARS")
Please wait...
geneExpFile1: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/GeneExpExample5000.txt
gene id column in geneExpFile1: 1
expression value column(s) in geneExpFile1: 7 9 12 15 18
total number of reads uniquely mapped to genome obtained from sample1: 345504 354981 334557
geneExpFile2: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/GeneExpExample5000.txt
gene id column in geneExpFile2: 1
expression value column(s) in geneExpFile2: 8 10 11 13 16
total number of reads uniquely mapped to genome obtained from sample2: 274430 274486 264999
method to identify differentially expressed genes: MARS
pValue threshold: 0.001
output directory: none
The outputDir is not specified!
Please wait …
Identifying differentially expressed genes ...
Please wait patiently ...
output …
The results can be observed in directory: none
描述
此函数用于识别来自 RNA 测序数据的差异表达基因。它将来自 RNA 测序数据的两个样本的唯一映射读取与基因注释作为输入。因此，用户应提前将读取（从样本的测序文库获得）映射到相应的基因组。
用法
DEGseq ( mapResultBatch1, mapResultBatch2, fileFormat="bed", readLength=32, strandInfo=FALSE, refFlat, groupLabel1="group1", groupLabel2="group2", method=c("LRT", "CTR", "FET", "MARS", "MATR", "FC"), pValue=1e-3, zScore=4, qValue=1e-3, foldChange=4, thresholdKind=1, outputDir="none", normalMethod=c("none", "loess", "median"), depthKind=1, replicate1="none", replicate2="none", replicateLabel1="replicate1", replicateLabel2="replicate2" )
参数
| 参数 | 描述 |
|---|---|
| mapResultBatch1 | 样本 1（或 method="CTR" 时的重复 1）技术重复的唯一映射结果文件。 |
| MapResultBatch2 | 样本 2（或 method="CTR" 时的重复 2）技术重复的唯一映射结果文件。 |
| fileFormat | 文件格式："bed" 或 "eland"。 "bed" 格式示例：chr12 7 38 readID 2 "eland" 格式示例：readID chr12.fa 7 U2 F 注意：字段分隔符字符为 TAB。并且文件必须遵循其中一个示例的格式。 |
| ReadLength | 读取的长度（仅在 fileFormat="eland" 时使用）。strandInfo 在 cDNA 克隆过程中是否保留了链信息。 "TRUE"：保留，"FALSE"：未保留。 |
| RefFlat | UCSC refFlat 格式的基因注释文件。 |
| GroupLabel1 | 图中组 1 的标签。 |
| GroupLabel2 | 图中组 2 的标签。 |
| Method | 识别差异表达基因的方法。可能的方法包括： "LRT"：似然比检验，"CTR"：检查两个技术重复之间的变异是否可以用随机抽样模型解释，"FET"：Fisher 精确检验，"MARS"：基于 MA 图的随机抽样模型方法，"MATR"：基于 MA 图的技术重复方法，"FC"：MA 图上的倍数变化阈值。 |
| pValue | p 值阈值（对于方法：LRT、FET、MARS、MATR）。仅在 thresholdKind=1 时使用。 |
| zScore | z 分数阈值（对于方法：MARS、MATR）。仅在 thresholdKind=2 时使用。 |
| qValue | q 值阈值（对于方法：LRT、FET、MARS、MATR）。仅在 thresholdKind=3 或 thresholdKind=4 时使用。 |
| ThresholdKind | 阈值类型。可能的类型包括：• 1：p 值阈值，• 2：z 分数阈值，• 3：q 值阈值，• 4：q 值阈值。 |
| foldChange | MA 图上的倍数变化阈值（对于方法：FC）。 |
| outputDir | 输出目录。 |
| normalMethod | 归一化方法："none"、"loess"、"median"。建议："none"。 |
| DepthKind | 1- 以每个复制样本中唯一比对到基因组的总读取数作为深度，0: 以每个复制样本中唯一比对到所有注释基因的总读取数作为深度。 我们建议使用 depthKind=1，尤其是在注释文件中的基因是所有基因的一部分时。 |
| replicate1 | 包含来自复制批次 1 的唯一比对读取的文件（仅在 method="MATR" 时使用）。 |
| replicate2 | 包含来自复制批次 2 的唯一比对读取的文件（仅在 method="MATR" 时使用）。 |
| ReplicateLabel1 | 图中重复批次 1 的标签（仅在 method="MATR" 时使用）。 |
| ReplicateLabel2 | 图中重复批次 2 的标签（仅在 method="MATR" 时使用）。 |
示例
> kidneyR1L1 <- system.file("data", "kidneyChr21.bed.txt", package = "DEGseq") > liverR1L2 <- system.file("data", "liverChr21.bed.txt", package = "DEGseq") > refFlat <- system.file("data", "refFlatChr21.txt", package = "DEGseq") > mapResultBatch1 <- c(kidneyR1L1) > mapResultBatch2 <- c(liverR1L2) > outputDir <- file.path(tempdir(), "DEGseqExample") > DEGseq(mapResultBatch1, mapResultBatch2, fileFormat = "bed",refFlat = refFlat, outputDir = outputDir, method = "LRT")
请稍候...
mapResultBatch1: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt mapResultBatch2: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/liverChr21.bed.txt 文件格式: bed refFlat: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/refFlatChr21.txt 在计算比对到基因的读取数时，忽略链信息！ 正在计算比对到每个基因的读取数 ... 这将需要几分钟，请耐心等待！
请稍候...
样本文件: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt 正在计算比对到每个基因的读取数。这将需要几分钟。
请稍候 …
总共 259 个独特基因已处理 0 个读取 (kidneyChr21.bed.txt) 已处理 10000 个读取 (kidneyChr21.bed.txt) 已处理 20000 个读取 (kidneyChr21.bed.txt) 已处理 30000 个读取 (kidneyChr21.bed.txt) 已处理 34304 个读取 (kidneyChr21.bed.txt) 总共耗时 0.328000 秒！
请稍候...
样本文件: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/liverChr21.bed.txt 正在计算比对到每个基因的读取数。这将需要几分钟。
请稍候 …
总共 259 个独特基因已处理 0 个读取 (liverChr21.bed.txt) 已处理 10000 个读取 (liverChr21.bed.txt) 已处理 20000 个读取 (liverChr21.bed.txt) 已处理 30000 个读取 (liverChr21.bed.txt) 已处理 30729 个读取 (liverChr21.bed.txt) 总共耗时 0.297000 秒！
请稍候...
geneExpFile1: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample/group1.exp geneExpFile1 中的基因 ID 列: 1 geneExpFile1 中的表达值列: 2 从样本 1 获取的唯一比对到基因组的总读取数: 34304 geneExpFile2: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample/group2.exp geneExpFile2 中的基因 ID 列: 1 geneExpFile2 中的表达值列: 2 从样本 2 获取的唯一比对到基因组的总读取数: 30729 用于识别差异表达基因的方法: LRT p 值阈值: 0.001 输出目录: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGseqExample
请稍候 …
正在识别差异表达基因 ... 请耐心等待 ... 输出 ...
已完成 ... 结果可在以下目录中查看: C:\DOCUME~1\wanglk\LOCALS~1\Temp\RtmpUvf7Fw/DEGs
getGeneExp
[edit | edit source]描述
此函数用于计算每个基因的读取数并计算 RPKM。它以来自样本的 RNA-seq 数据的唯一比对读取以及基因注释文件作为输入。因此，用户应提前将读取 (从样本的测序文库获得) 比对到相应的基因组。
用法
getGeneExp( mapResultBatch, fileFormat="bed", readLength=32, strandInfo=FALSE, refFlat, output=paste(mapResultBatch[1],".exp",sep=""), min.overlapPercent= 1 )
参数
| 参数 | 描述 |
|---|---|
| mapResultBatch | 包含样本的唯一比对结果文件的向量。注意: 样本可以具有多个技术重复。 |
| fileFormat | 文件格式："bed" 或 "eland"。 "bed" 格式示例：chr12 7 38 readID 2 "eland" 格式示例：readID chr12.fa 7 U2 F 注意：字段分隔符字符为 TAB。并且文件必须遵循其中一个示例的格式。 |
| readLength | 读取的长度 (仅在 fileFormat="eland" 时使用)。 strandInfo 链信息在 cDNA 克隆过程中是否保留。 • "TRUE" : 保留， • "FALSE": 未保留。 |
| refFlat | UCSC refFlat 格式的基因注释文件。 |
| output | 输出文件。 |
| min.overlapPercent | 读取与外显子重叠长度与读取自身长度的最小百分比，用于计算外显子的读取数。应 <=1。0: 读取与外显子之间至少重叠 1 个碱基。 |
示例
> kidneyR1L1 <- system.file("data", "kidneyChr21.bed.txt", package = "DEGseq") > refFlat <- system.file("data", "refFlatChr21.txt", package = "DEGseq") > mapResultBatch <- c(kidneyR1L1) > output <- file.path(tempdir(), "kidneyChr21.bed.exp") > getGeneExp(mapResultBatch, refFlat = refFlat, output = output) 请稍候... 样本文件: D:/myrpackage/DEGseq.Rcheck/DEGseq/data/kidneyChr21.bed.txt 正在计算比对到每个基因的读取数。这将需要几分钟。请稍候 ... 总共 259 个独特基因已处理 0 个读取 (kidneyChr21.bed.txt) 已处理 10000 个读取 (kidneyChr21.bed.txt) 已处理 20000 个读取 (kidneyChr21.bed.txt) 已处理 30000 个读取 (kidneyChr21.bed.txt) 已处理 34304 个读取 (kidneyChr21.bed.txt) 总共耗时 0.328000 秒！ > exp <- readGeneExp(file = output, geneCol = 1, valCol = c(2, + 3), label = c("raw count", "RPKM")) > exp[30:32, ] raw count RPKM C21orf131 0 0.000 C21orf15 0 0.000 C21orf2 51 665.789
readGeneExp
[edit | edit source]描述
此方法用于将基因表达值从文件读取到 R 工作空间中的矩阵。以便矩阵可以用作其他软件包的输入，例如 edgeR。该方法的输入是一个包含基因表达值的文件。
用法
readGeneExp(file, geneCol=1, valCol=2, label = NULL, header=TRUE, sep="")
参数
参数描述 file 包含基因表达值的文件。 GeneCol 文件中的基因 ID 列。 valCol 文件中要读取的表达值列。 label 列的标签。 Header 一个逻辑值，指示文件是否将其第一行作为变量名称。 sep 字段分隔符字符。如果 sep = "" (read.table 的默认值)，则分隔符是空白，即一个或多个空格、制表符、换行符或回车符。
示例
> geneExpFile <- system.file("data", "GeneExpExample1000.txt", + package = "DEGseq") > exp <- readGeneExp(file = geneExpFile, geneCol = 1, valCol = c(7, + 9, 12, 15, 18, 8, 10, 11, 13, 16)) > exp[30:32, ]
R1L1Kidney
R1L3Kidney
R1L7Kidney
R2L2Kidney
R2L6Kidney ENSG00000188976 73 77 68 70 82 ENSG00000187961 15 15 13 12 15 ENSG00000187583 1 1 3 0 3
R1L2Liver
R1L4Liver
R1L6Liver
R1L8Liver
R2L3Liver ENSG00000188976 34 56 45 55 42 ENSG00000187961 8 13 11 12 20 ENSG00000187583 0 1 0 0 2
samWrapper
[edit | edit source]描述
此函数是 samr 中函数的包装器。它用于识别两组具有多个重复的样本或来自不同个体 (例如，疾病样本与对照样本) 的两组样本的差异表达基因。
用法
samWrapper( geneExpFile1, geneCol1=1, expCol1=2, measure1=rep(1, length(expCol1) ), geneExpFile2, geneCol2=1, expCol2=2, measure2=rep(2, length(expCol2) ), header=TRUE, sep="", paired=FALSE, s0=NULL, s0.perc=NULL, nperms=100, testStatistic= c("standard","wilcoxon"), max.qValue=1e-3, in.foldchange=logged2=FALSE, output )
参数
参数描述 geneExpFile1 包含组 1 的基因表达值的文件。 geneCol1 geneExpFile1 中的基因 ID 列。 expCol1 geneExpFile1 中的表达值列。参见示例。 measure1
numeric vector of outcome measurements for group1. like c(1,1,1...) when paired=FALSE, or like c(-1,-2,-3,...) when paired=TRUE.
geneExpFile2 包含组 2 的基因表达值的文件。 geneCol2 geneExpFile2 中的基因 ID 列。 ExpCol2 geneExpFile2 中的表达值列。参见示例。 Measure2 组 2 的结果测量数值向量。如果 paired=FALSE，则类似 c(2,2,2...)，如果 paired=TRUE，则类似 c(1,2,3,...)。 header 一个逻辑值，指示 geneExpFile1 和 geneExpFile2 是否将其第一行作为变量名称。 sep 字段分隔符字符。如果 sep = "" (read.table 的默认值)，则分隔符是空白，即一个或多个空格、制表符、换行符或回车符。 paired 一个逻辑值，指示样本是否配对。 s0 测试统计量分母的可交换性因子；默认是自动选择。 s0.perc 用于 s0 的标准偏差值的百分位数；默认是自动选择；-1 表示 s0=0 (不同于 s0.perc=0，表示 s0=标准偏差值的第零百分位数= sd 值的最小值。 Nperms 用于估计错误发现率的置换次数。 TestStatistic 在两类非配对情况下使用的测试统计量。 "standard" (t 统计量) 或 "wilcoxon" (双样本 Wilcoxon 或 Mann-Whitney 检验)。建议使用 "standard"。 max.qValue 所需的最大 q 值；应 <1。 min.foldchange 所需的最小倍数变化；应 >1。默认值为零，表示不应用倍数变化标准。 logged2 一个逻辑值，指示表达值是否为 log2。 output 输出文件。
示例
> geneExpFile <- system.file("data", "GeneExpExample1000.txt", + package = "DEGseq") > set.seed(100) > geneExpFile1 <- geneExpFile > geneExpFile2 <- geneExpFile > output <- file.path(tempdir(), "samWrapperOut.txt") > expCol1 = c(7, 9, 12, 15, 18) > expCol2 = c(8, 10, 11, 13, 16) > measure1 = c(-1, -2, -3, -4, -5) > measure2 = c(1, 2, 3, 4, 5) > samWrapper(geneExpFile1 = geneExpFile1, geneCol1 = 1, expCol1 = expCol1, measure1 = measure1, geneExpFile2 = geneExpFile2, geneCol2 = 1, expCol2 = expCol2, measure2 = measure2, nperms = 100, min.foldchange = 2, max.qValue = 1e-04, output = output, paired = TRUE)
参考文献
[edit | edit source]Jiang,H. and Wong,W.H. (2009) Statistical inferences for isoform expression in RNASeq. Bioinformatics, 25, 1026–1032.
Likun Wang, Zhixing Feng, Xi Wang, Xiaowo Wang, and Xuegong Zhang. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data Bioinformatics Advance Access published on October 24, 2009, DOI 10.1093/bioinformatics/btp612.
王立昆; 冯志兴; 王曦; 王晓沃; 张雪峰. DEGseq 手册. 清华大学. 北京, 中国. 2009(B). 可从 <http://www.bioconductor.org/packages/devel/bioc/manuals/DEGseq/man/DEGseq.pdf> 访问, 访问日期: 2009 年 11 月 18 日
王立昆; 王曦. 如何使用 Degseq 包. 清华大学自动化系/ TNLIST 生物信息学实验室和生物信息学部门. 北京, 中国. 2009(A). 可从 <http://bioinfo.au.tsinghua.edu.cn/software/degseq/DEGseq.pdf> 访问, 访问日期: 2009 年 11 月 18 日