原創(chuàng) 小果 生信果?

以GSE70493為例，進(jìn)行R代碼差異分析結(jié)果可視化分析學(xué)習(xí)

代碼如下：

# 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")

#此段代碼是對所有基因進(jìn)行直方圖繪制，hist處是繪制直方圖的命令，col : 柱子的填充色 ，border：子的邊框的顏色，xlab：x軸名稱， ylab：y軸名稱，main：直方圖標(biāo)題名稱。

# summarize test results as "up", "down" or "not expressed"
dT <- decideTests(fit2, adjust.method="fdr", p.value=0.05)

#此段代碼是通過decideTests函數(shù)計算去統(tǒng)計差異基因的表達(dá)情況。

# Venn diagram of results
vennDiagram(dT, circle.col=palette())

#此段代碼是對差異結(jié)果進(jìn)行venn圖繪制，繪圖使用的函數(shù)是vennDiagram函數(shù)

# create Q-Q plot for t-statistic
t.good <- which(!is.na(fit2$F)) # filter out bad probes
qqt(fit2$t[t.good], fit2$df.total[t.good], main="Moderated t statistic")

#此段代碼是對數(shù)據(jù)進(jìn)行t檢驗，Q–Q (quantile-quantile) plot 即繪制散點圖

# volcano plot (log P-value vs log fold change)
colnames(fit2) # list contrast names
ct <- 1??????? # choose contrast of interest
volcanoplot(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值進(jìn)行火山圖繪制，首先列出了組名，然后選擇一點感興趣的組進(jìn)行火山圖的繪制，火山圖的繪圖使用的是volcanoplot函數(shù)。

# MD plot (log fold change vs mean log expression)
# highlight statistically significant (p-adj < 0.05) probes
plotMD(fit2, column=ct, status=dT[,ct], legend=F, pch=20, cex=1)
abline(h=0)

#此段代碼平均差異圖(MD-plot)是對數(shù)折疊變化(差異)與平均對數(shù)值(平均值)的曲線。在此處進(jìn)行強調(diào)統(tǒng)計顯著性(即p-adj 0.05)，繪圖所用函數(shù)是poltMD函數(shù) 。

# General expression data analysis
ex <- exprs(gset)

#此段代碼是提取了差異表達(dá)數(shù)據(jù)，進(jìn)行表達(dá)數(shù)據(jù)的分析

# box-and-whisker plot
dev.new(width=3+ncol(gset)/6, height=5)
ord <- order(gs)? # order samples by group
palette(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()

#此段代碼是進(jìn)行箱線圖的繪制，它是一種用作顯示一組數(shù)據(jù)分散情況資料的統(tǒng)計圖， 使用boxplot函數(shù)進(jìn)行繪制，dev.new打開繪圖模式，dev.off()關(guān)閉繪圖模式

# expression value distribution
par(mar=c(4,4,2,1))
title <- paste ("GSE70493", "/", annotation(gset), " value distribution", sep ="")
plotDensities(ex, group=gs, main=title, legend ="topright")

#此段代碼是統(tǒng)計表達(dá)式值并繪制分布密度圖 ，這里用到了plotDensities函數(shù)，包含表達(dá)矩陣，分組情況，標(biāo)題名稱等信息。

# mean-variance trend, helps to see if precision weights are needed
plotSA(fit2, main="Mean variance trend, GSE70493")

#此段代碼繪制是平均方差趨勢圖，使用plotSA函數(shù)，目的是為了明確是否需要去計算權(quán)重。