RNA-seq多Run合并、VST标准化、PCA、差异分析
ggplot2 PDF输出字体问题
下面是把代码包装成单独的函数，适用于分析多组差异分析
rm(list = ls())
alt3<-read.table("/Your_dir/alt3.events.quant.tsv",header = T,row.names = 1,sep = "\t")
alt5<-read.table("/Your_dir/alt5.events.quant.tsv",header = T,row.names = 1,sep = "\t")
es<-read.table("/Your_dir/es.events.quant.tsv",header = T,row.names = 1,sep = "\t")
ir<-read.table("/Your_dir/ir.events.quant.tsv",header = T,row.names = 1,sep = "\t")
#get the count expression matrix only
alt3<-alt3[,2:7]
alt5<-alt5[,2:7]
es<-es[,2:7]
ir<-ir[,2:7]
DEG<-function(x){
suppressMessages(library(DESeq2, quietly=T))
suppressMessages(library(ggplot2))
suppressMessages(library(ggrepel))
FC=1
pvalue=0.05
exprSet<-get(x)
group_list<-c(rep("Fed",3),rep("Fasting",3))
colData <- data.frame(row.names=colnames(exprSet), group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData=exprSet,
colData=colData,
design=~group_list)
suppressMessages(dds2 <- DESeq(dds))
res <-results(dds2, contrast=c("group_list","Fasting","Fed"))
resOrdered<-res[order(res$pvalue),]
resOrdered<-as.data.frame(resOrdered)
sigoutTab<-resOrdered[abs(resOrdered$log2FoldChange)>FC & resOrdered$pvalue<pvalue,]
all_deg<-paste0("/Your_dir/DEGs/Human_diffSplice_",as.character(x),"_DESeq2.csv")
sig_deg<-paste0("/Your_dir/DEGs/Human_diffSplice_",as.character(x),"_logFC_",as.character(FC),"_pvalue_",as.character(pvalue),"_DESeq2.csv")
write.csv(resOrdered,all_deg)
write.csv(sigoutTab,sig_deg)
vsd <- vst(dds2, blind=FALSE)
pcaData <- plotPCA(vsd, intgroup=c('group_list'), returnData=T)
xmin <- min(pcaData$PC1)
xmax <- max(pcaData$PC1)
xmin <- xmin - 0.3*(xmax - xmin)
g<-ggplot(pcaData, aes(PC1, PC2, color=group_list, shape=group_list)) +
geom_text_repel(aes(label=name)) +
geom_point(size=3) +
xlim(xmin, NA) +
theme_minimal()
ggsave(paste0("/Your_dir/DEGs/Human_diffSplice_",as.character(x),"_PCA_plot.pdf"),g,units = "in",width=5.5,height = 4.5)
}
analysis_list<-c("alt3","alt5","es","ir")
sapply(analysis_list,DEG)
#####Volcano plot
rm(list = ls())
library(ggplot2)
library(extrafont)
font_import()
fonts()
loadfonts()
alt3<-read.csv("/Your_dir/DEGs/Human_diffSplice_alt3_DESeq2.csv",header = T)
alt5<-read.csv("/Your_dir/DEGs/Human_diffSplice_alt5_DESeq2.csv",header = T)
es<-read.csv("/Your_dir/DEGs/Human_diffSplice_es_DESeq2.csv",header = T)
ir<-read.csv("/Your_dir/DEGs/Human_diffSplice_ir_DESeq2.csv",header = T)
analysis_list<-c("alt3","alt5","es","ir")
volcano<-function(x){
suppressMessages(library(DESeq2, quietly=T))
suppressMessages(library(ggplot2))
suppressMessages(library(ggrepel))
FC=1
pvalue=0.05
Degs<-get(x)
Degs$sigORnot = as.factor(ifelse(Degs$pvalue < pvalue & abs(Degs$log2FoldChange) > FC,
ifelse(Degs$log2FoldChange > FC ,'UP','DOWN'),'NOT'))
g <-ggplot(data=Degs, aes(x=log2FoldChange, y=-log10(Degs[,"pvalue"]), color=sigORnot)) +
geom_point(alpha=0.4, size=2) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("Log2(fold change)") + ylab(paste0("-Log10(ajust p value)")) +
scale_colour_manual(values = c('blue','grey','red'))+geom_hline(yintercept=(-log10(pvalue)),linetype=3)+geom_vline(xintercept=c(-(FC),FC),linetype=3)
g<-g+theme_classic()+
theme(axis.title.x = element_text(size=20,family="Arial",face="bold",colour = "black"),
axis.text.x = element_text(size=20,family="Arial",face="bold",colour = "black"),
axis.title.y = element_text(family="Arial",size=20,colour = "black",face = "bold"),
axis.text.y= element_text(size=20,family="Arial",face="bold",colour = "black"))
g<-g+theme(axis.line = element_line(size=1))+theme(axis.ticks = element_line(size=1))
g<-g+theme(legend.text=element_text(size=16,family="Arial",face="bold",colour = "black"),
legend.title=element_text(size=20,family="Arial",face="bold",colour = "black"))
g<-g+xlim(-max(abs(Degs$log2FoldChange)),max(abs(Degs$log2FoldChange)))
ggsave(paste0("/Your_dir/DEGs/volcano/Human_diffSplice_",as.character(x),"_logFC_",as.character(FC),"_pvalue_",as.character(pvalue),"_Volcano_plot.pdf"),g,units = "in",width=5.5,height = 4.5)
}
sapply(analysis_list,volcano)
网友评论