De novo リボソームプロファイリングデータによるアクティブに翻訳されたオープンリーディングフレームの同定

サマリー

リボソームを翻訳すると、コドンあたり3つのヌクレオチドがペプチドに解読される。リボソームプロファイリングによって捕捉されたmRNAに沿ったそれらの動きは、特徴的な三重項周期性を示すフットプリントを生成する。このプロトコルでは、RiboCodeを使用してリボソームプロファイリングデータからこの顕著な特徴を解読し、トランスクリプトーム全体でアクティブに翻訳されたオープンリーディングフレームを特定する方法について説明します。

要約

オープンリーディングフレーム(ORF)、特に小さなペプチドをコードし、特定の生理学的文脈の下で積極的に翻訳されるフレームの同定は、文脈依存性翻訳体の包括的な注釈付けにとって重要である。RNA上の翻訳リボソームの結合位置と密度を検出する技術であるリボソームプロファイリングは、翻訳がゲノムワイドスケールでどこで起こっているかを迅速に発見する手段を提供します。しかし、リボソームプロファイリングのために翻訳ORFを効率的かつ包括的に同定することは、バイオインフォマティクスにおいて些細な作業ではありません。ここでは、RiboCodeという名前の使いやすいパッケージで、リボソームプロファイリングデータ内の歪んだあいまいな信号から任意のサイズのORFを積極的に変換するように設計されています。この記事では、以前に公開したデータセットを例にとり、生データの前処理から最終的な出力結果ファイルの解釈まで、RiboCodeパイプライン全体について段階的な手順を説明します。さらに、注釈付きORFの翻訳速度を評価するための、各ORF上のリボソーム密度の可視化および定量化のための手順も詳細に説明される。要約すると、本稿は、翻訳、小型ORF、およびペプチドに関連する研究分野のための有用かつタイムリーな指示である。

概要

最近、ますます多くの研究により、コード遺伝子のORFおよび以前に注釈が付けられた遺伝子から翻訳されたペプチドが、長い非コードRNA(lncRNA)¹、²、³、⁴、⁵、⁶、^7、8などの非コードとして広く産生されることが明らかになりました。これらの翻訳ORFは、環境変化、ストレス、および細胞分化に応答するように細胞によって調節または誘導される¹、⁸、⁹、10^、11、^12、13。いくつかのORFの翻訳産物は、発生および生理学における多様な生物学的プロセスにおいて重要な調節的役割を果たすことが実証されている。例えば、Chngら¹⁴は、心血管の発達に重要なエラベラ(Ela、Apela/Ende/Toddlerとしても知られている)という名前のペプチドホルモンを発見しました。Pauliらは、エラが初期魚の胚における細胞移動を促進するマイトジェンとしても作用することを示唆した¹⁵。Magnyらは、カルシウム輸送を調節し、ショウジョウバエの心臓における規則的な筋肉収縮に影響を及ぼす30アミノ酸未満の2つのマイクロペプチドを報告した¹⁰。

そのようなペプチドがゲノムによってコードされている数と、それらが生物学的に関連性があるかどうかは不明のままである。したがって、これらの潜在的にコードORFの系統的同定は非常に望ましい。しかしながら、進化的保存16,17および質量分析^18,19などの伝統的なアプローチを用いて、これらのORFの産物(すなわち、タンパク質またはペプチド)を直接決定することは、両方のアプローチの検出効率が産生されたタンパク質またはペプチドの長さ、存在量、およびアミノ酸組成に依存するため、困難である。ヌクレオチド分解能でmRNA上のリボソーム占有率を同定する技術であるリボソームプロファイリングの出現は、その長さおよび組成に関係なく、異なる転写産物のコードポテンシャルを評価する正確な方法を提供しました³、^20、21。リボソームプロファイリングを用いて能動的に翻訳するORFを同定するために重要かつ頻繁に使用される特徴は、開始コドンから終止コドンまでのmRNA上のリボソームのフットプリントの3ヌクレオチド(3-nt)周期性である。しかし、リボソームプロファイリングデータには、ORFに沿った低シーケンシングリードとスパースシーケンシングリード、高いシーケンシングノイズ、リボソームRNA(rRNA)汚染など、いくつかの問題があることがよくあります。したがって、このようなデータによって生成される歪んだ曖昧なシグナルは、mRNA上のリボソームのフットプリントの3-nt周期性パターンを弱め、最終的に高信頼の翻訳ORFの同定を困難にする。

「RiboCode」という名前のパッケージは、修正されたウィルコクソン符号化検定とP値統合戦略を適応させ、ORFがオフフレームRPFよりも有意に多くのインフレームリボソーム保護フラグメント(RPF)を有するかどうかを調べた²²。これは、シミュレートされた実際のリボソームプロファイリングデータにおける翻訳体のde novo注釈に対して、非常に効率的で、敏感で、正確であることが実証された。ここでは、このツールを使用して、前の研究によって生成された生のリボソームプロファイリングシーケンシングデータセットから潜在的な翻訳ORFを検出する方法について説明します²³。これらのデータセットは、コントロール(si-Ctrl)およびEIF3E(si-eIF3e)小干渉RNA(siRNA)を導入したMCF-10A細胞のリボソーム占有プロファイルを比較することによって、翻訳におけるEIF3サブユニット「E」(EIF3E)の機能を探索するために使用されていた。これらのサンプルデータセットにRiboCodeを適用することにより、小さなペプチドまたはタンパク質をコードする可能性のある5,633の新規ORFを検出しました。これらのORFは、上流ORF(uORF)、下流ORF(dORF)、重複ORF、新規タンパク質コード遺伝子由来のORF(新規PCG)、新規非タンパク質コード遺伝子由来のORF(新規NonPCG)など、コード領域に対する位置に基づいて様々なタイプに分類された。uORFのRPF読み取り密度は、対照細胞と比較してEIF3E欠損細胞において有意に増加し、これは少なくとも部分的には、能動的に翻訳されるリボソームの濃縮によって引き起こされる可能性がある。EIF3E欠損細胞のコドン25番目から75^番目までの領域に局在するリボソーム蓄積は、初期段階における翻訳伸長の閉塞を示した。このプロトコルはまた、同定されたORF上のリボソームフットプリントの3-nt周期性パターンを調べるために、所望の領域のRPF密度を視覚化する方法も示している。これらの分析は、翻訳ORFを特定し、翻訳の規制を研究する上でのRiboCodeの強力な役割を示しています。

プロトコル

1. 環境設定とリボコードのインストール

1. Linux 端末ウィンドウを開き、conda 環境を作成します。
   conda create -n RiboCode python=3.8
2. 作成した環境に切り替えて、RiboCode と依存関係をインストールします。
   コンダはリボコードをアクティブにする
   conda install -c bioconda ribocode ribominer sra-tools fastx_toolkit cutadapt bowtie star samtools

2. データ準備

1. ゲノム参照ファイルを取得します。
   1. 参照シーケンスについては、https://www.ensembl.org/index.html のアンサンブルWebサイトにアクセスし、トップメニューの[ダウンロード]と左側のメニュー[FTPダウンロード]をクリックします。提示された表で、列 DNA (FASTA) の FASTA と、種が人間である行の FASTA をクリックします。開いたページで、Homo_sapiensのリンクをコピーします。GRCh38.dna.primary_assembly.fa.gz、ターミナルでダウンロードして解凍します。
      wget -c \
      http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
      gzip -d Homo_sapiens。GRCh38.dna.primary_assembly.fa.gz
   2. 参照注釈については、最後に開いた Web ページの [遺伝子セット] 列の GTF を右クリックします。Homo_sapiensのリンクをコピーします。GRCh38.104.gtf.gzを使用してダウンロードします。
      wget -c \
      http://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz
      gzip -d Homo_sapiens。GRCh38.104.gtf.gz
      注:GTFファイルは、3レベルの階層に編成されたゲノム注釈を含むため、アンサンブルのウェブサイトからGTFファイルを取得することをお勧めします、すなわち、各遺伝子にはエクソンとオプションの翻訳(例えば 、コード配列[CDS]、翻訳開始部位、翻訳終了部位)を含む転写産物が含まれています。UCSC または NCBI から取得した GTF ファイルなど、遺伝子または転写産物の注釈が欠落している場合は、 GTFupdate を使用して、完全な親子階層注釈 ( GTFupdate original.gtf > updated.gtf) を含む更新された GTF を生成します。.gff 形式のアノテーション・ファイルの場合は、AGAT ツールキット ²⁴ またはその他のツールを使用して .gtf 形式に変換します。
2. rRNA配列を取得します。
   1. UCSC Genome Browser を https://genome.ucsc.edu で開き、[ツール |] をクリックします。ドロップダウンリストのテーブルブラウザ。
   2. 開いたページで、クレードに 哺乳類 、ゲノムに ヒト 、グループに すべてのテーブル、テーブルに rmask 、領域に ゲノム を指定します。フィルターの場合は、[ 作成 ] をクリックして新しいページに移動し、 rRNA と一致するように repClass を設定します。
   3. [ 送信] をクリックし、出力形式を [シーケンス ] に設定し、ファイル名を hg38_rRNA.fa として出力します。最後に、[ 出力|の取得] をクリックします。rRNA配列を取得するには配列 を取得します。
3. シーケンス読み取りアーカイブ(SRA)からリボソームプロファイリングデータセットを取得します。
   1. si-eIF3e治療群の複製サンプルをダウンロードし、名前を変更します。
      ファスト Q ダンプ SRR9047190 SRR9047191 SRR9047192
      mv SRR9047190.fastq si-eIF3e-1.fastq
      mv SRR9047191.fastq si-eIF3e-2.fastq
      mv SRR9047192.fastq si-eIF3e-3.fastq
   2. 制御グループの複製サンプルをダウンロードし、名前を変更します。
      ファスト Q ダンプ SRR9047193 SRR9047194 SRR9047195
      mv SRR9047193.fastq si-Ctrl-1.fastq
      mv SRR9047194.fastq si-Ctrl-2.fastq
      mv SRR9047195.fastq si-Ctrl-3.fastq
      注: これらのサンプル データセットの SRA アクセッション ID は、GSE131074 を検索して Gene Expression Omnibus (GEO) Webサイト²⁵ から取得しました。

3. アダプターをトリムし、rRNA汚染を除去する

1. (オプション)シーケンス・データからアダプターを取り外します。この場合のように、アダプター・シーケンスが既にトリミングされている場合は、このステップをスキップしてください。それ以外の場合は、 cutadaptを使用して 、読み取りからアダプタをトリミングします。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   する
   cutadapt -m 15 --match-read-wildcards -a CTGTAGGCACCATCAAT \
   -o ${i}_trimmed.fastq ${i}.fastq
   完成です
   注: -a パラメーターの後のアダプター配列は、cDNA ライブラリーの準備によって異なります。15 より短い読み取り ( -m で指定) は、リボソームで保護されたフラグメントが通常このサイズより長いため、破棄されます。
2. 以下の手順でrRNA汚染を除去します。
   1. インデックスrRNA参照配列:
      蝶ネクタイビルド-f hg38_rRNA.fa hg38_rRNA
   2. 読み取りを rRNA 参照に合わせ、rRNA から発生した読み取りを除外します。
      for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
      する
      蝶ネクタイ -n 0 -y -a --norc --best --strata -S -p 4 -l 15 \
      --un=./${i}_noncontam.fastq hg38_rRNA -q ${i}.fastq ${i}.aln
      完成です
      -p は、タスクを並列に実行するためのスレッドの数を指定します。RPF 読み取りのサイズが比較的小さいことを考慮すると、報告されたアライメントが最適であることを保証するために、他の引数 (-n、-y、-a、-norc、--best、--strata、-l など) を指定する必要があります。詳細については、ボウタイのウェブサイト²⁶を参照してください。

4. クリーンな読み取り値をゲノムに合わせる

1. ゲノムインデックスを作成します。
   mkdir STAR_hg38_genome
   STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./STAR_hg38_genome --genomeFastaFiles Homo_sapiens.GRCh38.dna.primary_assembly.fa --sjdbGTFfile Homo_sapiens.GRCh38.104.gtf
2. クリーン読み取り(rRNA汚染なし)を作成した参照に合わせます。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   する
   STAR --runThreadN 8 --outFilterType Normal --outWigType wiggle --outWigStrand Stranded --outWigNorm RPM --outFilterMismatchNmax 1 --outFilterMultimapNmax 1 --genomeDir STAR_hg38_genome --readFilesIn ${i}_noncontam.fastq --outFileNamePrefix ${i}.--outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM GeneCounts --outSAMattributes All
   完成です
   注:テンプレート化されていないヌクレオチドは、逆転写酵素²⁷によって各読み取りの5'末端に頻繁に追加され、デフォルトでソフトクリッピングを実行するため、STARによって効率的にトリミングされます。STAR のパラメータについては、STAR マニュアル ²⁸ に説明があります。
3. 並べ替えとインデックスの配置ファイル。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   する
   samtools sort -T ${i}.Aligned.toTranscriptome.out.sorted \
   -o ${i}.Aligned.toTranscriptome.out.sorted.bam \
   ${i}.Aligned.toTranscriptome.out.bam
   samtools index ${i}.Aligned.toTranscriptome.out.sorted.bam
   samtools index ${i}.Aligned.sortedByCoord.out.bam
   完成です

5. RPFのサイズ選択とそのPサイトの識別

1. トランスクリプト注釈を準備します。
   -g prepare_transcriptsをHomo_sapiensします。GRCh38.104.gtf \
   -f Homo_sapiens。GRCh38.dna.primary_assembly.fa -o RiboCode_annot
   注: このコマンドは、GTF ファイルから mRNA 転写産物の必要な情報を収集し、FASTA ファイルからすべての mRNA 転写産物の配列を抽出します (各転写産物は、GTF ファイルで定義された構造に従ってエクソンをマージすることによってアセンブルされます)。
2. 特定の長さのRPFを選択し、そのPサイト位置を特定します。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   する
   メタプロット -a RiboCode_annot -r ${i}.Aligned.toTranscriptome.out.bam \
   -o ${i} -f0_percent 0.35 -pv1 0.001 -pv2 0.001
   完成です
   注: このコマンドは、注釈付き翻訳開始 (または停止) コドンの周囲に、各長さの整列された読み取り値の 5' 末端の集約プロファイルをプロットします。リード長依存性P部位は、主要リードの5'末端と開始コドンとの間のオフセット距離の分布プロット(例えば、図1B)を調べることによって手動で決定することができる。RiboCodeはまた、サンプルごとに構成ファイルを生成し、そこでは、有意な3nt周期性パターンを示す読み取りのPサイト位置が自動的に決定される。パラメーター -f0_percent、-pv1、および -pv2 は、読み取りフレームで強化された RPF 読み取りを選択するための比率しきい値と p 値カットオフを定義します。この例では、29、30、および 31 nt 読み取りの 5' 末端からの +12、+13、および +13 ヌクレオチドは、各構成ファイルで手動で定義されています。
3. 各サンプルの構成ファイルを編集し、それらをマージする
   メモ: 固有のORFのコンセンサスセットを生成し、後続の分析を実行するのに十分な読み取りカバレッジを確保するために、前のステップで選択したすべてのサンプルの読み取り値がマージされます。 merged_config.txtファイル (補足ファイル1)で定義された特定の長さの読み取りとそのPサイト情報は、次のステップでORFの翻訳可能性を評価するために使用されます。

6. 翻訳 ORFに注釈を付ける

1. リボコードを実行します。
   RiboCode -a RiboCode_annot -c merged_config.txt -l yes -g \
   -o RiboCode_ORFs_result -S ATG -M 5 -A CTG,GTG,TTG
   ここで、このコマンドの重要なパラメーターは次のとおりです。
   -c、入力ファイルのパス、選択した読み取りとそのPサイトの情報を含む設定ファイル。
   −lは、終止コドンの上流に複数の開始コドンを有する転写産物について、最長のORF(最も遠位の開始コドンから終止コドンまでの領域)がそれらの翻訳可能性を評価するために使用されるか否か。 noに設定すると、開始コドンが自動的に決定される。
   -sは、ORFの識別に使用される正準開始コドンである。
   −Aは、(任意選択で)ORF同定のために使用される非正準開始コドン(例えば、ヒトのためのCTG、GTG、およびTTG)と、他の種のミトコンドリアまたは核において異なる可能性がある²⁹。
   -mは、ORFの最小長(すなわち、アミノ酸)である。
   -o は、予測される ORF の詳細を含む出力ファイル名の接頭部です (補足ファイル 2)。
   -gと-bは、予測ORFをそれぞれgtfまたはベッド形式に出力します。

7. (オプション) ORF の定量化と統計

1. 各 ORF で RPF 読み取りをカウントします。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   する
   ORFcount -g RiboCode_ORFs_result_collapsed.gtf \
   -r ${i}.Aligned.sortedByCoord.out.bam -f 15 -l 5 -m 25 -M 35 \
   -o ${i}_ORF.counts -s yes -c intersection-strict
   完成です
   注: ORF の開始と終了の周囲に蓄積する可能性のあるリボソームを除外するために、最初の 15 個 ( -f で指定) と最後の 5 個のコドン ( -l で特定) に割り当てられた読み取り数はカウントされません。オプションで、カウントされる RPF の長さは、 25 ~ 35 nt (RPF の一般的なサイズ) の範囲に制限されます。
2. RiboCodeを使用して、検出されたORFの基本統計を計算します。
   Rscript RiboCode_utils。R
   注: RiboCode_utils。R (補足ファイル3)は、識別されたORFの数のカウント、ORFの長さの分布の表示、正規化されたRPF密度の計算(すなわち、RPKM、マップされた読み取り100万回あたりのキロベースあたりの読み取り)など、RiboCode出力に関する一連の統計を提供します。

8. (オプション) 予測ORFの可視化

1. ENSG00000100902_35292349_35292552_67からその転写産物上の所望のORFの開始コドンおよび終止コドン(例えば、 RiboCode_ORFs_result_collapsed.txt )の相対位置を取得する(補足ファイル3)。次に、RPF読み取りの密度をORFにプロットします。
   plot_orf_density -a RiboCode_annot -c merged_config.txt -t ENST00000622405 \
   -s 33 -e 236 --開始コドン ATG -o ENSG00000100902_35292349_35292552_67
   ここで 、-s と -e は、プロット ORF の平行移動の開始位置と終了位置を指定します。 --start-codonは ORFの開始コドンを定義し、図のタイトルに現れる。 -o は、出力ファイル名の接頭部を定義します。

9. (オプション)リボマイナーを用いたメタジーン解析

注: 以下の手順に従って、メタジーン解析を実行して、同定された注釈付きORFの翻訳に対する EIF3E ノックダウンの影響を評価します。

1. RiboMinerのトランスクリプトアノテーションを生成し、RiboCodeによって生成されたアノテーションファイルに基づいて各遺伝子の最も長いトランスクリプトを抽出します(ステップ5.1)。
   OutputTranscriptInfo -c RiboCode_annot/transcripts_cds.txt \
   -g Homo_sapiens。GRCh38.104.gtf -f RiboCode_annot/transcripts_sequence.fa \
   -o longest.transcripts.info.txt -O all.transcripts.info.txt
2. リボマイナーの構成ファイルを準備します。RiboCodeの metaplots コマンドで生成された設定ファイルをコピーし(ステップ5.4)、名前を「RiboMiner_config.txt」に変更します。次に、 補足ファイル 4 に示されている形式に従って変更します。
3. リボマイナーを用いたメタジーン解析
   1. MetageneAnalysisを使用して、転写産物全体のRPFの密度の集計および平均プロファイルを生成します。
      MetageneAnalysis -f RiboMiner_config.txt -c longest.transcripts.info.txt \
      -o MA_normed -U コドン -M RPKM -u 100 -d 400 -l 100 -n 10 -m 1 -e 5 --norm yes \
      -y 100 --タイプ UTR
      重要なパラメータがある場合: --type、 CDS または UTR 領域の分析。 --norm、読み取り密度を正規化したかどうか。 -y、各転写産物に用いたコドンの数; −U、 RPF密度をコドン レベルまたは nt レベルのいずれかでプロットする;− u および −dは、開始コドンまたは終止コドンに対する分析領域の範囲を定義する;−l、CDSの最小長(すなわち、コドンの数); -M は、トランスクリプトフィルタリングのモードであり、 カウント または RPKM のいずれかです。 分析 用の CDS の最小カウントまたは RPKM正規化された領域における CDS の -m 最小カウントまたは RPKM −eは、正規化領域から除外されたコドンの数である。
   2. 対照細胞とeIF3欠損細胞におけるmRNA上のリボソーム占有率を比較するための一連のpdfファイルを生成する。
      PlotMetageneAnalysis -i MA_normed_dataframe.txt -o MA_normed \
      -g si-Ctrl,si-eIF3e -r si-Ctrl-1,si-Ctrl-2,si-Ctrl-3__si-eIF3e-1,si-eIF3e-2,si-eIF3e-3 -u 100 -d 400 --mode mean
      注: PlotMetageneAnalysis は、一連の pdf ファイルを生成します。MetageneAnalysisとPlotMetageneAnalysisの使用法に関する詳細は、RiboMinerのウェブサイト³⁰で入手できます。

結果

リボソームプロファイリングデータセットの例は、アクセッション番号GSE131074でGEOデータベースに寄託されました。このプロトコルで使用されるすべてのファイルとコードは、補足ファイル 1 ~ 4 から入手できます。リボソームプロファイリングデータセット²³の公開セットにRiboCodeを適用することにより、コントロールおよびEIF3E siRNAで処理したMCF-10A細胞で能動的に翻訳された新規ORFを同定しました。翻訳リボソームによって結合される可能性が最も高いRPF読み取りを選択するために、シーケンシング読み取りの長さを調べ、既知の翻訳遺伝子上にマッピングされたRPFを使用してメタジーン解析を実施した。リードの長さの頻度分布は、ほとんどのRPFが25〜35nt(図1A)であり、予想通りリボソームによってカバーされるヌクレオチド配列に対応することを示した。異なる長さのRPFのPサイト位置は、それらの5'末...

ディスカッション

リボソームプロファイリングは、細胞におけるリボソームの作用をゲノムスケールで研究する前例のない機会を提供します。リボソームプロファイリングデータによって運ばれる情報を正確に解読することで、遺伝子または転写産物のどの領域が積極的に翻訳されているかについての洞察が得られる可能性があります。このステップバイステップのプロトコルは、パッケージのインストール、データ準備、コマンドの実行、結果の説明、データの視覚化など、RiboCodeを使用してリボソームプロファイリングデータを詳細に分析する方法に関するガイダンスを提供します。RiboCodeの解析結果は、翻訳が広範に行われ、コード遺伝子の注釈なしORFおよび以前は非コーディングであると考えられていた多くの転写産物で起こることを示した。下流の解析は、翻訳が起こるにつれて、リボソームが予測されたORFに沿って3ヌクレオチドステップで移動するという証拠を提供しました。しかし、翻訳の過程や産生されたペプチドが何らかの機能を果たすかどうかは不明である。それにもかかわらず、ゲノム上のORFを翻訳する正確な注釈は、これまで特徴付けられていなかった転写産物の機能を特定するエキサイティングな機会を生み出す可能性が...

材料

この記事で使用された材料の一覧

名前	会社	カタログ番号	コメント
Linux Anyを実行しているコンピュータ/サーバ		-
Anaconda または Miniconda	Anaconda	-	Anaconda: https://www.anaconda.comミニコンダ:https://docs.conda.io/en/latest/miniconda.html
R	R	Foundation-https://www.r-project.org/
Rstudio		Rstudio-https://www.rstudio.com/