【bulk RNA-seq】DESeq2差异分析 硬核
2024-04-16 07:42:09
心血来潮写了个bulk RNA-seq的代码，该有的基本都有了。省去了路径和具体的基因。以后更新一下用公共数据跑出来的图吧。最近比较忙。 包非原创，代码是原创。转载或引用请注明出处或与我联系。
版本：R 4.1
################################
## 导入需要用到的R package以及自定义函数
################################
library(DESeq2) #差异分析的包
library(clusterProfiler) #基因名转换&GO富集分析的包
library(dplyr)
library(RColorBrewer) #调色包
library(gplots)
library(corrplot)
MyName2Color <- function(name,palette = brewer.pal(9,"Set1"),is.row = FALSE,na.value = NULL){name.factor <- as.factor(name)if(!is.null(names(palette))){palette=palette[levels(name.factor)]}else{print("your palette order must adapt to you names level order")}name.color <- palette[name.factor]name.color <- as.matrix(name.color)if(is.row){name.color <- t(name.color)}if(!is.null(na.value)){name.color[is.na(name.color)]=na.value}return(name.color)
} #自定义函数
################################
## 导入原始数据 并标注样本信息
################################
setwd("/home/") #设置工作路径
read.count <- read.delim("read_count.tab") #读表达矩阵
rownames(read.count) <- read.count$ID
# sample.information <- read.csv() #读样品信息
read.count <- read.count[,-1] #去掉第一列
sample.information <- as.data.frame(colnames(read.count)) #变为dataframe格式
colnames(sample.information) <- "ID"
row.names(sample.information) <- sample.information$ID
sample.information$condition <- c("KO","KO","KO","WT","WT","WT") #设置样本信息################################
## DEseq2
################################
dds <- DESeqDataSetFromMatrix(countData = read.count, colData = sample.information, design = ~condition) #构建DEseq矩阵
dds.temporary <- rowSums(counts(dds) >= 10) >= 3 #过滤低表达基因，至少在3个样品中都有5个reads
dds <- dds[dds.temporary, ]
dds <- DESeq(dds) #运行DEseq2##### PCA 看看样品间的差异(比如有无明显批次效应)
vsdata <- rlog(dds, blind=FALSE) #归一化(样本数小于30用rlog函数)
plotPCA(vsdata, intgroup = "condition") #绘制样本的PCA图result <- results(dds,alpha = 0.1) #导出结果
summary(result) #查看结果的概览
sum(result$padj<0.1, na.rm = TRUE) #检查pvalue<0.1的基因个数
res_data <- merge(as.data.frame(result), #DEseq2的检验结果as.data.frame(counts(dds,normalize=TRUE)), #用DEseq2 normalize后的表达矩阵by="row.names",sort=FALSE) # 合并DEseq2的结果与标准化之后的表达矩阵Symbol <- bitr(res_data$Row.names,fromType = 'ENSEMBL',toType = 'SYMBOL',OrgDb = "org.Mm.eg.db") # ENSEMBL ID转基因名
res_data <- res_data[res_data$Row.names%in%Symbol$ENSEMBL,] #有部分ID转换失败，去除之
Symbol <- Symbol[!duplicated(Symbol$ENSEMBL),] #有部分ID转换失败，去除之
rownames(res_data) <- res_data$Row.names #改变行名称
rownames(Symbol) <- Symbol$ENSEMBL #改变行名称
res_data <- cbind(res_data, Symbol) #把基因名放到res_data里，列名为SYMBOL
res_data <- res_data[!duplicated(res_data$SYMBOL),] #去掉SYMBOL列里重复的行up_DEG <- subset(res_data, padj < 0.1 & log2FoldChange > 0) #筛选上调基因
up_DEG$cluster <- "up"
down_DEG <- subset(res_data, padj < 0.1 & log2FoldChange < 0) #筛选下调基因
down_DEG$cluster <- "down"
gene.info <- read.csv("/home/xxx.csv")#基因信息的表格
rownames(gene.info) <- gene.info$EnsemblGeneID
down_DEG_merge <- cbind(down_DEG,gene.info[rownames(gene.info)%in%rownames(down_DEG),])
up_DEG_merge <- cbind(up_DEG,gene.info[rownames(gene.info)%in%rownames(up_DEG),])####### 将全部基因的信息以及上下调基因的信息写入csv
write.csv(res_data, "all.csv") #写入csv表格，保存在工作目录下
write.csv(up_DEG_merge, "up.csv")
write.csv(down_DEG_merge, "down.csv")################################
up_and_down_gene <- rbind(down_DEG, up_DEG) #合并两个表格
rownames(up_and_down_gene) <- up_and_down_gene$SYMBOL##### 绘制heatmap
pdf("Normalized_count.pdf",20,25)
heatmap.2(x = as.matrix(log2(up_and_down_gene[,c("")]+1)),scale = "none", col=bluered(100), trace = "none", density.info = "none",RowSideColors = MyName2Color(up_and_down_gene$cluster, palette = c("#d53e4f", "#3288bd")))
dev.off()##### 样本间相关性分析
cor.sample <- cor(up_and_down_gene[,c("")])#相关性分析函数
pdf("correlation_afterNormalization.pdf",9,7) #做相关性分析图
corrplot(cor.sample, method="number", type="upper", order="hclust", tl.col="black", tl.srt=45)
dev.off()##### GO富集分析
GO_database <- 'org.Mm.eg.db' #设置GO数据集 这个数据集用BiocManager::install("org.Mm.eg.db")可以安装
universe <- bitr(Symbol$SYMBOL,fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = GO_database) #选取背景基因 并转换为ENTREZ ID
gene.UPregulate <- bitr(up_DEG$SYMBOL,fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = GO_database)
gene.DOWNregulate <- bitr(down_DEG$SYMBOL,fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = GO_database)GO.UP <- enrichGO(gene.UPregulate$ENTREZID,#GO富集分析OrgDb = GO_database,keyType = "ENTREZID",#设定读取的gene ID类型ont = "BP",#Biological Process）pvalueCutoff = 0.05,#设定p值阈值qvalueCutoff = 0.05,#设定q值阈值readable = T)
GO.DOWN <- enrichGO(gene.DOWNregulate$ENTREZID,#GO富集分析OrgDb = GO_database,keyType = "ENTREZID",#设定读取的gene ID类型ont = "BP",pvalueCutoff = 0.05,#设定p值阈值qvalueCutoff = 0.05,#设定q值阈值readable = T)pdf("GOplot.upregulate.pdf",12,5)
barplot(GO.UP)
dev.off()pdf("GOplot.downregulate.pdf",12,5)
barplot(GO.DOWN)
dev.off()##### xxx相关基因表达量的变化
#在kegg上找到xxx所在的xxx通路的下游基因xxx：即xxx，是受体基因 下游有xxx
rownames(res_data) <- res_data$SYMBOL
gene.count_normalized <- res_data[,c("")]
MyBoxplot_genecount <- function(read.count, gene){df <- as.data.frame(t(read.count[gene,]))df$sample2 <- c("KO","KO","KO","WT","WT","WT")colnames(df) <- c("Read_Count","Sample")# print(df)ggplot(df, aes(x = Sample, y = Read_Count)) + geom_boxplot() + labs(title = gene) + theme_bw() + #设置背景为无色theme(plot.title = element_text(face = "bold.italic", size = 19)) + #设置表格字体为斜体，以及字号大小theme(panel.grid = element_blank()) #设置标尺线为无色
}#写个函数方便批量做图pdf("xxx_associated_genes.pdf",4,4) #批量做图
MyBoxplot_genecount(gene.count_normalized, "xxx")
MyBoxplot_genecount(gene.count_normalized, "xxx")
MyBoxplot_genecount(gene.count_normalized, "xxx")
MyBoxplot_genecount(gene.count_normalized, "xxx")
dev.off()##### 火山图
data_VolcanoMap <- res_data[,c("log2FoldChange", "padj", "SYMBOL")]
data_VolcanoMap$UpDown <- "invariable" #加新的一列，列名为UpDown，方便做图的时候分类
data_VolcanoMap[data_VolcanoMap$SYMBOL%in%up_DEG$SYMBOL,]$UpDown <- "UP_regulated"
data_VolcanoMap[data_VolcanoMap$SYMBOL%in%down_DEG$SYMBOL,]$UpDown <- "DOWN_regulated"
data_VolcanoMap$nega_logPadj <- -log10(data_VolcanoMap$padj) #对Pvalue adj取负log10pdf("VolcanoMap.pdf",6.6,5)
ggplot(data = data_VolcanoMap) +geom_point(aes(x=log2FoldChange, y=nega_logPadj, color=UpDown)) +labs(title="Volcano Plot") +theme_bw() +theme(panel.grid = element_blank())
dev.off()【bulk RNA-seq】DESeq2差异分析 硬核相关推荐
- 【硬核】肝了一个月，Cisco网络工程师知识点总结
[硬核]肝了一个月,Cisco网络工程师知识点总结 高能预警,本文是我一个月前就开始写的,所以内容会非常长,当然也非常硬核,有实验,有命令,所以才写到了现在. 我相信90%的读者都不会一口气看完的,因 ...
- 超硬核！兔兔阿里p7学长给的面试知识库
一个阿里p7学长给的nosql面试知识库,绝对真实,学会了去面呀. 最近整理了一下超硬核系列的文章和面经系列的文章,可以持续关注下: 超硬核系列历史文章:(我保证每篇文章都有值得学习的地方,并且对小白 ...
- 超硬核！！！一篇文章搞定TCP、UDP、Socket、HTTP（详细网络编程内容+现实解释三次握手四次挥手+代码示例）【网络编程 1】
TCP.UDP.Socket 一天面试的经验: 什么是网络编程 网络编程中两个主要的问题 网络协议是什么 为什么要对网络协议分层 计算机网络体系结构 1 TCP / UDP 1.1 什么是TCP/IP ...
- 超硬核全套Java视频教程(学习路线+免费视频+配套资料)
文内福利,扫码免费领取 Hello,各位锋迷们,我是小千.很多学习Java的小伙伴都在找的全套免费java视频教程,这里全都有,资料齐全,拿来吧你! 零基础学Java的学习路线图是怎样的?! 曾经写过 ...
- 中科大硬核“毕业证”：“一生一芯”计划下，5位本科生带自研芯片毕业
作者 | 包云岗 编辑 | 伍杏玲 本文经作者授权转载自包云岗知乎 [CSDN编者按]近日,中国科学院大学五位本科生的硬核"毕业证"引发IT圈热议,在"一生一芯" ...
- AIの幕后人：探秘“硬核英雄”的超级武器
作者 | 云计算的阿晶 出品 | AI科技大本营(ID:rgznai100) 掐指一算八年之前,那时正是国内互联网卯足劲头起飞的一年,各行各业表现都很突出,尤其是与人们生活密切相关的手机,正大踏步地从 ...
- 硬核吃瓜！上万条数据撕开微博热搜真相
作者 | 徐麟 来源 | 转载自数据森麟(ID:shujusenlin) 吃瓜前言 关于新浪微博,向来都是各路吃瓜群众聚集之地,大家在微博中可以尽情吃瓜,各种类型的瓜应有尽有,只有你想不到的,没有你吃 ...
- 硬核！我的导师手写129页毕业论文，堪比打印！
点击上方"视学算法",选择加"星标"或"置顶" 重磅干货,第一时间送达 本文来源:浙江大学.中国青年报 "闲来无事翻了一下导师的毕 ...
- 全网刷爆！B站Up主何同学带火这只A股：最硬核桌子,苹果也做不到！
点击上方"视学算法",选择加"星标"或"置顶" 重磅干货,第一时间送达 B站UP主一条视频,直接让A股上市公司股价涨超13%? 失踪三个月之 ...
最新文章
- 2012年的这些经典书目你读了没？
- JAVA程序设计----关于字符串的一些基本问题处理
- 如何把Sql Server2005 数据库转换成Access
- [c/c++] c 操作mysql数据库
- paip兼容windows与linux的java类根目录路径的方法
- Hadoop 各组件介绍
- 测试方案和测试计划区别
- waf绕过—过360主机卫士sql注入
- PHP格式化 插件 vs code
- socket编程c++
- 郑州 - 天气总是灰蒙蒙的
- C语言的转义字符，八进制
- NVDIMM原理与应用之一：基本原理
- 宝塔面板怎么实名认证_云服务器安装宝塔面板完整教程
- 启动、停止elasticsearch的脚本（没有技术含量）
- Java socket服务端
- G4L-硬盘对拷工具
- 排查not eligible for getting processed by all BeanPostProcessors
- HUAWEI 擎云L420 折腾记 （搭建arm gcc、openocd 雅特力 MCU开发环境）
- 2017年大数据 云计算 物联网发展趋势