差异表达基因分析
分析结果输出
library(DESeq) dds<-DESeq(dds) res<-results(dds)
write.table(res,\,sep=\,row.names=TRUE) head(res)
log2foldchange(MAP):conditionCCvsCA Waldtestp-value:conditionCCvsCA DataFramewith6rowsand6columns
baseMeanlog2FoldChangelfcSEstatpvalue
gene1321173.2886810.262679590.20499831.28137422.000623e-01
gene13222.118367-0.052379520.4989589-0.10497769.163936e-01 gene132335.9737010.500545800.30380961.64756419.944215e-02 gene132488.4216610.176776050.24027270.73573094.618945e-01 gene132543.0018280.811431040.29193962.77944865.445127e-03 gene1326662.136259-1.053561050.1752230-6.01268801.824720e-09 padj
gene13213.790396e-01 gene13229.559679e-01 gene13232.337858e-01 gene13246.565731e-01 gene13252.447141e-02 gene13264.520861e-08
?
注： (1)rownames: 基因 ID (2)baseMean:所有样本矫正后的平均 reads 数 (3)log2FoldChange:取 log2 后的表达量差异 (4)pvalue:统计学差异显著性检验指标 (5)padj:校正后的 pvalue, padj 越小,表示基因表达差异越显著
查看整体分析结果
? summary summary(res)
out of 4190 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 595, 14% LFC < 0 (down) : 644, 15% outliers [1] : 0, 0% low counts [2] : 325, 7.8% (mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
MA 图
library(geneplotter)
plotMA(res,main=\,ylim=c(-2,2))
Heatmap 图
sum(res$padj<0.1,na.rm=TRUE)
library(\)
select<-order(rowMeans(counts(dds,normalized=TRUE)),decreasing=TRUE)[1:1000]
nt<-normTransform(dds)# defaults to log2(x+1) log2.norm.counts<-assay(nt)[select,]
df<-as.data.frame(colData(dds)[,c(\,\)]) pdf('heatmap1000.pdf',width=6,height=7)
pheatmap(log2.norm.counts,cluster_rows=TRUE,show_rownames=FALSE, cluster_cols=TRUE,annotation_col=df) dev.off()