Pythonバイオ/ツール/RNAseq-StringTie

サイトトップ

Pythonバイオ/ツール
├配列の切出し・加工
├RNAseq①リードデータ取得とクオリティチェック・トリミング
├RNAseq②リードの参照配列へのマッピング
├RNAseq③発現解析～カウントデータの解析
├RNAseq④発現解析（続き）
├RNAseq～GeneBankファイルからfasta+gffファイルへ変換
├RNAseq～VCFとBCF
├RNAseq～AP012030
└Variant-Calling

メニュー編集

RNAseqでStringTieを試す。 †

https://f1000research.com/articles/4-1521/v2:Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [version 2; peer review: 2 approved]

↑

インストール †

git clone https://github.com/gpertea/stringtie
cd stringtie
make release

which stringtie
/usr/local/bin/stringtie

ついでに

git clone https://github/gpertea/gffcomp
cd gffcomp
make
sudo ln -s ....../gffcomp/gffcomp /usr/local/bin/gffcomp  # makeではinstallしない

↑

実行 †

be sure to run HISAT2 with the --dta option for alignment, or your results will suffer.

cd ~/src/RNAseq-Saccha/Saccha
ls -lt *.bam
SRR453571.tagged.bam
 ...
SRR453566.tagged.bam

# 入力bamファイルはソートされていなければならない
samtools sort -@ 16 -o SRR453566.tagged.sorted.bam SRR453566.tagged.bam

stringtie SRR453566.tagged.sorted.bam -o SRR453566.gtf -p 16 -G s288c_e.gff -l SRR453566

stringtie SRR453566.tagged.sorted.bam -o SRR453566.gtf -p 16 -G s288c_e.gff -l SRR453566 -A SRR453566.abund.tab

Aを付けると gene abundance table を出力

↑

サンプルごとのGTFファイルをマージ †

stringtie --merge -p 8 -G s288c_e.gff -o SRR453566-71_merged.gtf SRR453566.gtf SRR453567.gtf SRR453568.gtf SRR453569.gtf SRR453570.gtf SRR453571.gtf

ついでに

gffcompare -r s288c_e.gff -G -o SRR453566-71_diff SRR453566-71_merged.gtf

↑

Salmonを試す †

Salmon provides accurate, fast, and bias-aware transcript expression estimates using dual-phase inference

Salmon provides accurate, fast, and bias-aware transcript expression estimates using dual-phase inference(PDF)

Salmon — Salmon 0.14.0 documentation

Salmon (A. thaliana, paired-end RNA-Seq) | Salmon を利用した転写産物の発現量定量
超高速にRNA seqのリードカウント（定量）を行う salmon - macでインフォマティクス]
salomonは単独で「アラインメント」と「定量」の両方を行う、すなわち、indexがついたtranscriptsのリファレンスとFASTQを入力として受け取り、中間のアラインメントファイルを生成せずに直接定量を実行する。結果かなりの時間とスペースを節約することができる。
RNA-seq 発現変動遺伝子解析：salmonとTCC-GUIを使って - quiita
salmonを用いることにより、alighnerを使用せず、転写物配列のインデックスを利用して、readsをカウントする。 TCC-GUIを用いて、PCA分析、クラスタリング、シルエット分析、これらの可視化、発現変動遺伝子(DEG)抽出を行う。どの程度、コンピュータ資源と解析時間を節約して、DEG抽出が可能かを確認する。

↑

Kallistoを試す †

Near-optimal RNA-Seq quantification(Nature) Near-optimal RNA-Seq quantification Near-optimal RNA-Seq quantification(PDF)

Kallisto About

Kallisto - macでインフォマティクス

kallisto (A. thliana, paired-end RNA-Seq) | kallisto を用いた A. thaliana paired-end リードの転写産物の定量

Kallistoを用いたRNA-seq解析パイプライン – 大阪大学医学部 Python会 <===

↑

Indexを作る †

Kallistoを用いたRNA-seq解析パイプライン
リファレンスのダウンロード
kallistoでは、transcriptにシュードアラインメントするので、リファレンスにはcDNAを用います。今回はGenCodeGenesのヒトtranscript sequencesのデータを用いました。

$ wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.transcripts.fa.gz

＞＞＞　（リファレンスの？）cDNAが必要

Salmon (A. thaliana, paired-end RNA-Seq)に曰く
転写産物（cDNA）のリファレンス配列が存在しない場合、全ゲノムの配列データ（FASTA）とアノテーションデータ（GFF/GTF）があれば、Cufflinks 中の gffread コマンドで転写産物のリファレンス配列（FASTA）作成できる。
というのだけれど、もしこれが単に配列データから切出してくるだけなら、mRNAがスプライシングやらエディティングやらを受けるので、本物ではないはず？？

gffread を使った transcripts fasta 転写物の配列取得 – バイオインフォ道場 [bioinfo-Dojo]に曰く
gffread を使った転写物配列の切り出しゲノム配列（genome.fa）と位置情報（GFF3 または GTFファイル）を準備して、gffreadを実行します。オプションを豊富に持ち、以下の例では-wを指定してexonを基にした転写物の切り出し（fastaファイルの生成）を行っています。
　　# gffread GFF3を使った転写物（transcript）配列の切り出し
　　$ gffread -w transcripts.fa -g genome.fa transcripts.gff3
ということで、これでいいのだろう。

gffread -w s288c_transcript.fa -g s288c.fna s288c_e.gff

インデックスを作る

kallisto index -i s288c.ix ../s288c_transcript.fa

これでインデックスファイルs288c.ixができたので、これを使って処理。

↑

Quant †

データはpairedなので、対象ファイルを２つ指定。singleの場合はオプションに--singleを指定する。またsingleの場合は-lでリード長を指定する必要がある。

#!/bin/bash
id=(SRR453566 SRR453567 SRR453568 SRR453569 SRR453570 SRR453571)
for item in ${id[@]}
do
  echo start mapping ${item} with Kallisto
  result_dir=${item}_exp_kallisto
  kallisto quant -i s288c.ix -o ${item} -l 101 -s 15 -b 100 ../paired_SRR453569_1.trim.fastq ../paired_SRR453569_2.trim.fastq
  #kallisto quant -i s288c.ix -o ${item} --single -l 101 -s 15 -b 100 ../${item}.fastq
done

結果のファイルは、ディレクトリSRR453566（又は...）の下に、abundance.h5 abundance.tsv, run_info.json。abundanceファイルの中にtpmが入っている。

target_id       length  eff_length      est_counts      tpm
rna0    363     263     1.08821 0.921395
rna1    228     128     0       0
rna2    1782    1682    0       0
rna3    387     287     0       0
rna4    381     281     0       0
rna5    381     281     0       0
rna6    285     185     5       6.01846
rna7    291     191     8       9.32704
rna8    3969    3869    78.6869 4.52888
rna9    1374    1274    732.109 127.966
rna10   1254    1154    589     113.657
rna11   1149    1049    967     205.276
rna12   639     539     123     50.8164
rna13   1509    1409    173     27.3415
rna14   2643    2543    396     34.6766
rna15   543     443     83      41.7217

RNA idだけなので、そこから対応するgeneを拾う必要あるか？

kallisto (A. thliana, paired-end RNA-Seq) | kallisto を用いた A. thaliana paired-end リードの転写産物の定量

Kallistoを用いたRNA-seq解析パイプライン – 大阪大学医学部 Python会

のいずれも、次ステップはTximport, DESeq2/edgeRと書いている。

tximport | RSEM/kallisto/Salmon の発現量データを edgeR/DESeq2 などに橋渡しする R パッケージに曰く
遺伝子発現量 kallisto が出力した発現量データは転写産物レベルでの発現量である。遺伝子レベルでの解析を行う場合は、これら転写産物レベルでの発現量を遺伝子レベルの発現量に換算する必要がある。この換算は tximport パッケージの summarizeToGene 関数で実行できる。

次のように summarizeToGene を実行するときに countsFromAbundance を指定してカウントデータを取得する必要がある。countsFromAbundance に指定できるスケーリング方法は scaledTPM と lengthScaledTPM の 2 種類があるが、どちらを指定してもいい。

↑