今天我們繼續(xù)介紹一款使用三代全長(zhǎng)轉(zhuǎn)錄本數(shù)據(jù)進(jìn)行轉(zhuǎn)錄本注釋和定量的工具 - Bambu亿笤。來(lái)自新加坡科技研究局 (A-STAR) 的Jonathan G?ke團(tuán)隊(duì)（圖1）開發(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)性分析嘿期。
大多數(shù)轉(zhuǎn)錄本定量方法依賴于固定的基因參考注釋文件品擎；然而，實(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ù)，通過基于機(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）篡悟。
一谜叹、軟件介紹
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ì)特征可見方法部分）來(lái)校正RNA連接（junction）區(qū)域或位點(diǎn)的比對(duì)董栽。使用相同剪切位點(diǎn)的校正序列被合并成一類序列 (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ò)展了原本的參考注釋）。第三步歧沪，利用擴(kuò)展參考基因注釋歹撒，將每一個(gè)Read class序列重新歸類標(biāo)注。在這個(gè)過程中诊胞，由于比對(duì)誤差造成的不精確的序列匹配可以被修正暖夭。第四步，使用擴(kuò)展參考基因注釋撵孤，對(duì)每一樣本中的轉(zhuǎn)錄本進(jìn)行最終的定量迈着，獲得基因/轉(zhuǎn)錄本的表達(dá)矩陣。
二邪码、軟件安裝
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")
建議：因?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 ---扫腺，則說明Bambu運(yùn)行和安裝成功岗照。
三、軟件使用
默認(rèn)模式下輸入文件:
- 經(jīng)過和參考基因組比對(duì)的序列文件，
.bam文件格式攒至。 - 參考基因組注釋文件厚者，
.gtf文件。 - 參考基因組文件迫吐，
.fa文件库菲。
1. 準(zhǔn)備參考基因組比對(duì)序列文件
因?yàn)锽ambu說明文檔沒有提供比對(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版本則沒有吨铸。
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. 不依賴于參考注釋的轉(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ì)象，可以通過以下方式訪問：
- assays(se) 返回轉(zhuǎn)錄本表達(dá)豐度矩陣抬纸，以count或CPM的形式咙俩。
- rowRanges(se) 返回GRangesList，包含所有已注釋和新發(fā)現(xiàn)的轉(zhuǎn)錄本湿故。
- rowData(se) 返回每個(gè)轉(zhuǎn)錄本的額外信息阿趁。
通過使用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)錄本感興趣悍募，可以過濾掉參考注釋的轉(zhuǎn)錄本蘑辑。
se.novel = se[mcols(se)$novelTranscript,]
writeBambuOutput(se.novel, path = "./bambu/")
五、可視化
通過plotBambu可以對(duì)基因/轉(zhuǎn)錄本進(jìn)行可視化（圖5）坠宴。
plotBambu(se, type = "annotation", gene_id = "ENSG00000107104")
轉(zhuǎn)錄本tx.9以及其對(duì)應(yīng)基因的其它轉(zhuǎn)錄本的結(jié)構(gòu)和表達(dá)量（圖6）以躯。
plotBambu(se, type = "annotation", transcript_id = "tx.9")
通過plotBambu展示樣本分組PCA聚類（圖7）。
參考文獻(xiàn)
- Chen, Ying, et al. "Context-aware transcript quantification from long-read RNA-seq data with Bambu." Nature methods. (2023)
- Bambu github： https://github.com/GoekeLab/bambu
