ゲノム DNA をシーケンスした、1_hakusai と 5_hiyokomame のリードをそれぞれハクサイとヒヨコ豆のリファレンスゲノム配列にマッピングしてみます。
TODO: 作成時は 1_hakusai_genome.sam だが、続きの説明が hakusai.sam なのを調整
リファレンス配列は GenomeBento/databases/genomes フォルダに用意した 3711-hakusai.fasta を、マッピングするファイルは MinION/1_hakusai のフォルダの中の fastq_pass にある全ての FASTQ ファイル *.fastq をつなげた samples/1_hakusai.fastq ファイルを指定しています。
cd ~/Desktop/GenomeBento
minimap2 -a databases/genomes/3711-hakusai.fasta samples/1_hakusai.fastq > mappings/1_hakusai_genome.samなお、minimap2 コマンドが見つからない場合、パスが通ったところにインストールされていない可能性があります。今回は GenomeBento フォルダの tools ディレクトリの中にも minimap2 コマンドを入れてありますので、./tools/minimap2 のように書き換えると実行できるはずです。
cd ~/Desktop/GenomeBento
./tools/minimap2 -a databases/genomes/3711-hakusai.fasta samples/1_hakusai.fastq > mappings/1_hakusai_genome.samヒヨコ豆のデータを使う場合は、ハクサイの例からデータベース名を 3827-chickpea.fasta に変更し、FASTQ ファイルは samples/5_hiyokomame.fastq を使います。
cd ~/Desktop/GenomeBento
minimap2 -a databases/genomes/3827-chickpea.fasta samples/5_hiyokomame.fastq > mapping/5_chickpea_genome.samヒヨコ豆の場合、以下の説明では hakusai を chickpea に読み替えてください。
今回ニンジンとトマトはゲノムではなく rbcL 領域を PCR で増幅した試料を使っていますが、これらも一応ゲノムを用意してありますので、ついでにマッピングしてみます。
cd ~/Desktop/GenomeBento
minimap2 -a databases/genomes/79200-carrot.fasta samples/2_ninjin.fastq > mapping/2_ninjin_genome.sam
minimap2 -a databases/genomes/4081-tomato.fasta samples/6_tomato.fastq > mapping/6_tomato_genome.samIGV は、米国 Broad Institute で開発されているゲノムブラウザです。リファレンスゲノムおよびマッピングされたリードを表示することができ、ヒトゲノムのように巨大なゲノムにも対応しています。
はじめて起動した場合、ヒトゲノムのバージョン hg19 が選ばれた状態になっています。
他の生物種のゲノムに切り替えるためには、IGV の Genomes メニューから「Load Genome from File...」を選択します。
ハクサイだったら 3711-hakusai.genome など、各グループで読んだ対象の食材の .genome ファイルを選択して Open をクリックします。
読み込めたらこのようにメニューからハクサイの各染色体に相当する DNA 配列がメニューに出てきます。
ここで読み込んだ 3711-hakusai.genome のような .genome ファイルはリファレンス配列の準備で作成したものです。
SAM ファイルは、できたままだとマッピングされたリードがゲノム座標順に並んでいないので、ゲノムブラウザで見る際に効率が悪くなります。このため、通常はソートしてインデックス作成したものを利用します。後述の samtools でもできますが、ここではIGVに内蔵されている igvtools を使ってみます。
Tools メニューから「Run igvtools...」を選んで igvtools を起動します。
Command を Sort に変更して、Input File の Browse ボタンをクリックして、作成した hakusai.sam を選択します。
Run ボタンでソートを実行すると hakusai.sorted.sam のように名前に .sorted がついた SAM ファイルができます。Done と表示されれば完了です（わりと一瞬で終わるかと思います）。
Command を Index に変更して、この hakusai.sorted.sam ファイルに高速化のためのインデックスを作成します（Sort から Index に変更しただけだと hakusai.sorted.sam ではなく元の hakusai.sam ファイルが選ばれている可能性があるので注意して sorted ファイルの方に変更します）。
結果、元の hakusai.sam ファイルに対して、ソートされた hakusai.sorted.sam ファイルと、それに対するインデックスファイル hakusai.sorted.sam.sai ができます。
「Load from File...」から、作成した hakusai.sorted.sam ファイルを選択します。
あとはズームしていくと、マッピングされたリードが出てきます。
リードの数が少ない場合、マップされた領域を探すのがちょっと大変ですが、SAM ファイルの中を見ると配列ごとにマッピングされた座標が分かるのでヒントになります。このことから、ゲノムを完全に読むには大量にシーケンスする必要があることが分かると思います。
BAM ファイルはSAMファイルをバイナリにして、省スペース化および高速化をはかるためのフォーマットです。SAMファイルをBAMファイルに変換するには samtools を使います。samtools を使うと、リードのマッピング率なども調べることができます。
samtools view -bS hakusai.sam > hakusai.bamsamtools sort hakusai.bam hakusai.sortedsamtools index hakusai.sorted.bamsamtools flagstat hakusai.sorted.bam今回のハクサイの例では下記のような結果になりました。
% samtools flagstat 1_hakusai_genome.sorted.bam
115592 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
105781 + 0 mapped (91.51%:nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (nan%:nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (nan%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
91%のリードがマッピングされたということで、なかなか良かったのではないでしょうか。
TODO: bwa から minimap2 に変更した部分を調整
bwa で MinION のロングリードをマッピングするには、サブコマンド mem を利用します。コマンドラインオプションは下記のようになります。
bwa mem リファレンス配列名 入力FASTQファイル > 出力SAMファイルマッピング対象とするリファレンス配列は、あらかじめ準備しておく必要があります。今回は、デスクトップの、GenomeBentoの中のdatabasesフォルダ以下に置いてあります。
GenomeBento/databases
|-- genomes
| |-- banana
| |-- cabbage
| |-- carrot
| |-- chickpea
| |-- hakusai
| |-- human
| |-- rice
| `-- tomato
`-- rbcL
|-- rbcL_all
`-- rbcL_bento
実際はそれぞれ、.amb, .ann, .bwt, .pac, .sa の拡張子がついたファイルの組になります。
これらのファイルの作り方に興味のある方は、リファレンス配列の準備を参照してください。