差异基因表达 火山图 python fpkm 基因差异表达分析软件

edgeR是一个研究重复计数数据差异表达的Bioconductor软件包。方法：基于负二项分布的统计方法，包括经验贝叶斯估计、精确检验、广义线性模型和准似然检验。应用：与RNA-seq一样，edgeR包也可用于其他测序数据，包括ChIP-seq、ATAC-seq、亚硫酸氢盐seq、SAGE和CAGE。

1.读入数据到DGEList对象

library(edgeR)
# 下载GSE63310 Supplementary file中的文件，解压

setwd(project_path)
files <- list.files(".","^GSM")  # 匹配以GSM开头的文件，mock


# 读入多个文件中的表达谱数据
dgelist<-readDGE(files,columns=c(1,3))
# class(dgelist) #  DGEList object 
# 也可以通过DGEList()得到DGEList object 

# 含有 dgelist$samples （"data.frame"）和 dgelist$counts （"matrix" "array"）
# dim(dgelist$counts)
# dgelist$samples group 需要重新定义 norm.factors 需要重新计算
# dgelist$counts 需要进一步处理：消除库大小，计算cpm等。
ncol(dgelist)  # 样品数
nrow(dgelist) # 基因数

2. 标准化

## 2.标准化。
# TMM:trimmed mean of M-values
# dgelist$sample$norm.factors值改变
dgelist <- calcNormFactors(dgelist, method = 'TMM')

# cpm <-cpm(dgelist)  # 没有log转化, cpm为表达矩阵。
# rpkm <- rpkm(dgelist)
# lcpm <-cpm(dgelist,log=T)  # log转化
# lrpkm <-rpkm(dgelist,log=T)  # log转化

3. 过滤count数低的基因

经验设置为cpm=1位为cutoff点。但是，这并不是最精准。因为随着测序深度增加，例如20million（2 千万），cpm=1 意味着 counts=20。阈值可能会有点高。测序深度低的话，例如2million（2百万），cpm=1 意味着counts=1。阈值可能会太低。此时可以使用自动过滤，或者根据cut.off.cpm=10/来计算。

# 过滤标准自己定。
# 至少2个样品的reads数大于1
keep <- rowSums(cpm(dgelist) > 1 ) >= 2 
dgelist <- dgelist[keep,]

4. 转录组数据聚类可视化

# 图形展示
lcpm <-cpm(dgelist,log=T)
plotMDS(lcpm,labels=group,col=c('red','blue'))

5. 建立分组矩阵

# 根据分组信息构建试验设计矩阵
# 分组信息中一定要是对照组在前，处理组在后？？？
group <- c(rep('control',5),rep('treat',6)) # mock

design <- model.matrix(~group)
class(design)  # [1] "matrix" "array"

6.模型拟合

edgrR涉及到差异表达分析的函数有很多： exactTest、glmFit、glmLRT、glmQLFit、glmQLFTest。 qCML估计离散度需要搭配 exact test 进行差异表达分析，对应 exactTest 函数。 而其他四个glm*都是与GLM模型搭配使用的函数。其中，glmFit 和 glmLRT 函数是配对使用的，用于 likelihood ratio test (似然比检验)，而 glmQLFit和 glmQLFTest则配对用于 quasi-likelihood F test (拟极大似然F检验)。 有两个方法glmQLFit()和glmQLFTest(),在两个检验方法中，首选QLFit，因为它反映了估计每个基因的离散度时的不确定性。当重复次数较少时，它提供了更强大和可靠的错误率控制

#估算基因表达值的离散度
dge <- estimateDisp(dgelist, design, robust = TRUE)

#负二项广义对数线性模型
fit <- glmFit(dge, design, robust = TRUE)
# Likelihood ratio tests for trend
results <- glmLRT(fit, coef=2)
# Extracts the most differentially expressed genes 
lrt <- topTags(results,n=nrow(dge$counts)) # FDR排序，所有基因
# lrt <- topTags(results)
# lrt <- topTags(results,n=100) # 前100基因
# 写入文件
write.table(lrt, 'control_treat.glmLRT.txt', sep = '\t', col.names = NA, quote = FALSE)

#拟似然负二项广义对数线性模型
fit <- glmQLFit(dge, design, robust = TRUE)
lrt <- topTags(glmQLFTest(fit), n = nrow(dge$counts))
write.table(lrt, 'control_treat.glmQLFit.txt', sep = '\t', col.names = NA, quote = FALSE)