データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ（PRJNA153493）をサンプルとして、kallisto で TAIR10 のゲノム配列へマッピングする。このデータセットは 8 サンプルを含み、4 サンプルが DEX-treated 35S で、他の 4 サンプルは mock-treated 35S である（Huang et al, 2012）。

RNA-Seq データ

データベースにアーカイブされた single-end RNA-Seq （PRJNA153493） のデータを ENA からダウンロードする。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR444/${seqlib}/${seqlib}.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)

for seqlib in ${SEQLIBS[@]}; do
    java -jar trimmomatic-0.38.jar SE -phred33 -threads 4  \
    ${seqlib}.fastq.gz ${seqlib}.clean.fastq.gz            \
    ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:60
done

ゲノム配列データ

次に Ensembl Plants から A. thaliana のゲノム配列（dna.toplevel）をダウンロードする。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz

マッピングには必要としないが、マッピングされたリードを遺伝子領域ごとに集計するときにアノテーションファイルも必要である。このアノテーションファイルも（GFF または GTF）も合わせてダウンロードしておく。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.40.gff3.gz
gunzip Arabidopsis_thaliana.TAIR10.40.gff3.gz

• Ensembl にはゲノムの全配列データの他に、cDNA（coding DNA）や CDS（cDNA から UTR などを取り除いたデータ）などの配列データも提供されている。ここでは、ゲノムの全配列をダウンロードしたいので、DNA (.toplevel.fa.gz) をダウンロードしてくる。
• ゲノムの配列データとして *_sm.toplevel.fa.gz も配布されているが、これを利用しないので、ダウンロードする必要がない。
• Ensembl にゲノム配列に対応したアノテーションデータ（GFF/GTF）も配布されている。このアノテーションデータは .toplevel.fa.gz の配列に対応している。

kallisto による転写産物の定量

RNA-Seq データ（*.fastq.gz）と転写産物の配列データ（Arabidopsis_thaliana.TAIR10.cdna.all.fa）を揃えたら、次に転写産物の配列データのインデックスを作成する。インデックスの名前は自由に付けられるが、ここでは TAIR10_kallisto_index と名付けている。この処理は数分で終わる。

kallisto index -i TAIR10_kallisto_index Arabidopsis_thaliana.TAIR10.cdna.all.fa
## [build] loading fasta file Arabidopsis_thaliana.TAIR10.cdna.all.fa
## [build] k-mer length: 31
## [build] warning: clipped off poly-A tail (longer than 10)
##         from 14 target sequences
## [build] warning: replaced 11 non-ACGUT characters in the input sequence
##         with pseudorandom nucleotides
## [build] counting k-mers ... done.
## [build] building target de Bruijn graph ...   done
## [build] creating equivalence classes ...  done
## [build] target de Bruijn graph has 173469 contigs and contains 45567218 k-mers

ls -l
## total 1047980
## -rw-r--r-- 1 ____ ________  99415585 Aug 18 07:45 Arabidopsis_thaliana.TAIR10.cdna.all.fa
## -rw-r--r-- 1 ____ ________ 973695034 Aug 18 08:27 TAIR10_kallisto_index

次に RNA-Seq データと転写産物のインデックスの照合を行い、転写産物の定量を行う。-i に上で作成したインデックスを指定する。-o には kallisto の出力結果の保存先を指定する。また、定量後に sleuth による発現変動遺伝子の検出を行う予定がある場合、-b オプションをつけて、ブートストラップ法による発現量の不確実性（「分散・ばらつき」のようなもの）を推定する。

single-end リードの場合は、シーケンシング用のライブラリー中のフラグメントの長さの平均値（-l）と偏差（-s）を与える必要がある。Illumina のライブラリーの場合は 180-200 bp と言われている。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    result_dir=${seqlib}_exp_kallisto
    kallisto quant -i TAIR10_kallisto_index -o ${result_dir} --single -l 200 -s 20 -b 100 ${seqlib}.clean.fastq.gz
done

kallisto はマッピングを伴わないので、BAM/SAM ファイルは生成されない。その代わり、推定された発現量はテキストファイル abundances.tsv に保存される。また、-b によるブートストラップの結果は abundances.h5 ファイルに保存される。abundances.h5 は HDF5 フォーマットで、sleuth に発現変動遺伝子検出時の入力するファイルである。また、kallisto を実行した時のオプションなどは run_info.json に保存される。

ls
## drwxr-xr-x 2 _____ ________  4096 Aug 25 13:46 SRR444595_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 13:56 SRR444596_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:05 SRR444597_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:13 SRR444598_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:21 SRR444599_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:32 SRR444600_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:42 SRR444601_exp_kallisto
## drwxr-xr-x 2 _____ ________  4096 Aug 25 14:52 SRR444602_exp_kallisto

ls SRR444595_exp_kallisto
## abundance.h5  abundance.tsv  run_info.json

cat SRR444595_exp_kallisto/run_info.json
## {
## 	"n_targets": 48359,
## 	"n_bootstraps": 100,
## 	"n_processed": 38088331,
## 	"n_pseudoaligned": 31834811,
## 	"n_unique": 15740131,
## 	"p_pseudoaligned": 83.6,
## 	"p_unique": 41.3,
## 	"kallisto_version": "0.44.0",
## 	"index_version": 10,
## 	"start_time": "Sat Aug 25 13:44:25 2018",
## 	"call": "kallisto quant -i TAIR10_kallisto_index -o SRR444595_exp_kallisto --single -l 100 -s 20 -b 100 SRR444595.clean.fastq.gz"
## }

head SRR444595_exp_kallisto/abundance.tsv
## target_id	length	eff_length	est_counts	tpm
## AT1G19090.1	2509	2410	111	2.0295
## AT1G18320.1	699	600	70.0648	5.14555
## AT5G11100.1	1709	1610	153	4.18745
## AT4G35335.1	1280	1181	1001	37.3481
## AT1G60930.1	3876	3777	92	1.07331
## AT1G17000.1	2351	2252	0	0
## AT3G55270.1	3130	3031	1758	25.5574
## AT4G16150.1	3181	3082	2124.63	30.3763
## AT3G18150.1	1370	1271	8	0.27735

sleuth で発現変動転写産物（遺伝子）を検出したい場合は、R の sleuth パッケージのマニュアルを参照して、kallisto の結果を入力として与えればいい。一方で、edgeR や DESEq2 などを使いたい場合は、R の tximport パッケージを介して kallisto の結果（abundances.tsv または abundances.h5）を edgeR や DESeq2 に入力するのが便利（実行例）。