STAR による single-end RNA-Seq リードの高速マッピング

STAR (A. thaliana, single-end RNA-Seq)

データベースにアーカイブされた 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%

ここで得られた SAM ファイルを htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。

References

  • Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013, 29(1):15-21. DOI: 10.1093/bioinformatics/bts635
  • Huang T, López-Giráldez F, Townsend JP, Irish VF. RBE controls microRNA164 expression to effect floral organogenesis. Development. 2012, 139(12):2161-9. DOI: 10.1242/dev.075069