データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ（PRJNA153493）をサンプルとして、Salmon で 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 の配列に対応している。

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 -r ${seqlib}.clean.fastq.gz -o ${result_dir}
done

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

ls -l
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:53 SRR444595_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:54 SRR444596_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:55 SRR444597_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:55 SRR444598_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:56 SRR444599_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:56 SRR444600_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:57 SRR444601_exp_salmon
## drwxr-xr-x 5 _____ _______  4096 Aug 25 13:57 SRR444602_exp_salmon

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

head SRR444595_exp_salmon/quant.sf
## Name	Length	EffectiveLength	TPM	NumReads
## AT1G19090.1	2509	2260.000	1.644046	111.000
## AT1G18320.1	699	450.000	7.260776	97.610
## AT5G11100.1	1709	1460.000	3.553679	155.000
## AT4G35335.1	1280	1031.000	32.661689	1006.000
## AT1G60930.1	3876	3627.000	0.849062	92.000
## AT1G17000.1	2351	2102.000	0.000000	0.000
## AT3G55270.1	3130	2881.000	20.402368	1756.000
## AT4G16150.1	3181	2932.000	24.530087	2148.640
## AT3G18150.1	1370	1121.000	0.268742	9.000

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