1扑庞、在R中安裝DESeq2軟件包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
2譬重、運(yùn)行
準(zhǔn)備工作
# read.table ,文件有header,第一行為row.name;
input_data <- read.table("deseq2_input.txt", header=TRUE, row.names=1)
# 取整,函數(shù)round
input_data <-round(input_data,digits = 0)
# 準(zhǔn)備文件
# as.matrix 將輸入文件轉(zhuǎn)換為表達(dá)矩陣；
input_data <- as.matrix(input_data)
# 控制條件：因子（對(duì)照2罐氨，實(shí)驗(yàn)2）
condition <- factor(c(rep("ct1", 2), rep("exp", 2)))
# input_data根據(jù)控制條件構(gòu)建data.frame
coldata <- data.frame(row.names=colnames(input_data), condition)
library(DESeq2)
# build deseq input matrix 構(gòu)建輸入矩陣
#countData作為矩陣的input_data臀规；colData Data.frame格式；控制條件design;
dds <- DESeqDataSetFromMatrix(countData=input_data,colData=coldata, design=~condition)
DESeq2 進(jìn)行差異分析
#構(gòu)建dds對(duì)象
dds <- DESeq(dds)
# 提取結(jié)果
#dds dataset格式轉(zhuǎn)換為DESeq2 中result數(shù)據(jù)格式栅隐，矯正值默認(rèn)0.1
res <- results(dds,alpha=0.05)
#查看res(DESeqresults格式),可以看到上下調(diào)基因
summary(res)
# 將進(jìn)過矯正后的表達(dá)量結(jié)果加進(jìn)去塔嬉；
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 輸出結(jié)果
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
#上調(diào)基因的提取;
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > up.gene.txt
#下調(diào)基因的提取;
awk '{if($3<-1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > down.gene.txt
根據(jù)第1列是Geneid，第7,8租悄，9谨究，10列是counts數(shù)，用awk提取出geneID和counts恰矩。
awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' out_counts.txt >subread_matrix.out
參考：
生物信息樹
