转录组测序火山图_转录组差异基因筛选标准

# 设置R语言镜像
# options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
# options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# 安装R包
# if(!require(c("ggthemes","ggpubr","ggthemes","ggrepel"))) install.packages(c("ggthemes","ggpubr","ggthemes","ggrepel"))
# BiocManager::install("DESeq2")

#载入R包
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(ggpubr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggrepel))
suppressPackageStartupMessages(library(ggthemes))

exprset <- read.table("RNA-seq_counts_matrix.csv",sep = ",",header = T,check.names = F)
rownames(exprset) <- exprset[,1]
exprset <- exprset[,-1]
exprset <- exprset[apply(exprset,1,function(x) sum(x > 1) > 1),] #先判断值是否为0，得到逻辑值，再取和，判断0的个数是否小于3
dim(exprset)
# 7428    4
head(exprset)

expcor <- cor(exprset, method = "spearman")
head(expcor)
pheatmap::pheatmap(expcor, clustering_method = "average",
                   treeheight_row = 0,treeheight_col = 0,
                   display_numbers = T)

# t_exprset <- t(exprset)
# t_exprset <- t_exprset[,names(sort(apply(t_exprset,2,mad),decreasing = T))[1:500]]
# out.dist <- dist(t_exprset,method = 'euclidean')
# out.hclust <- hclust(out.dist,method = 'complete')
# rect.hclust(out.hclust,k=3)
# plot(out.hclust,xlab = "",main = "")

group_list <- factor(c(rep("untrt",2),rep("trt",2))) #因子型变量
group_list
table(group_list)
## group_list
##   trt untrt 
##     2     2
colData <- data.frame(row.names = colnames(exprset),
                        group_list = group_list)
colData
dds <- DESeqDataSetFromMatrix(countData = exprset,
                               colData = colData,
                               design = ~group_list) #~在R里面用于构建公式对象，~左边为因变量，右边为自变量。
head(dds)
## class: DESeqDataSet 
## dim: 6 4 
## metadata(1): version
## assays(1): counts
## rownames(6): ENSMUSG00000000028 ENSMUSG00000000088 ...
##   ENSMUSG00000000134 ENSMUSG00000000142
## rowData names(0):
## colnames(4): control1 control2 treat1 treat2
## colData names(1): group_list
tem_f <- 'RNA-seq_DESeq2-dds.Rdata'

if (!file.exists(tem_f)) {
    dds <- DESeq(dds) # 标准化
    save(dds,file = tem_f)
  }
load(file = tem_f)
# 结果用`result()`函数提取
res <- results(dds,
              contrast = c("group_list","untrt","trt")) # 差异分析结果
resOrdered <- res[order(res$padj),] # 对结果按照调整后的p值进行排序
head(resOrdered)
summary(res)
## 
## out of 7428 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 465, 6.3%
## LFC < 0 (down)     : 507, 6.8%
## outliers [1]       : 0, 0%
## low counts [2]     : 2160, 29%
## (mean count < 4)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

DEG <- as.data.frame(resOrdered)
DESeq2_DEG <- na.omit(DEG)
diff <- subset(DESeq2_DEG,pvalue < 0.05) #先筛选P值
up <- subset(diff,log2FoldChange > 2) #上调
down <- subset(diff,log2FoldChange < -2) #下调
#可利用`write.csv()`函数保存文件

DEG_data <- DESeq2_DEG
DEG_data$logP <- -log10(DEG_data$padj) # 对差异基因矫正后p-value进行log10()转换
dim(DEG_data)
## [1] 5268    7
#将基因分为三类：not-siginficant，up，dowm
#将adj.P.value小于0.05，logFC大于2的基因设置为显著上调基因
#将adj.P.value小于0.05，logFC小于-2的基因设置为显著上调基因
DEG_data$Group <- "not-siginficant"
DEG_data$Group[which((DEG_data$padj < 0.05) & DEG_data$log2FoldChange > 2)] = "up-regulated"
DEG_data$Group[which((DEG_data$padj < 0.05) & DEG_data$log2FoldChange < -2)] = "down-regulated"
table(DEG_data$Group)
## 
##  down-regulated not-siginficant    up-regulated 
##             336            4659             273
DEG_data <- DEG_data[order(DEG_data$padj),]#对差异表达基因调整后的p值进行排序
#火山图中添加点(数据构建)
up_label <- head(DEG_data[DEG_data$Group == "up-regulated",],1)
down_label <- head(DEG_data[DEG_data$Group == "down-regulated",],1)
deg_label_gene <- data.frame(gene = c(rownames(up_label),rownames(down_label)),
                                label = c(rownames(up_label),rownames(down_label)))
DEG_data$gene <- rownames(DEG_data)
DEG_data <- merge(DEG_data,deg_label_gene,by = 'gene',all = T)
#不添加label
ggscatter(DEG_data,x = "log2FoldChange",y = "logP",
          color = "Group",
          palette = c("green","gray","red"),
          repel = T,
          ylab = "-log10(Padj)",
          size = 1) + 
  theme_base()+
  scale_y_continuous(limits = c(0,8))+
  scale_x_continuous(limits = c(-18,18))+
  geom_hline(yintercept = 1.3,linetype = "dashed")+
  geom_vline(xintercept = c(-2,2),linetype = "dashed")

#添加特定基因label
ggscatter(DEG_data,x = "log2FoldChange",y = "logP",
          color = "Group",
          palette = c("green","gray","red"),
          label = DEG_data$label,
          repel = T,
          ylab = "-log10(Padj)",
          size = 1) + 
  theme_base()+
  theme(element_line(size = 0),element_rect(size = 1.5))+ #坐标轴线条大小设置
  scale_y_continuous(limits = c(0,8))+
  scale_x_continuous(limits = c(-18,18))+
  geom_hline(yintercept = 1.3,linetype = "dashed")+
  geom_vline(xintercept = c(-2,2),linetype = "dashed")

  其中左侧菜单栏分别是Dataset--选择相关物种参考基因组；
  Filters--选择数据gene ID的类型，并输入gene ID，也存在其他类型的ID输入；
  Attributes--选择需要输出的ID类型；
  点击Result可以输出结果，并且支持文件下载。