山内のサイト

RNAseqでStringTieを試す。

インストール

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

サンプルごとの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)

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なので、対象ファイルを2つ指定。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 種類があるが、どちらを指定してもいい。

RSEM/kallisto/Salmon の発現量データを edgeR/DESeq2 などに橋渡しする R パッケージ tximport

RSEM/kallisto/Salmon の発現量データを edgeR/DESeq2 などに橋渡しする R パッケージ tximport

salmonを使う例

RNA-seq 発現変動遺伝子解析:salmonとTCC-GUIを使って - Qiita


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-06-29 (土) 08:52:53 (24d)