bambu: context-aware transcript quantification from long read
long read sequencing を用いたトランスクリプト検出、定量の手法は発表が続いている。これまでに発表されたものには、TALON、FLAIRなどがある。R packageとして発表されたbambuは比較的新しいものかもしれない。
bambu is a R package for multi-sample transcript discovery and quantification using long read RNA-Seq data. You can use bambu after read alignment to obtain expression estimates for known and novel transcripts and genes. The output from bambu can directly be used for visualization and downstream analysis such as differential gene expression or transcript usage.
https://github.com/GoekeLab/bambu#bambu-context-aware-transcript-quantification-from-long-read-rna-seq-data
information about bambu
- bambu on github
- documentation on bioconductor
- package on bioconductor
- package on bioconda
- bioRxiv article
- Nature method
bioconda channelでbambuをインストール
手早く試してみたく、condaによる導入を行なった。テスト用に環境を作った。
conda create -n py39 python=3.9
conda activate py39
conda search -c bioconda bioconductor-bambu最新バージョンはgithubのものと同じだった。
Loading channels: done
# Name Version Build Channel
bioconductor-bambu 1.0.0 r40h5f743cb_2 bioconda
bioconductor-bambu 1.0.2 r40h399db7b_1 bioconda
bioconductor-bambu 1.0.2 r40h5f743cb_0 bioconda
bioconductor-bambu 1.2.0 r41h399db7b_0 bioconda
bioconductor-bambu 2.0.0 r41h399db7b_0 bioconda
bioconductor-bambu 2.0.3 r41h399db7b_0 bioconda
bioconductor-bambu 2.0.5 r41h619a076_0 bioconda
bioconductor-bambu 2.0.6 r41h619a076_0 bioconda
bioconductor-bambu 2.0.6 r41hc247a5b_1 bioconda
bioconductor-bambu 3.0.1 r42hc247a5b_0 bioconda
bioconductor-bambu 3.0.5 r42hc247a5b_0 bioconda
bioconductor-bambu 3.0.6 r42hc247a5b_0 bioconda
bioconductor-bambu 3.0.8 r42hc247a5b_0 bioconda
bioconductor-bambu 3.0.8 r42hf17093f_1 bioconda
bioconductor-bambu 3.2.4 r43hf17093f_0 bioconda この環境にcondaで導入して確かめてみることにして、以下の通りに実行した。
conda install -c bioconda -c conda-forge bioconductor-bambu依存については、r関係のものをたくさん一緒にインストールしていた。コンソールから、py39の環境にインストールされたRを用いてライブラリをロードしてみた。
~/miniconda3/envs/py39/bin/RR version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-conda-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(bambu)
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: S4Vectors
Loading required package: BSgenome
Loading required package: Biostrings
Loading required package: XVector
Loading required package: rtracklayertutorialのテストスクリプトを実行
bambuのドキュメントには、チュートリアルが示されていたので、ライブラリが動くかどうか試してみた。
最初の実行の時、biocManagerがないというエラーを受け取ったので、この環境にcondaで導入して進めた。
conda install -c conda-forge r-biocmanagerテストスクリプトは以下の通り。
> test.bam <- system.file("extdata", "SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000.bam", package = "bambu")
> fa.file <- system.file("extdata", "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", package = "bambu")
> gtf.file <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu")
> bambuAnnotations <- prepareAnnotations(gtf.file)
> se <- bambu(reads = test.bam, annotations = bambuAnnotations, genome = fa.file)seを表示させてみると、確かに実行されて結果ができていることがわかった。
> se
class: RangedSummarizedExperiment
dim: 105 1
metadata(2): incompatibleCounts warnings
assays(4): counts CPM fullLengthCounts uniqueCounts
rownames(105): ENST00000190165 ENST00000305248 ... ENST00000622412
ENST00000628764
rowData names(11): TXNAME GENEID ... txid eqClassById
colnames(1): SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000
colData names(1): name
> se@metadata
$incompatibleCounts
GENEID SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000
1: ENSG00000064218 0
2: ENSG00000218839 0
3: ENSG00000172785 10
4: ENSG00000107104 41
5: ENSG00000202172 0
6: ENSG00000137090 0
7: ENSG00000107099 0
8: ENSG00000183784 0
9: ENSG00000170122 0
10: ENSG00000283921 0
11: ENSG00000227155 0
12: ENSG00000231808 0
13: ENSG00000229875 0
14: ENSG00000236875 0
15: ENSG00000227914 0
16: ENSG00000227518 0
17: ENSG00000227917 0
18: ENSG00000235880 0
19: ENSG00000277631 0
20: ENSG00000226403 0
21: ENSG00000228115 0
22: ENSG00000181404 8
23: ENSG00000235330 0
GENEID SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000
$warnings
$warnings$SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000
[1] "not all chromosomes present in reference annotations, annotations might be incomplete. Please compare objects on the same reference"
[2] "not all chromosomes from reads present in reference genome sequence, reads without reference chromosome sequence are dropped"
[3] "Bambu was unable to train a model on this sample, and is using a pretrained model" 結果を保存してみた。pathを指定して実行すると、以下の例ではproject下にbambuが作成されていて、結果が格納されていた。
writeBambuOutput(se, path = "./project/bambu/")ls -1 ./project/bambu/
counts_gene.txt
counts_transcript.txt
CPM_transcript.txt
extended_annotations.gtf
fullLengthCounts_transcript.txt
uniqueCounts_transcript.txthead extended_annotations.gtf
9 Bambu transcript 12134 13783 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620";
9 Bambu exon 12134 12190 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "1";
9 Bambu exon 12291 12340 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "2";
9 Bambu exon 12726 12834 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "3";
9 Bambu exon 13088 13157 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "4";
9 Bambu exon 13338 13487 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "5";
9 Bambu exon 13566 13783 . + . gene_id "ENSG00000236875"; transcript_id "ENST00000421620"; exon_number "6";
9 Bambu transcript 14521 29739 . - . gene_id "ENSG00000181404"; transcript_id "ENST00000442898";
9 Bambu exon 14521 14940 . - . gene_id "ENSG00000181404"; transcript_id "ENST00000442898"; exon_number "11";
9 Bambu exon 15081 15149 . - . gene_id "ENSG00000181404"; transcript_id "ENST00000442898"; exon_number "10";Rscriptで実行
環境下にインストールされたRscriptを用いて先程のテストスクリプトを実行した。
cat .test_script_bambu.R
library(bambu)
# Test run
test.bam <- system.file("extdata", "SGNex_A549_directRNA_replicate5_run1_chr9_1_1000000.bam", package = "bambu")
fa.file <- system.file("extdata", "Homo_sapiens.GRCh38.dna_sm.primary_assembly_chr9_1_1000000.fa", package = "bambu")
gtf.file <- system.file("extdata", "Homo_sapiens.GRCh38.91_chr9_1_1000000.gtf", package = "bambu")
bambuAnnotations <- prepareAnnotations(gtf.file)
se <- bambu(reads = test.bam, annotations = bambuAnnotations, genome = fa.file)nohup ~/miniconda3/envs/py39/bin/Rscript ~project/test-bambu.R > nohup.out.log &これによって無事実行されて、ファイルが生成された。