データベースにアーカイブされた A. thaliana の paired-end RNA-Seq データをサンプルとして、kallisto（Bray et al, 2016）による転写産物の定量を行う。ここで使用するデータは PRJNA480638 でアーカイブされている。nsra nsrb double mutant と対照群それぞれ 3 サンプルからなる計 6 サンプルからなる。

RNA-Seq データ

データベースにアーカイブされた PRJNA480638 のデータを利用する。このデータは A. thaliana のデータで、nsra nsrb double mutant と対照群それぞれ 3 サンプルからなる計 6 サンプルからなる。サンプルは paired-end でシーケンシングされているので、リード情報が 1 サンプルに付き forward と reversed の 2 つのファイルに保存されている。この 6 サンプル 12 ファイルをローカルにダウンロードする。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR750/00${seqlib:9:10}/${seqlib}/${seqlib}_1.fastq.gz
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR750/00${seqlib:9:10}/${seqlib}/${seqlib}_2.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。フィルタリング後のファイルを *_1.clean.fastq.gz および *_2.clean.fastq.gz のような名前で保存する。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

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

転写産物の配列データ

次に Ensembl Plants から A. thaliana の転写産物（cDNA）のリファレンス配列をダウンロードする。ダウンロードしたデータに含まれている転写産物を調べてみると、48,359 個である。

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

grep ">" Arabidopsis_thaliana.TAIR10.cdna.all.fa | cut -d' ' -f 1 | wc
##  48359   48359  628959

転写産物（cDNA）のリファレンス配列が存在しない場合、全ゲノムの配列データ（FASTA）とアノテーションデータ（GFF/GTF）があれば、Cufflinks 中の gffread コマンドで転写産物のリファレンス配列（FASTA）作成できる。

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 オプションをつけて、ブートストラップ法による発現量の不確実性（「分散・ばらつき」のようなもの）を推定する。

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} -b 100 ${seqlib}_1.clean.fastq.gz ${seqlib}_2.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 18 12:04 SRR7508939_exp_kallisto
## drwxr-xr-x 2 _____ ________       4096 Aug 18 12:14 SRR7508940_exp_kallisto
## drwxr-xr-x 2 _____ ________       4096 Aug 18 12:25 SRR7508941_exp_kallisto
## drwxr-xr-x 2 _____ ________       4096 Aug 18 12:35 SRR7508942_exp_kallisto
## drwxr-xr-x 2 _____ ________       4096 Aug 18 12:46 SRR7508943_exp_kallisto
## drwxr-xr-x 2 _____ ________       4096 Aug 18 12:56 SRR7508944_exp_kallisto

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

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