質問:
エクソンとCDSを含むGTFファイルからタンパク質コーディング遺伝子を含むGTFを導出する
Tom Kelly
2018-12-10 17:48:42 UTC
view on stackexchange narkive permalink

互換性のあるファイルが必要な理由

Rパッケージで velocytoを実行して、RNA速度（細胞軌道）を分析しようとしています。シングルセルRNASeqデータ。セルレンジャーを使用して10倍のゲノミクスデータから単一セル分析を実行しました。

読み取りを正常に調整して織機ファイルを取得し、これらをRにインポートしました。ビネット。ただし、「遺伝子構造」に基づくRNA速度分析を再現することはできません。

例とは異なる生物（人間やマウスではない）で作業しているため、例で提供されている注釈データは再現されません。作業。この生物の最新の注釈用のGTFファイルがあります。ただし、機能として「エクソン」と「CDS」のみが含まれています。これが問題の原因のようです。 「find.ip.sites」関数には、「features」=「gene」および「attributes」の1つが「protein_coding」であるGTFが必要です。これらの要件はvelocyto.R関数にハードコードされています。

AtRTD2データセットから次のGTFファイルがあります。染色体ラベルは、Rの織機ファイルと一致します。

  Chr1TAIR10エクソン36313913。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1TAIR10エクソン39964276。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1TAIR10エクソン44864605。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1TAIR10エクソン47065095。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1TAIR10エクソン51745326。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1TAIR10エクソン54395899。 +。 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 37603913。 + 0transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1 TAIR10 CDS 39964276。 + 2transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 44864605。 + 0transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 47065095。 + 0transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 51745326。 + 0transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 54395630。 + 0transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1Araport11エクソン67887069。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン71577450に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン75647649に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン77627835に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン79427987に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン82368325に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン84178464に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11エクソン85718737に注意してください。 -。 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11 CDS 73157450に注意してください。 -1つのtranscript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」に注意してください。
Chr1 Araport11 CDS 75647649。 -0transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11 CDS 77627835に注意してください。 --2transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11 CDS 79427987に注意してください。 -0transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; 「ARV1」; Chr1Araport11 CDS 82368325に注意してください。 -0transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020";注「ARV1」;

GTFまたはGFF2の仕様は次のとおりです。

フィールド

フィールドはタブで区切る必要があります。また、各計画線の最後のフィールドを除くすべてに値が含まれている必要があります。 「空の」列は「。」で示す必要があります。

  1. seqname-染色体または足場の名前。染色体名は、「chr」プレフィックスの有無にかかわらず指定できます。重要な注意：seqnameは、Ensembl内で使用されるもの、つまり、種やアセンブリなどの追加コンテンツなしで、標準の染色体名またはスキャフォールドIDなどのEnsembl識別子である必要があります。以下のGFF出力の例を参照してください。

  2. source-この機能を生成したプログラムの名前、またはデータソース（データベースまたはプロジェクト名）

  3. feature-機能タイプ名、例：遺伝子、バリエーション、類似性

  4. start-機能の開始位置。シーケンス番号は1から始まります。

  5. end-フィーチャの終了位置。シーケンス番号は1から始まります。

  6. スコア-浮動小数点値。

  7. ストランド- +（順方向）または-（逆方向）として定義されます。

  8. フレーム-「0」、「1」、または「2」のいずれか。 「0」は、特徴の最初の塩基がコドンの最初の塩基であることを示し、「1」は、2番目の塩基がコドンの最初の塩基であることを示します。

  9. 属性-タグと値のペアのセミコロンで区切られたリスト。各機能に関する追加情報を提供します。

    10. ol>

私が探しているものfor

エクソンとCDSからタンパク質コーディング遺伝子を含む互換性のあるGTFファイルを導出する方法はありますか？機能に遺伝子を含み、属性にprotein_codingを含むGTFを作成したいと思います。既存のツールまたはスクリプトでこれを行うことは可能ですか？

これまでに試したこと

「find.ip.sites」関数のソースコードを変更して実行できますこれらの機能が欠落している私のGTFファイル。ただし、これにはRcppで記述されたパッケージから内部関数を実行する必要があり、ワークフローがパッケージの将来の更新を上書きすることを意味します。十分な長さのイントロンが識別されていないため、残りのビネットを実行するとエラーが返されます（しきい値を低く設定することも、と呼ばれる一般線形モデルと互換性がありません）。したがって、関数のソースコードに警告するよりも、互換性のあるGTFまたはGFF3ファイルを生成する方が良いと思います。

GTFファイルを対象としていますが、パッケージ関数はGTFについて説明されているにもかかわらず、GFF3ファイルをインポートできます。ドキュメントで。 Cufflinksのgffreadを使用してGFF3ファイルを生成し、「feature」=「mRNA」を「gene」に置き換えて「protein_coding」を追加してみました。これは、速度アルゴリズムの実行時にもエラーを返します。セルレンジャーの入力として使用されるGTFファイル、またはセルレンジャーによって生成されるファイルのいずれに対しても機能しません。エクソンとCDSの注釈のみを含むGTFファイルに基づいてタンパク質コーディング遺伝子とイントロン/エクソンの境界に注釈を付ける方法はありますか？

カフリンク2.2.1を使用して、GFF3が生成されました

  gffread --E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gffsed -i '/ mRNA / s / gene_name = AT / gene_type = protein_coding; gene_name = AT / g' AtRTD2_19April2016.gff'sed -i gene / g'AtRTD2_19April2016.gff

これは私が試したGFF3ファイルです：

  ＃gffread --E AtRTD2_19April2016.gtf -o- ## gff-version3Chr1TAIR10遺伝子36315899。 +。 ID = AT1G01010.1; geneID = AT1G01010; gene_type = protein_coding; gene_name = AT1G01010Chr1TAIR10エクソン36313913。 +。親= AT1G01010.1Chr1TAIR10エクソン39964276。 +。親= AT1G01010.1Chr1TAIR10エクソン44864605。 +。親= AT1G01010.1Chr1TAIR10エクソン47065095。 +。親= AT1G01010.1Chr1TAIR10エクソン51745326。 +。親= AT1G01010.1Chr1TAIR10エクソン54395899。 +。親= AT1G01010.1Chr1 TAIR10 CDS 37603913。 +0親= AT1G01010.1Chr1 TAIR10 CDS 39964276。 +2親= AT1G01010.1Chr1 TAIR10 CDS 44864605。 +0親= AT1G01010.1Chr1 TAIR10 CDS 47065095。 +0親= AT1G01010.1Chr1 TAIR10 CDS 51745326。 +0親= AT1G01010.1Chr1 TAIR10 CDS 54395630。 + 0 Parent = AT1G01010.1Chr1Araport11遺伝子67888737。 +。 ID = AT1G01020_P2; geneID = AT1G01020; gene_type = protein_coding; gene_name = AT1G01020 Chr1Araport11エクソン67887069。 -。親= AT1G01020_P2Chr1Araport11エクソン71577450。 -。親= AT1G01020_P2Chr1Araport11エクソン75647649。 -。親= AT1G01020_P2Chr1Araport11エクソン77627835。 -。親= AT1G01020_P2Chr1Araport11エクソン79427987。 -。親= AT1G01020_P2Chr1Araport11エクソン82368325。 -。親= AT1G01020_P2Chr1Araport11エクソン84178464。 -。親= AT1G01020_P2Chr1Araport11エクソン85718737。 -。親= AT1G01020_P2Chr1 Araport11 CDS 73157450。 -1親= AT1G01020_P2Chr1 Araport11 CDS 75647649。 -0親= AT1G01020_P2Chr1 Araport11 CDS 77627835。 --2親= AT1G01020_P2
Chr1 Araport11 CDS 79427987。 -0親= AT1G01020_P2Chr1 Araport11 CDS 82368325。 --0 Parent = AT1G01020_P2

GTFまたはGFF2の仕様は次のとおりです。

GF3 aの仕様>は次のとおりです。

フィールド

フィールドはタブで区切る必要があります。また、各計画線の最後のフィールドを除くすべてに値が含まれている必要があります。 「空の」列は「。」で示す必要があります。

  1. seqid-染色体または足場の名前。染色体名は、「chr」プレフィックスの有無にかかわらず指定できます。重要な注意：seq IDは、Ensembl内で使用されるもの、つまり、種やアセンブリなどの追加コンテンツなしで、標準の染色体名またはスキャフォールドIDなどのEnsembl識別子である必要があります。以下のGFF出力の例を参照してください。

  2. source-この機能を生成したプログラムの名前、またはデータソース（データベースまたはプロジェクト名）

  3. type-機能のタイプ。 SOFAシーケンスオントロジーからの用語またはアクセッションである必要があります

  4. start-機能の開始位置で、シーケンス番号は1から始まります。

  5. end-機能の終了位置。シーケンス番号は1から始まります。

  6. score-浮動小数点値。

  7. ストランド-+（順方向）または-（逆方向）として定義されます。

  8. フェーズ-「0」、「1」、または「2」のいずれか。 「0」は、フィーチャの最初の塩基がコドンの最初の塩基であることを示し、「1」は、2番目の塩基がコドンの最初の塩基であることを示します。

  9. 属性-タグと値のペアのセミコロンで区切られたリスト。各機能に関する追加情報を提供します。これらのタグの一部は事前定義されています。 ID、名前、エイリアス、親-詳細については、GFFのドキュメントを参照してください。

    10. ol>
質問を[編集]して、2つのGTFファイルのそれぞれを数行表示して、意味がより明確にわかるようにすると、おそらく役に立ちます。
これにはおそらくAWKを使用できます。これは、「CDS」に似た「遺伝子」機能であるか、「最初のエクソンの開始」から「最後のエクソンの終了」までです。
質問と「私が探しているもの」は少し紛らわしいと思います。具体的には、「探しているもの」でGTFを生成したいが、タイトルのGTFからタンパク質のアミノ酸配列を取得したいようです。タイトルと「私が探しているもの」を明確にしたり、一貫性を持たせる方法はありますか？
[Ensembl植物のシロイヌナズナの完全なGTF]（ftp://ftp.ensemblgenomes.org/pub/plants/current/gtf/arabidopsis_thaliana/）があります。
@Emily_Ensembl,のおかげで、この参照配列に基づいて「グローバル」RNA速度を計算する速度を得ることができました。他の誰かがこの問題に遭遇した場合、Rコマンドと互換性を持たせるために、シェルユーティリティを使用してloomファイルを同じGTFファイルで再計算する必要があります。
1 回答:
Ian Sudbery
2018-12-11 16:23:38 UTC
view on stackexchange narkive permalink

あなたの場合、@ Emily_Ensemblのアドバイスに従い、EnsemblのArabidopsisGTFを使用することを強くお勧めします。しかし、将来の参考のために、Ensembl GTFが利用できない場合は、 cgat

cgatモジュールの gtf クラスを使用してこのようなものを構築できます。依存関係は、こちらの手順に従ってインストールできます。具体的には、インストールスクリプトをダウンロードして実行できます。これにより、必要なバージョンの依存関係を実行するための専用のconda環境がセットアップされます。詳細なドキュメントはここにあります。

 ＃ダウンロードインストールスクリプト：curl-O https://raw.githubusercontent.com/cgat-developers/cgat-apps/ master / install.sh＃は開発バージョンをインストールします（推奨、まだ製品バージョンはありません）：bash install.sh --devel

これをインストールするには、python3とanacondaが必要です。 conda環境で次のスクリプトを実行する必要がある場合があります。

  ＃show available environmentconda info --envs＃activate by pathsource activate / home / user / local / bin / conda-install / envs / cgat-a

次に、python3スクリプトとして以下を実行します。 convert_gtf.py などのファイルの内容をコピーします。次に、gtfファイルのターミナルで実行します。

  pythonconvert_gtf.pygenes.gtf >genes_new.gtf  
  #python 3.6.4import sysfrom cgat import GTFfrom cgatcore import iotoolsfor gene in GTF.flat_gene_iterator（GTF.iterator（iotools.open_file（sys.argv [1]））、strict = False）：gene_start = min（ e.start for e in gene）gene_end = max（e.end for e in gene）if any（e.feature == "CDS" for e in gene）：gene_type = "protein_coding" else：gene_type = "non_coding" for遺伝子内のエクソン：exon_line.setAttribute（ "gene_biotype"、gene_type）gene_line = GTF.Entry（）。fromGTF（gene [0]）gene_line.feature = "gene" gene_line.start = gene_start gene_line.end = gene_end
＃遺伝子系統のtranscript_idが何であるかは明確ではありません。 ＃技術的に言えば、GTFエントリにはtranscript_idsを含める必要があるため、Ensembl GTFはGTF標準に準拠していませんが、＃アンサンブルGTFには準拠していません。 if hasattr（gene_line、 "gene_name"）：gene_line.attributes ["gene_name"] = gene [0] .gene_name else：gene_line.attributes ["gene_name"] = gene [0] .gene_id gene_line.attributes ["gene_biotype"] = Gene_type print（str（gene_line）+ "\ n"）for exon in gene：print（str（exon_line）+ "\ n"）
それを知るのは素晴らしいことです。これはRでのツールの使用に関係するため、ライブラリをインストールしてPythonに慣れていない人のために実行する手順があると便利です。
このスクリプトをテストし、小さな変更を加えました。これが私にとってうまくいったことです（GTFを隠すために：それがベロサイトと互換性があるかどうかはまだわかりません）。
修正していただきありがとうございます。最近状況が変わったので、バージョンがわからなくなってしまいました。 CGATチームはインストールの簡素化に取り組んでいると思います。いくつかのポイント： `if hasattr（gene_line、" gene_name） `ビットで何を達成しようとしているのかよくわかりません：おそらく、すでに` gene_name`属性がある場合は、設定する必要はありません。また、結果の行には `transcript_id`が含まれないことに注意してください。これはEnsembl GTFにも当てはまりますが、これらのGTFがGTF仕様を意味するわけではないことを意味します。
遺伝子名が見つからないという問題があり、遺伝子記号の代わりにTAIRIDが使用されることがよくあります。したがって、これはシロイヌナズナに固有のものである可能性があります。最も重要なのは、これらの行が定義されていないため、エクソンのループを修正することです。
属性に `gene_name`を追加する必要があるかもしれないことを私は理解していません。それは、現在、` gene_name`の存在をテストし、すでに定義されている値に設定しているだけです。 1。これは冗長なようです。


ⓘ
このQ＆Aは英語から自動的に翻訳されました。オリジナルのコンテンツはstackexchangeで入手できます。これは、配布されているcc by-sa 4.0ライセンスに感謝します。
Loading...