今天来给大家分享的是：怎么一行命令完成RNAseq数据差异分析+火山图+散点图。
一、准备3个文件：
1、基因表达count文件：gene_count_edger.txt
格式如下：（行是基因、列是样本）
image
2、 基因表达FPKM文件：gene_exp_edger.txt
image
3、 样本分组文件：group.txt（第一列是样本、第二列是分组名）
image
二、运行代码
1、确保已经安装了R或RStudio，如果没有安装R包edgeR和ggplot2，可以先安装一下。这里使用edgeR包进行RNAseq的差异分析。
安装R包代码：
if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install(c('edgeR'))
install.packages("ggplot2")
下载代码edgeR.R
library(edgeR)
library(ggplot2)
data <- read.table("gene_count_edger.txt",sep="\t",row.names=1,header=TRUE)
group0 <- read.table("group.txt",sep="\t",row.names=1,header=TRUE,as.is = TRUE)
exp <- read.table('gene_exp_edger.txt', sep = '\t', row.names = 1,header=TRUE)
group <- group0[,'group']
dgelist <- DGEList(counts = data, group = group)
#过滤低表达,CPM标准化
keep <- rowSums(cpm(dgelist) > 1 ) >= 2
dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]
#TMM 标准化
dgelist<- calcNormFactors(dgelist, method = 'TMM')
#估算离散值
dgelist = estimateCommonDisp(dgelist)
dgelist = estimateTagwiseDisp(dgelist)
dgelist = exactTest(dgelist)
tTag = topTags(dgelist,n=NULL)
tTag <- as.data.frame(tTag)
g1 <- unique(group)[1]
g2 <- unique(group)[2]
diff <- tTag[((tTag$logFC >= 1 | tTag$logFC <= -1) & (tTag$FDR<0.05)),]
write.csv(tTag,file = paste(g2,"_vs_",g1,"_edgeR_all.csv",sep=""))
write.csv(diff,file = paste(g2,"_vs_",g1,"_edgeR_different.csv",sep=""))
#画散点图
g1_exp = exp[rownames(tTag),rownames(group0)[which(group0$group==g1)]]
g2_exp = exp[rownames(tTag),rownames(group0)[which(group0$group==g2)]]
g1_mean = apply(g1_exp,1,mean)
g2_mean = apply(g2_exp,1,mean)
type=rep('No',length(g1_mean))
type[which(tTag$logFC > 1 & tTag$FDR < 0.05)] = "Up"
type[which(tTag$logFC < -1 & tTag$FDR < 0.05)] = "Down"
datam = data.frame(g1_mean,g2_mean,logFC=tTag$logFC,FDR=tTag$FDR,type,stringsAsFactors=FALSE)
##散点图
ggplot(datam,aes(log2(g1_mean),log2(g2_mean),colour=type))+
geom_point(stat="identity",size=1)+theme(legend.title=element_blank())+scale_color_manual(values =c("Down"='blue',"No"='grey',"Up"='orange'))+
labs(x=paste(g1,' Log2(FPKM)'),y=paste(g2,' Log2(FPKM)'),title=paste(g2,' VS ',g1,sep=""))+
coord_cartesian(ylim=c(-10,10),xlim=c(-10,10))+geom_segment(aes(x = -10, y = -10, xend = 10, yend = 10),size=1,colour="#999999",linetype="dotted")+theme(plot.title = element_text(hjust = 0.5),title=element_text(face="bold",size=15,colour="black"),axis.title=element_text(face="bold",size=13,colour="black"),axis.text.x=element_text(face="bold",size=12,colour="black"),axis.text.y=element_text(face="bold",size=12,colour="black"),legend.text=element_text(face="bold",size=13,colour="black"))
ggsave("diff_gene_scatter.pdf", width=6, height=6)
ggplot(datam,aes(logFC,-log10(FDR),colour=type))+ geom_point(stat="identity",size=1.2)+theme(legend.title=element_blank())+scale_color_manual(values =c("Down"='blue',"No"='grey',"Up"='orange'))+
labs(x="Log2 (FC)",y="-Lg10 (FDR)",title=paste(g2,' VS ',g1,sep=""))+coord_cartesian(xlim=c(-12,12),ylim=c(0,15))+
geom_hline(aes(yintercept=1.3),colour="white",size=1.1)+
geom_vline(aes(xintercept =-1),colour="white",size=1.1)+geom_vline(aes(xintercept =1),colour="white",size=1.1)+
theme(axis.title.y = element_text(vjust=-0.1),axis.title.x = element_text( vjust=-0.6),title = element_text( vjust=0.8))+theme(plot.title = element_text(hjust = 0.5),title=element_text(face="bold",size=15,colour="black"),axis.title=element_text(face="bold",size=13,colour="black"),axis.text.x=element_text(face="bold",size=10,colour="black"),axis.text.y=element_text(face="bold",size=10,colour="black"),legend.text=element_text(face="bold",size=12,colour="black"))
ggsave("diff_gene_volcano.pdf", width=7, height=6)
2、将前面准备的三个文件和代码放到同一个文件夹，如果自己的文件名跟前面列的文件名不一致需要先改成一样。
3、代码运行
方法一：RStudio直接运行整个edgeR.R文件夹
方法二：打开DOS运行窗口或者linux服务器窗口，运行命令
Rscript edgeR.R
三、结果展示
1、差异分析结果文件：
image
注：文件夹里会生成2个结果文件，一个是所有基因的差异分析结果，一个是按FDR<0.05 & FC 2倍以上筛选的显著差异基因文件。
1、差异分析散点图
image
2、差异分析火山图
image
怎么样，是不是突然觉得原来数据分析可以这么简单！！！
今天的分享就到此吧，下次继续。
image
网友评论