今天来给大家分享的是:怎么一行命令完成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)