STAR（Bray et al, 2016）を使って、A. thaliana の paired-end RNA-Seq データを、ゲノム配列にマッピングしてみる。ここで使用するデータは 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 のゲノム配列（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.gzSTAR によるマッピング
RNA-Seq データ（*.fastq.gz）とゲノムの配列データ（Arabidopsis_thaliana.TAIR10.cdna.all.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=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)
for seqlib in ${SEQLIBS[@]}; do
result_dir=${seqlib}_exp_star
STAR --runThreadN 4 --genomeDir TAIR10_STAR_index --readFilesIn ${seqlib}_1.clean.fastq.gz ${seqlib}_2.clean.fastq.gz --readFilesCommand zcat --genomeLoad NoSharedMemory --outFilterMultimapNmax 1 --outFileNamePrefix ${result_dir}
doneマッピング結果は、あらかじめ設定したプレフィックス（result_dir）のついたファイルに保存される。上のコマンドを実行した場合 SRR7508939_exp_star から始まるファイルが生成される。
ls SRR7508939_exp_star*
## SRR7508939_exp_starAligned.out.sam SRR7508939_exp_starLog.final.out SRR7508939_exp_starLog.out SRR7508939_exp_starLog.progress.out SRR7508939_exp_starSJ.out.tab*.sam はマッピング結果が保存されている SAM ファイルである。*.final.out は、マッピング率などのメタ情報が書かれている。このデータ（SRR7508939）の Uniquely mapped reads が全体の 97.6% である。
samtools flagstat SRR7508939_exp_starAligned.out.sam
## 51783768 + 0 in total (QC-passed reads + QC-failed reads)
## 0 + 0 secondary
## 0 + 0 supplementary
## 0 + 0 duplicates
## 51783768 + 0 mapped (100.00% : N/A)
## 51783768 + 0 paired in sequencing
## 25891884 + 0 read1
## 25891884 + 0 read2
## 51783768 + 0 properly paired (100.00% : N/A)
## 51783768 + 0 with itself and mate mapped
## 0 + 0 singletons (0.00% : N/A)
## 0 + 0 with mate mapped to a different chr
## 0 + 0 with mate mapped to a different chr (mapQ>=5)
cat SRR7508939_exp_starLog.final.out
## Started job on | Aug 18 23:15:13
## Started mapping on | Aug 18 23:15:16
## Finished on | Aug 18 23:32:01
## Mapping speed, Million of reads per hour | 95.02
##
## Number of input reads | 26526105
## Average input read length | 199
## UNIQUE READS:
## Uniquely mapped reads number | 25891884
## Uniquely mapped reads % | 97.61%
## Average mapped length | 198.87
## Number of splices: Total | 14810474
## Number of splices: Annotated (sjdb) | 0
## Number of splices: GT/AG | 14670193
## Number of splices: GC/AG | 131385
## Number of splices: AT/AC | 4174
## Number of splices: Non-canonical | 4722
## Mismatch rate per base, % | 0.06%
## Deletion rate per base | 0.00%
## Deletion average length | 1.22
## Insertion rate per base | 0.00%
## Insertion average length | 1.06
## 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 | 285251
## % of reads mapped to too many loci | 1.08%
## UNMAPPED READS:
## % of reads unmapped: too many mismatches | 0.00%
## % of reads unmapped: too short | 1.32%
## % of reads unmapped: other | 0.00%
## CHIMERIC READS:
## Number of chimeric reads | 0
## % of chimeric reads | 0.00%ここで得られた SAM ファイルを htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。
References
- STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013, 29(1):15-21. DOI: 10.1093/bioinformatics/bts635