2020-02-26 代作图 Volcano plot火山图

今下午有个买家要求代做差异分析，顺便出个火山图

屏幕截图 2020-02-26 22.19.34.png

library(DESeq2)
library(limma)
mydata <- read.table("normalize.txt",head=T) ## 输入文件为基因count矩阵，列名为样本，行名为基因
rt=as.matrix(mydata)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
mydata <- avereps(data)
mydata <-round(mydata,0)
mysampleDesign <- data.frame(row.names = colnames(mydata),group_list = c(rep("BB",11),rep("SS",5))) ### 后面的为对照组
dds <- DESeqDataSetFromMatrix(countData=mydata, colData=mysampleDesign, design=~ group_list) ##构建dds对象
head(dds) 
dds <- DESeq(dds) ## 数据预处理，标准化等，计算差异
#estimating size factors
#estimating dispersions
#gene-wise dispersion estimates
#mean-dispersion relationship
#final dispersion estimates
#fitting model and testing

res_BBvSS <- results(dds,contrast=c("group_list", "BB","SS")) #差异结果
head(res_BBvSS) 
summary(res_BBvSS)

resOrdered_BBvSS <- res_BBvSS[order(res_BBvSS$padj), ]  ## 按照padj列排序
resSig_BBvSS <-subset(resOrdered_BBvSS, padj < 0.1 & abs(log2FoldChange) > 2) #筛选显著表达基因 筛选条件padj < 0.1 并且 log2FoldChange 绝对值大于2
summary(resSig_BBvSS) 

resSig <- as.data.frame(resSig_BBvSS) ## 转换为数据框类型，方便下面分析
head(resSig)
write.csv(resSig,file= "resSig_BBvSS_padj0.1.csv") ## 输出差异基因DEseq结果，csv格式
write.table(resSig, file="resSig_BBvSS_padj0.1.txt",sep="\t",quote=F)  ## 输出差异基因DEseq结果，txt格式

后面自己瞎改了一张

1_volcano.jpeg