这部分主要做一些数据可视化,富集分析暂时放下一部分,如果想跳过这里,请直接移步RNA-seq(9):富集分析
参考资料:
Analyzing RNA-seq data with DESeq2
Count-Based Differential Expression Analysis of RNA-seq Data
An MA plot is an application of a Bland–Altman plot for visual representation of genomic data. The plot visualizes the differences between measurements taken in two samples, by transforming the data onto M (log ratio) and A (mean average) scales, then plotting these values. Though originally applied in the context of two channel DNA microarray gene expression data, MA plots are also used to visualise high-throughput sequencing analysis.
MA这部分代码主要参考hoptop,并进行修改
In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down.
plotMA(res,ylim=c(-2,2))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
结果如下:
mean of normalized counts.jpeg
It is more useful visualize the MA-plot for the shrunken log2 fold changes, which remove the noise associated with log2 fold changes from low count genes without requiring arbitrary filtering thresholds.
res
结果已经按padj排序了,所以这次要按照行名升序再排列回来,否则和dds
不一致res_order<-res[order(row.names(res)),]
res = res_order
res.shrink <- lfcShrink(dds, contrast = c("condition","treat","control"), res=res)
plotMA(res.shrink, ylim = c(-5,5))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
mean of normalized count _shrinked.jpeg
DESeq2提供了一个plotCounts()
函数来查看某一个感兴趣的gene在组间的差别。counts会根据groups分组。更多的参数请输入命令?plotCounts
下面我们来看plot两个genes
plotCounts
命令# 不画图,只显示数据
plotCounts(dds, gene=which.min(res$padj), intgroup="condition", returnData=TRUE)
#只画图,不显示数据
plotCounts(dds, gene="ENSMUSG00000024045", intgroup="condition", returnData=FAULSE)
# Plot it
plotCounts(dds, gene="ENSMUSG00000024045", intgroup="condition", returnData=TRUE) %>%
ggplot(aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("Akap8")
boxplot_Akap8.jpeg
d <- plotCounts(dds, gene="ENSMUSG00000024045", intgroup="condition", returnData=TRUE)
ggplot(d, aes(x=condition, y=count)) +
geom_point(aes(color= condition),size= 4, position=position_jitter(w=0.5,h=0)) +
scale_y_log10(breaks=c(25,100,400))+ ggtitle("Akap8")
Rplot.jpeg
?plotPCA
查看帮助文件。
用法如下 vsdata <- vst(dds, blind=FALSE)
plotPCA(vsdata, intgroup="condition")
PCA.jpeg
根据不同的数据转换方式,可以产生不同类型的heatmap
library("pheatmap")
select<-order(rowMeans(counts(dds, normalized = TRUE)),
decreasing = TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition","sizeFactor")])
# this gives log2(n + 1)
ntd <- normTransform(dds)
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
cluster_cols=FALSE, annotation_col=df)
上面这两幅图看起来没什么区别,我暂且只放一张
heatmap_ntd.jpeg
dist
函数来获得sample-to-sample距离。距离矩阵热图中可以清楚看到samples之间的相似与否的总概。需要给heatmap函数基于sample距离提供等级聚类hc。#sample to sample heatmap
sampleDists <- dist(t(assay(vsd)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors)
sample-to-sample-heatmap.jpeg