我的R语言版本是3.6.1
BiocManager::install("DESeq2")
使用library(DESeq2)
加载的时候遇到报错 :载入了名字空间‘rlang’ 0.4.0,但需要的是>= 0.4.2
解决办法:将rlang包手动删除,rlang所在的路径是\R-3.6.1\library\rlang。然后使用命令install.packages("rlang")
重新安装就可以了
BiocManager::install("pasilla")
library(pasilla)
pasCts<-system.file("extdata",
"pasilla_gene_counts.tsv",
package="pasilla",
mustWork = T)
pasCts
pasAnno<-system.file("extdata",
"pasilla_sample_annotation.csv",
package="pasilla",
mustWork=T)
pasAnno
df<-read.csv(pasCts,sep="\t",row.names = "gene_id")
cts<-as.matrix(df)
head(cts)
coldata<-read.csv(pasAnno,row.names = 1)
head(coldata)
coldata<-coldata[,c("condition","type")]
rownames(coldata)<-sub("fb","",rownames(coldata))
all(rownames(coldata) %in% colnames(cts))
cts<-cts[,rownames(coldata)]
all(rownames(coldata) == colnames(cts))
获得了两个数据集
> coldata
condition type
treated1 treated single-read
treated2 treated paired-end
treated3 treated paired-end
untreated1 untreated single-read
untreated2 untreated single-read
untreated3 untreated paired-end
untreated4 untreated paired-end
> head(cts)
treated1 treated2 treated3 untreated1 untreated2 untreated3
FBgn0000003 0 0 1 0 0 0
FBgn0000008 140 88 70 92 161 76
FBgn0000014 4 0 0 5 1 0
FBgn0000015 1 0 0 0 2 1
FBgn0000017 6205 3072 3334 4664 8714 3564
FBgn0000018 722 299 308 583 761 245
untreated4
FBgn0000003 0
FBgn0000008 70
FBgn0000014 0
FBgn0000015 2
FBgn0000017 3150
FBgn0000018 310
library(DESeq2)
dds<-DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design= ~ condition)
dds
keep<-rowSums(counts(dds))>=10
dds<-dds[keep,]
dds<-DESeq(dds)
res<-results(dds,contrast=c("condition","treated","untreated"))
resdata<-merge(as.data.frame(res),
as.data.frame(counts(dds,normalized=T)),
by="row.names",sort=F)
DEG<-resdata
logFC_cutoff<-with(DEG,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
DEG$change<-as.factor(ifelse(DEG$pvalue<0.05&abs(DEG$log2FoldChange)>logFC_cutoff,
ifelse(DEG$log2FoldChange>logFC_cutoff,"UP","DOWN"),
"NOT"))
this_title <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
DEG<-na.omit(DEG)
library(ggplot2)
ggplot(data=DEG,aes(x=log2FoldChange,y=-log10(pvalue),color=change))+
geom_point(alpha=0.4,size=1.75)+
labs(x="log2 fold change")+ ylab("-log10 pvalue")+
ggtitle(this_title)+theme_bw(base_size = 20)+
theme(plot.title = element_text(size=15,hjust=0.5))+
scale_color_manual(values=c('blue','black','red'))
image.png
绘制火山图学到的新的ggplot2知识点
theme_bw(base_size=20)
改变了图片中整体字体的大小
theme(plot.title = element_text(size=15,hjust=0.5))
将绘图标题调整到了中间