使用的例子curatedAdipoRNA
动机
所有的样品在这个数据集来自3 t3-l1细胞系。的MDI感应媒体,被用来诱导脂肪细胞的分化。数据集的两个重要变量时间和阶段对应时间点和分化阶段,当样本被抓获。理想情况下,这个数据集应该被视为一次课程。然而,对于本例的目的,我们只使用两个时间点0样本和24小时和把他们作为独立的团体。这个例子的目的是展示如何应用于一个典型的微分表达式分析数据集。重点是解释的数据和元数据adipo_counts在每一个主要的分析。我们开始通过过滤低质量样品和低计数的基因。然后我们应用DESeq2方法用默认值。
过滤低质量的样品
首先,我们的子集adipo_counts对象的所有样品时间点0或24。样品的总数量是30;22在0小时和8个样品在24小时。特性/基因集的总数是23916。
#构造子集数为0和24小时se < -adipo_counts [, adipo_counts美元时间%,%c(0,24)]#显示功能的数量,时间和样品组昏暗的(se)30 # > [1]23916表(se美元时间)# ># > 0 24# > 22日8自质量量度报告/文件运行,我们需要SSR * id为每个样本。注意,一些样品会有多个文件。在这种情况下,因为有些样品paired-end,所以每个人都会有两个文件SRR \ * _1和SRR \ * _2。
#过滤低质量的样品#赤图书馆布局表(se美元library_layout)# ># >配对单# > 12日18#在qc检查文件的数量qc < -se美元质量控制表(长度(qc))# ># > 1 2# > 18 12#压扁qc列表qc < -unlist(质量控制,递归=假)长度(质量控制)42 # > [1]的质量控制对象的colData包含的输出fastqc在一个qc_read类。更多信息可以通过调用访问这个对象? fastqcr:: qc_read。在这里,我们只使用per_base_sequence_quality过滤掉低质量的样品。这绝不是足够的质量控制,但它应该驱动点回家。
#提取per_base_sequence_qualityper_base < -拉普兰人(质量控制,函数(x) {df < -x [[“per_base_sequence_quality”]]df% > %选择(意思是基础)% > %变换(基础=strsplit(as.character(基地)“- - -”))% > %unnest运算(基地)% > %变异(基础=as.numeric(基地)})% > %bind_rows(.id =“运行”)整理数据后,我们得到一个data.frame三列;运行,的意思是和基地运行ID,平均质量分数和碱值在每个阅读。fastqc提供详细的文档和其他质量控制模块。注意阅读运行之间的长度变化显著,平均评分的平均值是合适的。
#快速查看质量得分总结(per_base)# >运行基础的意思# >长度:3408分钟:1.00分钟。:10.74# >类:字符瞿。1:21.00瞿。1:34.45# >模式:字符值:42.00中值:36.44意思是:50.24 # >的意思是:35.28瞿瞿# > 3日:70.00第三:37.90# > Max。:368.00马克斯。:40.17识别低质量的样品,我们分类的运行长度和run_average是读的长度和平均每基平均分数。下图应该更容易明白为什么这些截止在这种情况下使用。
#发现低质量per_base < -per_base% > %group_by(运行)% > %变异(长度=马克斯(基地)>150年,run_average =的意思是(意味着)>34)#情节平均每基地质量per_base% > %ggplot(aes(x =基地,y =的意思是,组=运行时,颜色=run_average))+geom_line()+facet_wrap(~长度,尺度=“free_x”)
“坏”的运行IDs样本用于从数据集中删除它们。
#运行ids的低质量的样品bad_samples < -data.frame(样品=独特的(per_base美元运行[per_base美元run_average= =假)))bad_samples < -单独的(bad_samples坳=样品,在=c(“id”,“运行”),9月=”\ \”。)#子集数对象se2 < -se (,!se美元id%,%bad_samples美元id)表(se2美元时间)# ># > 0 2419 # > 6过滤低计数基因
识别低计数功能/基因(可能不是表达),我们只保留特性至少10读入2个或更多的样品。然后我们对象排除这些基因子集。
#过滤低计数的基因low_counts < -应用(分析(se2),1,函数(x)长度(x (x>10])> =2)表(low_counts)# > low_counts# >假真# > 10254 13662#构造子集数对象se3 < -se2 (low_counts,)应用微分表达式使用DESeq2
DESeq2是一个很好的记录和广泛应用为微分表达式分析R包。在这里我们使用默认值DESeq找到样品之间的基因是谦恭地表达在24小时和0小时。
#差异表达分析se3美元时间< -因素(se3美元时间)dds < -DESeqDataSet(se3~时间)# >重命名的第一个元素化验“计数”# >将计数转换为整数模式dds < -DESeq(dds)# >估计大小的因素# >估计分散# > gene-wise分散估计# >平均分散关系# >最终分散估计# >拟合模型和测试# >——取代异常值和改装为53基因# >——DESeq论点minReplicatesForReplace = 7# >——保留原始计数在计数(dds)# >估计分散# >拟合模型和测试res < -结果(dds)表(res美元padj<1。)# ># >假真# > 6194 7462下一个!
在这个例子中,我们没有试图纠正研究因素之间可能混淆的结果。这怎么可能,我们使用主成分分析情节有一些这些因素在以下图表。第一个使用时间因素是兴趣的因素。我们可以看到,DESeq转换做得很好分离的样品预计组。然而,这似乎也时间玩的不是唯一的因素。例如,我们在第二和第三个图显示其他两个因素library_layout和instrument_model这或许可以解释一些样品之间的差异。这是预期的,因为数据收集来自不同研究使用不同协议和不同的测序仪器。因此,有必要考虑这些差异获得可靠的结果。有多种方法等RUV)删除不必要的变化(政府和代理变量分析(上海广电)。
#解释variabceplotPCA(rlog(dds),intgroup =“时间”)
plotPCA(rlog(dds),intgroup =“library_layout”)
plotPCA(rlog(dds),intgroup =“instrument_model”)
援引的研究数据的子集
说到学习,如前所述研究对象包含完整信息的原始研究的引用发表的数据。请使用这个数据集时引用。
#键的研究数据的子集独特的(se3美元bibtexkey)# > [1]“Duteil2014”“zhao_fto-dependent_2014”# > [3]“siersbaek_molecular_2014”“Lim2015”# > [5]“brunmeir_comparative_2016”“Brier2017”# > [7]“park_distinct_2017”“chen_diabetes_2018”# > [9]“siersbaek_dynamic_2017”“ryu_metabolic_2018”