内容

1简介

这是Bioc2016研讨会“用R和Bioconductor分析单细胞RNA-seq数据”的第一部分。

在本部分中，我们将介绍单细胞RNA-Seq质量控制(QC)和标准化司康饼包中。该包可在Github并将在不久的将来提交给Bioconductor。

单细胞RNA测序(scRNA-Seq)技术为在不同和复杂的哺乳动物组织中进行转录组范围内的分析开辟了道路，促进了对新细胞亚群及其功能作用的公正鉴定。与其他高通量分析一样，在scRNA-Seq数据中观察到的异质性的一小部分是由批处理效应和其他技术人为因素造成的。特别是，这些协议对微量起始mRNA的依赖可能导致广泛的“退出效应”，其中表达的转录本被遗漏。由于这些分析固有的偏差，数据规范化是任何下游分析之前的必要步骤。此外，由于scRNA-Seq研究设计广泛应用于该领域，我们不能指望找到一个适用于所有问题的解决方案。

司康饼支持一个合理的、数据驱动的框架，用于评估各种标准化工作流的有效性，鼓励用户在最终确定数据标准化策略之前探索其数据集固有的权衡。我们为并行运行多个规范化工作流提供了一个接口。我们还提供了对工作流程进行排名和可视化权衡的工具。我们导入了传统批量测序中常用的规范化模块，并提供了集成用户指定规范化模块的支持。

1．1数据集示例

本研讨会的前两部分将分析一小组细胞沿发育轨迹。我们将从从标准转录组比对管道获得的原始数据对象开始。原始数据和重要的参考数据可以直接从车间包中加载。

##加载示例数据Data (ws_input) ls()
# #[1]“批处理”“生物”“休息”“氯”# #[5]“cl3”“cl7”“厘米”“cont_levels”# #[9]“contrMat”“数”“德”“defaultMar”# #[13]“适合”“genesDendro”“genesF”“genesOneAll”# #[17]“genesPairs”“香港”“手动”“mast_res”# #[21]“规范”“norm_logcounts”“主成分分析”“plotCMar”# #[25]“质量控制”“res1”“它”“res3”# #[29]“rs”“sca”“se”“wh_rm”

的计数数据帧包含从96个单细胞库的tophat比对到mm10参考基因组的特征级读取计数(Trapnell, Pachter, and Salzberg 2009)．质量控制包含从中获得的库对齐度量皮卡德．

colnames (qc)
##[1]“nreads”“naligned”“ralign”##[4]“total_dup”“primer”“pct_ribosomal_bases”##[7]“pct_coding_bases”“pct_utr_bases”“pct_intronic_bases”##[10]“pct_intergenic_bases”“pct_mrna_bases”“median_cv_coverage”##[13]“median_5prime_bias”“median_3prime_bias”

德而且香港阳性和阴性对照基因集来源于人群研究。批处理而且生物分别是因素标注批次和时间点。考虑这些因素的联合分布:

##批次及生物条件(诱导后时间)联合分配表(批次、生物)
# # # #生物批time0 time1 time2历史问题time4 time5 # # Y01 10 0 0 0 0 0 # # Y02 0 0 0 6 0 0 # # Y03 0 0 0 8 0 0 # # Y04 8 0 0 0 0 0 # # Y05A 0 10 0 0 0 0 # # Y05B 0 9 0 0 0 0 # # Y06 0 0 5 0 0 0 # # Y07A 0 0 9 0 0 0 # # Y07B 0 0 10 0 0 0 # # Y08 0 6 0 0 0 0 # #日元0 0 0 0 2 0 # # Y11 0 0 0 0 4 0 # # Y12A 0 0 0 0 0 2 # # Y12B 0 0 0 0 0 2 # # Y13 0 0 0 0 0 5

注意，每个时间点由多个批次组成。这种嵌套设计在scRNA-Seq研究中很常见，因为可以同时收集/测序的细胞数量受到限制。

1．2技术变异性和批量效应

我们将首先可视化细胞水平的质量读数，从与基因组对齐的读数比例开始:

#配色方案cc <- c(brewer。pal(9，“Set1”)，brewer。pal(8， "Set2")， brewery .pal(9，"Set3")) #读取比例映射到基因组Barplot (qc$RALIGN[order(batch)]， col=cc[batch][order(batch)]， border=cc[batch][order(batch)]， main="已映射读取百分比，按批次着色")legend("bottomleft"， legend=levels(batch)， fill=cc,cex=0.4)

我们看到批次之间没有重要的差异，而有少数单元的对齐效率特别低。这些细胞可以通过过滤程序去除。我们也可以考虑每个库中的读取数:

Barplot (qc$NREADS[order(batch)]， col=cc[batch][order(batch)]， border=cc[batch][order(batch)]， main=" total number of reads, colored by batch") legend("topright"， legend=levels(batch)， fill=cc, cex=0.4)

我们看到不同批次之间的覆盖率有很大差异，尽管一些技术上的可变性可以通过过滤来解决，但我们仍然需要仔细地规范化数据以使批次具有可比性。

可视化单个指标的分布非常有帮助，但重要的是要注意这些指标通常是相关的。有时考虑质量矩阵的主成分可能更有用，以确定库变化的主要因素:

qpc = prcomp(qc,center = TRUE，刻度。= TRUE) barplot(cumsum(qpc$sdev^2)/sum(qpc$sdev^2)， border="灰色"，xlab="PC"， ylab="方差比例"，main="质量PCA")

即使已经量化了14个质量指标，PCA告诉我们只需要少量的pc就可以描述大部分方差(例如3个来解释74%)。现在我们将在批处理上下文中可视化第一台PC的分布:

barplot(qpc$x[，1][order(batch)]， col=cc[batch][order(batch)]， border=cc[batch][order(batch)]， main="Quality PC1, colored by batch") legend("bottomleft"， legend=levels(batch)， fill=cc, cex=0.8)

这个PC似乎与批处理相关(特别是Y04)，但它也强调了5个与其他库显著不同的离群库。虽然我们不会在这个研讨会上讨论这些pc在过滤方面的应用，但我们将展示如何将它们应用于标准化。

1.3中途退学的特点

接下来我们考虑一个独特的单细胞问题:辍学。香港包含一个基因的列表，被认为是无处不在的，统一表达在目标组织。因为我们假设这些基因在所有细胞中都有表达，所以我们可以将所有0标记为退出事件。下面，我们将检测失败建模为平均表达式的逻辑函数，与该领域采用的标准退出逻辑模型一致。

# Mean log10(x+1)表达式mu_obs = rowMeans(log10(counts[hk，]+1)) # dropouts drop_outs = counts[hk，] == 0 #故障逻辑回归模型ref.glms = list() for (si in 1:dim(drop_outs)[2]){fit = glm(cbind(drop_outs[，si]，1 - drop_outs[，si]) ~ mu_obs,family=binomial(logit)) ref.glms[[si]] = fit$coefficients}

列表ref.glm包含每个拟合的截距和斜率。我们现在可以可视化拟合曲线和相应的曲线下面积(AUCs):

par(mfrow=c(1,2)) #绘制故障曲线并计算AUC图(NULL, main =“假阴性率曲线”，ylim =c(0,1)，xlim =c(0,6)， ylab =“故障概率”，xlab =“Mean log10表达式”)x = (0:60)/10 AUC = NULL for(si in 1:ncol(counts)){y = 1/(exp(-ref. 0)){y = 1/(exp(-ref. 0))glm[[如果]][1]——ref.glms[[如果]][2]* x) + 1) AUC (si) = (y)之和/ 10行(x, 1 / (exp (ref。glms[[si]][1] - ref.glms[[si]][2] * x) + 1)， type = 'l'， lwd = 2, col=cc[batch][si])} # FNR AUC Barplot (AUC[order(batch)]， col=cc[batch][order(batch)]， border=cc[batch][order(batch)]， main="FNR AUC, colored by batch") legend("topleft"， legend=levels(batch)， fill=cc, cex=0.4)

您可能已经注意到，这里看到的5个异常值单元与上面质量度量PCA中突出显示的5个异常值相同。像AUC这样的退出摘要对于评估单细胞库质量非常有用。

1．4的司康饼工作流

到目前为止，我们只描述了数据集的潜在问题。现在我们将采取措施解决这些问题!我们今天将使用的基本qc和标准化管道将允许我们:

SCONE的规范化工作流模板由3个模块组成:

  1. 数据imputation:用dropout模型下的期望值替换零丰度值。在本例中，我们只考虑身份imputation(无imputation)。
  2. 缩放或分位数归一化:i)将每个样本的转录组丰度按单个因素缩放的归一化或ii)更复杂的偏移量，以匹配样本中的分位数。例如:TMM或DESeq比例因子，上四分位数归一化，或全分位数归一化。
  3. 基于回归的方法，用于从数据中去除不需要的相关变化。例子:RUVg(Risso et al. 2014)或上述质量主成分的回归。

2步骤1:样本过滤metric_sample_filter

其中最基本的样本滤波函数司康饼是metric_sample_filter．它接受一个表达式矩阵作为输入。输出取决于提供的参数，但通常由4个逻辑组成，指定每个样本在4个指标上失败(TRUE)或通过(FALSE)基于阈值的过滤器

如果这4个参数中的任何一个缺少必需的参数，函数将简单地返回NA而不是相应的逻辑。

mfilt_report <- metric_sample_filter(expr = counts, nreads = qc$ nreads,ralign = qc$ ralign, suff_nreads = 10^5, suff_ralign = 90, pos_controls = hk, zcut = 3,mix = FALSE, plot = TRUE)

在上面的调用中，我们设置了以下参数:

2.1关于阈值选择

让我们仔细看看重新对齐过滤器背后的计算。在选择阈值90时，metric_sample_filter考虑了4个值:

  1. hard_ralign在美国，默认的“硬”门槛是15——相当宽容……
  2. ３ （zcut)乘以低于均值的标准差ralign价值。
  3. ３ （zcut)乘以中位数以下的MADralign价值。
  4. suff_ralign时，足够阈值设置为90。
hist(qc$RALIGN, breaks = 0:100) #硬阈值abline(v = 15, col = "yellow"， lwd = 2) # 3 (zcut)均值RALIGN值以下的标准差abline(v = mean(qc$RALIGN) - 3*sd(qc$RALIGN)， col = "green"， lwd = 2) # 3 (zcut)均值RALIGN值以下的MADs值abline(v = median(qc$RALIGN) - 3*mad(qc$RALIGN)， col = "red"， lwd = 2) #充分阈值abline(v = 90, col = "grey"，lwd = 2) #最终阈值是1)足够的最低阈值和2)所有其他人的马克斯打= min(90年,马克斯(c(15日意味着(qc RALIGN美元)- 3 * sd (qc RALIGN美元),中等(qc RALIGN美元)- 3 *疯了(qc RALIGN美元))))abline (v =打坳=“蓝色”,lwd = 2, lty = 2)传说(“topleft”,传说= c(“硬”、“sd”、“疯狂”、“足够”,“最终”),lwd = 2,坳= c(“黄”、“绿色”、“红色”、“灰色”、“蓝色”),lty = c (1, 1, 1, 1, 2), cex = 5)

我们在这里看到第三个阈值超过了充分阈值，因此metric_sample_filter确定足够的阈值。注意，如果混合= TRUE一个附加的准则包括:双模态度量分布适合于一个双组分混合模型，并定义了一个阈值，相对于“较好”组分的均值和标准差。作为metric_sample_filter依赖于候选阈值的最大值，我们建议用户将此函数视为自由筛选器。

2．2应用样本过滤器

与metric_sample_filter输出在手，过滤出差的样本是相当简单的:

缺少哪些阈值?(width) is_na_filt = unlist(lapply(is.na(mfilt_report)，any)) #识别不符合任何阈值的样本m_sampfilter = !apply(simplify2array(mfilt_report[!is_na_filt])，1,any) #过滤样本fcounts = counts[，m_sampfilter] fqc = qc[m_sampfilter，] fbatch = batch[m_sampfilter] fbio = bio[m_sampfilter] #简单基因过滤器filterCount <- function(counts, nRead=5, nCell=5){Filter <- apply(counts, 1，function(x) length(x[x>=nRead])>=nCell) return(filter)} genefilter <- filterCount(fcounts) #过滤基因fcounts = fcounts[genefilter，] fhk = hk[hk %in% rownames(fcounts)] fde = de[de %in% rownames(fcounts)]

3.步骤2:运行规范化工作流并对其进行评分司康饼

如前所述，不仅如此司康饼规范化表达式数据，但它也为评估各种规范化工作流的性能提供了一个框架。为了运行司康饼函数时，我们将需要决定要比较哪些工作流。

3．1使用运行= FALSE

的论点司康饼允许很大的灵活性，但有时您需要运行非常特定的模块组合。为此，司康饼可以磨合运行= FALSE模式，只返回要执行的工作流的数据帧。

params <- scone(expr = as.matrix(fcounts)，scaling = c(none = identity,deseq = DESEQ_FN, tmm = TMM_FN, uqp = UQ_FN_POS, FQT_FN)， ruv_negcon = fhk, k_ruv = 3, qc = as.matrix(fqc)， k_qc = 3, bio = fbio,adjust_bio = "yes"， batch = fbatch,adjust_batch = "yes"， run = FALSE) head(params)
# # imputation_method scaling_method # #没有,没有,no_uv, no_bio, no_batch没有没有# #没有,deseq, no_uv, no_bio, no_batch没有deseq # #, tmm, no_uv, no_bio, no_batch没有tmm # #, uqp, no_uv, no_bio, no_batch没有uqp # #, fq, no_uv, no_bio, no_batch fq # #没有,没有没有,ruv_k = 1, no_bio, no_batch没有没有# # uv_factors adjust_biology adjust_batch # #没有,没有,no_uv, no_bio, no_batch no_uv no_bio no_batch # #, deseq, no_uv, no_bio, no_batch no_uv no_bio no_batch # #, tmm, no_uv, no_bio, no_batch no_uvNo_bio no_batch ## none,uqp,no_uv, No_bio,no_batch no_uv No_bio no_batch ## none,fq,no_uv, No_bio,no_batch no_uv No_bio no_batch ## none,none,ruv_k=1, No_bio,no_batch ruv_k=1 No_bio no_batch

在上面的调用中，我们设置了以下参数:

这些参数转换为以下选项集:

应用(params, 2,独特的)
# # $ imputation_method # #[1]”没有“# # # # $ scaling_method # #[1]“没有”“deseq”“tmm”“uqp”“fq uv_factors美元# # # # # #[1]“no_uv”“ruv_k = 1”“ruv_k = 2”“ruv_k = 3”“qc_k = 1”“qc_k = 2”“qc_k = 3”adjust_biology美元# # # # # #[1]“no_bio”“生物”adjust_batch美元# # # # # #[1]“no_batch”“批处理”

虽然在下游分析之前对生物学进行调整通常是有问题的，但在单细胞环境中可能特别危险:scRNA-Seq数据中有更多的样本固有生物学变异性，这些变异性无法通过暴露来捕捉。在我们的嵌套设计场景中，我们只会在考虑批量效应时考虑生物调整。我们可以生产一个更新的参数个数反映这一选择的数据帧:

is_filtered = (params$adjust_biology == "bio") & (params$adjust_batch != "batch") params = params[!is_screened,)

3.2调用司康饼与运行= TRUE

现在我们已经选择了工作流，我们可以运行了司康饼在运行= TRUE模式。此模式提供了一些附加参数，包括一个可选参数参数个数参数来传递运行= FALSE模式。为了理解这些论点，我们必须首先理解用于评估每种标准化的8个指标。前6个指标依赖于通过PCA(默认)将规范化数据降至3维。每个指标都有一个正(越高越好)或负(越低越好)的特征。

res <- scone(expr = as.matrix(fcounts)，scaling = c(none = identity, deseq = DESEQ_FN, tmm = TMM_FN, uqp = UQ_FN_POS, fq = FQT_FN)， ruv_negcon = fhk, k_ruv = 3, qc = as.matrix(fqc)， k_qc = 3, bio = fbio,adjust_bio = "yes"， batch = fbatch,adjust_batch = "yes"， run = TRUE,params = params, eval_poscon = fde, eval_kclust = 2:3, conditional_pam = TRUE)

在上面的调用中，我们设置了以下参数:

运行该命令将花费几分钟时间。输出将包含一个包含四个元素的列表:

名(res)
## [1] "normalized_data" "metrics" "scores" "params"

normalized_data包含规范化表达式数据列表(对数级);中所见的相同约定对每个表达式矩阵进行命名参数个数论点。

指标包含每个标准化的8个原始指标。分数包含指标乘以它们的签名(或“分数”)，以及第九列，其中包含该标准化的平均分数。中的规范化工作流normalized_data，指标,分数按平均分降序排列。

头(res几十美元)
## BIO_SIL BATCH_SIL PAM_SIL ## none,fq,ruv_k=1,no_bio,no_batch 0.08653280 0.1508225 0.3741884 ## none,fq,ruv_k=2,no_bio,no_batch 0.09865276 0.1415136 0.3629964 ## none,deseq,ruv_k=2,no_bio,no_batch 0.09721160 0.1722529 0.3983736 ## none,deseq,ruv_k=1,no_bio,no_batch 0.08828919 0.1640219 0.4215769 ## none,fq,qc_k=1,no_bio,no_batch 0.08318459 0.1581687 0.4326760 ## EXP_QC_COR EXP_UV_COR ## # none,fq,ruv_k=1,no_bio,no_batch ## # none-0.2132825 -0.09864541 0.5405376 ## none,fq,ruv_k=2,no_bio,no_batch -0.2012192 -0.08994268 0.5313734 ## none,fq,ruv_k=3,no_bio,no_batch -0.2300302 -0.06228916 0.5779535 ## none,deseq,ruv_k=2,no_bio,no_batch -0.2984754 -0.16153256 0.5042606 ## none,deseq,ruv_k=1,no_bio,no_batch -0.3437892 -0.16409015 0.4916737 ## none,fq,qc_k=1,no_bio,no_batch -0.2336496 -0.28368051 0.5126642 ## RLE_MED RLE_IQR mean_score ## none,fq,ruv_k=1,no_bio,no_batch -0.002169450 -0.9977360 -0.01996901 ## none,fq,ruv_k=2,no_bio,no_batch -0.002856720 -1.0625194 -0.02704359 ## none,fq,ruv_k=3,no_bio,no_batch -0.002686065 -1.1112346 -0.02814046 ## none,deseq,ruv_k=2,no_bio,no_batch -0.007600576 -0.9843497 -0.03498245 ## none,deseq,ruv_k=1,no_bio,no_batch -0.007537978 -0.9375725 -0.03592851 ## none,fq,qc_k=1,no_bio,no_batch -0.003546027 -0.9857528 -0.03999193

4步骤3:为下游分析选择规范化

根据我们的排序标准，看起来没有,fq ruv_k = 1 no_bio no_batch与其他规范化工作流相比，性能更好。将此方法与其他方法可视化的一个有用方法是biplot_colored函数

pc_obj = prcomp(res$scores[，-ncol(res$scores)]，center = TRUE,scale = FALSE) bp_obj = biplot_colored(pc_obj,y = -res$scores[，ncol(res$scores)]，expand = .6)

我们根据相应方法的平均分给上面的每个点上色。我们可以看到工作流主要通过三个分数聚类:BIO_SIL、EXP_WV_COR和EXP_UV_COR。由于生物学的聚类作用很大，我们可以突出生物学调整对应的几点:

Bp_obj = biplot_colored(pc_obj,y = -res$scores[，ncol(res$scores)]，expand = .6) points(Bp_obj [grepl("，bio,batch"，rownames(Bp_obj))，]， PCH = 1, col = "red"， cex = 1) points(Bp_obj [grepl("，bio,batch"，rownames(Bp_obj))，]， PCH = 1, col = "red"， cex = 1.5)

我们看到，左下角的中等性能集群仅由这些方法组成。这个例子强调了调整生物学的危险，即使是在嵌套设计场景中。虽然生物学调整的数据可能更强烈地反映先前的差异，但其他指标将表明数据表现不佳。我们也可以考虑只考虑batch的情况:

Bp_obj = biplot_colored(pc_obj,y = -res$scores[，ncol(res$scores)]，expand = .6) points(Bp_obj [grepl("，no_bio,batch"，rownames(Bp_obj))，]， PCH = 1, col = "red"， cex = 1) points(Bp_obj [grepl("，no_bio,batch"，rownames(Bp_obj))，]， PCH = 1, col = "red"， cex = 1.5)

在这里，我们看到右边性能较差的集群完全由这些方法组成。这源于一个简单的事实，即批量在生物学中是嵌套的，对前者的调整可以消除关键的生物学变异。最后，我们将可视化性能最好的方法以及它与非标准化数据的关系:

bp_obj = biplot_colored (pc_obj y = - r几十美元[,ncol (res几十美元)],扩大= 6)点(t (bp_obj [1]), pch = 1 =“红色”,上校cex = 1)点(t (bp_obj [1]), pch = 1 =“红色”,上校cex = 1.5)点(t (bp_obj [rownames (bp_obj) = = rownames (params) [1]]), pch = 1 =“蓝色”,上校cex = 1)点(t (bp_obj [rownames (bp_obj) = = rownames (params) [1]]), pch = 1 =“蓝色”,上校cex = 1.5)箭头(bp_obj [rownames (bp_obj) = = rownames (params) [1]] [1], bp_obj [rownames (bp_obj) = = rownames (params) [1]] [2],Bp_obj [1，][1]， Bp_obj [1，][2]， lty = 2, LWD = 2)

箭头指向从“no-op”归一化到SCONE中最高级归一化的一条线。我们看到SCONE选择了一种介于两个极端之间的方法，减少了不需要的变化信号(定义为fhk)，同时保存生物信号(定义为fbio而且fde）.

5会话信息

sessionInfo ()
## R版本3.3.0(2016-05-03)##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 14.04.4 LTS ## ##语言环境:## [1]LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_US。UTF-8 LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:##[1]并行stats4统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:[7] SummarizedExperiment_1.3.5 Biobase_2.33.0 ## [9] GenomicRanges_1.25.8 GenomeInfoDb_1.9.1 ## [11] IRanges_2.7.11 S4Vectors_0.11.7 ## [13] BiocGenerics_0.19.1 bioc2016singlecell_0.0.1 ## [15] devtools_1.11.1 BiocStyle_2.1.10 ## ##通过命名空间加载(并且没有附加):# # # # [1] colorspace_1.2-6 hwriter_1.3.2 [3] class_7.3-14 modeltools_0.2-21 # # [5] mclust_5.2 XVector_0.13.2 # # [7] flexmix_2.3-13 AnnotationDbi_1.35.3 # # [9] mvtnorm_1.0-5 xml2_0.1.2 # # [11] codetools_0.2-14 splines_3.3.0 # # [13] R.methodsS3_1.7.1 doParallel_1.0.10 # # [15] bold_0.3.5 DESeq_1.25.0 # # [17] robustbase_0.92-6 geneplotter_1.51.0 # # [19] knitr_1.13 ade4_1.7-4 # # [21] jsonlite_0.9.22 locfdr_1.1-8 # # [23] Rsamtools_1.25.0 phylobase_0.8.2 # # [25] gridBase_0.4-7 annotate_1.51.0 # # [27]kernlab_0.9-24 R.oo_1.20.0 ## [29] rentrez_1.0.2 httr_1.2.0 Matrix_1.2-6 ## [37] htmltools_0.3.5 tools_3.3.0 ## [39] gtable_0.2.0 taxize_0.7.8 ## [41] reshape2_1.4.1 dplyr_0.4.3 ## [43] ShortRead_1.31.0 Rcpp_0.12.5 ## [45] NMF_0.20.6 trimcluster_0.1- 1 ## [47] Biostrings_2.41.4 gdata_2.17.0 ## [51] rtracklayer_1.33.7 iterators_1.0.8 ## [53] fpc_2.1-10 string_1 .0.0 ## [55]rngtools_1.2.4 gtools_3.5.0 ## [57] xml_1 .2.0 ## [59] edgeR_3.15.0 DEoptimR_1.0-4 ## [61] zlibbioc_1.19.0 MASS_7.3-45 ## [63] yaml_2.1.13 RUVSeq_1.7.2 ## [67] memoise_1.0.0 ggplot2_2.1.0 ## [69] pkgmaker_0.22 segmented_0.5-1.4 ## [71] biomaRt_2.29.2 latticeExtra_0.6-28 ## [73] stringi_1.1.1 RSQLite_1.0.0 ## [75] genefilter_1.55.2 foreach_1.4.3 ##[77]基因组features_1 .25.14 caTools_1.17.1 ## [79] boot_1.3-18 BiocParallel_1.7.4 ## [81]chron_2.3-47 prabclus_2.2-6 ## [83] matrixStats_0.50.2 bitops_1.0-6 ## [87] lattice_0.20-33 GenomicAlignments_1.9.4 ## [89] rredlist_0.1.0 plyr_1.8.4 ## [91] magrittr_1.5 R6_2.1.2 ## [93] gplots_3.0.1 DBI_0.4-1 ## [95] whisker_0.3-2 withr_1.0.2 ## [97] mixtools_1.0.4 RCurl_1.95-4.8 ## [101] nnet_7.3-12 tibble_1.0 ## [103] KernSmooth_2.23-15 howmany_0.3-1 ## [107] rmarkdown_0.9.6 RNeXML_2.0.6 ## [109] grid_3.3.0 data.table_1.9.6 ## [111] digest_0.6.9 diptest_0.75-7 ## [113] xtable_1.8-2 tidyr_0.5.1 ## [115] R.utils_2.3.0 munsell_0.4.3 ## [117] registry_0.3

瑞索，戴维，约翰·恩盖，特伦斯·p·斯毕德和桑德琳·杜多特，2014。“使用控制基因或样本的因子分析的Rna-Seq数据的标准化。”自然生物技术32: 896 - 902。

特拉普内尔，科尔，利奥·帕切特和史蒂文·l·萨尔茨伯格，2009。TopHat:用Rna-Seq发现剪接生物信息学25(9): 1105-11。