###DESeq2 DEseq2针对有生物学重复的样本。
###数据来源https://4va.github.io/biodatasci/r-rnaseq-airway.html#data_needed
metadata <- read.csv("airway_metadata.csv")
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData=counts,
colData=metadata,
design=~dex,
tidy=TRUE)
###其中countData存放counts数据，colData存放样本信息的数据，design就是实验设计。
#####数据预处理
dds1 <- dds [ rowSums(counts(dds)) > 1, ]
##方差齐性转换，rlog更适用于样本量小于30的数据，相应的vst负责的是样本量大于30的数据
rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)
sampleDists <- dist(t(assay(rld)))##样本距离，查看样本之间的相似度和差异程度
library("pheatmap")
library("RColorBrewer")#自动生成一组颜色
sampleDistMatrix <- as.matrix( sampleDists )
colors <- colorRampPalette( rev(brewer.pal(8, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors)#聚类的作用是判断样本之间的相似度
##PCA分析
#还有一种可视化样本-样本距离的方法就是主成分分析。
##DESeq2提供了专门的方法用于作图
plotPCA(rld,intgroup=c("dex"))
########差异分析
dds <- DESeq(dds)#包含了以下三个函数在内
#dds <- estimateSizeFactors(dds)
#计算归一化系数sizeFactor
##dds <- estimateDispersions(dds)
##估计基因的离散程度
##dds <- nbinomWaldTest(dds)
##统计检验，差异分析
##获得分析结果
res1 <- results(dds,tidy=TRUE)#tidy为TRUE是数据框，否则为S4对象
metadata$dex=c("a","a","a","b","b","c","c","c")#重新定义分组,3个组
dds <- DESeqDataSetFromMatrix(countData=counts,
colData=metadata,
design=~dex,
tidy=TRUE)
dds <- DESeq(dds)
resultsNames(dds)#两两对比，分别查看结果
res1 <- results(dds, name="dex_c_vs_a",tidy = TRUE)
res2 <- results(dds, name="dex_b_vs_a")
res3 <- results(dds, name="Intercept",tidy=TRUE)
##结果的简单统计
summary(res1)#summary函数对于数据框和S4对象的统计结果不同
summary(res2) #这种总结更直白
#results函数中可以修改阈值
res.05 <- results(dds, alpha=0.05)
#alpha为padj的阈值，默认padj=0.1。
resLFC <- results(dds, lfcThreshold=1)
#提升log2 fold change threshold，结果中不满足lfc阈值的gene的p值都是1。
#############################################################################
##值得一提的是DESeq2软件独有的normlization方法！相当于一种归一化处理，主要用于后续画图
dds1 <- rlogTransformation(dds) ## 得到经过DESeq2软件normlization的表达矩阵,counts表达矩阵不能直接做图
exprSet_new=assay(dds1)
par(cex = 0.7)#指定符号的大小
n.sample=ncol(counts)#样本数
cols <- rainbow(n.sample*1.2)#生成颜色
colnames(counts)
c1=column_to_rownames(counts,"ensgene")#把ensgene列作为行名
boxplot(c1, col = cols,main="expression value",las=2)#离散程度较高
boxplot(exprSet_new, col = cols,main="expression value",las=2)#归一化后结果才有意义
##res=arrange(res,padj)
anno <- read.csv("annotables_grch38.csv")#探针注释表格
##res_anno = res %>%inner_join(anno,by=c("row"="ensgene"))
res_anno=inner_join(res3,anno,by=c("row"="ensgene"))#注意res3为数据框结构
res_anno %>% filter(padj<0.05) #筛选padj小于0.05的值
res_anno %>%
filter(padj<0.05) %>%
write.csv("sigresults.csv")#通道函数省去了来回赋值，前一项结果传递到后一项
SPARCL1 <- res_anno %>% filter(symbol=="SPARCL1")
View(SPARCL1)#SPARCL1对的基因ID是ENSG00000103196
plotCounts(dds, gene="ENSG00000152583", intgroup="dex")
###值得注意的是plotCount返回的不仅仅是图，也可以是数据
plotCounts(dds, gene="ENSG00000152583", intgroup="dex", returnData=TRUE)
library(ggplot2)
plotCounts(dds, gene="ENSG00000152583", intgroup="dex", returnData=TRUE) %>%
ggplot(aes(dex, count)) +
geom_boxplot(aes(fill=dex)) +
scale_y_log10() +
ggtitle("SPARCL1")
library(edgeR)
##，使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。
colnames(counts)
count1=column_to_rownames(counts,"ensgene")#同样不能使用FPKM值，只能使用count表达矩阵
###排序
count2=count1[,c(1,3,5,7,2,4,6,8)]
group.list=c(rep("control",4),rep("treated",4))
group.list=factor(group.list)
group.list=relevel(group.list,ref = "control")
###第一步： 构建DGEList对象
dge <- DGEList(counts=count2,group=group.list)
####第二步： 过滤 low counts数据。与DESeq2的预过滤不同，DESeq2的预过滤只是为了改善后续运算性能，
#这一步可以省略不做，在结果处过滤也是可以的
# 简单粗暴的方法 基因至少在某一些文库的count超过10 ~ 15 才被认为是表达。
keep <- rowSums(dge$counts) > 50
# 利用CPM标准化 0.5(即阈值）等于 10/(最小的文库的 read count数 /1000000)
keep <- rowSums(cpm(dge) > 0.5 ) >=2
table(keep)
###第三步： 根据组成偏好(composition bias)标准化
#利用库的大小来进行标准化TMM
dge <- calcNormFactors(dge)
?calcNormFactors
##TMM <- calcNormFactors(d, method="TMM")
##TMM <- calcNormFactors(d, method="RLE")
##TMM <- calcNormFactors(d, method="upperquartile")
plotMDS(dge, col = rep(c('red', 'blue'), each = 4), dim = c(1, 2))
#第四步： 实验设计矩阵(Design matrix) 建立分组变量
design <- model.matrix(~0+group.list)
design1 <- model.matrix(~group.list)
##~0+group :不包括比较矩阵
##~group:包括了比较矩阵
rownames(design)<-colnames(dge)
colnames(design)<-levels(group.list)
#第五步： 估计离散值（Dispersion）。前面已经提到负二项分布（negative binomial，NB)需要均值和离散值两个参数。
#edgeR对每个基因都估测一个经验贝叶斯稳健离散值（mpirical Bayes moderated dispersion），
#还有一个公共离散值（common dispersion，所有基因的经验贝叶斯稳健离散值的均值）以及一个趋势离散值
##install.packages("statmod")
dge <- estimateDisp(dge, design, robust = TRUE)
## estimateDisp()实际上是个组合函数，可以一步得到多个计算结果
##dge <- estimateGLMCommonDisp(dge,design)
##dge <- estimateGLMTrendedDisp(dge, design)
##dge <- estimateGLMTagwiseDisp(dge, design)
plotBCV(dge)
####第六步 差异分析
##负二项式广义对数线性模型
fit <- glmQLFit(dge, design, robust=TRUE) #拟合模型
cntr.vs.T <- makeContrasts(control-treated, levels=design)
res <- glmQLFTest(fit, contrast=cntr.vs.T) #统计检验
##用有默认截距的方法
#fit <- glmFit(dge, design1, robust = TRUE)
#lrt <- glmLRT(fit)
edger.deg=topTags(res, n = nrow(res$counts))
dge_de <- decideTestsDGE(res, adjust.method = 'fdr', p.value = 0.05)
#查看默认方法获得的差异基因
summary(dge_de)
plotMD(res, status = dge_de, values = c(1, -1), col = c('blue', 'red'))
#作图观测
abline(h = c(-1, 1), col = 'gray', lty = 2)
##############类似然负二项式广义对数线性模型
###这里用的是glmQLFTest是因为前面用了glmQLTFit进行拟合，所以需要用QL F-test进行检验。
#ig.edger <- res$table[p.adjust(res$table$PValue, method = "BH") < 0.01, ]
##如果前面用的是glmFit，那么对应的就是glmLRT. 作者称QL F-test更加严格。多重试验矫正用的也是BH方法。
fit2 <- glmFit(dge, design)
#利用似然比检验来检验差异表达基因
fit2 <- glmLRT(fit2, contrast=c(-1,1))
#输出靠前的差异表达基因
DEG=topTags(fit2, n=nrow(exp))
DEG=as.data.frame(DEG)
########################################################
##使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。
##原始数据尺度转换，也是进行normalization操作
cpm <- cpm(count1)
lcpm <- cpm(count1, log=TRUE, prior.count=2)
library(limma)
##进行差异分析时常用limma。虽然它是针对芯片数据开发的，但也有limma-voom可以分析转录组数据
dge <- DGEList(counts=count2)
group.list=c(rep("control",4),rep("treated",4))
group.list=factor(group.list)
group.list=relevel(group.list,ref = "control")
##分组矩阵 (design) ：就是将表达矩阵的列（各个样本）分成几组（例如最简单的case - control，或者一些时间序列的样本day0, day1, day2 ...）【通过model.matrix()得到】
design <- model.matrix(~0+group.list)
rownames(design)<-colnames(dge)
colnames(design)<-levels(group.list)
### 注意： calcNormFactors并不会标准化数据，只是计算标准化因子
dge <- calcNormFactors(dge)
v <- voom(dge,design, normalize="quantile")
##voom到底做了什么转换？
首先原始counts转换成log2的CPM（counts per million reads ），
这里的per million reads是根据之前calcNormFactors计算的norm.factors进行规定的；
然后根据每个基因的log2CPM制作了线性模型，并计算了残差 ；
然后利用了平均表达量（红线）拟合了sqrt(residual standard deviation)；
最后得到的平滑曲线可以用来得到每个基因和样本的权重
##lmFit()：线性拟合模型构建
fit <- lmFit(v, design)
constrasts = paste(rev(levels(group.list)),collapse = "-")
##比较矩阵（contrast）：意思就是如何指定函数去进行组间比较【通过makeContrasts()得到】
cont.matrix <- makeContrasts(contrasts=constrasts,levels = design)
### 比较每个基因
fit2=contrasts.fit(fit,cont.matrix)
###eBayes()使用trend=TRUE对标准误差进行经验贝叶斯平滑，计算每个对比中每个基因的moderated t-statistic和log-odds。
fit2=eBayes(fit2)
##topTable()给出一个最有可能在给定对比下差异表达的基因列表。
###使用coef参数，这里设为1，也就是表示根据上面makeContrasts的第一个（group2-group1）来提取结果
DEG = topTable(fit2, coef=constrasts,sort.by = "P", n=Inf)
##adjust="BH"表示对p值的校正方法，包括了："none", "BH", "BY" and "holm"。
###那么为啥要对P值进行校正呢？
#p值是针对单次检验，假设采用的p值为小于0.05，我们通常认为这个基因在两组样本中的表达是有显著差异的，但是仍旧有5%的概率表示这个基因并不是差异基因。
#但是，当两组样本中有20000个基因采用同样的检验方式进行统计检验时，就会遇到一个问题：单次犯错的概率为0.05， 如果进行20000次检验，那么就会有0.05*20000=1000 个基因在组间的差异被错误估计
# topTable 列出差异显著基因
DEG2=topTable(fit2, coef=1, adjust="BH")
##去掉那些NA值
DEG = na.omit(DEG)
edgeR 使用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因。特别地，经验贝叶斯用于通过在基因之间来调节跨基因的过度离散程度。 使用类似于Fisher精确检验但适应过度分散数据的精确检验用于评估每个基因的差异表达。edgeR 在默认情况下，执行TMM标准化程序以考虑样本之间的不同测序深度，通过Benjamini-Hochberg用于控制FDR 。
Limma包基于线性模型建模。 它最初设计用于分析微阵列数据，但最近已扩展到RNA-seq数据。 根据limma用户指南的当前建议是使用edgeR包的TMM标准化和“voom”转换，其本质上将标准化数据取对数（基数2）并估计它们的均值 - 方差关系以确定在线性建模之前每次观察的权重。 默认情况下，Benjamini-Hochberg程序用于估计FDR 。
DESeq使用类似于edgeR的负二项式模型，与edgeR类似，执行缩放因子归一化以考虑不同样本的变化的测序深度，并且Benjamini-Hochberg用于控制FDR。 DESeq能够分析具有少量重复的实验。DESeq技术上可以在没有任何生物学重复的情况下进行实验。DESeq2是在DESeq基础上更新的软件。
这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
