データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ(PRJNA153493)をサンプルとして、Bowtie2 で TAIR10 のゲノム配列へマッピングする。このデータセットは 8 サンプルを含み、4 サンプルが DEX-treated 35S で、他の 4 サンプルは mock-treated 35S である(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 の配列に対応している。
TopHat によるマッピング
RNA-Seq データ(*.fastq.gz)とゲノムの配列データ(Arabidopsis_thaliana.TAIR10.dna.toplevel.fa)を揃えたら、次にゲノムの配列データのインデックスを作成する。TopHat2 では Bowtie2 のインデックスを利用する。ここで Bowtie2 のコマンド bowtie2-build
を使用してインデックスを作成する。インデックスの名前は自由に付けられるが、ここでは TAIR10 と名付けている。ゲノムが大きいとき、1 時間近くかかる場合がある。
bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna.toplevel.fa TAIR10
インデックスを作成し終えると、そのフォルダには TAIR10.1.bt2、TAIR10.2.bt2、TAIR10.3.bt2、TAIR10.4.bt2、TAIR10.rev.1.bt2、TAIR10.rev.2.bt2 の 6 つのファイルが生成される。これらのファイルのプレフィックスが、インデックスの名前に指定した TAIR10 となっている。また、ヒトゲノムなどのように、ゲノムが大きい場合、インデックスの拡張子は .bt2 ではなく .bt2l になる場合がある。
次に TopHat2 を使って RNA-Seq データをゲノム上にマッピングする。サンプルが 8 個あるので、ここでは bash の for
を使って、各サンプルに対して Bowtie2 マッピングを行う。
SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)
for seqlib in ${SEQLIBS[@]}; do
tophat2 -p 4 -o ${seqlib}_tophat2 TAIR10 ${seqlib}.clean.fastq.gz
done
TopHat2 によるマッピング結果は -o
で指定したフォルダ(例:SRR7508939_tophat2 など)に保存される。accepted_hits.bam がリードをゲノム配列にマッピングした結果である。マッピング率などは align_summary.txt に保存されている。
ls SRR444595_tophat2
## accepted_hits.bam align_summary.txt deletions.bed insertions.bed junctions.bed logs prep_reads.info unmapped.bam
cat SRR444595_tophat2/align_summary.txt
## Reads:
## Input : 38088331
## Mapped : 33200505 (87.2% of input)
## of these: 1377354 ( 4.1%) have multiple alignments (0 have >20)
## 87.2% overall read mapping rate.
ここで得られた BAM ファイル(accepted_hits.bam)に変換しソートした後に、 htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。SAM ファイルをソートしながら BAM ファイルに変換するには、samtools を利用する。
SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)
for seqlib in ${SEQLIBS[@]}; do
samtools sort -O bam -o ${seqlib}.sorted.bam ${seqlib}.bam
samtools index ${seqlib}.sorted.bam
done
References
- RBE controls microRNA164 expression to effect floral organogenesis. Development. 2012, 139(12):2161-9. DOI: 10.1242/dev.075069