PowerPoint プレゼンテーション

Size: px

Start display at page:

Download "PowerPoint プレゼンテーション"

えつまいせき
4 years ago
Views:

1 V1 次世代シークエンサ実習 II

2 本講義にあたって代表的な解析の流れを紹介します論文でよく使用されているツールを使用しますコマンドを沢山実行しますスペルミスが心配な方はコマンド例がありますのでコピーして実行してください /home/admin1409/amelieff/ngs/reseq_command.txt マークのコマンドは実行してください実行が遅れてもあせらずに応用や課題の間に追い付いてください 2

3 本講義の内容 Reseq解析 RNA-seq解析公開データ取得公開データ取得クオリティコントロールクオリティコントロールマッピングマッピング変異検出発現定量 SNVとIndel検出を行います FPKMを算出します 3

4 Reseq解析検出可能な変異ショートリードのシーケンスでも様々な変異を検出可能 SNV InDel Duplication Translocation Inversion CNV 検出アルゴリズムとソフトウェア Paired-end mapping Split-read mapping Others Complex : : : BreakDancer VariationHunter Pindel CREST DELLY 充分に精度が高いとは言えません 4

5 Reseq解析パイプライン解析パイプラインとはあるソフトの出力結果が次のソフトの入力ファイルとなる連続した解析処理の流れ 5

6 Reseq解析パイプライン今日は一部のコマンドを実行しますサンプル間比較遺伝モデルを使用した絞り込み Genotype imputation など様々 6

7 Reseq解析データ酵母のゲノムのリファレンス取得 7

8 Reseq解析データ酵母のゲノムのリファレンス取得リファレンスのfastaのみではなくマッピングソフトのインデックスファイルや遺伝子情報ファイルも一緒に圧縮されて公開しています 8

9 Reseq解析データ酵母のゲノムのリファレンス取得実行済みダウンロードして解凍します $ wget ftp://igenome:g3nom3s4u@ussdftp.illumina.com/saccharomyces_cerevisiae/ncbi /build3.1/saccharomyces_cerevisiae_ncbi_build3.1.tar.gz $ tar zxvf Saccharomyces_cerevisiae_NCBI_build3.1.tar.gz ファイル容量が大きいため使用しないデータを一部削除しています 9

10 Reseq解析データ酵母のゲノムのリファレンスを確認 $ cd /home/admin1409/genome/saccharomyces_cerevis iae/ncbi/build3.1/sequence $ ll 10

11 Reseq解析データ酵母のゲノムのリファレンスを確認 $ ll WholeGenomeFasta 11

12 Reseq解析データ酵母のゲノムのリファレンスを確認 $ less WholeGenomeFasta/genome.fa ヘッダにはコンティグ名が記載されます q で閲覧を終了します 12

13 Reseq解析データ酵母のゲノムのリファレンスを確認 $ less WholeGenomeFasta/genome.fa.fai インデックスファイルを開きます SamToolsで作成できます 1列目 2列目 3列目 4列名 5列目コンティグ名 fastaファイルのヘッダコンティグの長さファイルの先頭から見た染色体の第一塩基目の位置 fastaの1行の文字数各行のバイト数 13

14 応用 ) ヒトリファレンスの話 GRCh Build37 + デコイ配列 Version 5 ヒト Whole Genome Sequencing Clone をヒトゲノム + ヒトヘルペスウイルス HHV-4 にマッピングしてよくマップできなかったものを集めたものサイズ : 合計 35.4Mb N50=22.9kb 特徴 : 50% はサテライト配列またはシンプルリピート 20% はレトロトランスポゾン現在は 2013/12/24 にメジャーアップした GRCh38 が公開されています 14

15 応用ヒトリファレンスの話 GRCh Build37 デコイ配列 Version 5 Without Decoy 最大カバレージ 1112 With Decoy 最大カバレージ 817 Reseq解析はリファレンスに対して変異検出するのでリファレンス自体がどの程度確かなのかが非常に大切 15

16 Reseq解析データ酵母のゲノムのリファレンスを確認 $ ll BWAIndex/version0.6.0 リンクの l BWA バージョン0.6のインデックスファイルを開きますシンボリックリンク名 -> 実体のファイル 16

17 Reseq解析データリファレンスのインデックスを作成 BWA バージョン0.7のインデックスファイルを作成します BWAの使い方を確認します既にbwa0.7にPATHが通っている場合はbwaだけで実行できます $ /home/admin1409/amelieff/ngs/bwa /bwa index $ mkdir BWAIndex/version $ cd BWAIndex/version

18 Reseq解析データリファレンスのインデックスを作成シンボリックリンクを作成します $ ln s 実体のファイル $ ln -s../../wholegenomefasta/genome.fa $ ll 18

19 Reseq解析データリファレンスのインデックスを作成インデックスを作成します $ /home/admin1409/amelieff/ngs/bwa /bwa index genome.fa $ ll 19

Reseq解析データシーケンスデータ取得 http://trace.ddbj.nig.ac.jp/dra/index.

20 Reseq解析データシーケンスデータ取得 DDBJのSequence Read Archive Search 20

21 Reseq解析データシーケンスデータ取得 Accessionに ERR と入力 Search 21

22 Reseq解析データシーケンスデータ取得ここからダウンロード実験の詳細 NavigationエリアのExperiment ERX をクリック 22

23 Reseq解析データシーケンスデータ取得他にもシーケンサのプラットフォームやリード長などの情報も記載されています Whole Genome Sequencing 23

24 Reseq解析データシーケンスデータ取得実行済みダウンロードします $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/er A038218/ERX015989/ERR038793_1.fastq.bz2 $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/er A038218/ERX015989/ERR038793_2.fastq.bz2 24

25 Reseq解析データシーケンスデータ取得実行済み解凍して先頭1000リードを抽出します $ bunzip2 ERR038793_1.fastq.bz2 $ bunzip2 ERR038793_2.fastq.bz2 $ head ERR038793_1.fastq > 1K_ERR038793_1.fastq $ head ERR038793_2.fastq > 1K_ERR038793_2.fastq 25

26 Reseq解析データシーケンスデータを確認 $ cd /home/admin1409/amelieff/ngs $ ll : 行数を数えます 1リードは4行で表記されます $ wc -l 1K_ERR038793_1.fastq K_ERR038793_1.fastq 26

27 Reseq解析クオリティコントロールシーケンスデータのクオリティを確認インストールされているFastQCのバージョンと使い方を確認します $ fastqc -version FastQC v Fastqのみではなく bamとsamも入力可能 $ fastqc -help 複数のファイルも指定可能 27

28 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 28

29 応用とあるシーケンスデータの実例最初の１塩基の割合が不自然リード末端でクオリティが低下マッピング率が低下や変異の偽陽性が増加するなどの問題を引き起こすシーケンス技術が向上しクオリティの高いデータを目にする機会が増えましたが試料シーケンストリミングなどに問題がないか確認することをおすすめします 29

30 Reseq解析クオリティコントロールクオリティ30以上の塩基が90%未満のリードを削除インストールされているfastq_quality_filterの使い方を確認します $ fastq_quality_filter -h 30

31 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. 以降の解析は片側のリードのみ使用します 31

クオリティ20未満が80%以上のリードを除去クオリティ20未満の末端をトリム未知の塩基(N)が多いリード除去

32 応用クオリティコントロールの順番も大切 FASTQ形式にマッチするかチェックデータクオリティチェック FastQC Illumina CASAVA filter [Y] を除去ロングリードの場合リードの大半が除外されてしまう可能性クオリティ20未満が80%以上のリードを除去クオリティ20未満の末端をトリム未知の塩基(N)が多いリード除去配列長が短いリード除去片側のみのリードを除外ペアエンドリードの場合ペアが揃っていないとマッピングソフトが停止する可能性データクオリティチェック FastQC 32

33 Reseq解析マッピング Bwa memコマンドの使い方を確認 $ bwa /bwa mem RG read groups platform (PL) および sample (SM)が必要 PLの例 454, LS454, Illumina, Solid, ABI_Solid 33

34 Reseq解析マッピングマッピング $ cd reseq $../bwa /bwa mem -R "@RG tid:1k_err038793_1 tsm:err tpl:illumina" /home/admin1409/genome/saccharomyces_cerevisiae/ncbi /build3.1/sequence/bwaindex/version0.7.10/genome.fa 1K_ERR038793_1_qual.fastq > 1K_ERR038793_1_qual.sam $ ll 34

35 Reseq解析マッピング SAMをBAMに変換 $ samtools view -Sb 1K_ERR038793_1_qual.sam > 1K_ERR038793_1_qual.bam $ ll -h 1/4程度にファイルサイズが小さくなりました 35

36 Reseq解析マッピングソートとインデキシング $ samtools sort 1K_ERR038793_1_qual.bam 1K_ERR038793_1_qual_sorted $ samtools index 1K_ERR038793_1_qual_sorted.bam $ ll 36

37 Reseq解析マッピングマッピングされたリード数 $ samtools idxstats 1K_ERR038793_1_qual_sorted.bam コンティグ名コンティグの長さマッピングされたリードマッピングされなかったリードの順に表示されます 3列目を足し合わせるとマッピングされたリード数がわかります 37

38 応用列の合計を計算するコマンド $ samtools idxstats 1K_ERR038793_1_qual_sorted.bam > tmp $ awk '{a += $3} END {print a}' tmp １行読み込むたびに３列目を a に足す 803 マッピングされたリード $ awk '{a += $4} END {print a}' tmp 0 マッピングされなかったリード 802リードのfastqをマッピングしたはずが 1本増えていますマルチヒットしたリードがあると考えられます 38

$bam awk '{print $1}' sort uniq -c 登場した回数配列ID １回以上登場した配列を探します $ samtools view 1K_ERR038793_1_qual_sorted.$

39 応用マルチヒットしたリードを探索 $ samtools view 1K_ERR038793_1_qual_sorted.bam awk '{print $1}' sort uniq -c 登場した回数配列ID １回以上登場した配列を探します $ samtools view 1K_ERR038793_1_qual_sorted.bam awk '{print $1}' sort uniq -c awk '$1>1 {print}' 本当に２回マッピングされているか確認します 39

40 応用マルチヒットしたリードを探索 $ samtools view 1K_ERR038793_1_qual_sorted.bam grep ERR 異なるコンティグにマッピングされています 40

41 Reseq解析変異検出 GATK UnifiedGenotyperコマンドの使い方を確認 $ java -jar /usr/local/src/genomeanalysistk g91f02df/genomeanalysistk.jar -T UnifiedGenotyper -h SNP, INDEL, BOTH から選べますデフォルトはSNP 41

42 Reseq解析変異検出 SNV/Indel検出 $ java -jar /usr/local/src/genomeanalysistk g91f02df/genomeanalysistk.jar -T UnifiedGenotyper glm BOTH -R /home/admin1409/genome/saccharomyces_cerevisiae/ncbi /build3.1/sequence/bwaindex/version0.7.10/genome.fa -I 1K_ERR038793_1_qual_sorted.bam -o 1K_ERR038793_1_qual_sorted.vcf $ ll 42

43 Reseq解析変異検出検出したSNV/Indelの数を確認 $ less 1K_ERR038793_1_qual_sorted.vcf $ awk '!/^#/' 1K_ERR038793_1_qual_sorted.vcf wc -l 個の変異が検出されました 43

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

45 Reseq解析変異検出検出したSNV/Indelを可視化ジェノタイプがC/Tのヘテロカバレージが６ $ igv.sh 45

46 Genomes Load Genome from File /home/admin1409/genome/ Saccharomyces_cerevisiae/ NCBI/build3.1/Sequence/ WholeGenomeFasta まで移動 Genome.fa ファイルを選択 46

47 File Load from File /home/admin1409/amelieff/ ngs/reseq まで移動 47

48 I:111 と入力ジェノタイプが C/T のヘテロカバレージが 6 48

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

50 応用変異のフィルタリング GATKのVariantFiltrationコマンドでフィルタリングをします $ java -jar /usr/local/src/genomeanalysistk g91f02df/genomeanalysistk.jar -T VariantFiltration -R /home/admin1409/genome/saccharomyces_cerevisiae/ncbi/bui ld3.1/sequence/bwaindex/version0.7.10/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 が記入されます 50

51 本講義の内容 Reseq解析 RNA-seq解析公開データ取得公開データ取得クオリティコントロールマッピング変異検出クオリティコントロールマッピング発現定量 FPKMを算出します 51

V1 ゲノム R e s e q 変異解析 Copyright Amelieff Corporation All Rights Reserved.

V1 ゲノム R e s e q 変異解析 Copyright Amelieff Corporation All Rights Reserved. V1 ゲノム R e s e q 変異解析本講義にあたって代表的な解析の流れを紹介します論文でよく使用されているツールを使用しますコマンドを沢山実行しますスペルミスが心配な方はコマンド例がありますのでコピーして実行してくださいマークのコマンドは実行してください実行が遅れてもあせらずに応用や課題の間に追い付いてください 2 本講義の内容 Reseq解析 RNA-seq解析公開データ取得