データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ（PRJNA153493）をサンプルとして、STAR（Dobin et al, 2013）で 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 の配列に対応している。
STAR によるマッピング
RNA-Seq データ（*.fastq.gz）とゲノム配列データ（Arabidopsis_thaliana.TAIR10.dna.toplevel.fa）を揃えたら、次にゲノムの配列データのインデックスを作成する。インデックスの名前は自由に付けられるが、ここでは TAIR10_STAR_index と名付けている。ゲノムが大きいとき、1 時間近くかかる場合がある。また、ゲノムが大きいとき、より大きなメモリを確保できないとメモリ不足でエラーになる。
mkdir TAIR10_STAR_index
STAR --runMode genomeGenerate --genomeDir TAIR10_STAR_index --genomeFastaFiles Arabidopsis_thaliana.TAIR10.dna.toplevel.fa --limitGenomeGenerateRAM 3400000000
## Aug 18 22:42:50 ..... started STAR run
## Aug 18 22:42:50 ... starting to generate Genome files
## Aug 18 22:42:53 ... starting to sort Suffix Array. This may take a long time...
## Aug 18 22:42:54 ... sorting Suffix Array chunks and saving them to disk...
## Aug 18 22:45:29 ... loading chunks from disk, packing SA...
## Aug 18 22:45:32 ... finished generating suffix array
## Aug 18 22:45:32 ... generating Suffix Array index
## Aug 18 22:46:24 ... completed Suffix Array index
## Aug 18 22:46:24 ... writing Genome to disk ...
## Aug 18 22:46:24 ... writing Suffix Array to disk ...
## Aug 18 22:46:26 ... writing SAindex to disk
## Aug 18 22:46:28 ..... finished successfully
## Aug 18 22:46:28 ..... started STAR run
## Aug 18 22:46:28 ... starting to generate Genome files次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。マッピングもメモリを大量に必要とするので、STAR を実行するときにメモリを確保しておく。
SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)
for seqlib in ${SEQLIBS[@]}; do
result_dir=${seqlib}_exp_star
STAR --runThreadN 4 --genomeDir TAIR10_STAR_index --readFilesIn ${seqlib}.clean.fastq.gz --readFilesCommand zcat --genomeLoad NoSharedMemory --outFilterMultimapNmax 1 --outFileNamePrefix ${result_dir}
doneマッピング結果は、あらかじめ設定したプレフィックス（result_dir）のついたファイルに保存される。上のコマンドを実行した場合 SRR7508939_exp_star から始まるファイルが生成される。
ls SRR444595_exp_star*
## SRR444595_exp_starAligned.out.sam SRR444595_exp_starLog.final.out SRR444595_exp_starLog.out SRR444595_exp_starLog.progress.out SSRR444595_exp_starSJ.out.tab*.sam はマッピング結果が保存されている SAM ファイルである。*.final.out は、マッピング率などのメタ情報が書かれている。
samtools flagstat SRR444595_exp_starAligned.out.sam
## 35767352 + 0 in total (QC-passed reads + QC-failed reads)
## 0 + 0 secondary
## 0 + 0 supplementary
## 0 + 0 duplicates
## 35767352 + 0 mapped (100.00% : N/A)
## 0 + 0 paired in sequencing
## 0 + 0 read1
## 0 + 0 read2
## 0 + 0 properly paired (N/A : N/A)
## 0 + 0 with itself and mate mapped
## 0 + 0 singletons (N/A : N/A)
## 0 + 0 with mate mapped to a different chr
## 0 + 0 with mate mapped to a different chr (mapQ>=5)
cat SRR444595_exp_starLog.final.out
## Started job on | Aug 20 07:52:28
## Started mapping on | Aug 20 07:52:33
## Finished on | Aug 20 08:11:29
## Mapping speed, Million of reads per hour | 120.70
##
## Number of input reads | 38088331
## Average input read length | 35
## UNIQUE READS:
## Uniquely mapped reads number | 35767352
## Uniquely mapped reads % | 93.91%
## Average mapped length | 35.46
## Number of splices: Total | 2145703
## Number of splices: Annotated (sjdb) | 0
## Number of splices: GT/AG | 2121416
## Number of splices: GC/AG | 17617
## Number of splices: AT/AC | 828
## Number of splices: Non-canonical | 5842
## Mismatch rate per base, % | 0.42%
## Deletion rate per base | 0.00%
## Deletion average length | 1.42
## Insertion rate per base | 0.00%
## Insertion average length | 1.25
## MULTI-MAPPING READS:
## Number of reads mapped to multiple loci | 0
## % of reads mapped to multiple loci | 0.00%
## Number of reads mapped to too many loci | 1686440
## % of reads mapped to too many loci | 4.43%
## UNMAPPED READS:
## % of reads unmapped: too many mismatches | 0.00%
## % of reads unmapped: too short | 1.66%
## % of reads unmapped: other | 0.01%
## CHIMERIC READS:
## Number of chimeric reads | 0
## % of chimeric reads | 0.00%STAR のマッピング結果から uniquely mapped reads を取得するには、次のプションをつけて samtools 実行すれば良い。（paired-end の場合は、ペア情報を考慮する必要があるため、異なるオプションを用いる。）
samtools -q 255 SRR444595_exp_starAligned.out.sam | grep -v "@SQ" | wc -lReferences
- STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013, 29(1):15-21. DOI: 10.1093/bioinformatics/bts635
- RBE controls microRNA164 expression to effect floral organogenesis. Development. 2012, 139(12):2161-9. DOI: 10.1242/dev.075069