データベースにアーカイブされた A. thaliana の paired-end RNA-Seq データをサンプルとして、Salmon（Patro et al, 2017）による転写産物の定量を行う。ここで、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）作成できる。

Salmon による転写産物の定量

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

salmon index -t Arabidopsis_thaliana.TAIR10.cdna.all.fa -i TAIR10_salmon_index
## Version Info: This is the most recent version of Salmon.
## index ["TAIR10_salmon_index"] did not previously exist  . . . creating it
## [2018-08-18 10:52:57.633] [jLog] [info] building index
## RapMap Indexer
## 
## [Step 1 of 4] : counting k-mers
## [2018-08-18 10:52:57.912] [jointLog] [warning] Entry with header [AT5G40315.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:52:58.182] [jointLog] [warning] Entry with header [AT1G33355.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:52:58.458] [jointLog] [warning] Entry with header [ATMG00665.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:52:59.602] [jointLog] [warning] Entry with header [AT5G23115.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:52:59.720] [jointLog] [warning] Entry with header [AT3G23122.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:52:59.918] [jointLog] [warning] Entry with header [AT4G09355.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:53:00.060] [jointLog] [warning] Entry with header [AT5G24575.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:53:00.435] [jointLog] [warning] Entry with header [AT1G64633.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:53:00.574] [jointLog] [warning] Entry with header [AT4G33735.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## [2018-08-18 10:53:00.622] [jointLog] [warning] Entry with header [AT2G21105.1], had length less than the k-mer length of 31 (perhaps after poly-A clipping)
## Elapsed time: 3.01856s
## 
## [2018-08-18 10:53:00.654] [jointLog] [warning] Removed 91 transcripts that were sequence duplicates of indexed transcripts.
## [2018-08-18 10:53:00.654] [jointLog] [warning] If you wish to retain duplicate transcripts, please use the `--keepDuplicates` flag
## Replaced 11 non-ATCG nucleotides
## Clipped poly-A tails from 14 transcripts
## Building rank-select dictionary and saving to disk done
## Elapsed time: 0.0131726s
## Writing sequence data to file . . . done
## Elapsed time: 0.126709s
## [info] Building 32-bit suffix array (length of generalized text is 86407032)
## Building suffix array . . . success
## saving to disk . . . done
## Elapsed time: 5.83075s
## done
## Elapsed time: 16.1958s
## processed 86000000 positions
## khash had 46155873 keys
## saving hash to disk . . . done
## Elapsed time: 15.1383s
## [2018-08-18 10:54:17.623] [jLog] [info] done building index

ls -l
## total 1047984
## -rw-r--r-- 1 _____ ________  99415585 Aug 18 07:45 Arabidopsis_thaliana.TAIR10.cdna.all.fa
## drwxr-xr-x 2 _____ ________      4096 Aug 18 10:54 TAIR10_salmon_index

ls TAIR10_salmon_index
## duplicate_clusters.tsv  hash.bin  header.json  indexing.log  quasi_index.log  refInfo.json  rsd.bin  sa.bin  txpInfo.bin  versionInfo.json

次に RNA-Seq リードを、インデックスを介して、転写産物の配列上に擬似的にマッピングする。-i に上で作成したインデックスを指定する。-o には Salmon の出力結果の保存先を指定する。-l A （小文字エル）はリードの方向性（strand）を自動的に予測するオプションである。また、-1 と -2 にはそれぞれ paired-end の FASTQ ファイルを指定する。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    result_dir=${seqlib}_exp_salmon
    salmon quant -i TAIR10_salmon_index -l A -1 ${seqlib}_1.clean.fastq.gz -2 ${seqlib}_2.clean.fastq.gz -o ${result_dir}
done

Salmon による定量結果は -o で指定したフォルダに保存される。発現量は Salmon の出力に指定したフォルダの中に quant.sf ファイルに保存されている。

ls -l
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:04 SRR7508939_exp_salmon
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:05 SRR7508940_exp_salmon
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:06 SRR7508941_exp_salmon
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:07 SRR7508942_exp_salmon
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:09 SRR7508943_exp_salmon
## drwxr-xr-x 5 _____ ________       4096 Aug 18 12:10 SRR7508944_exp_salmon

ls SRR7508939_exp_salmon
## aux_info  cmd_info.json  lib_format_counts.json  libParams  logs  quant.sf

head SRR7508939_exp_salmon/quant.sf
## Name	Length	EffectiveLength	TPM	NumReads
## AT1G19090.1	2509	2335.730	0.000000	0.000
## AT1G18320.1	699	525.749	1.159374	12.232
## AT5G11100.1	1709	1535.730	4.997076	154.000
## AT4G35335.1	1280	1106.730	20.487071	455.000
## AT1G60930.1	3876	3702.730	2.126401	158.000
## AT1G17000.1	2351	2177.730	0.000000	0.000
## AT3G55270.1	3130	2956.730	25.853762	1534.000
## AT4G16150.1	3181	3007.730	39.663895	2394.000
## AT3G18150.1	1370	1196.730	0.000000	0.000

ここで定量した転写産物の発現量を利用して、発現変動転写産物（遺伝子）を検出したい場合、Salmon の出力結果を、R の tximport パッケージを介して、edgeR や DESeq2 に与えればいい（実行例）。