V1

ゲノムReseq、変異解析

Copyright © Amelieff Corporation All Rights Reserved.

本講義にあたって• 代表的な解析の流れを紹介します– 論文でよく使用されているツールを使用します• コマンドを沢山実行します– スペルミスが心配な方は、コマンド例がありますのでコピーして実行してください•

TRY!

マークのコマンドは実行してください。
– 実行が遅れてもあせらずに、応用や課題の間に追い付いてください

Copyright © Amelieff Corporation All Rights Reserved.

2

本講義の内容• Reseq解析• RNA-seq解析

公開データ取得

公開データ取得↓↓

クオリティコントロール

クオリティコントロール↓↓

マッピングマッピング

↓

↓

変異検出

発現定量
SNVとIndel検出を行います。

Copyright © Amelieff Corporation All Rights Reserved.

FPKMを算出します。

3Reseq解析：検出可能な変異• ショートリードのシーケンスでも様々な変異を検出可能
SNVInDelDuplication Translocation

Inversion
CNV• 検出アルゴリズムとソフトウェア
Paired-end mapping
Split-read mapping
Others、Complex:::

BreakDancer、VariationHunter
PindelCREST、DELLY

充分に精度が高いとは言えません。
Copyright © Amelieff Corporation All Rights Reserved.

4Reseq解析：パイプライン
データ取得 → クオリティコントロール → マッピング→変異検出

解析パイプラインとは「あるソフトの出力結果が、次のソフトの入力ファイルとなる」連続した解析処理の流れ。
Copyright © Amelieff Corporation All Rights Reserved.

5Reseq解析：パイプライン
データ取得 → クオリティコントロール → マッピング→変異検出

今日は一部のコマンドを実行します。

サンプル間比較、遺伝モデルを使用した絞り込み
Genotype imputation など様々。
Copyright © Amelieff Corporation All Rights Reserved.

6

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンス取得– http://support.illumina.com/sequencing/sequencing_software/igenome.html

リファレンスのfastaのみではなく、マッピングソフトのインデックスファイルや遺伝子情報ファイルも一緒に圧縮されて公開しています。
Copyright © Amelieff Corporation All Rights Reserved.

7

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンス取得（実行済み）
ダウンロードして、解凍します。
$ wget ftp://igenome:G3nom3s4u@ussdftp.illumina.com/Saccharomyces_cerevisiae/NCBI/build3.1/Saccharom
yces_cerevisiae_NCBI_build3.1.tar.gz$ tar zxvf Saccharomyces_cerevisiae_NCBI_build3.1.tar.gz※お手元のテストデータでは、使用しないデータを一部削除しています
Copyright © Amelieff Corporation All Rights Reserved.

8

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンスを確認$ cd /home/ユーザ名/Desktop/amelieff/Scerevisiae$ llCopyright © Amelieff Corporation All Rights Reserved.

9

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンスを確認$ ll WholeGenomeFasta

Copyright © Amelieff Corporation All Rights Reserved.

10

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンスを確認$ less WholeGenomeFasta/genome.fa
ヘッダには、コンティグ名が記載されます。

「q」で閲覧を終了します。
Copyright © Amelieff Corporation All Rights Reserved.

11

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンスを確認$ less WholeGenomeFasta/genome.fa.fai
インデックスファイルを開きます。
SamToolsで作成できます。

：
1列目：
2列目：
3列目：
4列名：
5列目：

コンティグ名（fastaファイルのヘッダ）
コンティグの長さ
ファイルの先頭から見た、染色体の第一塩基目の位置
fastaの1行の文字数
各行のバイト数
Copyright © Amelieff Corporation All Rights Reserved.

12応用）ヒトリファレンスの話

GRCh Build37 ＋ デコイ配列 Version 5
ヒトWhole Genome Sequencing Cloneを「ヒトゲノム＋ヒトヘルペスウイルスHHV-4 」にマッピングして、よくマップできなかったものを集めたもの。
サイズ：
特徴 ：

合計35.4Mb、N50=22.9kb
50%はサテライト配列またはシンプルリピート、20%はレトロトランスポゾン※現在は、2013/12/24にメジャーアップしたGRCh38が公開されています。
Copyright © Amelieff Corporation All Rights Reserved.

13応用）ヒトリファレンスの話

GRCh Build37 ＋ デコイ配列 Version 5
Without Decoy
最大カバレージ： 1112

With Decoy

最大カバレージ： 817

Reseq解析は、リファレンスに対して変異検出するので、リファレンス自体がどの程度確かなのかが非常に大切
Copyright © Amelieff Corporation All Rights Reserved.

14

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• 酵母のゲノムのリファレンスを確認$ ll Scerevisiae/BWAIndex/
リンクの「l」

BWAのインデックス
ファイルを開きます

シンボリックリンク名 -> 実体のファイル…

Copyright © Amelieff Corporation All Rights Reserved.

15

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• リファレンスのインデックスを作成
BWA バージョン0.7のインデックスファイルを作成します。
BWAの使い方を確認します。

$ bwa index$ mkdir BWAIndex/version0.7.12$ cd BWAIndex/version0.7.12

Copyright © Amelieff Corporation All Rights Reserved.

16

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• リファレンスのインデックスを作成
シンボリックリンクを作成します。
$ ln –s 実体のファイル$ ln -s ../../WholeGenomeFasta/genome.fa$ ll…

Copyright © Amelieff Corporation All Rights Reserved.

17

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• リファレンスのインデックスを作成
インデックスを作成します。
$ bwa index genome.fa$ ll…

Copyright © Amelieff Corporation All Rights Reserved.

18

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得

http://trace.ddbj.nig.ac.jp/dra/index.html

DDBJのSequence Read Archive →Copyright © Amelieff Corporation All Rights Reserved.

Search19

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得

Accessionに「ERR038793」と入力 → Search
Copyright © Amelieff Corporation All Rights Reserved.

20

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得
ここからダウンロード

実験の詳細

NavigationエリアのExperiment→「ERX015989」をクリック

Copyright © Amelieff Corporation All Rights Reserved.

21

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得
他にも、シーケンサのプラットフォームやリード長などの情報も記載されています。
Whole Genome Sequencing

Copyright © Amelieff Corporation All Rights Reserved.

22

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得（実行済み）
ダウンロードします。
$ wgetftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/ERA038/ERA038218/ERX015989/ERR038793_1.fastq.bz2$ wgetftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/ERA038/ERA038218/ERX015989/ERR038793_2.fastq.bz2

Copyright © Amelieff Corporation All Rights Reserved.

23Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータ取得（実行済み）
解凍して、先頭1000リードを抽出します。
$ bunzip2 ERR038793_1.fastq.bz2$ bunzip2 ERR038793_2.fastq.bz2$ head -4000 ERR038793_1.fastq > 1K_ERR038793_1.fastq$ head -4000 ERR038793_2.fastq > 1K_ERR038793_2.fastq

Copyright © Amelieff Corporation All Rights Reserved.

24

TRY!

Reseq解析：データ
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータを確認$ cd /home/ユーザ名/Desktop/amelieff/$ ll:

行数を数えます。 1リードは4行で表記されます。
$wc -l 1K_ERR038793_1.fastq

4000 1K_ERR038793_1.fastq
Copyright © Amelieff Corporation All Rights Reserved.

25

TRY!

Reseq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータのクオリティを確認
インストールされているFastQCの、バージョンと使い方を確認します。
$ fastqc -version
FastQC v0.10.1
Fastqのみではなく、bamとsamも入力可能$ fastqc -help：：

複数のファイルも指定可能

Copyright © Amelieff Corporation All Rights Reserved.

26

TRY!

Reseq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→変異検出• シーケンスデータのクオリティを確認
FastQCを実行します。
$ mkdir reseq$ fastqc -o reseq -f fastq 1K_ERR038793_1.fastq 1K_ERR038793_2.fastq

fastqc_report.htmlを、ウェブブラウザで開きます。
$ firefox reseq/1K_ERR038793_1_fastqc/fastqc_report.html$ firefox reseq/1K_ERR038793_2_fastqc/fastqc_report.html

Copyright © Amelieff Corporation All Rights Reserved.

27応用）とあるシーケンスデータの実例

最初の１塩基の割合が不自然

リード末端でクオリティが低下

マッピング率が低下や、変異の偽陽性が増加するなどの問題を引き起こす。

シーケンス技術が向上しクオリティの高いデータを目にする機会が増えましたが、試料・シーケンス・トリミングなどに、問題がないか確認することをおすすめします。
Copyright © Amelieff Corporation All Rights Reserved.

28

TRY!

Reseq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→変異検出• クオリティ30以上の塩基が90%未満のリードを削除
インストールされているfastq_quality_filterの使い方を確認します。
$ fastq_quality_filter -h

Copyright © Amelieff Corporation All Rights Reserved.

29

TRY!

Reseq解析：クオリティコントロール
データ取得 → クオリティコントロール → マッピング→変異検出• クオリティ30以上の塩基が90%未満のリードを削除$ fastq_quality_filter -i 1K_ERR038793_1.fastq-o reseq/1K_ERR038793_1_qual.fastq -q 30 -p 90 -Q 33 -v
Quality cut-off: 30
ターミナルに直接解析のサマリー
Minimum percentage: 90を出力するソフトもあります。
Input: 1000 reads.
Output: 802 reads.
discarded 198 (19%) low-quality reads.
以降の解析は、片側のリードのみ使用します。
Copyright © Amelieff Corporation All Rights Reserved.

30応用）クオリティコントロールの順番も大切
FASTQ形式にマッチするかチェック
データクオリティチェック（FastQC）
Illumina CASAVA filter [Y] を除去

ロングリードの場合、リードの大半が除外されてしまう可能性。

クオリティ20未満が80%以上のリードを除去
クオリティ20未満の末端をトリム
未知の塩基(N)が多いリード除去
配列長が短いリード除去
片側のみのリードを除外

ペアエンドリードの場合、ペアが揃っていないとマッピングソフトが停止する可能性。

データクオリティチェック（FastQC）

Copyright © Amelieff Corporation All Rights Reserved.

31

TRY!

Reseq解析：マッピング
データ取得 → クオリティコントロール → マッピング→変異検出• Bwa memコマンドの使い方を確認$ bwa mem※RG（read groups）
platform (PL) および sample (SM)が必要
PLの例：454, LS454, Illumina, Solid, ABI_SolidCopyright © Amelieff Corporation All Rights Reserved.

32

TRY!

Reseq解析：マッピング
データ取得 → クオリティコントロール → マッピング→変異検出• マッピング$ cd reseq$ bwa mem -R "@RG¥tID:1K_ERR038793_1¥tSM:ERR038793¥tPL:Illumina"/home/ユーザ名/Desktop/amelieff/Scerevisiae/BWAIndex/genome.fa
1K_ERR038793_1_qual.fastq > 1K_ERR038793_1_qual.sam$ llCopyright © Amelieff Corporation All Rights Reserved.

33

TRY!

Reseq解析：マッピング
データ取得 → クオリティコントロール → マッピング→変異検出• SAMをBAMに変換$ samtools view -Sb 1K_ERR038793_1_qual.sam > 1K_ERR038793_1_qual.bam$ ll -h1/4程度にファイルサイズが小さくなりました。

Copyright © Amelieff Corporation All Rights Reserved.

34

TRY!

Reseq解析：マッピング
データ取得 → クオリティコントロール → マッピング→変異検出• ソートとインデキシング$ samtools sort 1K_ERR038793_1_qual.bam 1K_ERR038793_1_qual_sorted$ samtools index 1K_ERR038793_1_qual_sorted.bam$ ll

Copyright © Amelieff Corporation All Rights Reserved.

35

TRY!

Reseq解析：マッピング
データ取得 → クオリティコントロール → マッピング→変異検出• マッピングされたリード数$ samtools idxstats 1K_ERR038793_1_qual_sorted.bam

コンティグ名、コンティグの長さ、マッピングされたリード、マッピングされなかったリードの順に表示されます。

3列目を足し合わせると、マッピングされたリード数がわかります。
Copyright © Amelieff Corporation All Rights Reserved.

36

応用）

列の合計を計算するコマンド$ samtools idxstats 1K_ERR038793_1_qual_sorted.bam > tmp$ awk '{a += $3} END {print a}' tmp
１行読み込むたびに、３列目を「a」に足す。
803

マッピングされたリード$ awk '{a += $4} END {print a}' tmp0

マッピングされなかったリード

802リードのfastqをマッピングしたはずが、1本増えています。
マルチヒットしたリードがあると考えられます。

Copyright © Amelieff Corporation All Rights Reserved.

37

TRY!

Reseq解析：変異検出
データ取得 → クオリティコントロール → マッピング→変異検出• GATK UnifiedGenotyperコマンドの使い方を確認$ java -jar /usr/local/src/GenomeAnalysisTK-1.6-13g91f02df/GenomeAnalysisTK.jar -T UnifiedGenotyper -h

SNP, INDEL, BOTH から選べます。デフォルトはSNPCopyright © Amelieff Corporation All Rights Reserved.

38

TRY!

Reseq解析：変異検出
データ取得 → クオリティコントロール → マッピング→変異検出• SNV/Indel検出$ java -jar /usr/local/src/GenomeAnalysisTK-1.6-13-g91f02df/GenomeAnalysisTK.jar-T UnifiedGenotyper -glm BOTH
-R /home/ユーザ名/Desktop/amelieff/Scerevisiae/WholeGenomeFasta/genome.fa-I 1K_ERR038793_1_qual_sorted.bam -o 1K_ERR038793_1_qual_sorted.vcf$ ll

Copyright © Amelieff Corporation All Rights Reserved.

39

TRY!

Reseq解析：変異検出
データ取得 → クオリティコントロール → マッピング→変異検出• 検出したSNV/Indelを可視化$ less 1K_ERR038793_1_qual_sorted.vcf…
ジェノタイプがC/Tのヘテロ

Copyright © Amelieff Corporation All Rights Reserved.

カバレージが６40

TRY!

Reseq解析：変異検出
データ取得 → クオリティコントロール → マッピング→変異検出• 検出したSNV/Indelの数を確認$ awk '!/^#/' 1K_ERR038793_1_qual_sorted.vcf | wc -l
100
100個の変異が検出されました
検出されるSNV/Indel数は、使用するソフトウェアのバージョンやパラメータにより変動します

Copyright © Amelieff Corporation All Rights Reserved.

41応用）リアライメント

リアライメントは必要？
BWAでは、 1本のリードに複数の変異が含まれる場合に、アライメントスコアの計算上、SNVやIndelの正確な位置を決めることが出来ません。
このような領域を対象領域として抜き出して、改めて丁寧にアライメントを行う。

$ igv.shCopyright © Amelieff Corporation All Rights Reserved.

42

Genomes →Load Genome from File…

TRY!

/home/ユーザ名/Desktop/amelieff/Scerevisiae/Wh
oleGenomeFasta/genome.faまで移動

Genome.faファイルを選択

Copyright © Amelieff Corporation All Rights Reserved.

43

File →

Load from File…

TRY!
/home/ユーザ名/Desktop/amelieff/reseqまで移動

Copyright © Amelieff Corporation All Rights Reserved.

44

TRY!

「I:111」と入力

ジェノタイプがC/Tのヘテロ

カバレージが６

Copyright © Amelieff Corporation All Rights Reserved.

45応用）Indelの見方
TTの欠失

TTTの欠失…
ジェノタイプは、AT/Aホモポリマーではシーケンスエラーによって、偽陽性のIndelが検出されやすい。

Copyright © Amelieff Corporation All Rights Reserved.

46応用）変異のフィルタリング•GATKのVariantFiltrationコマンドでフィルタリングをします$ java -jar /usr/local/src/GenomeAnalysisTK-1.6-13-g91f02df/GenomeAnalysisTK.jar-T VariantFiltration–R /home/ユーザ名/Desktop/amelieff/Scerevisiae/WholeGenomeFasta/genome.fa-V 1K_ERR038793_1_qual_sorted.vcf -o 1K_ERR038793_1_qual_sorted_fil.vcf--clusterWindowSize 10 --filterExpression "DP < 10" --filterName "LowCoverage"

VCFファイルのFILTER列に、条件を通過した場合“PASS”、そうでない場合は “filterName”が記入されます。

Copyright © Amelieff Corporation All Rights Reserved.

47応用）遺伝子情報のアノテーション• snpEff…変異に対して遺伝子名や転写産物の情報、変異の影響などを付与します
snpEffを実行するには、snpEffをインストールした後、対応するゲノムのデータベースをダウンロードしておきます。
例：ヒトhg19データベースをダウンロードする$ java -jar /usr/local/src/snpEff/snpEff.jar download hg19

対応する生物種のデータベースがない場合は、データベースを作成する必要があります。
$ mkdir data/sacCer$ cd data/sacCer$ wgethttp://downloads.yeastgenome.org/curation/chromosomal_feature/saccharomyces_cerevisiae.gff$ mv saccharomyces_cerevisiae.gff genes.gff$ echo "sacCer.genome : Yeast" >> snpEff.config$ java -Xmx1G -jar snpEff.jar build -gff3 sacCerCopyright © Amelieff Corporation All Rights Reserved.

48応用）遺伝子情報のアノテーション• snpEff…変異に対して遺伝子名や転写産物の情報、変異の影響などを付与します$ java –Xmx10G –jar /usr/local/src/snpEff/snpEff.jar eff-c /usr/local/src/snpEff/snpEff.config -i vcf sacCer -o vcf
1K_ERR038793_1_qual_sorted_fil.vcf 1>
1K_ERR038793_1_qual_sorted_fil_snpeff.vcf

Copyright © Amelieff Corporation All Rights Reserved.

49

本講義の内容• Reseq解析• RNA-seq解析

公開データ取得

公開データ取得↓↓

クオリティコントロール

クオリティコントロール

↓↓

マッピング

マッピング

↓↓

変異検出

発現定量

Copyright © Amelieff Corporation All Rights Reserved.

50