#首先升级R
install.packages("installr")
library(installr)
updateR()
library("BiocManager")
#Bioconductor version 3.12 (BiocManager 1.30.15), R 4.0.5 (2021-03-31) Bioconductor version '3.12' is out-of-date; the current release version '3.13' is available with R version '4.1'; see https://bioconductor.org/install
#重新更新Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.13")
#更新成功
library("BiocManager")
#通过BiocManager安装我们常用的R包
BiocManager::install("ggplot2")
BiocManager::install(c("ggplot2","ggtree","DESeq2"))
#一次安装多个包
BiocManager::install("DESeq2")
#错误: 无法载入程辑包‘GenomeInfoDb’
BiocManager::install("GenomeInfoDbData")
#测试软件DESeq2
library("DESeq2")
#开始，DESeq2包分析差异表达基因简单来说只有三步：构建dds矩阵，标准化，以及进行差异分析。
library(DESeq2) #加载包
library(apeglm)
#Error in library(apeglm) : 不存在叫‘apeglm’这个名字的程辑包
BiocManager::install("apeglm")
####测试
library(DESeq2) #加载包
library(apeglm)
安装成功。
以下为差异基因表达分析的相关脚本：
pvalue(pval): 统计学差异显著性检验指标。
qvalue(p-adjusted): 校正后的p值（qvalue=padj=FDR=Corrected p-Value=p-adjusted），是对p值进行了多重假设检验，能更好地控制假阳性率。
校正后的p值不同的几种表现形式，都是基于BH的方法进行多重假设检验得到的。校正后的p值不同的展现形式是因为不同的分析软件产生的。
作者：村长吃火锅
链接：https://www.jianshu.com/p/b1f6b50fde0e
来源：简书 著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0/
E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_48 vs 0/
E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_48 vs 24/
致病疫霉
###############################################################################################
library(DESeq2) #加载包
library(apeglm)
data.pi_Sae24 = read.table("E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0/pi_Sae24_count.csv",sep = ',',header = T,row.names = 1)
data.pi_Sly24 = read.table("E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0/pi_Sly24_count.csv",sep = ',',header = T,row.names = 1)
data.pi_Sme24 = read.table("E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0/pi_Sme24_count.csv",sep = ',',header = T,row.names = 1)
data.pi_Stu24 = read.table("E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0/pi_Stu24_count2.csv",sep = ',',header = T,row.names = 1)
condition.pi_Sae24 <- factor(c(rep('0',3),rep('24',5)),levels = c("0","24"))
condition.pi_Sly24 <- factor(c(rep('0',3),rep('24',5)),levels = c("0","24"))
condition.pi_Sme24 <- factor(c(rep('0',3),rep('24',5)),levels = c("0","24"))
condition.pi_Stu24 <- factor(c(rep('0',3),rep('24',5)),levels = c("0","24"))
colData.pi_Sae24 = data.frame(row.names= colnames(data.pi_Sae24),condition.pi_Sae24)
colData.pi_Sly24 = data.frame(row.names= colnames(data.pi_Sly24),condition.pi_Sly24)
colData.pi_Sme24 = data.frame(row.names= colnames(data.pi_Sme24),condition.pi_Sme24)
colData.pi_Stu24 = data.frame(row.names= colnames(data.pi_Stu24),condition.pi_Stu24)
dds.pi_Sae24 <- DESeqDataSetFromMatrix(data.pi_Sae24,colData.pi_Sae24,design = ~ condition.pi_Sae24)
dds.pi_Sly24 <- DESeqDataSetFromMatrix(data.pi_Sly24,colData.pi_Sly24,design = ~ condition.pi_Sly24)
dds.pi_Sme24 <- DESeqDataSetFromMatrix(data.pi_Sme24,colData.pi_Sme24,design = ~ condition.pi_Sme24)
dds.pi_Stu24 <- DESeqDataSetFromMatrix(data.pi_Stu24,colData.pi_Stu24,design = ~ condition.pi_Stu24)
dds.pi_Sae24 <-DESeq(dds.pi_Sae24)
resultsNames(dds.pi_Sae24)
dds.pi_Sly24 <-DESeq(dds.pi_Sly24)
resultsNames(dds.pi_Sly24)
dds.pi_Sme24 <-DESeq(dds.pi_Sme24)
resultsNames(dds.pi_Sme24)
dds.pi_Stu24 <-DESeq(dds.pi_Stu24)
resultsNames(dds.pi_Stu24)
res.pi_Sae24.24.0 = results(dds.pi_Sae24,contrast=c("condition.pi_Sae24","24","0"))
res.pi_Sly24.24.0 = results(dds.pi_Sly24,contrast=c("condition.pi_Sly24","24","0"))
res.pi_Sme24.24.0 = results(dds.pi_Sme24,contrast=c("condition.pi_Sme24","24","0"))
res.pi_Stu24.24.0 = results(dds.pi_Stu24,contrast=c("condition.pi_Stu24","24","0"))
order.res.pi_Sae24.24.0 <- res.pi_Sae24.24.0[order(res.pi_Sae24.24.0$pvalue),]
order.res.pi_Sly24.24.0 <- res.pi_Sly24.24.0[order(res.pi_Sly24.24.0$pvalue),]
order.res.pi_Sme24.24.0 <- res.pi_Sme24.24.0[order(res.pi_Sme24.24.0$pvalue),]
order.res.pi_Stu24.24.0 <- res.pi_Stu24.24.0[order(res.pi_Stu24.24.0$pvalue),]
setwd("E:/SolPanTP/RawData/诺禾转录组测序（310）/05 数据分析/DESeq2/log2FC+p value/pi/pi_24 vs 0")
write.csv(order.res.pi_Sae24.24.0,file = "order.res.pi_Sae24.24.0.csv",quote = F)
write.csv(order.res.pi_Sly24.24.0,file = "order.res.pi_Sly24.24.0.csv",quote = F)
write.csv(order.res.pi_Sme24.24.0,file = "order.res.pi_Sme24.24.0.csv",quote = F)
write.csv(order.res.pi_Stu24.24.0,file = "order.res.pi_Stu24.24.0.csv",quote = F)
diffgene_pi_Sae24.24.0_deseq2 <- subset(order.res.pi_Sae24.24.0,padj < 0.01 & abs(log2FoldChange) >2)
diffgene_pi_Sly24.24.0_deseq2 <- subset(order.res.pi_Sly24.24.0,padj < 0.01 & abs(log2FoldChange) >2)
diffgene_pi_Sme24.24.0_deseq2 <- subset(order.res.pi_Sme24.24.0,padj < 0.01 & abs(log2FoldChange) >2)
diffgene_pi_Stu24.24.0_deseq2 <- subset(order.res.pi_Stu24.24.0,padj < 0.01 & abs(log2FoldChange) >2)
write.csv(diffgene_pi_Sae24.24.0_deseq2,file = "diffgene_pi_Sae24.24.0_deseq2.csv",quote = F)
write.csv(diffgene_pi_Sly24.24.0_deseq2,file = "diffgene_pi_Sly24.24.0_deseq2.csv",quote = F)
write.csv(diffgene_pi_Sme24.24.0_deseq2,file = "diffgene_pi_Sme24.24.0_deseq2.csv",quote = F)
write.csv(diffgene_pi_Stu24.24.0_deseq2,file = "diffgene_pi_Stu24.24.0_deseq2.csv",quote = F)
网友评论