分析竹子的Long Read RNA-Seq数据

陈颖,沈安杰,温玉琪,乔纳森Göke

简介

Bambu是一种从长读RNA-Seq数据中发现和量化转录本的方法。Bambu使用对齐的reads和基因组参考注释作为输入,并将返回所有已知转录本和新发现转录本的丰度估计值。Bambu使用来自参考注释的信息来纠正拼接连接上的错位,然后减少对齐的读取以读取等效类,并使用此信息在所有感兴趣的样本中识别新的转录本。然后将读值分配给转录本,并使用期望最大化算法获得表达估计。在这里,我们提出了一个示例工作流,用于分析来自新加坡纳米孔表达项目(SG-NEx)的两个人类癌细胞系的纳米孔长读取rna测序数据。

快速入门:bambu的转录本发现和量化

安装

你可以从github安装bambu:

一般使用

默认运行模式bambu正在使用一组对齐的读取(bam文件)、参考基因组注释(gtf文件、TxDb对象或bambuAnnotation对象)和参考基因组序列(fasta文件或BSgenome)。bambu将返回一个总结的实验对象,其中包含注释的和新的转录本和转录本表达估计的基因组坐标。我们强烈建议使用与基因组比对相同的注释。如果你有一个gtf文件和fasta文件,你可以用以下选项运行bambu:

bambu返回一个summarizeexperiment对象,可以通过如下方式访问:

从纳米孔RNA-Seq数据中识别和量化转录本表达的完整工作流程

为了证明Bambu的使用,我们使用了NanoporeRNASeq包中的牛津纳米孔测序生成的长读RNA-Seq数据,该包由来自两个人类细胞系(K562和MCF7)的6个样本组成SG-NEx项目.每个细胞系有3个重复,1次直接RNA测序和2次cDNA测序。读取的数据与22号染色体(Grch38)对齐,并存储为bam文件。在这个工作流中,我们将演示如何应用bambu这些bam文件识别新的转录本和估计转录本表达水平,可视化结果,并识别差异表达的基因和转录本。

输入数据

基因组序列(fasta文件/ BSGenome对象)

bambu另外还需要一个基因组序列,用于纠正读对齐中的剪接连接。理想情况下,我们建议使用与bambu对齐相同的基因组序列文件。

作为一个选项,用户也可以选择使用BSgenome对象:

基因组注释(bambu注释对象/ gtf文件/ TxDb对象)

bambu还需要一个引用文本注释对象,该对象用于纠正读取对齐、训练用于文本发现和量化的模型。注释对象可以从一个gtf文件中创建:

注释对象也可以从TxDb对象中创建:

注释对象可以被存储并再次用于重新运行bambu。的注释对象NanoporeRNASeq使用by function in中的函数从GTF文件中准备的包bambu

转录本发现和定量

运行多个样本

如果一个样本有多个重复,或者计划在不同条件之间进行比较分析,那么将所有样本一起运行而不是单独运行可能会更有益。这可以通过提供您想要一起分析的所有bam文件的路径向量来实现。

同时运行样本的优势包括:在多个样本中识别的新转录本被分配统一的id,从而可以在不同样本之间进行比较分析。当观察新的转录本时,这对于下游差异表达分析尤其重要。运行多个示例可以是多线程的。当运行多个样本时,默认情况下,bambu将在每个样本上分别训练一个模型,并在每个样本中分别评分。

如果您需要在多个配置中组合样本(例如不同的成对比较),我们建议使用中间的rcFiles来节省处理时间。

调整发现的敏感性(前后分析)

在进行转录本发现时,需要在敏感性(检测到的真实小说转录本的数量)和准确性(有多少真实小说转录本)之间取得平衡。为了控制这种平衡,bambu采用新发现率(NDR)作为主要参数。NDR阈值近似于bambu输出的候选小说的比例,相对于它发现的已知转录本的数量,也就是说,NDR为0.1意味着通过阈值的所有转录本中有10%被归类为小说。

如果你使用的基因组期望有大量的新转录本,建议使用更高的NDR,这样就不会错过这些新转录本。相反,如果您使用的是注释良好的基因组,我们建议降低NDR阈值以减少假阳性的出现。默认情况下,与使用人类参考注释(Hg38)训练的默认模型相比,根据预测的注释完整性水平自动为用户选择NDR阈值。

中使用NDR参数手动选择NDR值bambu

或者,可以在没有阈值的情况下运行转录本发现,生成一个GRangesList注释对象,其中包含用NDR评分评分的所有转录本。注意,这意味着在运行时将quant = FALSEbambu.注释可以通过它们的NDR评分(见下面的例子)、发现和量化步骤之间的读取计数和基因读取比例进行过滤,或者用于其他类型的分析。

此外,高级用户在运行时可以通过opt.discovery访问其他阈值bambu(见参数)。

想象的结果

bambu提供可视化和探索结果的函数。当使用多个样本时,我们可以用热图可视化所有样本的相关性和聚类:

此外,我们还可以用二维PCA图来可视化相关关系。

除了可视化样本之间的相关性,bambu还提供了对单个基因的扩展注释和表达估计的可视化功能。在这里,我们研究了基因ENSG00000099968,并可视化了所有样本中注释的和新异构体的转录坐标以及这些异构体的表达水平。

##> [[1]] ##> TableGrob (3 * 1)“排列”:3 grobs ##> z cell name grob ##> 1 1 1(2,2,1 -1)排列gtable[布局]##> 2 2(3,3,1 -1)排列gtable[布局]##> 3 3(1,1,1 -1)排列文本[GRID.text.244]

从转录物表达中获得基因表达估计

对于下游分析,我们将向描述样本的对象添加感兴趣的条件。这里我们对两种细胞系的比较感兴趣:

要获得使用可分配给每个基因的所有读取(包括与所有现有注释不兼容的读取)的准确基因表达估计值,可以运行以下命令。看一下输出,我们可以看到也有新的基因被鉴定出来

我们可以再次使用该函数通过热图或PCA图来可视化6个样本的基因表达数据。正如预期的那样,来自同一细胞系的样本比跨细胞系的样本显示出更高的相关性。

保存数据(gtf/text)

bambu包括一个函数来编写扩展注释,转录本和基因表达估计,包括任何新发现的基因和转录本到文本文件。

bambu还包括一个函数,只导出扩展注释到GTF文件:

不同用例的高级用法

使用预训练的模型

默认情况下,bambu为了训练特定于样本的模型,需要在样本中检测到至少1000个注释副本。在用例中,这是不可能的bambu将使用默认的预训练模型来计算每个读类的转录概率分数(TPS)。如果用户认为他们的注释不足以进行特定样本的训练(例如,如果他们怀疑样本中存在很大比例的真实小说文本),则可以强制执行这种行为。当您希望NDR校准不受使用低质量注释训练的模型的影响时,这是有利的。

默认的预训练模型在SGNex_HepG2_directRNA_replicate5_run1上进行训练,具有以下特征:

基因组:Homo_sapiens.GRCh38.dna_sm.primary_assembly
注释:Homo_sapiens.GRCh38.91
读计数:7,861,846
技术:纳米孔(ONT)
库准备:directRNA
基地呼叫准确率:79%
平均读长度:1093

我们发现,预先训练的模型在跨物种边界(拟南芥)和不同技术(PacBio)上都能成功工作,与使用样本特定模型相比,性能仅略有下降。预先训练的模型在测序质量差异较大的样本中并不总是有效的,或者如果文库制备导致转录组整体结构的偏差。在这种情况下,我们建议使用来自具有高质量参考注释的不同样本的类似数据来训练一个新模型(参见在另一个物种/数据集上训练一个模型并应用它).

对于本文中的示例,使用的数据相当小,因此默认使用预训练的模型。因此,您不会看到输出中的差异,但是您可以使用更多的数据进行尝试(请参见SG-NEx数据

新转录本发现

在生物体还没有参考注释或不可靠注释的情况下,bambu可以在de-novo模式下运行。在从头开始模式下,bambu不训练模型,而是使用预先训练的模型来分类新的转录本(参见使用预训练的模型.要了解如何为更密切相关的生物/样本训练一个新模型,请参阅在另一个物种/数据集上训练一个模型并应用它.没有注释bambu无法校准NDR输出,也无法建议一个阈值,而是使用TPS作为阈值。因此,您应该提供一个手动NDR阈值(调整发现的敏感性(前后分析)),并注意输出的精度不太可能与应用的阈值线性匹配。TPS阈值为(> 1-NDR)。如果没有提供NDR,则使用默认的NDR阈值<0.1(有效TPS阈值为> 0.9)。就像在调整发现的敏感性(前后分析)NDR为1可以输出所有可能的读类及其TPS分数

存储和使用预处理文件(rcFiles)

第一步bambu涉及到读类的构造,这占了运行时间的很大一部分。如果我们想在同一个数据集上使用不同的配置多次执行转录本发现和量化,这可能会很耗时。NDR,或样本的组合),特别是当样本很大的时候。为了缓解这个问题,我们可以将读类信息存储为读类文件(rcFiles)bambu运行。类中的输入参数bambu主函数为后续bambu运行。

rcFiles可以通过两种方式生成,一种是在quant和discovery为FALSE时作为bambu()函数的直接输出,另一种是在为rcOutdir参数提供路径时作为写入输出。当rcFiles使用rcOutdir输出时,这是使用BiocFileCache完成的。有关如何访问、使用和识别这些文件的详细信息,请参见在这里.下面是一个简短的示例。

使用rcOutDir生成预处理文件

这将为每个提供读取的示例文件在提供的目录中存储一个预处理过的rcFile。要访问这些文件以供将来使用,我们建议使用BioCFileCache包,该包提供标识示例所需的元数据。

info对象是一个tibble,它将文件名(fpath)与示例(rname)相关联,以帮助您确定需要哪个.rds文件。

在以样本为索引的列表形式将quant和discovery设置为false时,也会生成此输出。

由于这是一个中间对象,因此不适合用于一般用例。对于可能出现的任何潜在的高级用例,我们将在下面记录该对象。

rowData(se [[1]])##>数据帧有5601行和23列# # >空空的。rc链。rc startSD endSD readCount.posStrand##> <因子> <因子> <数字> <数字> <整数># # > rc.122+NA NA1# # > rc.222* NA NA 1# # > rc.322+NA NA1# # > rc.422+NA NA1# # > rc.522* NA NA 1##> ... ... ... ... ... ...# # > rcunsplicedNew.9522* NA NA 0# # > rcunsplicedNew.9622* NA NA 1# # > rcunsplicedNew.9722* 9.89949 0.0 0# # > rcunsplicedNew.9822* 1282.19444 1481.3 72# # > rcunsplicedNew.9922* NA NA 0##> intronStarts intronEnds##> <字符> <字符># # > rc.110716989,10718464,10.. 10717993,10719774,10..# # > rc.210717385 10721323# # > rc.310717714 10720835# # > rc.410718330,10718568 10718490,10723770# # > rc.510718577,10723906 10723868,10729611##> ... ... ...# # > rcunsplicedNew.95NA NA# # > rcunsplicedNew.96NA NA# # > rcunsplicedNew.97NA NA# # > rcunsplicedNew.98NA NA# # > rcunsplicedNew.99NA NA##> confidenceType readCount读取##>   # # > rc.1lowConfidenceJunctionReads 1 1# # > rc.2lowConfidenceJunctionReads 1 2# # > rc.3lowConfidenceJunctionReads 1 3# # > rc.4lowConfidenceJunctionReads 1 4# # > rc.5lowConfidenceJunctionReads 1 5##> ... ... ... ...# # > rcunsplicedNew.95unsplicedNew 1 1020# # > rcunsplicedNew.96unsplicedNew 1 1021# # > rcunsplicedNew.97unsplicedNew 2 1022,1023# # > rcunsplicedNew.98unsplicedNew 96 1024,1025,1031,...# # > rcunsplicedNew.99unsplicedNew 1 1176##> GENEID novelGene numExons generedprop##> <字符> <逻辑> <整数> <数字># # > rc.1BambuGene1 TRUE 8 0.142857# # > rc.2BambuGene328 TRUE 2 0.333333# # > rc.3BambuGene2 TRUE 2 1.000000# # > rc.4BambuGene1 TRUE 3 0.142857# # > rc.5BambuGene1 TRUE 3 0.142857##> ... ... ... ... ...# # > rcunsplicedNew.95unstranded.Gene69 TRUE 1 1.0000000# # > rcunsplicedNew.96unstranded.Gene70 TRUE 1 1.0000000# # > rcunsplicedNew.97Ensg00000172967 false 1 0.0124224# # > rcunsplicedNew.98Ensg00000172967 false 1 0.5962733# # > rcunsplicedNew.99BambuGene445 TRUE 1 0.5000000##> geneReadCount等于兼容numstart numAend##>     # # > rc.17 NA NA NA NA# # > rc.23.呐呐呐呐呐# # > rc.31呐呐呐呐呐# # > rc.47 NA NA NA NA# # > rc.57 NA NA NA NA##> ... ... ... ... ... ...# # > rcunsplicedNew.951呐呐呐呐呐# # > rcunsplicedNew.961呐呐呐呐呐# # > rcunsplicedNew.97161 FALSE 0 3 5# # > rcunsplicedNew.98161 FALSE 0 1 6# # > rcunsplicedNew.992呐呐呐呐呐##> numTstart numTend txScore。noFit txScore##>    # # > rc.1呐呐呐呐呐# # > rc.2呐呐呐呐呐# # > rc.3呐呐呐呐呐# # > rc.4呐呐呐呐呐# # > rc.5呐呐呐呐呐##> ... ... ... ... ...# # > rcunsplicedNew.95呐呐呐呐呐# # > rcunsplicedNew.96呐呐呐呐呐# # > rcunsplicedNew.9724 0.0256271 0.00047955# # > rcunsplicedNew.983.10.2085004 0.00814320# # > rcunsplicedNew.99呐呐呐呐呐

下游分析

鉴别差异表达基因

在分析RNA-Seq数据时,最常见的任务之一是在互测条件下分析差异基因表达。这里我们使用DESeq2寻找MCF7和K562细胞系之间的差异表达基因。类似于使用来自的结果大马哈鱼,估计来自bambu首先是四舍五入。

图书馆(DESeq2)dds < -DESeqDataSetFromMatrix化验(seGene.multiSample)数),colData =colData(se.multiSample),设计=条件)dds.deseq < -DESeq(dds)德杰尼勒斯< -DESeq2::结果(dds.deseqindependentFiltering =(德杰尼勒斯订单(德杰尼勒斯padj)))##> log2 fold change (MLE):条件MCF7 vs K562Wald测试p值:条件MCF7 vs K562##>数据帧6行6列##> baseMean log2FoldChange lfcSE stat##> <数字> <数字> <数字> <数字>##> ensg00000185686 582.8059 -7.29920 0.491019 -14.86541##> ensg00000283633 103.4946 -9.26144 1.248518 -7.41794##> ensg00000197077 27.1259 9.15620 1.336144 6.85270##> ensg00000240972 2407.1317 2.39864 0.357799 6.70386##> ensg00000099977 243.9014 1.88449 0.312224 6.03571##> ensg00000169635 46.3432 -3.38693 0.572718 -5.91379# # > pvalue# # > <数字>##> ensg00000185686 0.0000000000000000000000000000000000000000000000000000000000000552742##> ensg00000283633 0.0000000000001189547320822303753823584211316525715478130##> ensg00000197077 0.0000000000072467181792469874010513474977371485693730668##> ensg00000240972 0.0000000000202984969843757815079117861902072194116297688##> ensg00000099977 0.0000000015826591938522844434687589343477460979148219167##> ensg00000169635 0.0000000033432534225140065505430468612664204886009144957# # > padj# # > <数字>##> ensg00000185686 0.000000000000000000000000000000000000000000000000000000158084##> ensg00000283633 0.0000000000170105266877589434019982155840301622264088##> ensg00000197077 0.0000000006908537997548794176378101559695202316113694##> ensg00000240972 0.0000000014513425343828683551974728381029657031664470##> ensg00000099977 0.0000000905281058883506705800033173212049142364321597##> ensg00000169635 0.0000001593617464731676388205768978706400318401392724

差异表达基因的快速总结

我们还可以将不同异构体的ma图可视化。然而,在不需要任意滤波阈值的情况下,使用原始对数倍变化结果可视化ma图将受到低计数基因log2倍变化相关噪声的影响。如在DESeq2教程.我们对效果大小应用了相同的收缩来改善可视化。

鉴别差异转录本用法

我们使用DEXSeq检测替代使用的异构体。

我们可以看到ma图

得到帮助

问题和问题可以在Bioconductor支持站点提出(一旦bambu可以通过Bioconductor获得):https://support.bioconductor.org.请在你的帖子上贴上bambu。

或者,可以在bambu Github存储库中提出问题:https://github.com/GoekeLab/bambu

援引bambu

一份描述竹子的手稿目前正在准备中。如果您使用bambu进行研究,请使用以下doi: 10.5281/zenodo.3900025引用。

会话信息

sessionInfo()##> R版本4.2.2 (2022-10-31)##>平台:x86_64-pc-linux-gnu运行在Ubuntu 20.04.5 LTS下# # >矩阵产品:默认##> BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so##> LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so# # ># # >语言环境:##> [1] LC_CTYPE=en_US。utf - 8 LC_NUMERIC = C##> [3] LC_TIME=en_GB 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]基地# # >##>其他附加包:##> [1] DEXSeq_1.44.0##> [2] RColorBrewer_1.1-3注释dbi_1 .60.0##> [4] BiocParallel_1.32.5##> [5] apeglm_1.20.0##> [6] DESeq2_1.38.3##> [7] ggplot2_3.4.0BSgenome.Hsapiens.NCBI.GRCh38_1.3.1000##> [9] Rsamtools_2.14.0纳米orernaseq_1 .8.0##>[11]实验bbb_2 .6.0##> [12] AnnotationHub_3.6.0##> [13] BiocFileCache_2.6.0##> [14] dbplyr_2.3.0##> [15] bambu_3.0.7##> [16] BSgenome_1.66.2##> [17] rtracklayer_1.58.0##> [18] Biostrings_2.66.0##> [19] XVector_0.38.0总结实验_1.28.0##> [21] Biobase_2.58.0##>[22]基因组范围_1.50.2GenomeInfoDb_1.34.7IRanges_2.32.0##> [25] S4Vectors_0.36.1##> [26] BiocGenerics_0.44.0##> [27] MatrixGenerics_1.10.0##> [28] matrixStats_0.63.0# # >##>通过命名空间加载(并且没有附加):##> [1] utf8_1.2.2 tidyselect_1.2.0##> [3] RSQLite_2.2.20 htmlwidgets_1.6.1##> [5] grid_4.2.2 munsell_0.5.0codetools_0.2-18 statmod_1.5.0##> [9] interp_1.1-3 xgboost_1.7.3.1##>[11]与thr_2.5.0 colorspace_2.0-3##> [13] filelock_1.0.2 OrganismDbi_1.40.0##> [15] highr_0.10 knitr_1.41##> [17] rstudioapi_0.14 labeling_0.4.2##> [19] bbmle_1.0.25 GenomeInfoDbData_1.2.9##> [21] hwriter_1.3.2.1 bit64_4.0.5farver_2.1.1 coda_0.19-4##> [25] vctrs_0.5.1 generics_0.1.3##> [27] xfun_0.36 biovizBase_1.46.0##> [29] R6_2.5.1 doParallel_1.0.17##> [31] clue_0.3-63 locfit_1.5-9.7##> [33] AnnotationFilter_1.22.0 bitops_1.0-7##> [35] cachem_1.0.6 reshape_0.8.9##> [37] DelayedArray_0.24.0 assertthat_0.2.1##>[39]许诺_1.2.0.1 BiocIO_1.8.0##> [41] scales_1.2.1 nnet_7.3-18##> [43] gtable_0.3.1 cairo 1.6-0##> [45] ggbio_1.46.0 ensembldb_2.22.0##> [47] rlang_1.0.6 genefilter_1.80.3##> [49] GlobalOptions_0.1.2 splines_4.2.2##> [51] lazyeval_0.2.2 dichromat_2.0-0.1##> [53] checkmate_2.1.0 BiocManager_1.30.19##> [55] yaml_2.3.6 reshape2_1.4.4##>[57]基因组特征_1.5 .3 backports_1.4.1##> [59] httpuv_1.6.8 Hmisc_4.7-2##> [61] RBGL_1.74.0 tools_4.2.2##> [63] ellipsis_0.3.2 jquerylib_0.1.4##> [65] Rcpp_1.0.9 plyr_1.8.8##> [67] base64enc_0.1-3 progress_1.2.2##> [69] zlibbioc_1.44.0 purrr_1.0.1##> [71] RCurl_1.98-1.9 prettyunits_1.1.1##> [73] rpart_4.1.19 deldir_1.0-6##> [75] GetoptLong_1.0.5 cluster_2.1.4##> [77] magrittr_2.0.3 data.table_1.14.6##> [79] magick_2.7.3 circlize_0.4.15##> [81] mvtnorm_1.1-3 ProtGenerics_1.30.0##> [83] hms_1.1.2 mime_0.12##> [85] evaluate_0.20 xtable_1.8-4##> [87] XML_3.99-0.13 emdbook_1.3.12##> [89] jpeg_0.1-10 gridExtra_2.3 .##> [91] shape_1.4.6 bdsmatrix_1.3-6##> [93] compiler_4.2.2 biomaRt_2.54.0##> [95] tibble_3.1.8 crayon_1.5.2##> [97] htmltools_0.5.4 later_1.3.0[99] Formula_1.2-4 tidyr_1.2.1##> [101] geneplotter_1.76.0 DBI_1.1.3##> [103] formatR_1.14 ComplexHeatmap_2.14.0##> [105] MASS_7.3-58.1 rappdirs_0.3.3##> [107] Matrix_1.5-3 cli_3.6.0##> [109] parallel_4.2.2 pkgconfig_2.0.3##> [111] genome alignments_1.34.0 numderivatives _2016.8-1.1##> [113] foreign_0.8-84 xml2_1.3.3##> [115] foreach_1.5.2 annotate_1.76.0##> [117] bslib_0.4.2 string_1 .5.0##> [119] VariantAnnotation_1.44.0 digest_0.6.31##> [121] graph_1.76.0 rmarkdown_2.19##> [123] htmltable_name .4.1 restfulr_0.0.15##> [125] curl_5.0.0 shiny_1.7.4##> [127] rjson_0.2.21 lifecycle_1.0.3##> [129] jsonlite_1.8.4 fansi_1.0.3##> [131] pillar_1.8.1 lattice_0.20-45##> [133] GGally_2.1.2 KEGGREST_1.38.0##> [135] fastmap_1.1.0 httr_1.4.4##>[137]幸存者3.5-0 interactiveDisplayBase_1.36.0##> [139] glue_1.6.2 png_0.1-8##> [141] iterators_1.0.14 BiocVersion_3.16.0##> [143] bit_4.0.5 stringi_1.7.12##> [145] sass_0.4.4 blob_1.2.3##> [147] latticeExtra_0.6-30 memoise_2.0.1##> [149] dplyr_1.0.10