一行命令完成RNAseq差异分析+火山图+散点图

2020-08-06 本文已影响0人邓老师呦

今天来给大家分享的是：怎么一行命令完成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服务器窗口，运行命令