一行命令完成RNAseq差异分析+画图

今天来给大家分享的是：怎么一行命令完成RNAseq数据差异分析+火山图+散点图。

一、准备3个文件：

1、基因表达count文件：gene_count.txt 格式如下：（行是基因、列是样本）

2、基因表达FPKM文件：gene_exp.txt

3、样本分组文件：group.txt（第一列是样本、第二列是分组名）

二、运行代码

1、确保已经安装了R或RStudio，如果没有安装R包edgeR和ggplot2，可以先安装一下。这里使用edgeR包的exactTest函数进行RNAseq的差异分析。

安装R包代码：

if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install(c('edgeR'))
install.packages("ggplot2")

R包已经装好的忽略上面步骤，直接复制下面代码，另存为 edgeR.R。

library(edgeR)library(ggplot2)#读取文件data <- read.table("gene_count.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.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]