NGSデータ解析まとめ

サカナ研究者の手探りNGS解析（おもに進化生物学）

IGVの使い方(1): ヒトゲノムannotation情報の取り込み

NGS勉強会 IGV Genome

IGVにRefSeqなど、ヒトゲノム上の遺伝子アノテーションを取り込む方法についてです。

IGVでは、"File" メニューの"Load from Server"から、ゲノム上の位置に関連付けられたアノテーション情報を読み込むことができます。

(1) "File" -> "Load from Server"の中身について

このメニューを選択すると、以下のような画面が出ます。

f:id:hashiyuki:20160107165903j:plain

このうち、Available Datasetの一番上にある"Annotations"の中に、いくつかの呼び出し可能なデータがあります。"Annotations"の左側の三角をクリックすると、下図のようにデータの各項目が出てきます。

f:id:hashiyuki:20160107182254j:plain

さらに各項目の左側の三角をクリックすると、各データ内の各項目が出てきます（下図）。たとえば"Genes"の項目の中にも、"Ensembl Genes"や"UCSC Genes"など、異なる由来のアノテーションが含まれます。

f:id:hashiyuki:20160107183657j:plain

(2) "Annotations"の各項目について

  • Genes: 遺伝子のアノテーション。EnsemblやUCSCなど、複数の異なる機関によるアノテーションが存在する。通常はどれか一つでOK。
  • Phenotype and Disease Associations: 表現型および病気との関連（直訳）、"OMIM Genes"は遺伝子、"GWAS catalog"はGenome wide association studiesの情報で、SNPsです。
  • Sequence and Regulation: CpG islandの位置及び領域ごとのGC%を表示します。"GC percentage"がわりと重いです。要注意。
  • Variation and Repeats: "dbSNP"はSNPsのデータベースで、"Repeat Masker"はそのまま、RepeatMaskerで予測した繰り返し配列の位置を表示します。
  • Comparative Genomics: これについては別のエントリで詳細を見てみます（予定）


(3) 一般的な注意

"File" -> "Load from Server"でデータを取り込む際に、データ項目のtreeの上流の方でチェックを入れてしまうと（たとえば、"Annotations"そのものにチェックを入れる）、膨大なデータを取り込もうとするので、IGV自体が止まってしまいます。ただ、これはちょっとわかりにくい仕様のような気がします。IGVのdocumentsを見ると、ちゃんと以下のような注意書きがあるのですが。

Warning: Selecting a folder selects all of its subfolders and all of the datasets in those folders. This can potentially be a huge amount of data. To be sure you are loading only the data you want, open collapsed folders and select only the datasets of interest.
和訳- 注意：一つのフォルダを選択すると、そのフォルダ内のすべてのサブフォルダ及びそれらのフォルダ内のすべてのデータセットを選択してしまいます。この場合、潜在的に膨大な量のデータになることがあります。すべてのフォルダを開いて、必要なデータセットだけを選択して、必要なデータのみを読み込むよう、気を付けて下さい。

IGV (Integrative Genomic Viewer)を使ってみる

NGS勉強会 Genome IGV

次に、IGVで実際にリファレンスゲノムにreadsがマップされた結果を見てみます。

まず、IGVを開きます。デスクトップ上のIGVのエイリアスをダブルクリックします。このとき、IGVでは"Loading genome"という表示が一瞬出ますが、このときリファレンスゲノムとしてhg19を読み込んでいるようです。

IGVが起動したら、自分のサンプルを読み込みます。

(1) BAM fileを読み込む

"File"メニューから"Load from File..."を選択し、開きたいBAM file（indexあり）を選択します。
Mappingされた状態の見方ですが、まず画面左上の2番目の選択する場所（図1）で、"All"から"chr1"などに選択を変更します。

f:id:hashiyuki:20151207155057j:plain
図1. IGVの染色体番号を選択する部分

次に、IGV画面右上のZoomを選択する部分（図2）で、倍率を適当に上げていくと（+側に動かす）、BAMファイルに記録されている個々のreadの位置が表示されます。

f:id:hashiyuki:20151207155110j:plain
図2. IGVの右上（Zoom）

ただし、最初は見ている染色体上の位置が、何もないところ（セントロメア上）になっているので、赤で示されているカーソルを適当な位置に動かします（図3）。

f:id:hashiyuki:20151207160016j:plain
図3. 染色体の位置とカーソル（赤い四角）

(2) VCF file（SNP等のアノテーション）を読み込む

"File"メニューから"Load from File..."を選択し、開きたいVCF fileを選択します。VCFファイルは圧縮された状態（.gz）をそのまま開けます（解凍する必要なし）。

うまく開けると図4のような状態になります。

f:id:hashiyuki:20151207160533j:plain
図4. IGV window (BAM fileとVCF fileを読み込んだ状態）

"Sample"の行に示されているのがreferenceのgenotypeの状態で、水色はreferenceと異なるalleleのhomozygote, 青はheterozygoteになっています。色の部分にカーソルを合わせると、各SNPの情報を見ることができます（おそらくVCF fileに記載されているもの）。

IGVの具体的な使い方については次回以降のエントリで説明します。

IGV (Integrative Genomic Viewer)のインストール（Macの場合）

NGS勉強会 Genome メモ Tools IGV

Broad Instituteが提供しているIGV (Integrative Genomic Viewer)のインストールについての覚え書きです。

まず最初に、以下のIGVのPageを開きます。

Home | Integrative Genomics Viewer

"Download"のところをクリックします。最初は、ユーザ登録を求められるので、登録のページに移動してメールアドレスを入力します。メールアドレスを送ると、ダウンロードのページに移動できます。

IGVのインストール方法は3つあります。
(1) Mac applicationとしてIGVをインストールする
(2) "Use the Java Web Start buttons"を使用する
(3) バイナリをダウンロードして、コマンドラインからIGVを立ち上げる

このうち(1)ですが、このAppを使用するためには、MacにJava 7がインストールされている必要があります。Java 6 (JRE1.6)を使用している場合は、(2)の方法でインストールする必要があります。自分のMacでJavaのversionを確認するためには、ターミナルを起動して

java -verison

と入力します。

たとえば私のMac (OS10.8.5) <--ちょっと古い　では、こうでした。Version 6です。

java version "1.6.0_65"
Java(TM) SE Runtime Environment (build 1.6.0_65-b14-462-11M4609)
Java HotSpot(TM) 64-Bit Server VM (build 20.65-b04-462, mixed mode)

そこで、(2)の方法で行きます（もし最初からJava 7がインストールされている場合は、(1)でダウンロードしたアプリをダブルクリックするだけでIGVが開くと思います）。

まず、Dowanloadページから"Java Web Start"をインストールします。うまくいけば、ダウンロードされたファイル"igv_hm.jnlp"（注）をダブルクリックするとIGVのインストールが始まり、起動します。他のソフト（テキストエディタ等）で開いてしまう場合、「右クリック」-> 「このアプリケーションで開く」-> 「Java Web Start.app」とします。

うまく行かない場合、Java Runtime Environmentをupdateする必要があり、この場合はOSの指示に従って最新のJREをインストールします。

インストールが完了すると、IGVのエイリアスがデスクトップに生成されます。

(3) の方法については試していません（Linuxなどではこちらを使う？）


注：igv_hm.jnlpがダウンロードされる位置

/private/var/folders/nz/cgn91zkn2vg1lp8fbnvmbgvr0000gn/T/igv_hm.jnlp

BAM fileからunmapped readsを抽出する（samtoolsを使用する）

Genome samtools メモ

BAM file (.bam)から、reference sequencesにmapされなかったreadsを抽出する方法について。以下、ちょっとした覚え書き。

BAM fileはbwaなどでreference sequence(s)にNGS readsをマッピングした結果のoutput file であるSAM file (.sam)を圧縮したファイル。基本的には、samtools view コマンドで作成します。

samtools view -bS example.sam > example.bam

マッピング結果から多型サイトを抽出したり、contigをexportしたり、などの作業は、BAM fileをもとに行うことが多いです。また、BAM indexを作成して、samtools tviewやTabletなどでmapping結果を見る際にも使用します。

samtools sort example.bam example.sort
samtools index example.sort.bam
samtools view example.sort.bam reference.fasta

以下は、BAM fileから"mappingされなかった" readsをFASTQ formatで抽出する方法について。

具体的には、novoalignの解説サイト

Extracting Unmapped Reads from a BAM File Produced by NovoAlign | Novocraft

を参考にしています。ていうか、以下はこのサイトの（ほぼ）日本語訳です。

          • -

以下ここから

Novoalignで作られたBAM fileからunmapped readsを抽出する (extracting)方法について

Introduction

Referenceにmapされなかった配列をBAM fileから抽出することがが必要となるケースはしばしば生じる。Samtoolsを使うことで、この作業は容易に行うことができる。
（中略）

この作業は以下の2つのステップで行われる。

1. Unmapped readsを、readnameでsortされたBAM fileから抽出する
2. BAM fileをFASTQ read fileに変換する

FASTQ fileの生成をskipして、(unmapped) readsを直接、BAM fileから再度アラインメントすることも可能である。BAMをFASTQに変換するために、Hydra packageの"bamToFastq"を使用する。
（中略）

Novoalignによるアラインメント

# Illumina paired-end: mean=500, sd=100
novoalign -d database.nix -f reads1.fastq -reads2.fastq -oSAM  -i 500 100 | samtools view -bS - alignments

Unmapped readsに対するフィルタリング

次に、以下3つのカテゴリーでreadsのフィルタリングを行う。（注：ここではIllumina NGSによるpaired end readsを想定している）

1. Paired endの一方はマッピングされているが、それ自体はマッピングされていない
2. それ自体はマッピングされているが、paired endのもう一方はマッピングされていない
3. Paired endの両方ともマッピングされていない

1-3の各カテゴリは、以下のsamtools filtering commandで抽出することができる。

samtools view -u  -f 4 -F264 alignments.bam  > tmps1.bam
samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam
samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam

上のコマンドは、sortされたBAM fileでも、sortされていないBAM fileでも可能である。

次に、これら各カテゴリのBAM fileを結合 (merge)し、readsの名前順にsortする。ここでのsortはpaired endの相手を正しく得るために必要である。

samtools merge -u - tmps[123].bam | samtools sort -n - unmapped

ReadsをBAM fileから直接再アラインメントする

novoalign  -F BAMPE -f unmapped.bam -d …

または

ReadsをFASTQ formatで抽出する

次に、unmapped read pairsをそれぞれのFASTQ filesに抽出する（注：paired end readsを想定）。ここではHydra-SV packageのbamToFastqプログラムを利用する。

bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq

Single endの場合はもっと簡単に、以下の1行コマンドで抽出できる。

samtools view -uf 4 alignments.bam | bamToFastq -bam stdin -fq1 unmapped.fastq -fq2  empty.fastq

ここまで

          • -

わからないのは、3つのカテゴリをsamtoolsで抽出するところの

samtools view -u  -f 4 -F 264 alignments.bam  > tmps1.bam
samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam
samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam

の部分ですが、samtoolsにおけるoption -fや-Fを理解するためには、SAM fileの形式を理解する必要がありそうです。SAM fileがどのように記述されているか、については、日本語で以下のサイトが参考になりました。

http://crusade1096.web.fc2.com/sam.html

SAM fileの2列目は、FLAGという形式で、readの状況を表わしているようです（2進数の10進法表示！？）。いずれにしても、samtoolsを使いこなすには、まずSAM fileの意味を理解する必要がありそうです・・・

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

NGS勉強会 多型解析 VCF format Genome samtools

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:"のところでアップロードして、解析を実行します。解析結果の見方はこちらを参照して下さい。