【RNA-seq自学09】数据分析之差异基因分析DEseq2、可视化

1、在R中安装DESeq2软件包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

2、运行

准备工作

# read.table ,文件有header,第一行为row.name;
input_data <- read.table("deseq2_input.txt", header=TRUE, row.names=1)

# 取整,函数round
input_data <-round(input_data,digits = 0)

# 准备文件
# as.matrix 将输入文件转换为表达矩阵；
input_data <- as.matrix(input_data)

# 控制条件：因子（对照2，实验2）
condition <- factor(c(rep("ct1", 2), rep("exp", 2)))

# input_data根据控制条件构建data.frame
coldata <- data.frame(row.names=colnames(input_data), condition)

library(DESeq2)
# build deseq input matrix 构建输入矩阵
#countData作为矩阵的input_data；colData Data.frame格式；控制条件design;
dds <- DESeqDataSetFromMatrix(countData=input_data,colData=coldata, design=~condition)

DESeq2 进行差异分析

#构建dds对象
dds <- DESeq(dds)
# 提取结果
#dds dataset格式转换为DESeq2 中result数据格式，矫正值默认0.1
res <- results(dds,alpha=0.05)
#查看res(DESeqresults格式),可以看到上下调基因
summary(res)

# 将进过矫正后的表达量结果加进去；
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
names(resdata)[1] <- "Gene"
#查看(resdata)
# output result 输出结果
write.table(resdata,file="diffexpr-results.txt",sep = "\t",quote = F, row.names = F)

diffexpr-results.txt

可视化展示

plotMA(res)
maplot <- function (res, thresh=0.05, labelsig=TRUE,...){
with(res,plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x", ...))
with(subset(res, padj<thresh), points(baseMean, log2FoldChange, col="red", pch=20,cex=1.5))
}
png("diffexpr-malot.png",1500,1000,pointsize=20)
maplot(resdata, main="MA Plot")
dev.off()

画火山图

install.packages("ggrepel")
library(ggplot2)
library(ggrepel)

resdata$significant <- as.factor(resdata$padj<0.05 & abs(resdata$log2FoldChange) > 1)
ggplot(data=resdata, aes(x=log2FoldChange, y =-log10(padj),color =significant)) +
geom_point() +
ylim(0, 8)+
scale_color_manual(values =c("black","red"))+
geom_hline(yintercept = log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))+
labs(title="Volcanoplot_biotrainee (by sunyi)", x="log2 (fold change)",y="-log10 (padj)")+

theme(plot.title = element_text(hjust = 0.5))+
geom_text_repel(data=subset(resdata, -log10(padj) > 6), aes(label=Gene),col="black",alpha =0.8)

火山图

差异基因的提取

#查看符合条件的差异基因
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |head
#查看差异基因有多少行
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |wc -l
#提取某几列
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 |head
#上调基因的提取;
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > up.gene.txt
#下调基因的提取;
awk '{if($3<-1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > down.gene.txt

根据第1列是Geneid，第7,8，9，10列是counts数，用awk提取出geneID和counts。

awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' out_counts.txt >subread_matrix.out

参考：
生物信息树

最后编辑于：2020.07.08 17:33:17