我该分享点什么？分享点什么能让大家有兴趣？分享点什么能帮到大家？这是最近一直困扰我的三个问题。
昨晚在翻腾公众号的时候，发现2022年01月20日公众号发表了《转录组数据分析 RNA-seq》，NGS分析相关的推文就一直没有更新过，而且该文仅分享了上游分析的众多pipeline之一，并未涉及下游分析以及绘图，于是乎有了这篇推文。
数据简介
当我们拿到了基因的表达矩阵，接下来就是差异分析。接下来我将以2020年发表在BMC的一篇鹌鹑的文章中部分RNAseq数据为例进行分析。（文末附基因表达矩阵的百度云链接，大家可自行下载练手，不必下载原始测序文件重新跑一遍上游分析）
Morris KM, Hindle MM, Boitard S, et al. The quail genome: insights into social behaviour, seasonal biology and infectious disease response. BMC Biol. 2020;18(1):14. Published 2020 Feb 12. doi:10.1186/s12915-020-0743-4
基因定量矩阵文件格式如下（第一列为基因名，第二列为基因所在染色体，第三列为基因的起始位置，第四列为基因的终止位置，后八列为8个文库的基因定量信息，前4个文库为对照组，后四个文库为处理组）：
DESeq2的介绍
DESeq2是一个为高维计量数据的归一化、可视化和差异表达分析而设计的一个R语言包。它通过经验贝叶斯方法(empirical Bayes techniques)来估计对数倍数变化(log2foldchange）和离差的先验值，并计算这些统计量的后验值。
它由美国北卡罗莱纳大学教授Michael Love(michaelisaiahlove@gmail.com)于2014年发布，目前仍在更新与维护中，是目前差异表达分析方面最常用的R包。
当然像limma，edgeR也是常用的计算差异的R包，至于三者想用哪个，随你个人喜好。
DESeq2的安装
#这一步请在R或Rstudio中运行
if (!requireNamespace("BiocManager")) install.packages("BiocManager") BiocManager::install("DEseq2")
DESeq2差异分析
# 加载DESseq2 R包 library(DESeq2) # 读取基因表达矩阵数据 countData <- read.table("C:\\Users\\Jeff\\OneDrive\\桌面\\online_count.txt",row.names=1,header=F)) ''' BTW 如果你不想输入太长的绝对路径，你可以先设置工作路径，就像python中os模块中os.chdir()一样 setwd("C:\\Users\\Jeff\\OneDrive\\桌面") countData <- read.table("online_count.txt",row.names=1,header=F)) ''' countData <- countData[,c(4:11)] # 查看数据的头部 head(countData) # 删除表达量均值小于1的基因 countData <- countData[rowMeans(countData)>1,] # 构造一个处理组-列交联矩阵 condition <- factor(c(rep("control",4),rep("treatment",4))) colData <- data.frame(row.names=colnames(countData), condition=condition) # 构造DESeq2计算矩阵 dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~condition) # 差异分析 dds <- DESeq(dds) # 提取差异分析结果 res <- results(dds, contrast=c("condition","treatment","control")) # 按照矫正后的p值进行结果的排序 res <- res[order(res$padj),] # 提取所有的差异基因，其中阈值设置为p-adjust value < 0.01，log2FoldChange的绝对值大于1 # 另外如若想单独提取上下调基因的相关信息，修改subset的筛选条件即可： #p-adjust value < 0.01，log2FoldChange > 1 #p-adjust value < 0.01，log2FoldChange <- 1 diff_gene_deseq2 <-subset(res,padj < 0.01 & (log2FoldChange > 1 | log2FoldChange < -1)) diff_gene_deseq2 <- row.names(diff_gene_deseq2) # 整合差异结果以及各文库标准化后的定量数值至同一数据框 resdata <- merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE) #得到csv格式的差异表达分析结果 write.csv(resdata,file= "DESEQ_result.csv",row.names = F)
看到此处，有人想问数据有了，火山图不是最经典的转录组图吗？怎么画图？有没有快捷作图的代码？QAQ
Yep 复制~ 粘贴~
#导入R包，如果未安装请按DESeq2的方法安装 library(RColorBrewer) library(ggpubr) library(ggthemes) #在这里我们定义一个函数对接DESeq2产出的结果 #函数需要传入res和pdf_path两个参数 VolcanoPlot=function(res,pdf_path) { res$logP <- -log10(res$padj) #padj取-log值生成新列 res <- as.data.frame(res) #转为数据框形式 res$group <- 'not-significant' #根据选定的阈值对基因添加label（上调？下调？还是没差异？） res$group[which((res$padj <=0.01) & (res$log2FoldChange >=1) )] = 'up-regulated' res$group[which((res$padj <=0.01) & (res$log2FoldChange <=-1) )] = 'down-regulated' table(res$group) #简单看下多少个上调基因多少个下调基因 ggscatter(res, #画图咯 x = 'log2FoldChange', #x轴数据 y ='logP' , #y轴数据 color = 'group', #用标签区分散点颜色为三个颜色 size = 1, #散点大小 palette = c("#2f5688","#BBBBBB","#CC0000"), #手动定义三个颜色应为什么色 xlab = "log2FoldChange", #x轴标签 ylab = "-log10(padj)" #y轴标签 )+theme(text = element_text(size = 10),legend.position = "right")+ #theme中定义全局字体大小和图例的位置（在右边） geom_hline(yintercept = 2,linetype='dashed')+ #绘制y轴阈值线 即padj =0.01处 geom_vline(xintercept = c(-1,1),linetype='dashed') #绘制x轴阈值线 即log2FoldChange分别为1和-1处 ggsave(pdf_path, width = 6.4, height = 3.6 ,dpi = 300) #保存图片，指定pdf高度和宽短，dpi分辨率为300 return(TRUE) } #在这里 我们运行一下我们定义的函数 VolcanoPlot(res,"tmp.pdf")
## PCA library("FactoMineR") library("factoextra") count.pca <- counts(dds) %>% as.matrix %>% t %>% PCA(.,graph = F) pdf("PCA.pdf") fviz_pca_ind(count.pca, col.ind=condition, mean.point=T, addEllipses = F, legend.title="Groups")+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid")) dev.off() ## Heatmap library(pheatmap) pdf("Heatmap.pdf") normalized_counts <- counts(dds, normalized=TRUE) normalized_counts_mad <- apply(normalized_counts, 1, mad) normalized_counts <- normalized_counts[order(normalized_counts_mad, decreasing=T), ] Heatmap_Data=normalized_counts[rownames(diff_gene_deseq2),] annotation_col = data.frame( Conditions = factor(c(rep("control", 4),rep("treatment", 4)) )) colnames(Heatmap_Data)=factor(colnames(Heatmap_Data)) pheatmap(Heatmap_Data, scale = "row", show_rownames = F, cluster_cols = F, main = "DAR Heatmap", cutree_rows = 2, gaps_col = 4, color = colorRampPalette(c("navy", "white", "firebrick3"))(50), ColSideColors = rep(c("pink", "purple"), each = 4)) dev.off() ### pearson library(amap) library(gplots) rld <- rlog(dds, blind=TRUE) rlogMat <- assay(rld) rlogMat <- rlogMat[order(normalized_counts_mad, decreasing=T), ] hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100) pearson_cor <- as.matrix(cor(rlogMat, method="pearson")) hc <- hcluster(t(rlogMat), method="pearson") heatmap.2(pearson_cor, Colv = as.dendrogram(hc), Rowv=as.dendrogram(hc), symm=T, trace="none", col=hmcol, main="The pearson correlation of each sample") dev.off()
小结
好了，今天摸鱼时间有限，实在不能再摸了，课题还亟待处理，只能先写到这里了，接下来我将不定时的继续更新转录组相关的下游个性化分析，转录组系列写完之后我们就开始Chip-seq，ATAC-seq，非模式物种的单细胞ATAC-seq等多组学的上下游数据分析，欢迎讨论~
互相学习，加油！
