全長(zhǎng)轉(zhuǎn)錄組 | 三代全長(zhǎng)轉(zhuǎn)錄組分析流程（PacBio & ONT ）-- Bambu

今天我們繼續(xù)介紹一款使用三代全長(zhǎng)轉(zhuǎn)錄本數(shù)據(jù)進(jìn)行轉(zhuǎn)錄本注釋和定量的工具 - Bambu。來(lái)自新加坡科技研究局 (A-STAR) 的Jonathan G?ke團(tuán)隊(duì)（圖1）開(kāi)發(fā)的長(zhǎng)度長(zhǎng)RNA-seq轉(zhuǎn)錄組分析工具Bambu，于2023年6月12日發(fā)表在《Nature Methods》雜志上，題目為Context-aware transcript quantification from long-read RNA-seq data with Bambu。該工具基于機(jī)器學(xué)習(xí)來(lái)識(shí)別和表征新轉(zhuǎn)錄本，從而能夠?qū)Σ煌锓N和樣本進(jìn)行適應(yīng)性分析。

圖1. Jonathan G?ke團(tuán)隊(duì)

大多數(shù)轉(zhuǎn)錄本定量方法依賴(lài)于固定的基因參考注釋文件；然而，實(shí)際情況下的轉(zhuǎn)錄組是動(dòng)態(tài)的，需要根據(jù)當(dāng)時(shí)的情境情況來(lái)判斷，靜態(tài)轉(zhuǎn)錄組注釋文件包含了一些基因的非活躍轉(zhuǎn)錄本（isoforms)，而對(duì)于另一些基因存在注釋不完整的情況。Bambu，一種利用長(zhǎng)讀長(zhǎng)RNA-seq測(cè)序數(shù)據(jù)，通過(guò)基于機(jī)器學(xué)習(xí)的轉(zhuǎn)錄本鑒定方法，實(shí)現(xiàn)對(duì)實(shí)際測(cè)序數(shù)據(jù)里轉(zhuǎn)錄本的定量。為了識(shí)別新的轉(zhuǎn)錄本，Bambu 估計(jì)了新發(fā)現(xiàn)率（novel discovery rate），用可解釋的、精度校準(zhǔn)的單一參數(shù)替代了任意的單樣本閾值。Bambu 保留了全長(zhǎng)和獨(dú)特的轉(zhuǎn)錄本序列，使其在存在非活躍轉(zhuǎn)錄本（isoforms) 的情況下能夠進(jìn)行準(zhǔn)確的定量。與現(xiàn)有的轉(zhuǎn)錄本鑒定方法相比，Bambu 在不損失靈敏性的情況下實(shí)現(xiàn)了更高的精度。依據(jù)實(shí)際數(shù)據(jù)情況的動(dòng)態(tài)注釋改善了新的和已知的轉(zhuǎn)錄本的定量。利用長(zhǎng)讀長(zhǎng)RNA測(cè)序數(shù)據(jù)和機(jī)器學(xué)習(xí)，Bambu 促進(jìn)了準(zhǔn)確的轉(zhuǎn)錄本鑒定和定量 （圖2）。

圖2. Bambu文章摘要

一、軟件介紹

Bambu 是一個(gè)利用長(zhǎng)讀長(zhǎng)RNA-Seq數(shù)據(jù)進(jìn)行多樣本轉(zhuǎn)錄本鑒定和定量的 R包。在進(jìn)行序列比對(duì)后，可以使用Bambu 獲取已知和新轉(zhuǎn)錄本以及基因的表達(dá)量。Bambu 的輸出可以直接用于可視化和下游分析，例如差異基因表達(dá)或轉(zhuǎn)錄本使用情況等。

Bambu 分析流程可以分成四步（圖3. a-d）：首先，采用概率模型，利用參考注釋、基因組序列以及從數(shù)據(jù)中提取的特征（詳細(xì)特征可見(jiàn)方法部分）來(lái)校正RNA連接（junction）區(qū)域或位點(diǎn)的比對(duì)。使用相同剪切位點(diǎn)的校正序列被合并成一類(lèi)序列 (Read class, RCs）。第二步，整合來(lái)自所有樣本的合并序列（Read classes），并對(duì)跨樣本間的新發(fā)現(xiàn)率（novel discovery rate，NDR）進(jìn)行計(jì)算。小于特定NDR閾值的Read class序列被歸為新的轉(zhuǎn)錄本，并將新注釋的轉(zhuǎn)錄本加入?yún)⒖嫁D(zhuǎn)錄本注釋中，形成基于實(shí)際樣本的注釋庫(kù)（擴(kuò)展了原本的參考注釋?zhuān)?strong>第三步，利用擴(kuò)展參考基因注釋?zhuān)瑢⒚恳粋€(gè)Read class序列重新歸類(lèi)標(biāo)注。在這個(gè)過(guò)程中，由于比對(duì)誤差造成的不精確的序列匹配可以被修正。第四步，使用擴(kuò)展參考基因注釋?zhuān)瑢?duì)每一樣本中的轉(zhuǎn)錄本進(jìn)行最終的定量，獲得基因/轉(zhuǎn)錄本的表達(dá)矩陣。

圖3. Bambu分析流程

二、軟件安裝

Bambu: https://github.com/GoekeLab/bambu
版本：v3.2.4

1. 使用Bioconductor安裝Bambu

在R中，或Rstudio，或Rstudio-server環(huán)境中運(yùn)行以下命令進(jìn)行安裝，安裝完成如圖4所示。

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")  #如果已經(jīng)安裝BiocManager可忽略此行命令
BiocManager::install("bambu")
圖4. Bambu的R包安裝

建議：因?yàn)槿L(zhǎng)轉(zhuǎn)錄組數(shù)據(jù)量一般較大以及計(jì)算量大，一般臺(tái)式機(jī)或者自己的筆記本電腦部署的本地R軟件，可能導(dǎo)致占用大量?jī)?nèi)存和計(jì)算資源。建議使用服務(wù)器，部署Rstudio-server，這樣就可以在服務(wù)器上進(jìn)行數(shù)據(jù)存儲(chǔ)和計(jì)算了。具體安裝部署步驟可參考官網(wǎng)：https://posit.co/download/rstudio-server/。

2. 測(cè)試Bambu的安裝

用自帶的測(cè)試數(shù)據(jù)確認(rèn)Bambu安裝成功。

#加載bambu軟件包
library(bambu)

#運(yùn)行以下命令進(jìn)行測(cè)序數(shù)據(jù)，參考基因組，參考基因組注釋文件的導(dǎo)入，并進(jìn)行全長(zhǎng)轉(zhuǎn)錄組分析
test.bam <- system.file("extdata", "SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000.bam", package = "bambu") #讀入測(cè)序數(shù)據(jù)

fa.file <- system.file("extdata", "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", package = "bambu") #讀入?yún)⒖蓟蚪M

gtf.file <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu") #讀入?yún)⒖蓟蚪M注釋文件

bambuAnnotations <- prepareAnnotations(gtf.file) #參考基因組注釋預(yù)處理

se <- bambu(reads = test.bam, annotations = bambuAnnotations, genome = fa.file) #全長(zhǎng)轉(zhuǎn)錄組分析

運(yùn)行完最后一行，得到 --- Start isoform quantification --- 和--- Finished running Bambu ---，則說(shuō)明Bambu運(yùn)行和安裝成功。

三、軟件使用

默認(rèn)模式下輸入文件:

  • 經(jīng)過(guò)和參考基因組比對(duì)的序列文件，.bam文件格式。
  • 參考基因組注釋文件，.gtf文件。
  • 參考基因組文件，.fa文件。
1. 準(zhǔn)備參考基因組比對(duì)序列文件

因?yàn)锽ambu說(shuō)明文檔沒(méi)有提供比對(duì)方法，所以這里采用常用的長(zhǎng)度長(zhǎng)（long-read）比對(duì)軟件minimap2。

#根據(jù)三代測(cè)序平臺(tái)和建庫(kù)方法選擇合適的運(yùn)行命令，一步法
$ minimap2 -ax splice:hq -uf ref.fa iso-seq.fq | samtools sort -@ 12 -o align.bam --write-index -        # PacBio Iso-seq/traditional cDNA
$ minimap2 -ax splice ref.fa nanopore-cdna.fa | samtools sort -@ 12 -o align.bam --write-index -      # Nanopore 2D cDNA-seq
$ minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq | samtools sort -@ 12 -o align.bam --write-index -  # Nanopore Direct RNA-seq

#分步進(jìn)行示例
$ minimap2 -ax splice:hq -uf ref.fa iso-seq.fq | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam        # PacBio Iso-seq/traditional cDNA
$ minimap2 -ax splice ref.fa nanopore-cdna.fa | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam      # Nanopore 2D cDNA-seq
$ minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq | samtools view -@ 12 -bS | samtools sort -@ 12 -o align.bam  # Nanopore Direct RNA-seq

#對(duì)bam文件進(jìn)行索引
$ samtools index align.bam

注意:

  • PacBio平臺(tái)和ONT平臺(tái)產(chǎn)生的數(shù)據(jù)，具體參數(shù)有所不同，詳細(xì)請(qǐng)參考minimap2使用文檔。
  • 這里使用的參考基因組和Bambu使用的要統(tǒng)一起來(lái)。
  • 需要提前安裝samtools，留意安裝版本。
  • 因?yàn)閙inimap2輸出的是.sam文件格式，所以使用samtools將.sam轉(zhuǎn)換成.bam，并且使用samtools sort 對(duì) .bam 進(jìn)行排序，然后用samtools index進(jìn)行索引。上面展示了一步轉(zhuǎn).bam和分步轉(zhuǎn).bam的示例。
  • samtools v1.18版本有--write-index - 選項(xiàng); samtools v1.9版本則沒(méi)有。
2. 讀取樣本序列.bam文件
#可同時(shí)讀入多個(gè)樣本文件
samples <- c(S1.bam, S1.bam, S1.bam)
3. 準(zhǔn)備參考基因組文件

如果人和小鼠可以在Gencode上下載，如遇到其它物種可以從Ensembl上下載。

4. 準(zhǔn)備參考基因組注釋
annotations <- prepareAnnotations( *.gft ）

#保存注釋文件為rds文件
saveRDS(annotations, ”/path/to/annotations.rds” )

#下次讀取注釋文件
annotations <- readRDS("/path/to/annotations.rds")
5. 運(yùn)行Bambu
se <- bambu(reads = samples, annotations = annotations, genome = genome.fa)
6. 只鑒定新轉(zhuǎn)錄本，不定量
se.discoveryOnly <- bambu(reads = samples, annotations = annotations, genome = genome.fa, quant = FALSE)
7. 只對(duì)已知的轉(zhuǎn)綠本和基因進(jìn)行定量，不做新的轉(zhuǎn)錄本鑒定
se.quantOnly <- bambu(reads = samples, annotations = annotations, genome = genome.fa, discovery = FALSE)
8. 不依賴(lài)于參考注釋的轉(zhuǎn)錄本組裝
novelAnnotations <- bambu(reads = test.bam, annotations = NULL, genome = fa.file, NDR = 0.5, quant = FALSE)

具體參數(shù)設(shè)置和調(diào)試參考官方github。

四、輸出文件

Bambu 運(yùn)行完后會(huì)返回一個(gè)SummarizedExperiment對(duì)象，可以通過(guò)以下方式訪問(wèn)：

  • assays(se) 返回轉(zhuǎn)錄本表達(dá)豐度矩陣，以count或CPM的形式。
  • rowRanges(se) 返回GRangesList，包含所有已注釋和新發(fā)現(xiàn)的轉(zhuǎn)錄本。
  • rowData(se) 返回每個(gè)轉(zhuǎn)錄本的額外信息。

通過(guò)使用assays()中的變量（如counts或CPM）提取轉(zhuǎn)錄本表達(dá)量：

  • assays(se)$counts - count表達(dá)量。
  • assays(se)$CPM - 測(cè)序深度歸一化后的表達(dá)量。
  • assays(se)$fullLengthCounts - 每個(gè)轉(zhuǎn)錄本的全長(zhǎng)序列count表達(dá)量。
  • assays(se)$uniqueCounts - 每個(gè)轉(zhuǎn)錄本唯一回貼序列的count表達(dá)量。

可以使用writeBambuOutput()輸出完整的結(jié)果文件。這個(gè)命令可以產(chǎn)生三個(gè)文件：1. 擴(kuò)展的.gtf文件。2.轉(zhuǎn)錄本的count表達(dá)矩陣。3.基因的count表達(dá)矩陣。

writeBambuOutput(se, path = "./bambu/")

如果只對(duì)新的轉(zhuǎn)錄本感興趣，可以過(guò)濾掉參考注釋的轉(zhuǎn)錄本。

se.novel = se[mcols(se)$novelTranscript,]
writeBambuOutput(se.novel, path = "./bambu/")

五、可視化

通過(guò)plotBambu可以對(duì)基因/轉(zhuǎn)錄本進(jìn)行可視化（圖5）。

plotBambu(se, type = "annotation", gene_id = "ENSG00000107104")
圖5. 轉(zhuǎn)錄本tx.9的外顯子組成和在樣本中的表達(dá)量

轉(zhuǎn)錄本tx.9以及其對(duì)應(yīng)基因的其它轉(zhuǎn)錄本的結(jié)構(gòu)和表達(dá)量（圖6）。

plotBambu(se, type = "annotation", transcript_id = "tx.9")
圖6. 轉(zhuǎn)錄本tx.9以及其對(duì)應(yīng)基因的其它轉(zhuǎn)錄本的結(jié)構(gòu)和表達(dá)量

通過(guò)plotBambu展示樣本分組PCA聚類(lèi)（圖7）。

圖7. 基于基因/轉(zhuǎn)錄本表達(dá)量的PCA聚類(lèi)

參考文獻(xiàn)

  1. Chen, Ying, et al. "Context-aware transcript quantification from long-read RNA-seq data with Bambu." Nature methods. (2023)
  2. Bambu github： https://github.com/GoekeLab/bambu
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容