实验室:基因级RNA-seq差异表达
背景
本实验室来源于:RNA-Seq工作流:基因水平的探索性分析和差异表达，作者:Michael Love, Simon Anders, Wolfgang Huber;马丁·摩根，2015年10月修改。
本实验室将引导您通过端到端RNA-Seq差异表达工作流程，使用DESeq2与其他Bioconductor包。完整的工作流程从FASTQ文件开始，但我们将在读取与参考基因组对齐并计算已知基因重叠的读取之后开始。我们将进行探索性数据分析(EDA)，差异基因表达分析DESeq2，并直观地探索结果。
许多其他的Bioconductor包在基因水平差异表达的统计推断中是重要的，包括Rsubread，刨边机，limma，BaySeq等。
实验数据
本工作流程中使用的数据是用地塞米松(一种具有抗炎作用的合成糖皮质激素类固醇)处理的气道平滑肌细胞的RNA-Seq实验。例如，糖皮质激素被用于哮喘患者，以预防或减轻呼吸道炎症。在实验中，四株原代人气道平滑肌细胞系用1微摩尔地塞米松处理18小时。对于每一种细胞系，我们都有一个处理过的和一个未处理过的样本。实验参考如下:
Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q.“RNA-Seq转录组分析发现CRISPLD2是调节气道平滑肌细胞细胞因子功能的糖皮质激素响应基因。”公共科学图书馆，2014年6月13日;9(6):e99625。PMID:24926665．地理:GSE52778．
准备计数矩阵
作为输入,DESeq2包期望获得的计数数据，例如，从RNA-Seq或其他高通量测序实验，以整数值矩阵的形式。的值。我-第一行和j矩阵的第1列表示有多少个reads被映射到该基因上我在示例j．类似地，对于其他类型的分析，基质的行可能对应于例如结合区域(使用ChIP-Seq)或肽序列(使用定量质谱)。
计数值必须是测序读取的原始计数。这对于DESeq2的统计模型，因为只有实际的计数才能正确地评估测量精度。因此，请不要提供其他数量，如(四舍五入)标准化计数，或覆盖碱基对计数-这只会导致无意义的结果。
我们将在后面的课程中讨论如何将BAM文件中的数据汇总到计数表中。在这里，我们将“直接进入”，从一个准备好的SummarizedExperiment．
从SummarizedExperiment
我们现在用R数据()命令加载已准备好的SummarizedExperiment这是由与上面描述的Himes等人的论文相关的公开测序数据文件生成的。我们用来制造这个物体的步骤和前面几节中所讲的步骤是相同的，只是我们使用了所有的读取和所有的基因。有关用于创建此对象类型的确切步骤的详细信息装饰图案(“气道”)进入你的R会话。
库(气道)数据(“气道”)se <-气道
的信息SummarizedExperiment对象可以通过访问器函数访问。例如，要查看实际数据，即这里的读取计数，我们使用分析()函数。(头()函数将输出限制在前几行。)
(试验(se))
# # SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 # # ENSG00000000003 679 448 873 408 1138 1047 770 # # ENSG00000000005 0 0 0 0 0 0 0 # # ENSG00000000419 467 515 621 365 587 799 417 # # ENSG00000000457 260 211 263 164 245 331 233 55 # # ENSG00000000460 60 40 35 78 63 76 # # ENSG00000000938 0 0 2 0 1 0 0 # # SRR1039521 # # 572 # # ENSG00000000005 ENSG00000000003 0 # # ENSG00000000419 508 60 # # # # 229 # # ENSG00000000460 ENSG00000000457 ENSG00000000938 0
在这个计数矩阵中，每行代表一个Ensembl基因，每列代表一个测序的RNA文库，值给出了每个文库中对应基因的测序reads的原始数。我们还有每个样本的元数据(计数矩阵的列)。如果使用其他软件对读取进行计数，那么在这一步需要检查计数矩阵的列是否与列元数据的行对应。
我们可以快速检查数以百万计的与基因唯一排列的片段。
colSums(化验(se))
## srr1039508 srr1039509 srr1039512 srr1039513 srr1039516 srr1039517 srr1039520 srr1039521 ## 20637971 18809481 25348649 15163415 24448408 30818215 19126151 21164133
假设我们已经构造了aSummarizedExperiment使用前一节中描述的方法之一，我们现在需要确保对象包含关于样本的所有必要信息，也就是说，一个包含count矩阵列的元数据的表存储在colData槽:
colData (se)
运行avgLength实验样本## SRR1039508 SRR1039508 126 SRX384345 srr10395068 # SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346 srr508567 # SRR1039512 GSM1275866 N052611 trt untrt SRR1039512 126 SRX384349 srr508571 # SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350 SRS508572 ##SRR1039516 GSM1275870 N080611不被接收SRR1039516 120 SRX384353 srr508575 ## SRR1039517 GSM1275871 N080611不被接收SRR1039517 126 SRX384354 SRS508576 ## SRR1039520 GSM1275874 n06101011不被接收srr1039579 ## SRR1039521 GSM1275875 n06101011不被接收SRR1039508 SAMN02422669 ## SRR1039509 SAMN02422675 ## SRR1039512 SAMN02422678 ## SRR1039516 SAMN02422682 ## SRR1039517 SAMN02422673 ##生物样本## ### srr1039520 samn02422683 ## srr1039521 samn02422677
在这里，我们看到这个对象已经包含了一个翔实的colData槽-因为我们已经为您准备好了，如气道装饰图案。然而，当你使用你自己的数据时，你必须在这个阶段为实验添加相关的样本/表型信息。我们强烈建议将此信息保存在逗号分隔值(CSV)或制表符分隔值(TSV)文件中，该文件可以从Excel电子表格导出，并将其分配给colData槽，确保行与列相对应SummarizedExperiment．通过使用示例表的一列指定BAM文件，我们确保了这种通信。
请查看rowRanges ()总结实验;这些是计数发生的基因组范围。
rowRanges (se)
## GRanges对象，长度为64102:## $ENSG00000000003## seqnames ranges | exon_id exon_name ## | ## [1] X [99883667, 99884983] - | 667145 ENSE00001459322 ## [2] X [99885756, 99885863] - | 667146 ENSE00000868868 ## [3] X [99887482, 99887565] - | 667147 ENSE00000401072 ## [4] X [99887538, 99887565] - | 667148 ENSE00001849132 ## [5] X [99888402, 99888536] - | 667149 ENSE00003554016 ## # ... ... ... ... . ... ...## [13] x [99890555, 99890743] - | 667156 ense00003512331 ## [14] x [99891188, 99891686] - | 667158 ense00001886883 ## [15] x [99891605, 99891803] - | 667159 ense00001855382 ## [16] x [99891790,99892101] - | 667160 ense00001863395 ## [17] x [99894942, 99894988] - | 667156 ense00003512331 ## [14] x[99891188, 99891686]…## <64101更多元素> ## ------- ## seqinfo: 722个序列(1个循环)来自一个未指定的基因组
让我们看看数据的基本属性，特别是与RNAseq差分表达式分析中已知的重要统计因素相关的属性。
的图书馆的规模是每个样本映射的读的总数。使用colSums ()在分析()数据汇总库大小。
colSums(化验(气管))
## srr1039508 srr1039509 srr1039512 srr1039513 srr1039516 srr1039517 srr1039520 srr1039521 ## 20637971 18809481 25348649 15163415 24448408 30818215 19126151 21164133
- 样本之间的库大小如何变化?
- 为什么在评估差异表达时纳入库的大小是很重要的?
- 虽然这很容易理解，但为什么从统计学的角度来看，简单地按库的总大小进行缩放不能令人满意呢?
- 可以采用哪些不同的方法来估计库的大小?
- 可以采取什么统计方法来考虑图书馆的大小?
- (完成DESeq2工作流后回答)DESeq2采用什么方法来(a)估计库的大小;(b)将库大小纳入差分表达式分析?
使用rowMeans ()在分析()数据和任意一个嘘()或情节(密度())显示的分布平均基因表达在所有基因中。这可能有助于转换数据，并排除在所有样本中计数很少的基因。
均值<- row均值(化验(气道))xlim <-范围(log(1 +均值))plot(密度(log(1 +均值))，xlim=xlim)
Plot(密度(log(1 +表示[表示> 1]))，xlim=xlim)
很明显，在任何样本中没有任何表达的基因都不可能是差异表达的，这与所研究的假设无关。这个标准排除了多少基因?排除这些基因的好处是什么先天的，在评估任何假设之前?
(在完成DESeq2工作流后回答)通过扩展，似乎直观地知道存在一个表达阈值水平，低于该水平的差异基因表达不能被检测到，与所研究的假设无关。DESeq2如何解决这个问题?
MDS图试图表示n维点投影到2维或3维之间的距离。
D <- dist(t(log(1 + assay(气道))))MDS <- cmdscale(D) plot(MDS, pch=20, asp=1, cex=2)
图(mds, pch=20, asp=1, cex=2, col=气道$cell)
Plot (mds, pch=20, asp=1, cex=2, col=气道$dex)
- 计算样本之间的(欧几里得)距离，并使用多维缩放来在二维中表示这些距离。
- 以探索性的方式，颜色点基于细胞系(
气道细胞美元)和实验处理(气道美元敏捷）.
- 解释这些图，并提出这些图如何为后续的统计分析提供信息。
如果计数是泊松分布的，则计数的均值和方差之间存在线性关系。你能以一种直截了当的方式证明事实并非如此吗?
rowVars <- matrixStats::rowVars plot(rowVars(1 +化验(气道))~ rowMeans(1 +化验(气道))，log="xy")
##警告在xy。坐标(x, y, xlabel, ylabel, log): 30633个<= 0的y值从对数图中省略
从SummarizedExperiment来DESeqDataSet
我们将使用DESeq2评估差异表达包。的扩展版本SummarizedExperiment类,叫做DESeqDataSet．很容易从SummarizedExperiment来DESeqDataSet：
库("DESeq2") dds <- DESeqDataSet(se, design = ~ cell + dex)
“设计”参数是一个公式，它表达了每个基因的计数如何取决于基因中的变量colData．记住，你总是可以获得方法参数的信息?,如DESeqDataSet ?．
微分表达式分析
这样做很方便untrt第一层是在敏捷因子，这样默认的log2倍的变化会被计算为处理后的变化(默认R将选择第一个字母级别，记住:计算机不知道该做什么，除非你告诉它们)。这个函数relevel ()实现:
Dds $dex <- relevel(Dds $dex， "untrt")
此外，如果在任意点上有列的子集DESeqDataSet你应该类似地调用droplevels ()如果子集设置导致某些水平的样本为0。
运行管道
最后，我们准备运行差分表达式管道。准备好数据对象后，DESeq2现在可以通过对函数的一次调用来运行分析DESeq ()：
dds <- DESeq(dds)
估算尺寸因素
估算离散度
##基因分散估计
均值-色散关系
最终的离散度估计
装配模型和测试
这个函数将为它执行的各个步骤打印一条消息。手册页中有更详细的描述DESeq ?．简单地说，这些是:大小因素的估计(控制测序实验文库大小的差异)，每个基因离散度的估计，以及拟合广义线性模型。
一个DESeqDataSet返回，其中包含所有拟合信息，下一节描述如何从该对象中提取感兴趣的结果表。
构建结果表
调用结果()没有任何参数将提取估计的log2倍变化和p设计公式中最后一个变量的值。如果这个变量有2个以上的级别，结果()将提取结果表，以进行最后一级与第一级的比较。
(res <- results(dds))
## log2 fold change (MAP): dex trt vs untrt Wald测试p-value:eng00000000003 708.60217 -0.37415246 0.09884435 -3.7852692 0.0001535422 0.001289269 NA NA NA NA NA NA ## eng00000000419 520.29790 0.20206175 0.10974241 1.8412367 0.0655868755 0.197066699 ## eng00000000457 237.16304 0.13834540 0.2614244 0.7937652378 0.913855995 ## eng0000000046057.93263 -0.08445399 0.24990710 -0.3379415 0.7354072485 0.884141561 ## ... ... ... ... ... ... ...## lrg_94 0 na na na na na na ## lrg_96 0 na na na na na na ## lrg_97 0 na na na na na na ## lrg_98 0 na na na na na ## lrg_99 0 na na na na na na na
作为res是一个DataFrame对象时，它携带包含列含义信息的元数据:
mcols (res use.names = TRUE)
## 6行2列的数据帧##类型描述## <字符> <字符> ## baseMean所有样本归一化计数的中间平均值## log2FoldChange结果log2 fold change (MAP): dex trt vs untrt# # lfcSE结果标准错误:dex trt vs untrt# #统计结果Wald统计:dex trt vs untrt# # pvalue结果Wald测试p-value: dex trt vs untrt# # padj结果BH调整p值
第一列，baseMean，是对所有样本进行归一化计数值的平均值，除以大小因子。其余四列是指具体的对比，即比较泰爱泰党水平高于untrt水平为因素变量敏捷．参见帮助页结果()(通过输入结果呢?)，以获取其他对比资料。
列log2FoldChange是效应量估计。它告诉我们，与未处理的样本相比，由于使用地塞米松治疗，基因的表达似乎发生了多大的变化。该值以2为基数的对数刻度报告:例如，1.5的log2倍变化意味着基因的表达增加了一个乘法因子\(2^{1.5} \约2.82\)．
当然，这个估计有一个与之相关的不确定性，可以在列中找到lfcSE， log2倍变化估计的标准误差估计。我们还可以将特定效应量估计的不确定性表示为统计检验的结果。微分表达式检验的目的是检验数据是否提供了足够的证据来得出这个值确实不同于零的结论。DESeq2为每个基因执行a假设检验看是否有足够的证据来判定零假设处理对基因没有影响，并且观察到的处理与对照组之间的差异仅仅是由实验变异性引起的(即，你同样可以在同一处理组的不同样本之间期望的变异性类型)。与统计学中通常的情况一样，这个检验的结果报告为p值，并在列中找到它pvalue．(记住pValue表示在零假设所描述的情况下，看到与观测到的一样强，甚至更强的折叠变化的概率。)
我们还可以用下面这行代码总结结果，它报告了一些附加信息。
总结(res)
## ##出33469非零总读计数##调整p值< 0.1 ## LFC > 0(上):2617,7.8% ## LFC < 0(下):2203,6.6% ##异常值[1]:0,0% ##低计数[2]:15441,46% ##(平均计数< 5)##[1]见'cooksCutoff'参数?results ##[2]见'independentFiltering'参数?results
注意，在FDR水平为10%的地塞米松治疗下，有许多基因表达差异。这是有道理的，因为已知气道的平滑肌细胞对糖皮质激素有反应。然而，有两种方法可以更严格地确定哪些基因被认为是重要的:
- 降低错误发现率阈值
padj在结果表中)
- 函数从0提高log2倍变化阈值
lfcThreshold的观点结果()．看到DESeq2用小插图来演示这个论点的使用。
的子集p值res将NA(“不可用”)。这是DESeq ()该公司报告该基因所有计数为零的方式，因此没有应用测试。此外,p可以赋值NA如果该基因被排除在分析之外，因为它包含一个极端计数异常值。有关更多信息，请参阅小插图的离群值检测部分。
多个测试
高通量生物学的新手经常假设这些阈值p数值较低，比如0.05，就像在其他设置中经常做的那样，是合适的——但事实并非如此。我们简要解释为什么:
有5648个基因含有ap在33469个基因中，值低于0.05，测试成功报告ap值:
(res$pvalue < 0.05, na.rm=TRUE)
## [1] 5648
总和(! is.na (res pvalue美元))
## [1] 33469
现在，暂时假设原假设对所有基因都成立，即没有基因受到地塞米松治疗的影响。那么，根据的定义p价值，我们预计有5%的基因有p值低于0.05。这相当于1673个基因。如果我们只考虑带有a的基因列表p差异表达的值低于0.05，因此该列表应该包含高达1673 / 5648 = 30%的假阳性。
DESeq2使用基准R中描述的Benjamini-Hochberg (BH)调整p.adjust函数;简而言之，该方法为每个基因计算一个调整后的p值，它回答了以下问题:如果一个称为显著性的所有基因p值小于或等于该基因的值p值阈值，将是什么比例的假阳性(错误发现率(在上述计算的意义上)?这些值称为bh调整值p值，在列中给出padj的res对象。
因此，如果我们认为10%的假阳性是可以接受的，我们可以认为所有的基因都有一个调整p下面的值\(10% = 0.1\)是重要的。有多少这样的基因?
(res$padj < 0.1, na.rm=TRUE)
## [1] 4820
我们将结果表划分为这些基因的子集，然后按log2倍变化估计进行排序，以得到下调最强的显著基因。
resSig <-子集(res, padj < 0.1) head(resSig[order(resSig$log2FoldChange)，])
## log2 fold change (MAP): dex trt vs untrt Wald测试p-value:dex trt vs untrt ## DataFrame与6行6列## baseMean log2FoldChange lfcSE stat pvalue padj ## <数字> <数字> <数字> <数字> <数字> <数字> > ## ENSG00000162692 508.17023 -3.449475 0.1767138 -19.520120 7.406383e-85 9.537305e-82 ## ENSG00000105989 333.21469 -2.847364 0.1763098 -16.149771 1.139834e-58 5.871121e- 15 ## ENSG00000214814 243.27698 -2.753559 0.2235537 -12.317215 7.317772e-35 ##ENSG00000267339 26.23357 -2.704476 0.3519722 -7.683776 1.544666e-14 5.924944e-13 ## ENSG00000013293 244.49733 -2.641050 0.1992872 -13.252481 4.365309e-40 8.842448e-38
而且是最强烈的上调的订单()函数以递增的顺序给出指标，所以求递减顺序的一个简单方法是加a-的迹象。或者，你也可以使用参数减少= TRUE．
head(resSig[order(-resSig$log2FoldChange)，])
## log2 fold change (MAP): dex trt vs untrt Wald测试p-value:dex trt vs untrt ## DataFrame与6行6列## baseMean log2FoldChange lfcSE stat pvalue padj ## <数字> <数字> <数字> <数字> <数字> <数字> ## ENSG00000109906 385.07103 4.847164 0.3313657 14.62784 1.866158e-48 5.902298e-46 ## ensg00000159593 67.24305 4.830815 0.3314192 14.57615 3.983670e-48 1.196960e-45 ## ENSG00000163884 561.10717 4.074297 0.2104708 19.35802 1.744626e-83ENSG00000168309 159.52692 3.977191 0.2558532 15.54481 1.725211e-54 7.775524e-52
诊断的情节
一种快速可视化特定基因计数的方法是使用plotCounts ()函数作为参数DESeqDataSet，一个基因名称，以及用来绘制计数图的组。
topGene <- rownames(res)[it .min(res$padj)] data <- plotCounts(dds, gene=topGene, intgroup=c("dex")， returnData=TRUE)
方法还可以制作更多可自定义的图形ggplot ()函数从ggplot2包:
Library (ggplot2) ggplot(data, aes(x=dex, y=count, fill=dex)) + scale_y_log10() + geom_dotplot(binaxis="y"， stackdir="center")
## stat_bindot() '使用' bin = 30 '。用' binwidth '选择更好的值。
“ma图”为两组比较的实验提供了有用的概述。特定比较的log2倍变化被绘制在y轴上，由大小因子标准化的计数的平均值显示在x轴上(“M”表示负，因为对数比等于log - log，“a”表示平均)。
plotMA (res ylim = c (5,5))
每个基因用一个点表示。具有调节的基因\ (p \)低于阈值(此处为默认值0.1)的值以红色显示。的DESeq2包包含了log2倍的先验变化，导致低计数和高可变计数基因的log2倍变化被调节，这可以从图左侧的点分布缩小中看出。该图表明，只有具有较大平均归一化计数的基因包含足够的信息来产生重要的呼叫。
我们也可以在MA图上标记单个点。这里我们使用与()函数为结果对象的选定行绘制圆圈和文本。在与()函数，只有baseMean而且log2FoldChange的选定行的值res使用。
plotMA(res, ylim=c(-5,5)) with(res[topGene，]， {points(baseMean, log2FoldChange, col="dodgerblue"， cex=2, lwd=2) text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")})
一个基因是否被称为显著性不仅取决于它的LFC，还取决于它的组内变异性DESeq2量化为分散．对于强表达的基因，分散度可以理解为变异系数的平方:分散度值为0.01意味着基因的表达倾向于由典型差异\(\根号{0.01}= 10\%\)在同一处理组的样本之间。对于弱基因，泊松噪声是另一种噪声源。
这个函数plotDispEsts ()可视化DESeq2的离散度估计:
plotDispEsts (dds)
黑点是通过分别考虑每个基因的信息得到的每个基因的离散度估计值。除非有很多样本，否则这些值会在真实值附近剧烈波动。因此，我们拟合红色的趋势线，它显示了分散对均值的依赖，然后将每个基因的估计缩小到红线，以获得最终的估计(蓝点)，然后用于假设检验。在主要“云”点上方的蓝色圆圈是具有高基因分散估计的基因，标记为分散离群值。因此，这些估计值并没有向拟合的趋势线收缩。
另一个有用的诊断图是的直方图p值。
Hist (res$pvalue, breaks=20, col="grey50"， border="white")
通过排除数量非常小的基因，这个图变得更加平滑:
hist(res$pvalue[res$baseMean > 1]， break =20, col="grey50"， border="white")
独立的过滤
MA图突出了RNA-Seq数据的一个重要属性。对于弱表达的基因，我们没有机会看到差异表达，因为低读计数受到很高的泊松噪声的影响，任何生物效应都淹没在读计数的不确定性中。我们也可以通过检查小的比率来证明这一点p按平均归一化计数分类的基因的值(例如小于0.01):
#使用分位数函数qs <- c(0, quantile(res$baseMean[res$baseMean > 0]， 0:7/7)) #将基因切割到垃圾箱垃圾箱<- cut(res$baseMean, qs) #重命名的级别的垃圾箱使用中间点水平(垃圾箱)<- paste0("~"，圆(。5*qs[-1] + .5*qs[-length(qs)])) #计算每个bin比率<- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE)) #绘制这些比率barplot(比率，xlab=“平均归一化计数”，ylab=“小p值比率”)
乍一看，过滤掉这些基因似乎没有什么好处。毕竟，测试发现它们并不显著。但是，这些基因对多重测试调整有影响，如果去除这些基因，多重测试调整的性能会提高。通过从FDR程序的输入中删除弱表达基因，我们可以在保留的基因中发现更多重要的基因，从而提高了测试的能力。这种方法被称为独立的过滤．
这个词独立的强调一个重要的警告。只有当过滤准则独立于实际测试统计量时，这种过滤是允许的。否则，过滤将使测试无效，从而使BH过程的假设无效。这就是为什么我们要过滤平均值所有样本:该过滤器对处理组和对照组的样本分配是盲目的，因此是独立的。内部采用独立过滤软件DESeq2来自于genefilter包，其中包含对一篇描述独立过滤的统计基础的论文的参考。
注释:添加基因名
我们的结果表只使用Ensembl基因id，但基因名称可能更有信息。Bioconductor的注释包有助于将各种ID方案相互映射。
我们装载AnnotationDbi包和注释包org.Hs.eg.db：
库(org.Hs.eg.db)
这是用于的有机体注释包(“org”)智人(“Hs”)，组织形式为AnnotationDbi数据库包(“db”)，使用Entrez Gene id(“eg”)作为主键。要获得所有可用键类型的列表，使用:
列(org.Hs.eg.db)
##[19]“accnum”“alias”“ensembl”“ensemblprot”“ensembltrans”“entrezid”##[7]“酶”“证据”“证据all”“genename”“go”“goall”##[13]“ipi”“地图”“omim”“本体”“ontologyall”“路径”##[19]“pfam”“pmid”“prosite”“refseq”“符号”“ucsckg”##[25]“unigene”“uniprot”
res$hgnc_symbol <- unname(mapIds(org. hs . exe .db, rownames(res)， "SYMBOL"， "ENSEMBL")))
## 'select()'返回1:多个键和列之间的映射
res$entrezgene <- unname(mapIds(org. hs . exe .db, rownames(res)， "ENTREZID"， "ENSEMBL"))
## 'select()'返回1:多个键和列之间的映射
现在结果有了想要的外部基因id:
resOrdered <- res[order(res$pvalue)，] head(resOrdered)
## log2 fold change (MAP): dex trt vs untrt Wald测试p-value:dex trt vs untrt ## DataFrame与6行8列## baseMean log2FoldChange lfcSE stat pvalue padj ## <数字> <数字> <数字> <数字> <数字> <数字> ## ENSG00000152583 997.4398 4.313961 0.1721375 25.06114 1.320082e-138 2.379844e-134 ## ENSG00000165995 495.0929 3.186823 0.1281565 24.86665 1.708459e-136 1.540005e-132 ## ENSG00000120129 3409.0294 2.871488 0.1182491 24.28337 2.937761e-1301.324049e-126 ## ENSG00000189221 2341.7673 0.1366746 23.63567 1.975359e -120 ## ENSG00000211445 12285.6151 3.553360 0.1579821 22.49216 4.95359e -112 1.488067e-108 ## hgnc_symbol entrezgene ## <字符> <字符> ## ENSG00000152583 SPARCL1 8404 ## ENSG00000165995 CACNB2 783 ## ENSG00000101347 SAMHD1 25939 ## ENSG00000189221 MAOA 4128 ## ENSG00000211445 GPX3 2878
输出结果
您可以轻松地将结果表保存在CSV文件中，然后可以使用电子表格程序(如Excel)加载该文件。呼唤as.data.frame是否需要转换DataFrame对象(IRanges打包)到data.frame对象可以处理的write.csv．
write.csv (as.data.frame (resOrdered),文件=“results.csv”)
HapMap样品来自一个工厂，按加工日期订购。
基因组拷贝数与mRNA表达之间的相关性在TCGA卵巢癌Affymetrix微阵列数据集中鉴定了38个错误标记的样本。