macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

RNA seqのワークフロー タキシードプロトコル

 

2012年にnatureでtophatとcufflinksを使ったRNA-seq解析の手法が発表されている。

http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html

tophat、bowtie、cufflinksなどを使ったいわゆるタキシードプロトコルRNA seqの1つのワークフローとして知られており、検索すると多くのプロトコルが見つかる(検索例)。

その1つが以下になる。

上記の手法にのっとり、RNA seq解析を行ってみよう。使えそうなRのプログラムも試してみることにする。

生き物はなんでもよいが、今回は高精度なgtfファイルが手に入るシロイヌナズナRNA seqデータを使うことにする。SRAからアラビドプシスRNA-seqデータを探すと、以下の2群比較のデータが見つかった。これを使ってみる。 

Transcriptome analysis of arabidopsis kamchatica under normal and stress conditions

アラビドプシス属ではあるが、シロイヌナズナではない。Arabidopsis kamchatica(画像)を実験材料に使っている。このシーケンスデータをナズナゲノムに当てて定量していく。

 

シロイヌナズナfastaファイルとgtfファイルはEnsembl genomeの

http://plants.ensembl.org/info/website/ftp/index.html

からダウンロード。

 

nature protocolに載ったワークフローは以下になる。

f:id:kazumaxneo:20170703213533j:plain

 

Prestep、Quality check

 まずデータのクオリティを確認しておく。末端でクオリティが極端に悪くなる部位があればこの段階でクオリティトリムする。

ペアリードのR1を調べる。

fastqc -t 12 -o fastx_test_R1/ -f fastq DRR013381_1.fastq 
  • -t Specifies the number of files which can be processed simultaneously.
  • -o Create all output files in the specified output directory.
  • -f Bypasses the normal sequence file format detection and forces the program to use the specified format. Valid formats are bam,sam,bam_mapped,sam_mapped and fastq
  • -k Specifies the length of Kmer to look for in the Kmer content module. Specified Kmer length must be between 2 and 10. Default length is 7 if not specified.

出力されたhtmlを開く。 

f:id:kazumaxneo:20170630180215j:plain

トリミングは必要なさそうである。マッピングに移ろう。

 

step1、Mapping

cufflinksまでは、biological replicatesがあっても独立したサンプルとして扱っていく。今回はテストなので、chr1だけに限定してマッピングを行う(本当は良くない。アライメントの仕方次第でアーティファクトが出る可能性もゼロではない)。

bowtie2のidnex作成

bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa chr1

最後のchr1はユーザーが決めた名前である。この名前を使いtophatを走らせる。

 

tophat2でpaired-endリードをmapping

tophat2 -p 4 -I 50000 -i 20 --library-type fr-firststrand  
-o tophat_control1 -G Arabidopsis_thaliana.TAIR10.36.gtf chr1
DRR013381_1.fastq DRR013381_2.fastq
  • -p  num-threads  
  • -I  --max-intron-length (default: 500000)
  • -i  --min-intron-length  (default: 50)
  • -G --GTF
  • -o --output-dir (default: ./tophat_out)

(gtfをtophatの時点で指定してmappingする方法もある。その場合-Gをつける)。

アライメントファイルが出力される。今回はリファレンスとデータで種が違うので、アライメントは非常に困難になる。エラーは覚悟の上、--read-mismatches 4 --read-gap-length 4 --read-edit-dist 4も上のtohat2につけてランした(How to control the alignment of reads in terms of number of mismatches, gap length etc. ? を参考。edit distはwikiに説明がある)。ナズナはエコタイプの違いでもたくさん塩基置換が起きている。このパラメータでマッピング感度は上がるが、ランダムマッチも増える。どれくらい有効 or 無謀な方法なのかは分からない。

 

追記;マッピング率は40%ほどだった(default parameterで20%ほど)。

 

 tophatの出力はaccepted_hits.bamという定型名のため、サンプルの区別がつかなくなる恐れがある。出力されるbamをサンプル名が分かる名前にリネームする。

mv tophat_output/accepted_hits.bam tophat_output/control1.bam

 

のちにIGVで可視化するために、bamに indexをつけておく。

samtools index tophat_output/control1.bam

 gtfファイルを指定して解析した場合、次の2-4は飛ばしてstep5に進む。gtfを指定していない場合、cufflinksにアライメント状況からfeatureファイル(gtf)を作ってもらう必要がある。step2-4を実行する。(詳細については、topにリンクを貼ったnature protocolの論文box1. Cを参照)。

 今回はcuffmergeしてアライメント領域を1から再定義した方が賢明と思われるが、時間があまりない。でナズナのgtfをそのまま使う。

 

2、Cufflinks

アライメント部位の情報に基づき、orfを予測してgtfファイルを出力する。

A、 mRNAを再構成して、gftファイルを出力する。

cufflinks --no-update-check --overlap-radius 1 --library-type fr-firststrand -o cufflinks control1.bam
  • –overlap-radius <int> Transfrags that are separated by less than this distance get merged together, and the gap is filled. Default: 50 (in bp).
  • –no-update-check Turns off the automatic routine that contacts the Cufflinks server to check for a more recent version.

B、出力されるgtfファイルをサンプル名がわかる名前にリネームする。

mv cufflinks/transcripts.gtf cufflinks/control1.gtf 

 これをサンプル全てに対して行う。サンプルが多いならシェルスクリプトを書いて自動処理を考える。

 

3-4、Cuffmerge

このstepもgtfファイルをアライメント時に指定している場合は必要ない。cufflinksに読み込んで全サンプルのgtfを統合し、1つのgtfファイルを出力する。最初に読み込むファイルのリストを作る(サンプルが少ないならコマンド実行時に明示してもよい)。

echo cufflinks/control1.gtf > assemblies.txt
echo cufflinks/control2.gtf >> assemblies.txt
echo cufflinks/control3.gtf >> assemblies.txt
echo cufflinks/stress1.gtf >> assemblies.txt
echo cufflinks/stress2.gtf >> assemblies.txt
echo cufflinks/stress3.gtf >> assemblies.txt

 catで中身を確認しておく。

 

gtfのマージを実行する。

cuffmerge -s Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa  assemblies.txt

merged_asm/merged.gtfができる。このmeggeしたgtfファイルを使って以下の計算を進める。

   

5、Cuffdiff

cuffdiffコマンドで サンプル間の発現量差を計算する。

使い方は

cuffdiff [options] <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> [... sampleN_hits.sam]

実際のランは以下のように行った。

cuffdiff -o cuffdiff_output -p 12 --library-type fr-firststrand 
-u Arabidopsis_thaliana.TAIR10.36.gtf -b Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa 
-N --compatible-hits-norm -L control,stress
tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam
tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam
  • -o --output-dir
  • -p --num-threads
  • -b use bias correction - reference fasta required
  • -u use 'rescue method' for multi-reads
  • --library-type Library prep used for input reads 
  • -N --compatible-hits-norm count hits compatible with reference RNAs only 
  • -L comma-separated list of condition labels

-Lでつけた名前を今後使うことになる。区別がつくシンプルな名前が良いn > 2の実験でreplicatesがある場合、replicatesは","でつなぐ。異なるサンプルはスペースでつなぐmanual上のデータはcontrolとstressのtriplicateのデータなので、上のboxの様に記述している。そのほか、末端にある""はコマンドを複数行にわたって記述する時のマークである。

 

出力のcuffdiff_outputにはたくさんのファイルができる。

f:id:kazumaxneo:20170704170316j:plain

発現量の差を調べたファイルの中を見てみる。

gene_exp.diff

f:id:kazumaxneo:20170704150718j:plain

最初の10行だけ表示。 value(FPKMではない)とp_valueが表示されている。これだけでも十分解析できそうである。

FPKMはgenes.fplm.trackingの中にある。

 

 

おまけ

igvでbamを開いて、データを確認してみる。igvがない人は

brew install igv

で入る。

fastaとgtf、bamファイルをコマンドから指定して開く。

igv -g Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa 
Arabidopsis_thaliana.TAIR10.36.gtf,tophat_control1/control1.bam,tophat_control2/control2.bam,tophat_control3/control3.bam,tophat_stress1/stress1.bam,tophat_stress2/stress2.bam,tophat_stress3/stress3.bam

gtfファイル以降は ","でつないで記載する(gtf,bam1,bam2,bam3...)。

gtfがソートされてないと怒られるが、そのまま開く。(gtfのソートについてはバイオインフォ道場さんの記事が参考になります http://bioinfo-dojo.net/2017/04/01/igv_gtf_bed/)。

f:id:kazumaxneo:20170704190723j:plain

bedも選択した場合、gtfの下に追加表示される。

 

pop upが邪魔なら上のボタンからnever showをチェック。

f:id:kazumaxneo:20170704190809j:plain

 

 

6-18、CummeRbund

Rにデータを読み込む。ターミナルでRとタイプ

uesaka$ R

 

R version 3.3.0 (2016-05-03) -- "Supposedly Educational"

Copyright (C) 2016 The R Foundation for Statistical Computing

Platform: x86_64-apple-darwin13.4.0 (64-bit)

 

R は、自由なソフトウェアであり、「完全に無保証」です。 

一定の条件に従えば、自由にこれを再配布することができます。 

配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 

 

R は多くの貢献者による共同プロジェクトです。 

詳しくは 'contributors()' と入力してください。 

また、R や R のパッケージを出版物で引用する際の形式については 

'citation()' と入力してください。 

 

'demo()' と入力すればデモをみることができます。 

'help()' とすればオンラインヘルプが出ます。 

'help.start()' で HTML ブラウザによるヘルプがみられます。 

'q()' と入力すれば R を終了します。 

 

ターミナル内でRを起動する。3.3が起動している。

 

 

library("cummeRbund")を打ってcummeRbundを読み込む。cummeRbundを入れてない人はBiocunducterのパッケージをインストールする。

source("http://bioconductor.org/biocLite.R") #時間がかかります 
biocLite("cummeRbund")
biocLite("org.At.tair.db") #GO解析用
library("cummeRbund")
library("org.At.tair.db")

 cuffdiff出力ディレクトリを読み込む。

cuff = readCufflinks('diff_output')

 

cummeRbundのグラフコマンドを一通り実行してみる。

csDensity(genes(cuff))

 

csScatter(genes(cuff), 'stress1', 'control1') 

 

csScatterMatrix(genes(cuff))

 

csVolcanoMatrix(genes(cuff))

 

 

 

 cufflinks2.2.1にバグがあるのか、cummeRbundで正しく読み込めない。というかsqliteのデータのファイルサイズが妙に小さいので、正しくデータベースが構築できていない気がする。以前cuffldiffで解析したデータは正常にcummeRbundに読み込めたため、cummeRbundは問題ない可能性が高い。正しくcummeRbundに読み込め次第追記します。

 

cufflinksを使わず、featurecountsでリードをカウントし、edgeRで正規化、検定してDEGを出す流れは以下の2記事を参考にしてください。

 

リードのカウント

edgeRで検定


引用

 

 

 

 

 

 

 

 

 

メモ

補正前のリードカウントしただけのrawデータとFPKM補正値を比較。

featureCountsでraw count dataを計算(参考)。

featureCounts -T 12 -p -t exon -g gene_id -a Arabidopsis_thaliana.TAIR10.36.gtf -o featureCounts_counts.txt 
tophat_control1/control1.bam
tophat_control2/control2 bam
tophat_control3/control3.bam
tophat_stress1/stress1.bam
tophat_stress2/stress2.bam
tophat_stress3/stress3.bam

 

cufflinksでFPKM補正値を算出。

cufflinks -o output_FPKM -G Arabidopsis_thaliana.TAIR10.36.gtf 
-N --compatible-hits-norm --library-type fr-firststrand
tophat_control1/control1.bam 
tophat_control2/control2.bam 
tophat_control3/control3.bam
tophat_stress1/stress1.bam 
tophat_stress2/stress2.bam 
tophat_stress3/stress3.bam 
  • -G 既存のアイソフォームのみの計算の指示
  • -N --compatible-hits-norm 既存のアイソフォームにマッピングされるものだけでFPKMを計算
  •  -o FPKMの出力ディレクトリ名

 

 

featureCountsの出力は余計な行があるので削っておく。

cut -f 1,7-12 featureCounts_counts.txt > featureCounts_counts_trimmed.txt

 featureCounts_counts_trimmed.txtはこのような内容になっている。リードをカウントしているだけなので整数である。

Geneid control1 control2 control3 stress1 stress2 stress3
AT1G01010 0 0 0 0 0 0
AT1G01020 19 31 26 24 31 17
AT1G03987 0 0 1 0 0 0
AT1G01030 7 20 10 11 20 4
AT1G01040 44 77 44 91 77 86
AT1G03993 0 0 0 0 0 0
AT1G01046 0 0 0 0 0 0
AT1G01050 265 472 339 222 472 226
AT1G03997 0 0 0 0 0 0

 

 Rを起動してデータフレームに featureCounts_counts_trimmed.txtを読み込む。

>R
library("ggplot2") #ggplot2読み込み
library("reshape2") #reshape読み込み

rawCount = read.table("featureCounts_counts_trimmed.txt", header = TRUE, row.names = 1)

 

 ggplot2でヒストグラムを書く。

artificialCount = log2(rawCount + 1) #log2の前に0を消すため1を強制的に足す。
ggplot(artificialCount, aes(x = control1)) + ylab(expression(log[2](count + 1))) + geom_histogram(colour = "white", fill = "#525252", binwidth = 0.8)

fillでバーの色を指定(RGB表記)。

control1

f:id:kazumaxneo:20170715125723j:plain

今回は15GBx2くらいシーケンスしてるデータ使っている。リード数がゼロのORFが多いのは、シーケンス不足が理由ではないと推測される。アライメントができていないか、gtfの定義が正しくないorfが多い可能性がある。

残りの条件もcontrol1のところを変えてグラフ化する。

 

 

念のため、箱ひげ図も書く。

data <- read.table("featureCounts_counts_trimmed.txt", header = TRUE) #読み込み
logdata = log2(data + 1)
head(logdata) #確認
df <- melt(logdata) #ggplot2が認識できるdata.frame形式にmeltコマンドで変換。
head(df)

実行結果はこうなる。

> head(df)

     Geneid                          variable value

1 AT1G01010 tophat_control1.accepted_hits.bam     0

2 AT1G01020 tophat_control1.accepted_hits.bam    19

3 AT1G03987 tophat_control1.accepted_hits.bam     0

4 AT1G01030 tophat_control1.accepted_hits.bam     7

5 AT1G01040 tophat_control1.accepted_hits.bam    44

6 AT1G03993 tophat_control1.accepted_hits.bam     0

グラフを書く。

g <- ggplot(df,  aes (x = variable, y = value )) + geom_boxplot()  
plot(g)

 

raw reads

f:id:kazumaxneo:20170715131106j:plain

 

 FPKM補正

 

 

 

 

 TMM補正f:id:kazumaxneo:20170715130615j:plain