pythonってすごいね

RNAseqを用いた遺伝子発現量解析、機械学習を用いた回帰、分類などの解析を中心に記事を書いていきたいです！

StringTieを用いたリード数の計算 (RNAseq)

RNAseq

マッピング後のSAM、BAMから、各遺伝子ごとにリード数を算出します。今回は、StringTieについてまとめてみました。

SAMファイルは重たいので、BMAファイルに変換した方が解析がしやすいと思います。

The main input of the program is a BAM file with RNA-Seq read mappings which must be sorted by their genomic location (for example the accepted_hits.bam file produced by TopHat or the output of HISAT2

準備

マッピング後のファイルは事前にソートを行なっておく

samtools sort - input.sorted

実行コマンド

stringtie [input_sort.bam] [options]

使いそうな設定をいくつか抜粋

オプション 説明
-o 出力するGTFファイルパスの設定
-p スレッド数　デフォルト１
-G GTF、GFF3
-f 最も多く存在するアイソフォームに対して、最小のアイソフォーム転写物の存在量を設定　デフォルト0.01
-m 新規転写物の最小長を設定　デフォルトは２００
-c 出力の転写物名をレファレンスを基に表示する -Gオプションが必要
-s エキソンが一つしかない転写物の最小カバレッジ量を指定する
-A タブ区切りの出力を行う際の保存場所の指定（出力形式は下記を参照）

GTF、GFFファイルを指定する場合は、マッピングに使用したゲノムの染色体名と一致しているかを確認する。

小文字、大文字の違いで上手く処理できない場合がある。

例文

stringtie -p 8 -G  TAIR10_for_stringtie.gff -o count.txt -A table.txt sorted.bam

出力形式

f:id:pytyonttesugoine:20210222182740p:plain

  1. seqname: マップされたクロモソーム名
  2. source: 基となったGTFファイル
  3. feature: RNAの特徴 (e.g., exon, transcript, mRNA, 5'UTR).
  4. start: 開始位置 1-based index
  5. end: 終了位置 1-based index
  6. score: 現在は使用されていない。デフォルト表記は１０００？
  7. strand: ストランド
  8. frame: 表記なし（本来はCDSのフレーム）
  9. attributes: その他の特徴を記載
    • gene_id: 参照ファイルを基にした名前
    • transcript_id: 参照ファイルを基にした名前
    • exon_number: 何番目のエキソンかを示す
    • reference_id
    • ref_gene_id
    • ref_gene_name
    • cov: 各塩基単位のリード数
    • FPKM
    • TPM

FPKM、TPMに関しては下記のサイトを参照
RNA-seqブログ

-Aを用いた出力

f:id:pytyonttesugoine:20210222182736p:plain

  • Column 1 / 遺伝子ID
  • Column 2 / 遺伝子名
  • Column 3 / 使用したReferenceの名前
  • Column 4 / ストランド
  • Column 5 / 開始位置 (1-based index).
  • Column 6 / 終了位置 (1-based index).
  • Column 7 / Coverage 各塩基単位のリード数
  • Column 8 / FPKM
  • Column 9 / TPM

参考

StringTie Manual