表达差异基因分析

> # <差异基因分析>
> # 1.判断是否有BiocManager包，若不存在则安装
> options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))) #设置清华镜像，加速下载
> 
> if (!requireNamespace("BiocManager", quietly = TRUE))
+   install.packages("BiocManager")
> if (!requireNamespace("DESeq2", quietly = TRUE))
+   BiocManager::install('DESeq2')  #通过BiocManager安装DESeq2 
> library(DESeq2) #加载library

> #输入数据要求
> # DEseq2要求输入数据是由整数组成的矩阵
> # DESeq2要求矩阵是没有标准化的，一定记住用readcount
> 
> ##2.读入所有基因原始readscount表达矩阵，行为基因，列为样品
> A <- read.table(p, header = T, row.names = 1)
> B <- as.matrix(A) #转换成矩阵格式，保证都是数值

> coldata <- read.table("/home/shijm/Rlearning/R-Online-learning/data/sample_info.txt",header = T,row.names = 1)
> coldata <- coldata[, c("condition", "type")]

> dds <- DESeqDataSetFromMatrix(countData = B, colData = coldata, design = ~ condition);
> 
> # countData = B，readscount矩阵
> # colData = coldata,分组信息，根据这个才能在2组之间比较
> # design = ~ condition，公式，表示按照condition进行分析

> dds <- DESeq(dds)

res <- results(dds, contrast = c("condition", "treated", "untreated")) 
> resordered <- res[order(res$padj),]	 #按照矫正后的P-value padj从小到大排序
> sum(res$padj < 0.05, na.rm = TRUE)	#总结padj小于0.05显著差异的基因
[1] 356

plotMA(res)	#画火山图，横轴是标准化后的平均readscount，纵轴是差异倍数，大于0是上调，小于0是下调,蓝色点表示显著差异的基因

> pvalue<-res$padj
> pvalue<-as.matrix(pvalue)
> fc<-res$log2FoldChange
> fc<-as.matrix(fc)
> df<-data.frame(pvalue,fc)
> df<-na.omit(df,cols=c("pvalue","fc"))   这里去掉pvalue或fc列有NA值得行
> df$threshold[df$pvalue < 0.05 & df$fc>0 ] = "up"
> df$threshold[df$pvalue < 0.05 & df$fc<0 ] = "down"
> df$threshold[df$pvalue > 0.05 & df$fc>=0|df$fc<=0 ] = "non"
>

> ggplot(df,aes(x=fc,y=-log10(pvalue)))+
+   xlab("log2 Fold Change")+ylab("-log10P-Value")+
+   
+   geom_point(aes(color=threshold),size=1,alpha=0.6,)+
+   
+   scale_color_manual(values=c("red","black","blue")) +
+   
+   geom_hline(aes(yintercept=-log10(0.05)),colour="grey",size=1.2 ,linetype=2) + #增加水平间隔线
+   
+   geom_vline(aes(xintercept=0), colour="grey",size=1.2 ,linetype=2)+ #增加垂直间隔线
+   
+   theme_few()+theme(legend.title = element_blank()) #去掉网格背景和图注标签