BWA（Li et al, 2009, Li et al, 2010）はロングリードのマッピングを得意としている。ここで、トウモロコシのアセンブリゲノムの改良の目的でシーケンスされた DNA の 250 bp の paired-end リード（Jiao et al, 2017）を、BWA を使ってトウモロコシ（Zea mays）のゲノムにマッピングしてみる。

DNA-Seq データ

DNA-Seq のデータは PRJNA10769 の番号でアーカイブされている。この番号でアーカイブされているデータは HiSeq 2500 と PacBio RS II の 2 種類がある。今回は、HiSeq 2500 でシーケンシングされた paired-end リードを利用する。

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR296/001/SRR2960981/SRR2960981_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR296/001/SRR2960981/SRR2960981_2.fastq.gz

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、Trimmomatic を使って低クオリティリードの除去や短いリードの除去などを行う。フィルタリング後のファイルの名前を SRR2960981_1.clean.fastq.gz および SRR2960981_2.clean.fastq.gz として保存する。このリードデータは 2 x 80 GB 前後あるので、コンピュータの性能により、フィルタリング処理に数日かかる場合がある。BWA の実行テストであれば、最初の 10 万リードだけにを取り出して処理した方がいい。

java -jar trimmomatic-0.38.jar PE -phred33 -threads 4          \
    SRR2960981_1.fastq.gz SRR2960981_2.fastq.gz                \
    SRR2960981_1.clean.fastq.gz SRR2960981_1.unpaired.fastq.gz \
    SRR2960981_2.clean.fastq.gz SRR2960981_2.unpaired.fastq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:100

## gzip -dc SRR2960981_1.fastq.gz | head -n 400000 > SRR2960981_1.subset.fastq
## gzip -dc SRR2960981_2.fastq.gz | head -n 400000 > SRR2960981_2.subset.fastq
## java -jar trimmomatic-0.38.jar PE -phred33 -threads 4    \
##     SRR2960981_1.subset.fastq SRR2960981_2.subset.fastq  \
##     SRR2960981_1.clean.fastq SRR2960981_1.unpaired.fastq \
##     SRR2960981_2.clean.fastq SRR2960981_2.unpaired.fastq \
##     ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:100

ゲノムの配列データ

トウモロコシ（Zea mays）のゲノムは Ensembl Plants で公開されている。このページから、ゲノム配列のデータをダウンロードする。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/zea_mays/dna/Zea_mays.AGPv4.dna.toplevel.fa.gz
gunzip Zea_mays.AGPv4.dna.toplevel.fa.gz

マッピングには必要としないが、マッピングされたリードを遺伝子領域ごとに集計するときにアノテーションファイルも必要である。このアノテーションファイルも（GFF または GTF）も合わせてダウンロードしておく。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/zea_mays/Zea_mays.AGPv4.40.gff3.gz
unzip Zea_mays.AGPv4.40.gff3.gz

BWA によるマッピング

フィルタリングしたあとのリードデータ（SRR2960981_1.clean.fastq.gz と SRR2960981_2.clean.fastq.gz）とゲノムの配列データ（Zea_mays.AGPv4.dna.toplevel.fa）を準備できたので、次に BWA を使ってリードをゲノム配列にマッピングしていく。

BWA ではマッピングするときに、ゲノム配列から作成したインデックスを使用する。まず、このインデックスを作成する。インデックスの名前は任意につけることができる。ここではインデックスの名前を AGPv4 とした。

bwa index -p AGPv4 Zea_mays.AGPv4.dna.toplevel.fa

インデックスを作成し終えると、AGPv4.amb、AGPv4.ann、AGPv4.bwt、AGPv4.pac、AGPv4.sa の 5 ファイルが生成される。これらのファイルが BWA が使用するインデックスとなる。

次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。BWA によるマッピング結果は SAM ファイルに保存される。このリードデータは 2 x 80 GB 前後あるので、マッピングするのに非常に時間がかかる場合があるので、BWA を試すのが目的であれば、リードデータのサブセットを作成して実行するとよい。

bwa mem -t 4 TAIR10 SRR2960981_1.clean.fastq.gz SRR2960981_2.clean.fastq.gz > SRR2960981.sam

## gzip -dc SRR2960981_1.clean.fastq.gz | head -n 4000000 > SRR2960981_1.clean.sub.fastq
## gzip -dc SRR2960981_2.clean.fastq.gz | head -n 4000000 > SRR2960981_2.clean.sub.fastq
## bwa mem -t 4 TAIR10 SRR2960981_1.clean.sub.fastq SRR2960981_2.clean.sub.fastq > SRR2960981.sub.sam

ここで生成した SAM は、必要に応じて samtools で BAM に変換することができる。

samtools sort -O bam -o SRR2960981.sorted.bam SRR2960981.sam
samtools index SRR2960981.sorted.bam