- 生信入门（五）——使用DESeq2进行RNA-seq数据分析
- 四、探索性数据分析
- 1、简单EDA
- 2、EDA 的数据转换
- 3、主成分图
- 五、差异数据分析
- 1、标准DE步骤
- 2、最小效应量
- 六、AnnotationHub
- 1、查询AnnotationHub
- 2、映射id
- 七、Building Report
- 1、报告工具
- 2、Glimma
- 八、与ZINB-WaVE的集成
- 1、模拟飞溅
- 2、用飞溅模拟单细胞计数数据
本篇接上一篇，本篇做探索性数据分析，差异表达分析以及后面步骤 四、探索性数据分析 1、简单EDA
#探索性数据分析
#BiocManager::install("airway")
library("airway")
data("airway")
#指定untrt是dex 变量的参考级别
airway$dex<-relevel(airway$dex,"untrt")
airway$dex
#快速检查与基因唯一对齐的数百万个片段（round的第二个参数告诉要保留多少小数点）
round(colSums(assay(airway))/1e6,1)
#通过拉出SummarizedExperiment的colData插槽来检查有关样本的信息：
colData(airway)
table(airway$cell)
table(airway$dex)
#加载DESeq2，创建DESeqDataSet，测试地塞米松治疗之间的差异
library("DESeq2")
dds<-DESeqDataSet(airway,design = ~cell+dex)
#dds
#执行最小过滤以减小数据集的大小
#（如果4个或者更多样本的计数不为5或更多，我们不需要保留基因，因为这些基因将没有检测差异的统计能力，也没有计算样本之间距离的信息）
keep<-rowSums(counts(dds)>=5)>=4
table(keep)
dds<-dds[keep,]
#基本的探索性分析是检查每个样本计数的箱线图（对数据取对数，方便大基数主导箱线图）
boxplot(log10(counts(dds)+1))
dds<-estimateSizeFactors(dds)
#归一化处理
boxplot(log10(counts(dds,normalized=TRUE)+1))
运行结果
取计数的对数加上伪计数1是一种常见的变换，但是他往往会夸大低计数的采取方差，以至于它甚至大于样本组之间的生物变异。因此在DESeq2中，提供了产生对数尺度数据的转换，从而消除了系统趋势。一般推荐的转换是方差稳定转换或VST，他可以用VST函数调用：
VST函数不返回DESeqDataSet，因为其不返回计数，而是返回连续值（在log2尺度上），可以使用分析访问转换后的数据。assay函数转换
#EDA数据转换 （使用VST函数调用转换） vsd<-vst(dds) class(vsd) assay(vsd)[1:3,1:8]#使用VST函数返回的不是计数，而是返回连续值
运行结果
VST数据适用于计算样本之间的距离或执行PCA（principal components analysis）。简而言之，PCA图使我们能够将数据中最主要的变异轴可视化，这对于质量控制和了解样本见差异在不同条件下和条件内的大小都有用。可以看到PC1（数据中变化的主轴）将处理过的和未处理过的样本分开：
#生成主成分图 plotPCA(vsd,"dex")
运行结果
通过额外的ggplot2代码，可以指出哪些样本属于哪个细胞系。
注意：不建议将转换后的数据用于主要差异表达分析。相反，我们将使用原始计数和广义线性模型（GLM），该模型考虑低计数和高计数的预期方差。
# 增加额外GGPLOT2代码，指出哪些样本属于哪个细胞系
library("ggplot2")
pcaData<-plotPCA(vsd,intgroup=c("dex","cell"),returnData=TRUE)
# pcaData
percentVar<-round(100*attr(pcaData,"percentVar"))
ggplot(pcaData,aes(x=PC1,y=PC2,color=dex,shape=cell))+geom_point(size=3)+
xlab(paste0("PC1:",percentVar[1],"%variance"))+
ylab(paste0("PC2:",percentVar[2],"%variance"))+
coord_fixed()
运行结果
# 差异表达分析 # 标准DE步骤 dds<-DESeq(dds) res<-results(dds) # res包含每个基因的结果（与DESeqDataSet中的顺序相同） # 如果想查看顶级基因，可以这样排序 head(res[order(res$pvalue),]) # 绘制顶部基因的计数 plotCounts(dds,which.min(res$pvalue),"dex")
运行结果
使用以下方法检查用于地塞米松治疗导致的所有log2倍数变化——LFC超过计数平均值plotMA
# 使用以下方法检查用于地塞米松治疗导致的所有log2倍数变化（LFC）超过计数平均值plotMA plotMA(res,ylim=c(-5,5))
运行结果
上面的MA绘图左侧有很多不显着的大型LFC（灰色点）。由于对数计数的不精确，这些点获得了很大的LFC 。为了获得更多信息可视化和精确的基因按效应大小排序（对数倍数变化有时称为效应大小）建议使用DESeq2的功能来缩小LFC。——方法apeglm收缩估计器，在DESeq2的lfcShrink函数中可用：
# BiocManager::install("apeglm")
library("apeglm")
resultsNames(dds)
res2<-lfcShrink(dds,coef="dex_trt_vs_untrt",type="apeglm")
par(mfrow=c(1,2))
plotMA(res,ylim=c(-3,3),main="No shrinkage")
plotMA(res2,ylim=c(-3,3),main="apeglm")
运行结果
如果不想报告具有小LFC的重要基因，可以通过选择LFC并对此进行测试来指定最小的生物学意义效应大小。我们可以使用没有shrink的LFC或lfcShrink提供的LFC使用apeglm方法执行这样的阈值测试
注意：测试针对LFC阈 ≠ 测试针对0假设然后过滤上LFC值
apeglm方法提供S值时svalue=TRUE，或者当我们提供最小的影响大小（如上），这些类似于q值或者调整后的p值，因为s值 小于α应该有一个总的误号率或绝对值小于我们会给定的LFC阈值，其界限为α
res.lfc<-results(dds,lfcThreshold = 1) res.lfc2<-lfcShrink(dds,coef = "dex_trt_vs_untrt",type="apeglm",lfcThreshold = 1) par(mfrow=c(1,2)) plotMA(res.lfc,ylim=c(-5,5),main="NO shringkage ,LFC test") plotMA(res.lfc2,ylim=c(-5,5),main="apeglm,LFC test",alpha=0.01)六、AnnotationHub 1、查询AnnotationHub
使用AnnotationHub包将附加信息附加到结果表中。AnnotationHub为超过40000条注释记录提供了一个易于使用的界面。一条记录可以来自ENCODE，人类基因组的序列，一个TXDB含有约转录物和基因，或信息OrgDb含有有关特定生物体的生物标识符的一般信息
# annnotationhub
# BiocManager::install("AnnotationHub",force = TRUE)
library("AnnotationHub")
ah<-AnnotationHub()
使用以下代码启动浏览器，用于导航通过AnnotationHub可用的所有记录
display(ah)
运行结果
还可以使用带有查询功能的关键字进行查询
query(ah,c("OrgDb","Homo sapiens"))
运行结果
# 如果要下拉特定记录，可以使用双括号和记录名称 hs<-ah[['AH92581']] hs
运行结果
# rowbames结果表是ENSEMEL IDs，大部分在OrgDb中（尽管有几千都不在） columns(hs) table(rownames(res)%in%keys(hs,"ENSEMBL")) #可以使用maplds函数添加基因符号，keytype=ENSEMBL,column="SYMBOL" res$symbol<-mapIds(hs,rownames(res),column = "SYMBOL",keytype = "ENSEMBL") head(res)
运行结果
有许多从Bioconductor构建交互报告的软件包。其中两个是Reporting Tools和Glimma，他们都提供HTML报告，允许合作者检查基因组分析中的顶级基因（或任何感兴趣的特征）
编译Reporting Tools报告的代码：
# building report
BiocManager::install("ReportingTools")
library("ReportingTools")
tmp<-tempdir()
rep<-HTMLReport(shortName = "airway",title = "Airway DGE",
basePath = tmp,reportDirectory = "report")
publish(res,rep,dds,n=20,make.plots=TRUE,factor=dds$dex)
finish(rep)
#在浏览器中启动
browseURL(file.path(tmp,"report","airway.html"))
运行结果
另一个可以生成交互式报告的包是Glimma。gIMDPlot 构建了一个交互式MA绘图，其中将鼠标悬停在左侧MA绘图中的基因上将显示右侧样本的计数。但即将在工具提示和屏幕底部的列表中显示基因信息。将鼠标悬停在右侧的样本上将在工具提示中提供样本ID
# 还有Glimma生成交互式报告
library("Glimma")
status<-as.numeric(res$padj<.1)
anno<-data.frame(GeneID=rownames(res),symbol=res$symbol)
glMDPlot(res2,status=status,counts=counts(dds,normalized=TRUE),
groups=dds$dex,transform=FALSE,
samples=colnames(dds),anno=anno,
path=tmp,folder="glimma",launch=TRUE)#当launch=FALSE时，不显示，可以结合下面语句显示，若为TRUE 则不需要下面语句
#browseURL(file.path(tmp,"glimma","MD-Plot.html"))
运行结果
在最后一小节，展示DESeq2与另一个Bioconductor包zinbwave集成，以便对额外的零进行建模和解释（超过负二项式模型的预期）——对于单细胞RNA-seq实验很有用
此处，使用splatter包来模拟单细胞RNA-seq数据（Zappia,Phipson和Oshlack）
**注：虽然 ZINB-WaVE 和 ZINGER 等方法可以成功识别多余的零点，但它们不能轻易区分其根本原因，即技术（例如，丢失）和生物（例如，爆裂）零点.
**
当感兴趣的信号不在零分量中时，可以使用下面概述的零通胀加权方法。也就是说，如果您想找到跨细胞群转录爆发的生物学差异，以下方法将无法帮助您找到这些差异。相反，它有助于发现除零分量之外的计数差异（这些零是生物的还是技术的）。
# BiocManager::install("splatter")
library("splatter")
params<-newSplatParams()
params<-setParam(params,"de.facLoc",1)
params<-setParam(params,"de.facScale",.25)
params<-setParam(params,"dropout.type","experiment")
params<-setParam(params,"dropout.mid",3)
set.seed(1)
sim<-splatSimulate(params,group.prob=c(.5,.5),method = "groups")
# sim
plot(log10(rowMeans(assays(sim)[["TrueCounts"]])),
rowMeans(assays(sim)[["Dropout"]]))
运行结果
# 将存储真是的log2倍数变化进行比较
rowData(sim)$log2FC<-with(rowData(sim),log2(DEFacGroup2/DEFacGroup1))
# 负二项式分量在均值上的真实离差
rowData(sim)$trueDisp<-rowMeans(assays(sim)[["BCV"]])^2
gridlines<-c(1e-2,1e-1,1);cols<-c("blue","red","darkgreen")
with(rowData(sim)[rowData(sim)$GeneMean>1,],
plot(GeneMean,trueDisp,log="xy",xlim=c(1,300),ylim=c(.01,5)))
abline(h=gridlines,col=cols)
text(300,gridlines,labels=gridlines,col=cols,pos=3)
运行结果
今日告一段落，回见。嘻嘻嘻。国庆放假了。祝大家假期快乐！
