写在前面：

这部分主要做一些数据可视化，富集分析暂时放下一部分，如果想跳过这里，请直接移步RNA-seq(9):富集分析

---------------------------------------------------

参考资料：
Analyzing RNA-seq data with DESeq2
[Count-based Differential expression Analysis of RNA-seq Data]

1 MA plot

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.

• 没有经过 statistical moderation平缓log2 fold changes的情况

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

• 经过lfcShrink 收缩log2 fold change

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

2 Plot counts

DESeq2提供了一个plotCounts()函数来查看某一个感兴趣的gene在组间的差别。counts会根据groups分组。更多的参数请输入命令?plotCounts下面我们来看plot两个genes

• 一个是padj最小的gene
• 一个是
  直接用plotCounts命令

# 不画图，只显示数据 plotCounts(dds, gene=which.min(res$padj), intgroup="condition", returnData=TRUE) #只画图，不显示数据 plotCounts(dds, gene="ENSMUSG00000024045", intgroup="condition", returnData=FAULSE)

下面用ggplot来画Akap8的box图和point图

• boxplot

# 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

• point plot

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

3 PCA(principal components analysis)

• 上面的分析，我们使用的原始的counts数据。但是又一些下游其他分析比如热图(heatmap), PCA或聚类(clustering)我们需要data的转换后的格式，因为如何最好的计算未转换的counts的距离测度仍然不清楚。一个选择是进行log变换。但是因为很多samples的count为0（这意味着 log(0)=−∞，当然也可以使用家counts，比如y=log(n+1)或更普遍使用的y=log(n+n0 ），n代表count值，n0是某个正常数。
  但是也有一些其他的方法提供更好的理论矫正，其中有一个称为variance stabilizing transformation(VST),它消除了方差对mean均值的依赖，尤其是低均值时的高log counts的变异。
• DESeq2提供了plotPCA函数进行PCA分析。?plotPCA查看帮助文件。
  用法如下

vsdata <- vst(dds, blind=FALSE) plotPCA(vsdata, intgroup="condition")


PCA.jpeg

4热图：两部分

4.1 count matrix 热图

根据不同的数据转换方式，可以产生不同类型的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

4.2 sample-to-sample distances热图

• 转换数据还可以做出样本聚类热图。用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