GEO2R分析R代码学习之差异分析结果可视化
收录于话题
GEO2R
以GSE70493为例，进行R代码差异分析结果可视化分析学习：
# Visualize and quality control test results.# Build histogram of P-values for all genes. Normal test# assumption is that most genes are not differentially expressed.tT2 <- topTable(fit2, adjust="fdr", sort.by="B", number=Inf)hist(tT2$adj.P.Val, col = "grey", border = "white", xlab = "P-adj",ylab = "Number of genes", main = "P-adj value distribution")
#此段代码是对所有基因进行直方图绘制，hist处是绘制直方图的命令，col : 柱子的填充色 ，border：子的边框的颜色，xlab：x轴名称， ylab：y轴名称，main：直方图标题名称。
# summarize test results as "up", "down" or "not expressed"dT <- decideTests(fit2, adjust.method="fdr", p.value=0.05)
#此段代码是通过decideTests函数计算去统计差异基因的表达情况。
# Venn diagram of resultsvennDiagram(dT, circle.col=palette())
#此段代码是对差异结果进行venn图绘制，绘图使用的函数是vennDiagram函数
# create Q-Q plot for t-statistict.good <- which(!is.na(fit2$F)) # filter out bad probesqqt(fit2$t[t.good], fit2$df.total[t.good], main="Moderated t statistic")
#此段代码是对数据进行t检验，Q–Q (quantile-quantile) plot 即绘制散点图
# volcano plot (log P-value vs log fold change)colnames(fit2) # list contrast namesct <- 1 # choose contrast of interestvolcanoplot(fit2, coef=ct, main=colnames(fit2)[ct], pch=20,highlight=length(which(dT[,ct]!=0)), names=rep('+', nrow(fit2)))
#此段代码是使用log P-value和，log fold change值进行火山图绘制，首先列出了组名，然后选择一点感兴趣的组进行火山图的绘制，火山图的绘图使用的是volcanoplot函数。
# MD plot (log fold change vs mean log expression)# highlight statistically significant (p-adj < 0.05) probesplotMD(fit2, column=ct, status=dT[,ct], legend=F, pch=20, cex=1)abline(h=0)
#此段代码平均差异图(MD-plot)是对数折叠变化(差异)与平均对数值(平均值)的曲线。在此处进行强调统计显著性(即p-adj <0.05)，绘图所用函数是poltMD函数 。
# General expression data analysisex <- exprs(gset)
#此段代码是提取了差异表达数据，进行表达数据的分析
# box-and-whisker plotdev.new(width=3+ncol(gset)/6, height=5)ord <- order(gs) # order samples by grouppalette(c("#1B9E77", "#7570B3", "#E7298A", "#E6AB02", "#D95F02","#66A61E", "#A6761D", "#B32424", "#B324B3", "#666666"))par(mar=c(7,4,2,1))title <- paste ("GSE70493", "/", annotation(gset), sep ="")boxplot(ex[,ord], boxwex=0.6, notch=T, main=title, outline=FALSE, las=2, col=gs[ord])legend("topleft", groups, fill=palette(), bty="n")dev.off()
#此段代码是进行箱线图的绘制，它是一种用作显示一组数据分散情况资料的统计图， 使用boxplot函数进行绘制，dev.new打开绘图模式，dev.offff()关闭绘图模式
# expression value distributionpar(mar=c(4,4,2,1))title <- paste ("GSE70493", "/", annotation(gset), " value distribution", sep ="")plotDensities(ex, group=gs, main=title, legend ="topright")
#此段代码是统计表达式值并绘制分布密度图 ，这里用到了plotDensities函数，包含表达矩阵，分组情况，标题名称等信息。
# mean-variance trend, helps to see if precision weights are neededplotSA(fit2, main="Mean variance trend, GSE70493")
#此段代码绘制是平均方差趋势图，使用plotSA函数，目的是为了明确是否需要去计算权重。
扫码关注我们
shengxinguoer
生信果
生信硬核知识解答
和小果一起学生信