RSEM は、内部で Bowtie、Bowtie2 および STAR を呼び出してマッピングを行い、マッピングを解釈して転写産物の発現量推定を行うことができる。ここで、RSEM を使って、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

RSEM では、転写産物の発現量推定を行うために、遺伝子と転写産物の対応関係を知る必要がある。そのため、ゲノム配列データに加えて、アノテーションファイルも（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

マッピング

RSEM はマッピング結果を利用して発現量推定を行う。そのため、まずリードをリファレンス配列にマッピングすることから始める。マッピングは Bowtie2、HISAT2 や STAR などのプログラムで行う。ここでは、STAR を用いた例を示す。

リファレンス配列に対して、STAR 用のインデックスを作成する。STAR は多くのメモリー容量を必要とするので、大きいゲノムを扱うときに十分にメモリ容量を確保しておく必要がある。

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

続けて、リードをインデックスにマッピングする。このとき、あとで RSEM で解析できるように STAR に --quantMode TranscriptomeSAM オプションをつける。

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 \
         --quantMode TranscriptomeSAM \
         --outSAMtype BAM SortedByCoordinate \
         --genomeLoad NoSharedMemory \
         --outFilterMultimapNmax 1 \
         --outFileNamePrefix ${result_dir}
done

STAR が正しく実行されると、SRR7508939_exp_star などの prefix のついたファイルが生成され、その中で SRR7508939_exp_star_Aligned.sortedByCoord.out.bam と SRR7508939_exp_star_Aligned.toTranscriptome.out.bam がマッピング結果である。後者の方は RSEM で使用する。

RSEM による転写産物の発現量定量

RSEM を利用するにあたって、リファンレンス配列のインデックスを作成する。インデックスを作成する際に、遺伝子と転写産物の対応関係を必要とするので、ゲノム配列の他にアノテーションも指定する。アノテーションが GFF ならば --gff3、GTF ならば --gtf オプションを利用する。

rsem-prepare-reference --gff3 Arabidopsis_thaliana.TAIR10.40.gff3 -p 4 \
         Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
         TAIR10_RSEM_index

次に toTranscriptome BAM ファイルを RSEM に処理させて発現量の定量を行う。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    rsem-calculate-expression --alignments --paired-end -p 4 \
             ${seqlib}_exp_star_Aligned.toTranscriptome.out.bam \
             TAIR10_RSEM_index \
             ${seqlib}_rsem
done

定量プロセスが正しく終了すると、SRR7508939_rsem.genes.results、SRR7508939_rsem.isoforms.results および SRR7508939_rsem.transcript.bam などのファイルが生成される。genes.results および isoforms.results ファイルにはリードカウント、TPM、および FPKM などの発現量データが保存されている。これらのファイルを R に読み込んで、続きの解析を行う。