2022/02/15 追記、コマンド修正
遺伝子やアイソフォームの発現変化を正確に定量することは、細胞の機能、分化、疾病の理解に不可欠である。ロングリードDirect RNA Sequencing (DRS) を用いた完全長ネイティブRNAのシーケンスは、RNAの断片化、cDNA合成、PCRを必要とするショートリードおよびロングリードのシーケンス法の多くの制限を克服できる可能性を持っている。しかし、DRSのために特別に設計されたツールは不足しており、複雑な生物における発現差を同定する能力は十分に特徴付けられていない。本著者らは、DRSにおける転写産物のアイソフォームを高速かつ正確に定量するためにNanoCountを開発し、類似の手法よりも優れていることを実証する。合成コントロールとヒトSH-SY5Y細胞のニューロン様細胞への分化を用い、DRSが正確にRNA発現を定量し、遺伝子とアイソフォームの発現差を同定することを示す。未分化と分化したSH-SY5Y細胞間で231の遺伝子、333のアイソフォーム、さらに27のアイソフォームスイッチの発現差が検出され、遺伝子とアイソフォームレベルで分化状態ごとにサンプルがクラスタ化された。ニューロン様細胞で発現が上昇した遺伝子は、ニューロン新生と関連していた。また、DRSで発見された数千の新規アイソフォームをNanoCountで定量したところ、同様にその発現量の差異を同定することができた。この結果は、NanoCountによるDRSアイソフォームの定量の向上を実証し、DRSが生物学的に関連した遺伝子とアイソフォームの発現差を同定する能力を確立したことを示すものである。
NanoCountはminimap2と互換性があるように設計されている。NanoCountは、リードをトランスクリプトームにアライメントして得たソートされたBAMファイルを必要とする。NanoCountの最初のステップでは、最適化された(カスタマイズ可能な)基準に基づいてリードをフィルタリングし、リードごとのベストアライメントが定義される。フィルタリング後、期待値最大化(EM)アルゴリズムを使用して、転写産物の存在量が推定される。まず、初期リード/トランスクリプトスコアハッシュテーブルを生成する(互換性インデックス)。各アラインメントに対して、リードごとのアラインメント数の割合に対応する初期スコアを定義する。つまり、マルチマッピングリードからのアライメントが重み付けされる。次に、EMループが開始され、転写産物の存在量を計算する。各転写物の存在量は、そのリファレンスにマッピングされたアラインメントの現在の互換性インデックスから抽出されたスコアの合計を、すべてのアラインメントの合計スコアで割ったものとして定義される。各アラインメントについて、現在の互換性スコアを、対応するリードがアラインメントされる転写物の存在量スコアの合計で割る。これは、真陽性である可能性がより高い、豊富な転写物にマッピングされたアラインメントのスコアを増加させる。これらの2つのステップは、収束に達するまで繰り返される。コンバージェンスは、連続する2つのEMラウンド間の累積転写産物量の差が小さいと定義される(カスタマイズ可能)。最後に、NanoCountは最終的な転写産物アバンダンステーブルに基づき、推定カウントと正規化TPM値を計算する。
最初のフィルタリングのステップでは、まず基本的なDRSフィルタリングを行い、アライメント率が低いリード(<0.5、カスタマイズ可能)、短いアライメント(<50 nt、カスタマイズ可能)、ネガティブストランドアライメントが破棄される。DRSのメカニズム上、シーケンスはRNAの3′末端(polyA tail)から始まるため、アライメントの3′末端は転写物の3′末端に近いことが望ましい。NanoCountは、アラインメント末端と転写産物の3′末端との最大塩基間距離を3′フィルタリングしきい値(MAX_DIST_3_PRIME、デフォルト:50 nt)で指定する。また、セカンダリーアライメントが有効とみなされるために満たすべきアライメントスコア(AS)の閾値(SEC_SCORING_THRESHOLD、デフォルト:0.95)を含んでいる。デフォルトでは、二次アライメントスコアがそのリードの最高アライメントスコアの95%以上であることを要求する。完全なドキュメントは、https://a-slide.github.io/NanoCount/ にある。
Documentation
https://a-slide.github.io/NanoCount/
https://a-slide.github.io/NanoCount/NanoCount_API_usage/
Excited to share our paper using NanoCount for accurate @nanopore direct RNA sequencing (DRS) quantification https://t.co/htMweJxYd5
— Josie Gleeson (@gleeson_josie) November 29, 2021
インストール
ubuntu18.04にて、condaで環境を作って導入した(python3.9)。
#conda(link)#aleg チャネルだけでは依存の一部が入らない。conda-forgeかbiocondaを指定
mamba create -n nanocount -y
conda activate nanocount
mamba install -c aleg -c bioconda nanocount -y
#samtoolsとminimap2
mamba install -c bioconda samtools minimap2 -y
#pip (link)
pip install NanoCount
> NanoCount -h
usage: NanoCount [-h] [--version] -i ALIGNMENT_FILE [-o COUNT_FILE] [-b FILTER_BAM_OUT] [-l MIN_ALIGNMENT_LENGTH] [-f MIN_QUERY_FRACTION_ALIGNED] [-s SEC_SCORING_VALUE] [-t SEC_SCORING_THRESHOLD]
[-c CONVERGENCE_TARGET] [-e MAX_EM_ROUNDS] [-x] [-p PRIMARY_SCORE] [-a] [-d MAX_DIST_3_PRIME] [-u MAX_DIST_5_PRIME] [-v] [-q]
NanoCount estimates transcripts abundance from Oxford Nanopore *direct-RNA sequencing* datasets, using an expectation-maximization approach like RSEM, Kallisto, salmon, etc to handle the uncertainty of multi-
mapping reads
optional arguments:
-h, --help show this help message and exit
--version show program's version number and exit
Input/Output options:
-i ALIGNMENT_FILE, --alignment_file ALIGNMENT_FILE
Sorted and indexed BAM or SAM file containing aligned ONT dRNA-Seq reads including secondary alignments (required) [str]
-o COUNT_FILE, --count_file COUNT_FILE
Output file path where to write estimated counts (TSV format) (default: None) [str]
-b FILTER_BAM_OUT, --filter_bam_out FILTER_BAM_OUT
Optional output file path where to write filtered reads selected by NanoCount to perform quantification estimation (BAM format) (default: None) [str]
Misc options:
-l MIN_ALIGNMENT_LENGTH, --min_alignment_length MIN_ALIGNMENT_LENGTH
Minimal length of the alignment to be considered valid (default: 50) [int]
-f MIN_QUERY_FRACTION_ALIGNED, --min_query_fraction_aligned MIN_QUERY_FRACTION_ALIGNED
Minimal fraction of the primary alignment query aligned to consider the read valid (default: 0.5) [float]
-s SEC_SCORING_VALUE, --sec_scoring_value SEC_SCORING_VALUE
Value to use for score thresholding of secondary alignments either "alignment_score" or "alignment_length" (default: alignment_score) [str]
-t SEC_SCORING_THRESHOLD, --sec_scoring_threshold SEC_SCORING_THRESHOLD
Fraction of the alignment score or the alignment length of secondary alignments compared to the primary alignment to be considered valid alignments (default: 0.95) [float]
-c CONVERGENCE_TARGET, --convergence_target CONVERGENCE_TARGET
Convergence target value of the cummulative difference between abundance values of successive EM round to trigger the end of the EM loop. (default: 0.005) [float]
-e MAX_EM_ROUNDS, --max_em_rounds MAX_EM_ROUNDS
Maximum number of EM rounds before triggering stop (default: 100) [int]
-x, --extra_tx_info Add transcripts length and zero coverage transcripts to the output file (required valid bam/sam header) (default: False) [boolean]
-p PRIMARY_SCORE, --primary_score PRIMARY_SCORE
Method to pick the best alignment for each read. By default ("alignment_score") uses the best alignment score (AS optional field), but it can be changed to use either the primary
alignment defined by the aligner ("primary") or the longest alignment ("alignment_length"). choices = [primary, alignment_score, alignment_length] (default: alignment_score) [str]
-a, --keep_suplementary
Retain any supplementary alignments and considered them like secondary alignments. Discarded by default. (default: False) [boolean]
-d MAX_DIST_3_PRIME, --max_dist_3_prime MAX_DIST_3_PRIME
Maximum distance of alignment end to 3 prime of transcript. In ONT dRNA-Seq reads are assumed to start from the polyA tail (-1 to deactivate) (default: 50) [int]
-u MAX_DIST_5_PRIME, --max_dist_5_prime MAX_DIST_5_PRIME
Maximum distance of alignment start to 5 prime of transcript. In conjunction with max_dist_3_prime it can be used to select near full transcript reads only (-1 to deactivate). (default:
-1) [int]
Verbosity options:
-v, --verbose Increase verbosity for QC and debugging (default: False) [boolean]
-q, --quiet Reduce verbosity (default: False) [boolean]
実行方法
NanoCountは、Oxford NanoporeのダイレクトRNAシーケンスデータセットにのみ使用されることを意図している。
1、minimap2を使って、リードをトランスクリプトームリファレンスにアライメントして、ソートされたbamを出力する。少なくとも10個のセカンダリーマッピングを保持するために、"-N 10"オプションを使用することが推奨されている(繰り返しの多いトランスクリプトームでは、この値を大きくすることも可能)。"-p 0"もつける。間違えてゲノム配列を使わないように注意する。
minimap2 -t 12 -ax map-ont -N 10 -p 0 ref.fa.gz DRS_reads.fq.gz | samtools view -@ 4 -bh > aligned.bam
- -p FLOAT min secondary-to-primary score ratio [0.8]
- -N INT retain at most INT secondary alignments [5]
トランスクリプトームリファレンスを使用するため、アライメントアルゴリズムはスプライスを認識する必要はない。
2、NanoCountのラン。転写産物の存在量の推定。
NanoCount -i aligned.bam -o transcript_counts.tsv -b output.bam
- -i Sorted and indexed BAM or SAM file containing aligned ONT dRNA-Seq reads including secondary alignments (required) [str]
- -o Output file path where to write estimated counts (TSV format) (default: None) [str]
- -b Optional output file path where to write filtered reads selected by NanoCount to perform quantification estimation (BAM format) (default: None) [str]
- -x Add transcripts length and zero coverage transcripts to the output file (required valid bam/sam header) (default: False) [boolean]
転写産物ごとのカウントデータを含むファイルが出力される。デフォルトでは、少なくとも1つのリードがマップされたTranscriptsのみ報告される。フィルタリング後に残ったアラインメントをBAMファイルとして出力するには”-b”オプションを使う。
出力例
transcript_counts.tsv
引用
Accurate expression quantification from nanopore direct RNA sequencing with NanoCount
Josie Gleeson, Adrien Leger, Yair D J Prawer, Tracy A Lane, Paul J Harrison, Wilfried Haerty, Michael B Clark
Nucleic Acids Research, Published: 25 November 2021
関連