Bambu是一种从长读RNA-Seq数据中发现和量化转录本的方法。Bambu使用对齐的reads和基因组参考注释作为输入,并将返回所有已知转录本和新发现转录本的丰度估计值。Bambu使用来自参考注释的信息来纠正拼接连接上的错位,然后减少对齐的读取以读取等效类,并使用此信息在所有感兴趣的样本中识别新的转录本。然后将读值分配给转录本,并使用期望最大化算法获得表达估计。在这里,我们提出了一个示例工作流,用于分析来自新加坡纳米孔表达项目(SG-NEx)的两个人类癌细胞系的纳米孔长读取rna测序数据。
你可以从github安装bambu:
默认运行模式bambu正在使用一组对齐的读取(bam文件)、参考基因组注释(gtf文件、TxDb对象或bambuAnnotation对象)和参考基因组序列(fasta文件或BSgenome)。bambu将返回一个总结的实验对象,其中包含注释的和新的转录本和转录本表达估计的基因组坐标。我们强烈建议使用与基因组比对相同的注释。如果你有一个gtf文件和fasta文件,你可以用以下选项运行bambu:
图书馆(bambu)测试。bam < -执行(“extdata”,“SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000.bam”,包=“bambu”)足总。文件< -执行(“extdata”,“Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa”,包=“bambu”)gtf。文件< -执行(“extdata”,“Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf”,包=“bambu”)bambuAnnotations < -prepareAnnotations(gtf.file)se < -bambu(读取=test.bam,注释=bambuAnnotations,基因组=fa.file)
bambu返回一个summarizeexperiment对象,可以通过如下方式访问:
为了证明Bambu的使用,我们使用了NanoporeRNASeq包中的牛津纳米孔测序生成的长读RNA-Seq数据,该包由来自两个人类细胞系(K562和MCF7)的6个样本组成SG-NEx项目.每个细胞系有3个重复,1次直接RNA测序和2次cDNA测序。读取的数据与22号染色体(Grch38)对齐,并存储为bam文件。在这个工作流中,我们将演示如何应用bambu这些bam文件识别新的转录本和估计转录本表达水平,可视化结果,并识别差异表达的基因和转录本。
bambu获取保存在bam文件中的基因组对齐。这里我们使用的是NanoporeRNASeq包,其中包含与人类第22号染色体前半部分对齐的读取minimap2.
图书馆(bambu)图书馆(NanoporeRNASeq)数据(“SGNexSamples”)SGNexSamples##>数据帧6行6列##> sample_id平台cellLine协议##> <字符> <字符> <字符> <字符> <字符>1 K562_directcDNA_repl..MinION K562直接dna白细胞##> 2 K562_directcDNA_repl..GridION K562直接dna白细胞##> 3 K562_directRNA_repli..GridION K562 directRNA白细胞##> 4 MCF7_directcDNA_repl..MinION MCF7 directcDNA乳房##> 5 MCF7_directcDNA_repl..GridION MCF7 directcDNA Breast##> 6 MCF7_directRNA_repli..GridION MCF7 directRNA乳房# # >文件名# # > <人物>##> 1 NanoporeRNASeq/versi..##> 2 NanoporeRNASeq/versi..##> 3 NanoporeRNASeq/versi..##> 4 NanoporeRNASeq/versi..##> 5 NanoporeRNASeq/versi..##> 6 NanoporeRNASeq/versi..图书馆(ExperimentHub)NanoporeData < -查询(ExperimentHub(),c(“NanoporeRNA”,“GRCh38”,“砰”))bamFile < -Rsamtools::BamFileList(NanoporeData [[“EH3808”]])
bambu另外还需要一个基因组序列,用于纠正读对齐中的剪接连接。理想情况下,我们建议使用与bambu对齐相同的基因组序列文件。
获取fasta文件的路径genomeSequenceData < -查询(ExperimentHub(),c(“NanoporeRNA”,“GRCh38”,“FASTA”))genomeSequence < -genomeSequenceData [[“EH7260”]]
作为一个选项,用户也可以选择使用BSgenome对象:
bambu还需要一个引用文本注释对象,该对象用于纠正读取对齐、训练用于文本发现和量化的模型。注释对象可以从一个gtf文件中创建:
gtf。文件< -执行(“extdata”,“Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf”,包=“bambu”)注释< -prepareAnnotations(gtf.file)
注释对象也可以从TxDb对象中创建:
txdb < -执行(“extdata”,“Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf”,包=“bambu”)注释< -prepareAnnotations(txdb)
注释对象可以被存储并再次用于重新运行bambu。的注释对象NanoporeRNASeq使用by function in中的函数从GTF文件中准备的包bambu.
数据(“HsChr22BambuAnnotation”)HsChr22BambuAnnotation长度为1500的GRangesList对象# # > ENST00000043402美元有2个范围和2个元数据列的GRanges对象:##> seqnames范围链| exon_rank exon_endrink##> | ##> [1] 22 20241415-20243110 - | 2 1##> [2] 22 20268071-20268531 - | 1 2##> -------##> seqinfo: 1个来自未指定基因组的序列;没有seqlengths# # ># # > ENST00000086933美元GRanges对象,包含3个范围和2个元数据列:##> seqnames范围链| exon_rank exon_endrink##> | > [1] 22 19148576-19149095 - | 3 1##> [2] 22 19149663-19149916 - | 2##> [3] 22 19150025-19150283 - | 1 3##> -------##> seqinfo: 1个来自未指定基因组的序列;没有seqlengths# # ># # > ENST00000155674美元有8个范围和2个元数据列的GRanges对象:##> seqnames范围链| exon_rank exon_endrink##> | ##> [1] 22 17137511-17138357 - | 8 1##> [2] 22 17138550-17138738 - | 7 2##> [3] 22 17141059-17141233 - | 6 3##> [4] 22 17143098-17143131 - | 5 4##> [5] 22 17145024-17145117 - | 45##> [6] 22 17148448-17148560 - | 3 6##> [7] 22 17149542-17149745 - | 2##> [8] 22 17165209-17165287 - | 1 8##> -------##> seqinfo: 1个来自未指定基因组的序列;没有seqlengths# # ># # >…##> <1497个元素>
接下来我们应用bambu输入数据(bam文件,注释,基因组序列)。Bambu将执行异构体发现来扩展所提供的注释,然后使用期望-最大化算法量化这些扩展注释中的转录本表达。在这里,我们将使用一个核心,它可以改变为并行处理多个文件。
se##>类:rangedsummarizeexperiment##> dim: 1501##>元数据(2):incompatibleCounts警告##> assays(4):计数CPM fullLengthCounts uniqueCounts##> rownames(1501): BambuTx1 ENST00000043402…ENST00000641933# # > ENST00000641967##> rowData name (11): TXNAME GENEID…txid eqClassById##> colnames(1): 21ce5934bf7470_3844##> colData names(1): name
用户还可以选择应用bambu只做抄本发现
用户还可以选择应用bambu只做定量(不发现异构体)
seUnextended##>类:rangedsummarizeexperiment##> dim: 1500##>元数据(2):incompatibleCounts警告##> assays(4):计数CPM fullLengthCounts uniqueCounts##> rownames(1500): ENST00000043402 ENST00000086933…ENST00000641933# # > ENST00000641967##> rowData名称(4):TXNAME GENEID txid eqClassById##> colnames(1): 21ce5934bf7470_3844##> colData names(1): name
如果一个样本有多个重复,或者计划在不同条件之间进行比较分析,那么将所有样本一起运行而不是单独运行可能会更有益。这可以通过提供您想要一起分析的所有bam文件的路径向量来实现。
bamFiles < -Rsamtools::BamFileList(NanoporeData [[“EH3808”NanoporeData]], [[“EH3809”]],NanoporeData [[“EH3810”NanoporeData]], [[“EH3811”NanoporeData]], [[“EH3812”]],NanoporeData [[“EH3813”]])se。multiSample < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence)
同时运行样本的优势包括:在多个样本中识别的新转录本被分配统一的id,从而可以在不同样本之间进行比较分析。当观察新的转录本时,这对于下游差异表达分析尤其重要。运行多个示例可以是多线程的。当运行多个样本时,默认情况下,bambu将在每个样本上分别训练一个模型,并在每个样本中分别评分。
如果您需要在多个配置中组合样本(例如不同的成对比较),我们建议使用中间的rcFiles来节省处理时间。
在进行转录本发现时,需要在敏感性(检测到的真实小说转录本的数量)和准确性(有多少真实小说转录本)之间取得平衡。为了控制这种平衡,bambu采用新发现率(NDR)作为主要参数。NDR阈值近似于bambu输出的候选小说的比例,相对于它发现的已知转录本的数量,也就是说,NDR为0.1意味着通过阈值的所有转录本中有10%被归类为小说。
如果你使用的基因组期望有大量的新转录本,建议使用更高的NDR,这样就不会错过这些新转录本。相反,如果您使用的是注释良好的基因组,我们建议降低NDR阈值以减少假阳性的出现。默认情况下,与使用人类参考注释(Hg38)训练的默认模型相比,根据预测的注释完整性水平自动为用户选择NDR阈值。
中使用NDR参数手动选择NDR值bambu:
或者,可以在没有阈值的情况下运行转录本发现,生成一个GRangesList注释对象,其中包含用NDR评分评分的所有转录本。注意,这意味着在运行时将quant = FALSEbambu.注释可以通过它们的NDR评分(见下面的例子)、发现和量化步骤之间的读取计数和基因读取比例进行过滤,或者用于其他类型的分析。
newAnnotations < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence,NDR =1,定量=假)注释。< -过滤newAnnotations [(!is.na(mcols(newAnnotations)$NDR)&mcols(newAnnotations)$NDR<0.1)|is.na(mcols(newAnnotations)$NDR)]se。NDR_1<-bambu(读取=bamFiles,注释=annotations.filtered,基因组=genomeSequence,NDR =1,发现=假)
此外,高级用户在运行时可以通过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]
对于下游分析,我们将向描述样本的对象添加感兴趣的条件。这里我们对两种细胞系的比较感兴趣:
要获得使用可分配给每个基因的所有读取(包括与所有现有注释不兼容的读取)的准确基因表达估计值,可以运行以下命令。看一下输出,我们可以看到也有新的基因被鉴定出来
seGene.multiSample<-transcriptToGeneExpression(se.multiSample)seGene.multiSample##>类:rangedsummarizeexperiment##> dim: 579 6# # >元数据(0):##> assays(2):计数CPM##> rownames(579): BambuGene1 BambuGene10…ENSG00000284654# # > ENSG00000284665##> rowData names(2): GENEID newGeneClass##> colnames(6): 21ce5934bf7470_3844 21ce597a36eb_3846…21ce591fdf695a_3854##> colData names(2): name条件
我们可以再次使用该函数通过热图或PCA图来可视化6个样本的基因表达数据。正如预期的那样,来自同一细胞系的样本比跨细胞系的样本显示出更高的相关性。
bambu包括一个函数来编写扩展注释,转录本和基因表达估计,包括任何新发现的基因和转录本到文本文件。
bambu还包括一个函数,只导出扩展注释到GTF文件:
默认情况下,bambu为了训练特定于样本的模型,需要在样本中检测到至少1000个注释副本。在用例中,这是不可能的bambu将使用默认的预训练模型来计算每个读类的转录概率分数(TPS)。如果用户认为他们的注释不足以进行特定样本的训练(例如,如果他们怀疑样本中存在很大比例的真实小说文本),则可以强制执行这种行为。当您希望NDR校准不受使用低质量注释训练的模型的影响时,这是有利的。
se < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence,opt.discovery =列表(fitReadClassModel =假))
默认的预训练模型在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分数
第一步bambu涉及到读类的构造,这占了运行时间的很大一部分。如果我们想在同一个数据集上使用不同的配置多次执行转录本发现和量化,这可能会很耗时。NDR,或样本的组合),特别是当样本很大的时候。为了缓解这个问题,我们可以将读类信息存储为读类文件(rcFiles)bambu运行。类中的输入参数bambu主函数为后续bambu运行。
rcFiles可以通过两种方式生成,一种是在quant和discovery为FALSE时作为bambu()函数的直接输出,另一种是在为rcOutdir参数提供路径时作为写入输出。当rcFiles使用rcOutdir输出时,这是使用BiocFileCache完成的。有关如何访问、使用和识别这些文件的详细信息,请参见在这里.下面是一个简短的示例。
使用rcOutDir生成预处理文件
这将为每个提供读取的示例文件在提供的目录中存储一个预处理过的rcFile。要访问这些文件以供将来使用,我们建议使用BioCFileCache包,该包提供标识示例所需的元数据。
info对象是一个tibble,它将文件名(fpath)与示例(rname)相关联,以帮助您确定需要哪个.rds文件。
信息##> #一个小游戏:6 × 10##> rid rname create…¹accat…²rpath rtype fpath last_…³etag过期# # > <空空的> <空空的> <空空的> <空空的> <空空的> <空空的> <空空的> <双> <空空的> <双># # > 1 BFC1 21 ce5934bf7470_…2023 - 0 2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA# # > 2 BFC2 21 ce597a36eb_38…2023 - 0…2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA# # > 3 BFC3 21 ce59128cfe40_……2023 - 0 2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA# # > 4 BFC4 21 ce594e23cf46_……2023 - 0 2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA# # > 5 BFC5 21 ce596d274518_……2023 - 0 2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA# # > 6 BFC6 21 ce591fdf695a_……2023 - 0 2023 - 0…/ tmp…rela…3 d4d…NA < NA > NA¹create_time,²access_time,##> #³last_modified_time#使用第一个文件运行bambuse < -bambu(读取=信息$rpath [1:3.),注释=HsChr22BambuAnnotation,基因组=genomeSequence)
在以样本为索引的列表形式将quant和discovery设置为false时,也会生成此输出。
se < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence,发现=假,定量=假)# # >|||0%||============|17%||=======================|33%||===================================|50%||===============================================|67%||==========================================================|83%||======================================================================|One hundred.%
由于这是一个中间对象,因此不适合用于一般用例。对于可能出现的任何潜在的高级用例,我们将在下面记录该对象。
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呐呐呐呐呐
有些用例需要知道哪些个人阅读支持特定的文本(小说和注释)。默认情况下,由于它会带来内存开销,该特性是关闭的,但可以使用trackReads参数打开。输出有三列:read_id,相等匹配的索引列表,兼容匹配的索引列表。这些索引匹配在rowRanges(se)中找到的注释
se < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence,trackReads =真正的)# # >|||0%||============|17%||=======================|33%||===================================|50%||===============================================|67%||==========================================================|83%||======================================================================|One hundred.%元数据(se)$readToTranscriptMaps [[1]]##> #一个tibble: 22,069 × 3##> readId equalMatches compatibleMatches##>
##> 1 9db30f53-2e3a-45bc-a01a-586782f00797 ##> 2 dec47cfc-becc-4ff1-9344-0275c2518a41 ##> 3 d544f2e7-a8cc-4a4c-94da-161213ba767b ##> 4 3a49d2c1-92e8-409c- bbc -0901cf6efa6c ##> 5 881bcf9d-9d04-44cb-8d69-a2c03eb9f417 ##> 6 4f41b1d0-69bf-4eeb-8eab-fd552ec1e991 ##> 7 a6d7ed99-683c-4f23-a138-d387059fc65e ##> 8 af2ffb1f-98e3-4e0b-befc-5571fd6af49e ##> 9 9a0ed6a4-6d9b- 4263 -ba97-858f93eca3a7 . ##> 9 9a0ed6a4-6d9b- 4269 -ba97-858f93eca3a7 . ##> ##> 10 1e6750b9-ed24-4626-96f3-c25625ffbfa3 ##> # + 22,059行
在没有或不能进行训练的情况下,默认模型也不适合样本(样本是由不同的技术、物种、条件等生成的),bambu提供训练新模型的选项,如果有注释良好的类似数据可用的话。例如,可以在拟南芥上训练一个模型,将其应用于未注释的植物样本。
#首先使用来自.bam的相关注释数据集训练模型se =bambu(读取=sample1.bam,注释=注释,基因组=fa.file,发现=假,定量=假,opt.discovery =列表(returnModel =真正的))#注意,discovery和quant需要设置为FALSE,或者你可以将它们设置为TRUE,并从rcFile中检索模型,只要returnModel = TRUE。# newDefaultModel = metadata(se[[1]])$model #[[1]]将选择在第一个样本上训练的模型#或者使用rcFile训练模型rcFile < -readRDS(pathToRcFile)newDefaultModel =trainBambu(rcFile)#在另一个样本上使用训练过的模型# sample2。嘭和嘭。File2表示注释不佳的样本的对齐reads和基因组se < -bambu(读取=sample2.bam,注释=零,基因组=fa.file2,opt.discovery =列表(defaultModels =newDefaultModel,fitReadClassModel =假))# trainBambu参数rcFile < -零, min.readCount =2, nrounds =50NDR。阈值=0.1, verbose =真正的
默认情况下bambu没有报告单外显子转录,因为已知它们有很高的假阳性频率,并且没有被使用的剪接连接bambu区分读类。不过bambu在单外显子转录本上训练一个单独的模型,这些预测可以访问并包含在注释中。
se < -bambu(读取=bamFiles,注释=HsChr22BambuAnnotation,基因组=genomeSequence,opt.discovery =列表(min.txScore.singleExon =0))# # >|||0%||============|17%||=======================|33%||===================================|50%||===============================================|67%||==========================================================|83%||======================================================================|One hundred.%
在分析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
差异表达基因的快速总结
总结(德杰尼勒斯)# # >##>从286个总读计数非零##>调整p值< 0.1##> LFC > 0: 19, 6.6%##> LFC < 0(下降):25,8.7%##>离群值[1]:0,0%##> low counts [2]: 0,0%##> (mean count < 0)##>[1]见'cooksCutoff'参数?results##>[2]见'independentFiltering'参数?results
我们还可以将不同异构体的ma图可视化。然而,在不需要任意滤波阈值的情况下,使用原始对数倍变化结果可视化ma图将受到低计数基因log2倍变化相关噪声的影响。如在DESeq2教程.我们对效果大小应用了相同的收缩来改善可视化。
图书馆(apeglm)resLFC < -lfcShrink(dds.deseq系数=“condition_MCF7_vs_K562”,类型=“apeglm”)plotMA(resLFCylim =c(-3.,3.))
我们使用DEXSeq检测替代使用的异构体。
图书馆(DEXSeq)dxd片剂< -DEXSeqDataSet(countData =轮(化验(se.multiSample)$数),sampleData =as.data.frame(colData(se.multiSample)),设计=~样本+外显子+条件:外显子,featureID =rowData(se.multiSample)$TXNAME,groupID =rowData(se.multiSample)$GENEID)dxr < -DEXSeq(dxd)头(dxr)# # >##>轻轨p值:满vs减# # >##>数据帧6行12列##> groupID featureID exonBaseMean分散##> <字符> <字符> <数字> <数字>##> BambuGene1:BambuTx1 BambuGene1 BambuTx1 0.546889 NA##> BambuGene28:BambuTx2 BambuGene28 BambuTx2 1.409260 NA##> BambuGene33:BambuTx3 BambuGene33 BambuTx3 1.819575 NA##> BambuGene10:BambuTx4 BambuGene10 BambuTx4 0.387383 NA##> BambuGene12:BambuTx5 BambuGene12 BambuTx5 6.927309 NA##> ENSG00000283023:BambuTx6 ENSG00000283023 BambuTx6 0.811498 2.56168##> stat pvalue padj K562 MCF7##> <数字> <数字> <数字> <数字> <数字>##> BambuGene1:BambuTx1 NA NA NA NA NA##> BambuGene28:BambuTx2 NA NA NA NA NA##> BambuGene33:BambuTx3 NA NA NA NA NA##> BambuGene10:BambuTx4 NA NA NA NA NA##> BambuGene12:BambuTx5 NA NA NA NA NA##> ENSG00000283023:BambuTx6 0.574596 0.448438 1 0.259637 0.426249##> log2fold_MCF7_K562基因组数据计数数据##> <数字> <矩阵> ##> BambuGene1:BambuTx1 NA 4:3:0:…##> BambuGene28:BambuTx2 NA 4:9:3:…##> BambuGene33:BambuTx3 NA 14:6:3:…##> BambuGene10:BambuTx4 NA 3:2:0:…##> BambuGene12:BambuTx5 NA 32:37:13:…##> ENSG00000283023:BambuTx6 1.4439 1:7:1:…
我们可以看到ma图
问题和问题可以在Bioconductor支持站点提出(一旦bambu可以通过Bioconductor获得):https://support.bioconductor.org.请在你的帖子上贴上bambu。
或者,可以在bambu Github存储库中提出问题:https://github.com/GoekeLab/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