DESeq2代码-详细注释

东北林业大学_LXG_曾可豪

已于 2023-07-04 10:01:40 修改

阅读量574

点赞数

文章标签： r语言

于 2023-07-03 21:33:58 首次发布

本文链接：https://blog.csdn.net/buchongmingle/article/details/131524108

版权

该文展示了如何利用R语言的DESeq2包对TPM数据进行差异表达分析。首先，读取并处理样本计数矩阵和样本信息，然后定义控制组和处理组，筛选表达量较高的基因。接着，创建DESeqDataSet对象并执行DESeq计算，获取基因的差异倍数和显著性p值。最后，对结果进行排序，标记上调、下调基因，并输出分析结果。

摘要由CSDN通过智能技术生成

1、准备数据

1.1、count数据矩阵

1.2、分组数据矩阵

6、总代码展示

#重开
rm(list = ls())

#基础准备
setwd("D:/great_work/R_wd/DESeq2")

#读取目的样本的TPM数据
library(data.table)
library(dplyr)
#用read.table读取，以设置第一列为行名，第一行为列名
#不检查非法名称，保证名称不变
count_data1 <- read.table("all-Count.txt", row.names=1, header = TRUE, check.names = FALSE)
col_data1 <- fread("col_data.txt")
control <- col_data1$sample[which(col_data1$condition=="control")]#匹配condition列的“control”字符并取出对应的样本名称
treatment <- col_data1$sample[which(col_data1$condition=="treatment")]
#合成目的序列的索引，顺序为“control、treatment”
Target_sample_index_value <- c(which((colnames(count_data1) %in% control) == 1),which((colnames(count_data1) %in% treatment) == 1))
#根据索引提取目的TPM列
count_data2 <- count_data1[,Target_sample_index_value,drop = FALSE]#drop=FALSE参数允许以变量的形式读取索引列表
#整理TPM数据框结构
count_data3 <- count_data2 %>%
  mutate_if(is.numeric, as.integer)#数值转为整数
count_data4 <- count_data3[rowMeans(count_data3)>1,]#筛选行均值大于1的行，去除极低表达基因

#col样本顺序需与TPM样本顺序相同
#检查样本顺序是否相同
col_data2 <- read.table("col_data.txt", row.names=1, header = TRUE, check.names = FALSE)
all(row.names(col_data2) == colnames(count_data4))
#构建样本和分组的对应数据框
col_data2$condition <- as.factor(col_data2$condition)#转换为因子
condition <- as.factor(col_data2$condition)#先对照后实验
condition = relevel(condition, "control")#以“control”组设置为基准

#检查有无缺失值NA
#若TRUE，将缺失值替换为0
any(is.na(count_data4))
#TPM_data5[is.na(TPM_data5)] <- 0

#构建DEseqDataSet对象

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_data4, colData = col_data2, design = ~ condition)

#计算差异倍数
dds <- DESeq(dds)

#查看基因的差异倍数，和显著性p值
#alt前，ref在后，意为alt相较于ref中哪些基因上调/下调
#也可设置其他比较组
#默认使用BH法修正P，alpha值为0.1
#alpha值常用0.05和0.01
res <- results(dds, contrast = c('condition', 'treatment', 'control'), alpha = 0.05)
#整理padj值：排序和修正NA
res$pvalue[is.na(res$pvalue)] <- 1
res$padj[is.na(res$padj)] <- 1
res = res[order(res$padj),]

#标记上调、下调和不显著
#便于下游绘制火山图
res[which(res$log2FoldChange >= 1 & res$padj < 0.05),'change'] <- 'up'
res[which(res$log2FoldChange <= -1 & res$padj < 0.05),'change'] <- 'down'
res [which(abs(res$log2FoldChange) <= 1 | res$padj >= 0.05),'change'] <- 'none'

#输出差异分析结果
write.csv(res,file="DEseq2_allresults.csv",quote = FALSE)