一 DESeq
DESeq的原理是什么以后會慢慢補(bǔ)坑，這次先上代碼.
DESeq必須用count數(shù)據(jù)＜妇Α７拷巍粤攒！
整理出表達(dá)矩陣所森，第一列為基因名囱持，后面為count值。
library(DESeq2)
##要求對象必須是矩陣
liver_gly<-as.matrix(liver_gly)
rownames(liver_gly)<-liver_gly[,1]
#表達(dá)值矩陣
exp=liver_gly[,2:ncol(liver_gly)]
##將基因名和樣本名做成一個列表
dimnames=list(rownames(exp),colnames(liver_gly)[-1])
data=matrix(as.numeric(as.matrix(exp)),nrow = nrow(exp),dimnames = dimnames)
#相同的行取平均值
data=avereps(data)
#這個包要求表達(dá)值必須為整數(shù)
data=round(data,0)
##設(shè)計分組信息
design=as.factor(colnames(liver_gly)[-1])
##構(gòu)建主程序
dds<-DESeqDataSetFromMatrix(data,DataFrame(design),design = ~design)
dds<-DESeq(dds,fitType = "local") ## or mean
res<-as.data.frame(results(dds))
res
接下來焕济，把fold change >2 <0.5 同時 P<0.05的篩選出來纷妆，一般建議選擇p adjust值，這樣會降低假陽性晴弃。
library(dplyr)
res<-cbind(rownames(res),res)
res_liver<-res %>% dplyr::filter((log2FoldChange>1 | log2FoldChange < (-1)) & padj < 0.05)
二 火山圖
library("labeling")
library(ggplot2)
#將符合要求的篩出來
threshold<-as.factor(
(res$log2FoldChange>1|res$log2FoldChange<(-1) &
res$pvalue<0.05 ))
pdf("/Users/baiyunfan/desktop/liver.pdf")
ggplot(res,aes(x= log2FoldChange,
y= -1*log10(res$pval),colour=threshold))+xlab("log2 fold-change")+ylab("-log10 p-value")+geom_point()
dev.off()
