WGSデータの参照ゲノム配列へのマッピング (5): VCF fileの生成、アノテーション

2015-08-06

Mapping結果の.bam fileから、SNPなどの多型を抽出して、VCF formatのfileを生成します。その後、Variant Effect Predictorを使ってSNPのアノテーションを行います。方法は以前紹介したこのエントリ（VCF fileの生成）や、このエントリ（Variant Effect Predictorの使い方その1）、このエントリ（Variant Effect Predictorの使い方その2）とほとんど同じです。

(1) 多型情報の抽出（VCF fileの生成）

samtools (bcftools)を使って以下のコマンドを実行します。

# mpileup // やや時間がかかる（30分-1時間）
samtools mpileup -g -f Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa ERR251633.chr6.sort.bam > ERR251633.chr6.bcf
# bcftools view
bcftools view -bvcg ERR251633.chr6.bcf > ERR251633.chr6.var.bcf
# bcftools vcfutils.pl
bcftools view ERR251633.chr6.var.bcf | vcfutils.pl varFilter -> ERR251633.chr6.var.vcf

(2) Variant Effect PredictorによるSNPsのアノテーション

Variant Effect Predictorのページに移動し、"ERR251633.chr6.var.vcf"を"Or upload file:"のところでアップロードして、解析を実行します。解析結果の見方はこちらを参照して下さい。

WGSデータの参照ゲノム配列へのマッピング (4): BWAによるマッピング

2015-08-06

ヒトゲノム参照配列に対して、NGSのリード（ERR251633）をbwaを使ってマッピングします。この辺りの解析は、以前のエントリと大体同じです。

ただし、今回はmappingの前にNGSのraw dataに前処理をします。MappingするFASTQ fileに対して、qualityの低いリードを除き、長さの短い（<50bp）リードも除去します。前処理を行うことで、bwaの解析にかかる時間の短縮と、精度の向上が期待できます（150805 --> 本当に？　どれくらい？）。具体的な方法は、井上さんによる以下のサイトを参考にしています。

参考URL: Trinity - 井上 潤

必要な解析ソフト（SolexaQA, DynamicTrim, LengthSort）
SolexaQA

(1) Index fileの作成

まずリファレンス配列のindexを作成します。ここでは時間の都合上、ヒトゲノム参照配列の6番染色体を使用します。データの入手はこちらを参照。

# bwa index実行 // 3分くらいかかる
bwa index -a bwtsw Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa

(2) NGS reads (.fastq) ファイルの前処理

.sra fileからFASTQ fileを抽出

以下のコマンドを実行します。ここではpaired endの右側・左側のデータをそれぞれ別のfileに分けます。

fastq-dump --split-files ERR251633.sra

DynamicTrimによるlow quality readsの除去

以下のようにコマンドを実行します。

DynamicTrim.pl ERR251633_1.fastq -h 20 &
DynamicTrim.pl ERR251633_2.fastq -h 20 &

LengthSortによる短すぎるreadsの除去

以下のようにコマンドを実行します。ここでは50bpより短いreadsを除きます。

LengthSort.pl ERR251633_1.fastq.trimmed ERR251633_2.fastq.trimmed -length 50

ここで生成される"ERR251633_1.fastq.trimmed.paired1"と"ERR251633_1.fastq.trimmed.paired2"が、以下のbwaでのmappingに使用するFASTQ fileになります。

(3) BWAによるmappingの実行

ここではbwa memを使用します。以下のコマンドを実行します。

# bwa mem実行 // 約21分 (MacPro 2 x 3.06GHz 6 core Intel Xeon, 32GB RAM)
bwa mem -M -t 12 Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa ../ERR251633_1.fastq.trimmed.paired1 ../ERR251633_1.fastq.trimmed.paired2 > ERR251633.chr6.sam

(4) samtoolsによるSAM->BAM変換、sorting, indexing

以下のコマンドを順次実行します。

samtools view -bS ERR251633.chr6.sam > ERR251633.chr6.bam    # .sam -> .bam変換
samtools sort ERR251633.chr6.bam ERR251633.chr6.sort    # sorting
samtools index ERR251633.chr6.sort.bam

.bam fileはsamtools tviewコマンドや、Tabletなどのviewerを使って見ることができます。

WGSデータの参照ゲノム配列へのマッピング (3): ヒトゲノム-リファレンス配列（GRCh37 or GRCh38）のダウンロード

2015-08-05

ヒトゲノムのリファレンス配列は、GRCh37またはhg19と呼ばれるデータセットがよく使われています（2013年12月以降は、GRCh38 (hg38)が最新のようです）。

GRCh37（またはGRCh38）は、UCSC Genome Browserのページからダウンロードできます。

hgdownload.soe.ucsc.edu

UCSC Genome Browser -> Genome -> Downloads とメニューをたどっていくと、このページに行きます。このページの下の方に"Human genome"というカテゴリがあるので、そこからFull data setをクリックすると、ダウンロードのページに移動できます。

ここでは、"chromFaMasked.tar.gz"をダウンロードします（GRCh37の場合）。展開すると各染色体ごとのfasta fileの入ったフォルダができます。

（2015-08-06追記）あとでVariant Effect Predictorを使う場合は、対応するリファレンスゲノム配列がGRCh38.p3なので、こちらを使用する必要があります。

WGSデータの参照ゲノム配列へのマッピング (2): データのダウンロード

2015-08-06

まず最初に、マッピングを行うNGSのサンプルデータを入手します。

このデータを使用します。
HapMapで使用されている、東京に住む日本人のデータになります。Illumina HiSeq 2000で読まれたデータ（ペアエンド）で、.sraフォーマットで8.9GBあります。詳しくはこのURLも参照して下さい。

データはこのFTPサイトに置かれています。Firefoxなどのブラウザなら、直接FTPサイトからダウンロードできます。Safariを使用している場合は、直接FTPサイトに接続できないので、ターミナルからftpコマンドを使って接続します。

以下のようにします。

(1) まず、データをダウンロードしたい（自分のMacの）ディレクトリに移動します。

たとえば、今日の日付フォルダ（"150806"など）を適当な場所に作成して、cdコマンドでそのフォルダ内に移動します。

(2) ターミナルのシェルに、以下のように入力します。

ftp ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR251/ERR251633/

うまくつながれば、一連の文字が表示されたあとに

…
250 CWD command successful
ftp>

と、"ftp>"コマンド入力を受け付けるモードになります。

このとき、"ERR251633/"の最後の"/"を忘れると

221 goodbye

と表示されて、FTPサイトから抜けてしまうので要注意です。

(3) FTPサイトのディレクトリの位置を確認します。"dir"コマンドを使います。

dir

と入力して

229 Entering Extended Passive Mode (|||50373|)
150 Opening ASCII mode data connection for file list
-r--r--r--   1 ftp      anonymous 8862493705 Apr 10  2013 ERR251633.sra
226 Transfer complete

のように表示されればOKです。

データのダウンロードは、"get"コマンドを使って以下のように行います。

get ERR251633.sra

8.9GBあるので、データのダウンロードにはそれなりに時間がかかります（一晩くらいかかる）。

ちなみに、データのダウンロードが終わってFTPサイトを抜けるときは

bye

と入力します。

追記（LFTPを使う）："LFTP"が使えれば、こちらのほうが早いと思います。LFTPはこのサイトから入手できます。Mac OSの場合、packageを展開するだけでインストールできました。使い方は、コマンドをftpからlftpにするだけで、あとは同じです。

WGSデータの参照ゲノム配列へのマッピング (1): 解析の概要

2015-08-06

今回は、NGS (Illumina HiSeq) によって配列決定されたヒト全ゲノムのデータを、ヒトの参照ゲノム配列にマッピングします。ここでの目標は、以下の3つです。

1. BWAを用いて、NGS readsをヒトゲノム参照配列にマッピングする
2. SNP情報を.vcf fileとして抽出する。Variant Effect Predictorを使ってアノテーションを付ける
3. 遺伝子の塩基配列を得る。各遺伝子にアノテーションを付ける
4. NGSで決定したヒト個体のゲノムの塩基配列をFASTA形式で得る

以下のエントリで、それぞれの作業について考えていきます（次回へ続く・・）。