ggplot2绘制富集分析柱状图

library("ggplot2")
pathway=read.table("COS3_kegg_total.xls",header=T,sep="\t")
pathbar = ggplot(pathway,aes(x=Pathway,y=-1*log10(Pvalue)))
pathbar + geom_bar(stat="identity")
pathbar + geom_bar(stat="identity") + coord_flip()
pathbar+geom_bar(stat="identity",color="red",fill="blue")+coord_flip()
pathbar+geom_bar(stat="identity",aes(fill=-1*log10(Pvalue)))+coord_flip()
porder=factor(as.integer(rownames(pathway)),labels=pathway$Pathway)
pathbar+geom_bar(stat="identity",aes(x=porder,fill=-1*log10(Pvalue)))+coord_flip()
pathbar+geom_bar(stat="identity",aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip()
porder=factor(as.integer(rownames(pathway)),labels=rev(pathway$Pathway))
pathbar+geom_bar(stat="identity",aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip()
pqbar=pathbar+geom_bar(stat="identity",aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip()+theme(legend.position="none")+labs(title="Top20 of pathway enrichment",y=expression(-log[10](Pvalue)))
pqbar
pqbar+geom_hline(yintercept=2,color=c("red"),linetype=4)

pie图绘制

x <- c(344,1619,1002,74,882,117,1292,58,68,158,73,303,139,579,92,93,2110,997)
label <-c("Black","Blue","Brown","Cyan","Green","Greenyellow","Grey","Grey60","Lightcyan","Magenta","Midnightblue","Pink","Purple","Red","Salmon","Tan","Turquoise","Yellow")
pie(x, labels = label,radius =1.8,col =c("black","blue","brown","cyan","green","greenyellow","grey","grey60","lightcyan","magenta","midnightblue","pink","purple","red","salmon","tan","turquoise","yellow"),clockwise = TRUE)
percent<-round(100*x/sum(x),2)
percent <-paste(percent, "%", sep = "")
pie(x, labels = percent,radius =1.8,col =c("black","blue","brown","cyan","green","greenyellow","grey","grey60","lightcyan","magenta","midnightblue","pink","purple","red","salmon","tan","turquoise","yellow"),clockwise = TRUE)
pie(x, labels = percent,main="饼图",radius =2.0,col =c("black","blue","brown","cyan","green","greenyellow","grey","grey60","lightcyan","magenta","midnightblue","pink","purple","red","salmon","tan","turquoise","yellow"),clockwise = TRUE)
pie(x, labels = percent,main="饼图",radius =2.0,col =c("black","blue","brown","cyan","green","greenyellow","grey","grey60","lightcyan","magenta","midnightblue","pink","purple","red","salmon","tan","turquoise","yellow"),clockwise = TRUE)
legend(locator(1),label, cex=1.2, fill=c("black","blue","brown","cyan","green","greenyellow","grey","grey60","lightcyan","magenta","midnightblue","pink","purple","red","salmon","tan","turquoise","yellow"))

WGCNA(my project)

setwd("D:/Ma_transcription/WGCNA/")
database <- read.table(file = "COS_deseq_counts_normalized.txt", sep = "\t", header = T, row.names = 1, stringsAsFactors = F)
library(reshape2)
library('WGCNA')
enableWGCNAThreads()#打开多线程
WGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],])
subname=sapply(colnames(database),function(x) strsplit(x,"_")[[1]][1])
datTraits = data.frame(gsm=names(database),
subtype=subname)
rownames(datTraits)=datTraits[,1]
head(datTraits)
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
datExpr <- WGCNA_matrix
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="green")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
sft$powerEstimate
net = blockwiseModules(datExpr, power = sft$powerEstimate,
                       maxBlockSize = 5000,TOMType = "unsigned", 
                       minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)

# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
table(mergedColors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
## assign all of the gene to their corresponding module 
## hclust for the genes.
#明确样本数和基因数
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#首先针对样本做个系统聚类树
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)
## 如果这个时候样本是有性状,或者临床表型的,可以加进去看看是否聚类合理
#针对前面构造的样品矩阵添加对应颜色
sample_colors <- numbers2colors(as.numeric(factor(datTraits$subtype)), 
                                colors = c("grey","blue","red","green"),signed = FALSE)
## 这个给样品添加对应颜色的代码需要自行修改以适应自己的数据分析项目。
#  sample_colors <- numbers2colors( datTraits ,signed = FALSE)
## 如果样品有多种分类情况,而且 datTraits 里面都是分类信息,那么可以直接用上面代码,
##当然,这样给的颜色不明显,意义不大。
#构造10个样品的系统聚类树及性状热图
par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
                    groupLabels = colnames(sample),
                    cex.dendroLabels = 0.8,
                    marAll = c(1, 4, 3, 1),
                    cex.rowText = 0.01,
design=model.matrix(~0+ datTraits$subtype)
colnames(design)=levels(datTraits$subtype)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同颜色的模块的ME值矩阵(样本vs模块)
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = greenWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))                 main = "Sample dendrogram and trait heatmap")
head(design)
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
## 算出每个模块跟基因的皮尔森相关系数矩阵
## MEs是每个模块在每个样本里面的值
## datExpr是每个基因在每个样本的表达量
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
## 只有连续型性状才能只有计算
## 这里把是否属于 CB 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
geneTraitSignificance = as.data.frame(cor(datExpr, CB, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(CB), sep="");
names(GSPvalue) = paste("p.GS.", names(CB), sep="")
module = "turquoise"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for CB",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

#首先针对所有基因画热图
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
geneTree = net$dendrograms[[1]]; 
#生成全基因不相似TOM矩阵
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6); 
plotTOM = dissTOM^7; 
diag(plotTOM) = NA; 
#TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

#然后随机选取部分基因作图
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

#最后画模块和性状的关系
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有连续型性状才能只有计算
## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, CB))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8,                         xLabelsAngle= 90)
# Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
## 模块的聚类图
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
                      plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
## 性状与模块热图
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
                      plotDendrograms = FALSE, xLabelsAngle = 90)
# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = 6); 
# Select module
module = "blue";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
nTop = 200;
IMConn = softConnectivity(datExpr[, modProbes]);
top = (rank(-IMConn) <= nTop)
vis = exportNetworkToCytoscape(modTOM[top, top], edgeFile = paste("CytoscapeInput-edges-top200", paste(module, collapse="-"), ".txt", sep=""), nodeFile = paste("CytoscapeInput-nodes-top200", paste(module, collapse="-"), ".txt", sep=""), weighted = TRUE, threshold = 0);

参考:lncRNA实战项目-第六步-WGCNA相关性分析

 

MA图绘制

# 导入ggplot2包
library(ggplot2)

# 设置好工作目录(到数据所在目录)
# 读取输入数据  "R0-vs-R3.isoforms.filter.tsv"
data = read.table("R0-vs-R3.isoforms.filter.tsv",header=T,row.names=1)

# 计算M值和A值,并将M作为y轴,A作为x轴
aes = aes(x=(log2(R0_fpkm)+log2(R3_fpkm))/2,y=log2(R0_fpkm)-log2(R3_fpkm))

# 绘制MAplot
ggplot(data=data,aes) + geom_point(aes(color=significant))

# 添加辅助线
ggplot(data=data,aes) + geom_hline(yintersect=0,linetype=4,color="blue") + 
geom_point(aes(color=significant))

# 改变点的大小
maplot = ggplot(data=data,aes) + geom_hline(yintersect=0,linetype=4,color="blue") + 
geom_point(aes(color=significant),size=1)

# 设置自定义染色
maplot + scale_color_manual(values=c("green","black","red"))

# 设置标题
maplot + scale_color_manual(values=c("green","black","red")) + 
labs(title="MAplot of R0-vs-R3",x="A",y="M")

# 保存MAplot
ggsave("R0-vs-R3.MAplot.png",width=8,height=6)

DESeq2使用流程

library("DESeq2")
database <- read.table(file = "COS_genes.count_table.xls", sep = "\t", header = T, row.names = 1)
countData <- database[,1:15]
condition <- factor(c(rep("COS0",3),rep("COS1",3),rep("COS3",3),rep("COS6",3),rep("COS12",3)), levels = c("COS0", "COS1","COS3","COS6","COS12"))
countData <- round(as.matrix(countData))
coldata <- data.frame(row.names = colnames(countData), condition)
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
dds <- dds[ rowSums(counts(dds)) > 1, ]
dds <- DESeq(dds)
res0vs1 <- results(dds, contrast = c("condition","COS1","COS0"))
table(res0vs1$pvalue <0.001)
res0vs1 <- res0vs1[order(res0vs1$pvalue),]
write.table(as.data.frame(res0vs1), file="res0vs1.xls", sep="\t", quote = F)

PCA图绘制

输入表达矩阵数据文件

library(ggplot2)
library(gmodels)
inname = "COS_deseq_counts_normalized.txt"
outname = "COS_PCA.png"
group <- factor(c(rep("COS0",3),rep("COS1",3),rep("COS3",3),rep("COS6",3),rep("COS12",3)), levels = c("COS0", "COS1","COS3","COS6","COS12"))
## step 1: 数据的读取和处理
# read the expr data
expr <- read.table(inname, header=T, row.names=1)
# transpose the data
data <- t(expr)
## step2:PCA分析
# do PCA 
data.pca <- fast.prcomp(data)
## step3: PCA结果解析
# fetch the proportion of PC1 and PC2
# 一般情况下PC1 + PC2 > 70% 二维PCA散点图才有效
a <- summary(data.pca)
tmp <- a[4]$importance
pro1 <- as.numeric(sprintf("%.3f",tmp[2,1]))*100
pro2 <- as.numeric(sprintf("%.3f",tmp[2,2]))*100

# 将成分矩阵转换为数据框
pc = as.data.frame(a$x)

# 给pc的数据框添加名称列和分组列(用来画图)
pc$group = group
pc$names = rownames(pc)

## step 4: 绘图
# draw PCA plot figure
xlab=paste("PC1(",pro1,"%)",sep="") 
ylab=paste("PC2(",pro2,"%)",sep="")
pca=ggplot(pc,aes(PC1,PC2)) + 
geom_point(size=3,aes(shape=group,color=group)) + 
geom_text(aes(label=names),size=4)+labs(x=xlab,y=ylab,title="PCA") + 
geom_hline(yintercept=0,linetype=4,color="grey") + 
geom_vline(xintercept=0,linetype=4,color="grey") + 
theme_bw()

# 保存结果
ggsave(outname,pca,width=10,height=8)

DESeq分析流程

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq")
library("DESeq")
database <- read.table(file = "macrophage_genes.count_table.matrix", sep = "\t", header = T, row.names = 1)
countData <- database[,1:3]
condition <- factor(c("A","B","C"))#无生物学重复
#type <- factor(c(rep("A",3), rep("B",3))) 有生物学重复
database <- round(as.matrix(countData))#取整数型
dds <- newCountDataSet(database,condition)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds, method="blind", sharingMode="fit-only" )#无生物学重复
dds <- estimateDispersions(dds) # 有生物学重复
resAvsB <- nbinomTest(dds,"A","B")
resAvsC <- nbinomTest(dds,"A","C")
table(resAvsB$pval <0.05)
table(resAvsC$pval <0.05)
resAvsB <- resAvsB[order(resAvsB$pval),]#排序
resAvsC <- resAvsC[order(resAvsC$pval),] 
plotMA(resAvsB, ylim = c(-5,5), col = ifelse(resAvsB$pval>=0.05, "gray32", "red3"),linecol = "#ff000080" )#画MA图
norcounts <- counts(dds, normalized=T)#提取标准化后的counts
write.table(as.data.frame(norcounts),
file="normalized_counts_ABC.txt",
sep="\t",
quote = F)
resAvsBup <- subset(resAvsB, pval < 0.05 & log2FoldChange >0 )
resAvsBdown <- subset(resAvsC, pval < 0.05 & log2FoldChange <0 )

ggplot2:气泡图及散点图小结

气泡图也属于散点图的一种,在散点图的基础上改变点的形状,大小和颜色

1.如何改变点的形状

用自带的mtcars演示

p = ggplot(mtcars,aes(wt,mpg))
p + geom_point(shape=17,size=4)
p + geom_point(aes(shape = factor(cyl)),size=4)

2.改变点的大小:气泡图

加载包、数据的读入及数据格式

library(ggplot2)
pathway = read.table("R0-vs-R3.path.richFactor.head20.tsv",header=T,sep="\t")

3.初始化数据

pp = ggplot(pathway,aes(richFactor,Pathway))
pp + geom_point()

气泡图(三维数据),size 可以是一个值,也可以是数据中的一列

pp + geom_point(aes(size=R0vsR3))

添加颜色变化——四维数据

scale_colour_gradient:自定义一个连续型的配色,常用于热图;

pbubble = pp + geom_point(aes(size=R0vsR3,color=-1*log10(Qvalue)))
pbubble
pbubble + scale_colour_gradient(low="green",high="red")


expression函数改变样式,[]是用来添加下标,^是用来添加上标

pr = pbubble + scale_colour_gradient(low="green",high="red") + labs(color=expression(-log[10](Qvalue)),size="Gene number",x="Rich factor",y="Pathway name",title="Top20 of pathway enrichment")
pr
pr + theme_bw()

散点图小结

• 颜色: color/colour (可渐变)
• 形状:shape
• 大小:size (可渐变)
• 透明度:alpha (练习)
• 改变坐标轴的范围: xlim()和ylim()函数
• 改变标题,坐标轴标题和legend标题:labs()函数
• 保存图片:ggsave()函数
支持多种图片格式:png,pdf,svg等

ggsave("pathway_enrichment.png")
ggsave("pathway_enrichment_gray.png",pr)
ggsave("volcano4.pdf",volcano,width=4,height=4)
ggsave("volcano8.pdf",volcano,width=8,height=8)

绘制直线小结

• geom_abline(slope=xx,intercept=xx): 绘
制直线
• geom_hline(yintercept=xx):绘制水平线
• geom_vline(xintercept=xx):绘制垂直线
• size: 线条粗细
• color: 线条颜色
• linetype: 线条类型

详细讲解参考ggplot2说明文档基迪奥在线视频

ggplot2:火山图

火山图是散点图的一种,也是描点,用 geom_point()绘图

1.加载包,读入数据及数据格式

library(ggplot2)
data = read.table("R0-vsR3.isoforms.filter.tsv",header=T,row.names=1)

2.绘图

r03 = ggplot(data,aes(log2FC,-1*log10(FDR)))
r03 + geom_point()

3.改变点的颜色

r03 + geom_point(color="red")
r03 + geom_point(aes(color="red"))
r03 + geom_point(aes(color=significant)) # 按照“ significant”这一列定义点的颜色;

4.设置坐标轴范围和标题

函数xlim(),ylim(),规定X、 Y轴范围;labs(title=“..”,x=“..”,y=“..”),确定标题和坐标轴标签;expression函数改变样式,[]是用来添加下标,^是用来添加上标

r03xy = r03 + geom_point(aes(color=significant)) + xlim(-4,4) + ylim(0,30)
r03xy + labs(title="Volcano plot",x="log2(FC)")
r03xy + labs(title="Volcano plot",x=expression(log[2](FC)),y=expression(-log[10](FDR)))

5.自定义颜色

scale_color_manual(): 自定义颜色配色;

r03xyp = r03xy + labs(title="Volcano plot",x=expression(log[2](FC)),y=expression(-log[10](FDR)))
r03xyp + scale_color_manual(values = c("green","black", "red"))
volcano = r03xyp + scale_color_manual(values = c("#00ba38","#619cff", "#f8766d"))

6.添加阈值线

geom_hline()添加水平线; geom_vline()添加垂直线

volcano+geom_hline(yintercept=1.3)+geom_vline(xintercept=c(-1,1))

7.改变线条类型

linetype参数0 = blank, 1 = solid, 2 = dashed, 3 = dotted, 4 = dotdash, 5 = longdash, 6 = twodash

volcano+geom_hline(yintercept=1.3,linetype=4)+geom_vline(xintercept=c(-1,1),linetype=4)

详细讲解查看帮助文档基迪奥在线视频

ggplot2:散点图

1.加载ggplot2包,读取数据及数据格式

library(ggplot2)
fpkm = read.table("all.fpkm",header=T,row.names=1)
head(fpkm,10)

2.绘图

mp = ggplot(fpkm,aes(C1_FPKM,C2_FPKM))
mp + geom_point()

3.取log10后重新作图

mp = ggplot(fpkm,aes(log10(C1_FPKM),log10(C2_FPKM)))
mp + geom_point()

4.添加直线

slope设置斜率,intercept设置截距,color设置线条颜色,size设置线条粗细

mp + geom_abline(slope=1,intercept=0,color="red",size=2) + geom_point()

ggplot2的核心思想就是图层,可以像PS一样一层层的添加图层,以加号连接,更详细参数可参考ggplot2的帮助文档http://docs.ggplot2.org/current/

视频教程:http://www.omicshare.com/class/home/index/singlev?id=34