BWA (Zea mays, paired-end DNA sequences)

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

References

  • Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009, 25(14):1754-60. DOI: 10.1093/bioinformatics/btp324
  • Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010, 26(5):589-95. DOI: 10.1093/bioinformatics/btp698
  • Jiao Y, Peluso P, Shi J, Liang T, Stitzer MC, Wang B, Campbell MS, Stein JC, Wei X, Chin CS, Guill K, Regulski M, Kumari S, Olson A, Gent J, Schneider KL, Wolfgruber TK, May MR, Springer NM, Antoniou E, McCombie WR, Presting GG, McMullen M, Ross-Ibarra J, Dawe RK, Hastie A, Rank DR, Ware D. Improved maize reference genome with single-molecule technologies. Nature. 2017, 546(7659):524-7. DOI: 10.1038/nature22971