PowerPoint プレゼンテーション

Size: px

Start display at page:

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

みいかやまがた
5 years ago
Views:

1 平成 28 年度 NGS ハンズオン講習会 Reseq 解析 2016 年 7 月 26 日

3 本講義の内容前半パート ( 講義 ) 後半パート ( 実習 ) Reseqとは検出可能な遺伝子変異解析パイプラインとは公開データの取得と利用クオリティコントロールとはマッピングとは変異検出とはアノテーションとはより高精度な分析のために後半パート ( 実習 ) で行うこと公開データの確認クオリティコントロールマッピング変異検出アノテーション解析結果の可視化まとめ最後に Amelieff Corporation All Rights Reserved 3

5 Reseq とは whole genome sequence whole exome sequence amplicon sequence target sequence Reseq DNA の変異検出を目的としたワークフローの総称 Reseq 1 公開データ取得 2 クオリティコントロール 3 マッピング 4 変異検出 SNVとIndelの検出を行います RNA-seq ( 明日実施 ) 1 公開データ取得 2 クオリティコントロール 3 マッピング 4 発現定量 FPKMを算出します Amelieff Corporation All Rights Reserved 5

reference insertion deletion SV (Structural Variation)

6 検出可能な遺伝子変異ショートリードのシーケンスでも様々な変異を検出可能 SNV (Single Nucleotide Variant) InDel (Insertion & Deletion) reference insertion deletion SV (Structural Variation) inversion translocation duplication Amelieff Corporation All Rights Reserved 6

7 検出可能な遺伝子変異各変異による影響例 SNV (Single Nucleotide Variant) InDel (Insertion & Deletion) Thr Tyr Thr (STOP) GAA ( グルタミン ) GTA ( バリン ) Chr. 9 Chr. 22 SV (Structural Variation) c-abl bcr translocation Philadelphia chromosome 慢性骨髄性白血病で見られる Amelieff Corporation All Rights Reserved 7

8 解析パイプラインとはあるソフトの出力結果が次のソフトの入力ファイルとなる連続した解析処理の流れ FASTQ BAM FastQC Trimmomatic 1. クオリティコントロール GATK 4. SNV / InDel 検出およびフィルタリング BWA 2. マッピング SnpEff 5. アノテーション GATK BAM 3. アライメントおよびベースクオリティのリキャリブレーション VCF Amelieff Corporation All Rights Reserved 8

公開データの取得と利用今回の解析に必要なデータ ( ダウンロード済み ) リファレンスゲノム http://support.illumina.com/sequencing/sequencing_software/ igenome.html 解析対象のシーケンスデータ ftp://ftp.ddbj.nig.ac.

10 公開データの取得と利用今回の解析に必要なデータ ( ダウンロード済み ) リファレンスゲノム igenome.html 解析対象のシーケンスデータ ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/era /ERX015989/ERR038793_1.fastq.bz2 ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/era /ERX015989/ERR038793_2.fastq.bz2 Amelieff Corporation All Rights Reserved 10

11 公開データの取得と利用酵母のリファレンスゲノムデータの取得方法 $ wget ftp://igenome:g3nom3s4u@ussdftp.illumina.com/saccharomyces_cerevisiae/ncbi/build3.1/saccha romyces_cerevisiae_ncbi_build3.1.tar.gz $ tar zxvf Saccharomyces_cerevisiae_NCBI_build3.1.tar.gz llumina の Web ページからリファレンスゲノムを取得し解凍 ( 実行済み ) $ ls -l /home/iu/genome/saccer3 : -rwxrwxr-x 1 iu iu 月 4 19:50 genome.fa -rwxrwxr-x 1 iu iu 14 7 月 4 19:50 genome.fa.amb -rwxrwxr-x 1 iu iu 月 4 19:50 genome.fa.ann : /home/iu/genome/saccer3 に配置してあります Amelieff Corporation All Rights Reserved 11

公開データの取得と利用解析対象のシーケンスデータの取得方法 1 ( 実行済み ) http://trace.ddbj.nig.ac.jp/dra/index.

公開データの取得と利用解析対象のシーケンスデータの取得方法 3 ( 実行済み ) 実験詳細を確認

公開データの取得と利用解析対象のシーケンスデータの取得方法 4 ( 実行済み )

16 公開データの取得と利用解析対象のシーケンスデータの取得方法 5 ( 実行済み ) CUI でダウンロード $ cd /home/iu/reseq/data $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/era /ERX015989/ERR038793_1.fastq.bz2 $ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/era038/era /ERX015989/ERR038793_2.fastq.bz2 /home/iu/reseq/data にダウンロードしてあります Amelieff Corporation All Rights Reserved 16

17 公開データの取得と利用解析対象のシーケンスデータの取得方法 6 ( 実行済み ) 圧縮ファイルの解凍 $ ls ERR038793_1.fastq.bz2 ERR038793_2.fastq.bz2 今回用いるデータは bz2 形式で圧縮されていました $ bzip2 -d ERR038793_1.fastq.bz2 $ bzip2 -d ERR038793_2.fastq.bz2 /home/iu/reseq/data には解凍済のファイルが配置してありますソフトウェアによっては圧縮されたままのファイルを扱えるものもあります Amelieff Corporation All Rights Reserved 17

read quality クオリティコントロールとは公開データをそのまま使うのは危険公開データ測定環境の違いシーケンス結果のクオリティアダプター配列の有無タグの有無

19 クオリティコントロールとはゲノム解析で用いられる主なクオリティコントロールツール FastQC クオリティチェック用ソフトウェア FASTX-toolkit C で書かれた多機能クオリティコントロールツール PRINSEQ Perl で書かれた多機能クオリティコントロールツール Trimmomatic Java で書かれたトリミングツール etc... Amelieff Corporation All Rights Reserved 19

マッピングとは各リードはリファレンスゲノムのどこに位置するか調べる sequence read reference genome Reseq 解析は

21 マッピングとはゲノム解析で用いられる主なマッピングツール BWA Indel に強いギャップ許容型のマッピングツール Bowtie2 ショートリード用のマッピングツール SOAP2 大量ショートリード高速マッピングツール Indel 不可 Novoalign NovoCraft 社の製品ギャップ許容型のマッピングツール etc... Amelieff Corporation All Rights Reserved 21

23 変異検出とは変異検出の前に 1 ~Realignment~ リアライメントとは? 1 本のリードに複数の変異が含まれる場合にアライメントスコアの計算上 SNV や Indel の正確な位置を決定できないことがありますこのような領域を対象領域として抜き出して改めて丁寧にアライメントを行うことで変異検出の信頼性を高めることが出来ます Amelieff Corporation All Rights Reserved 23

24 報告されているクオリティ値との差変異検出とは変異検出の前に 2 ~Base Recalibration~ 10 ベースリキャリブレーションとは? 変異検出のアルゴリズムはクオリティスコアに大きく影響されますこの行程では既知の SNP 情報を用いて測定環境により異なるクオリティスコアをノーマライズすることで測定環境に依存しない変異検出が可能となります AA AG AC AT GC Amelieff Corporation All Rights Reserved 24

25 変異検出とはゲノム解析で用いられる変異検出ツール A) GATK HaplotypeCaller GATK UnifiedGenotyper VarScan etc... B) A) Charles D. Warden et al., Peer J, 2014 B) Hasan et al., Human Genomics, 2015 Amelieff Corporation All Rights Reserved 25

変異検出とは最新版 GATK Realignmentは GATK HaplotypeCallerや MuTect2を使用すれば必要ないということで GATK3.

29 後半パート ( 実習 ) で行うこと本日実際に行う解析フロー FastQCによるクオリティチェック Trimmomaticによるクオリティコントロール BWA memによるマッピング GATK HaplotypeCallerによる変異検出 GATK VariantFiltrationによる変異のフィルタリング snpeffによるアノテーション IGVによる解析結果の可視化 Amelieff Corporation All Rights Reserved 29

32 公開データの確認 fasta ファイルの中身を見てみる $ less /home/iu/genome/saccer3/genome.fa >chri CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACC CACACACACACATCCTAACACTACCCTAACACAGCCCTAATCTAACCCTG GCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTAC CCTGTCCCATTCAACCATACCACTCCGAACCACCATCCATCCCTCTACTT ACTACCACTCACCCACCGTTACCCTCCAATTACCCATATCCAACCCACTG : 1 行目 : コンティグ名 2 行目以降 : 実際の配列情報 q で閲覧を終了 Amelieff Corporation All Rights Reserved 32

33 公開データの確認解析対象の fastq ファイルを確認 $ less 1 length=100 GGACAAGGTTACTTCCTAGATGCTATATGTCCCTACGGCCTTGTCTAACACCATCCAGCATGCA ATAAGGTGACATAGATATACCCACACACCACACCCT +ERR length=100 D/DDBD@B>DFFEEEEEEEEF@FDEEEBEDBBDDD:AEEE<>CB?FCFF@F?FBFF@?:EEE:E EBEEEB=EEE.>>?=AD=8CDFFFFFEFEF@C?;DC fastqファイルを見てみる 1 行目配列 IDと付加情報 2 行目 : 塩基配列 3 行目 : + 配列 IDと付加情報 4 行目 : クオリティ fastqファイルは1リードあたり4 行で表記されます Amelieff Corporation All Rights Reserved 33

34 公開データの確認解析対象データのリード数を確認 $ wc -l data/err038793_1.fastq data/err038793_1.fastq 行あるのでリード数は / 4 = となる $ wc -l data/err038793_2.fastq data/err038793_2.fastq ペアエンドなのでERR038793_2.fastqはもちろん同じリード数を持つ Amelieff Corporation All Rights Reserved 34

35 クオリティコントロールシーケンスクオリティチェックソフトウェア FastQC の紹介 $ fastqc -v FastQC v バージョンを確認 (2016 年 7 月現在最新版は v0.11.5) $ fastqc -h SYNOPSIS FastQC - A high throughput sequence QC analysis tool fastqc seqfile1 seqfile2.. seqfilen fastqc [-o output dir] [--(no)extract] [-f fastq bam sam] [-c contaminant file] seqfile1.. seqfilen :.fastq 以外に.sam や.bam も可能複数ファイルの指定も可能 Amelieff Corporation All Rights Reserved 35

36 クオリティコントロールシーケンスクオリティチェックソフトウェア FastQC の実行 $ mkdir fastqc_before $ fastqc -o fastqc_before -f fastq data/err038793_1.fastq data/err038793_2.fastq $ ls fastqc_before ERR038793_1_fastqc ERR038793_1_fastqc.zip ERR038793_2_fastqc ERR038793_2_fastqc.zip 解析結果の html ファイルができているのでこれをブラウザ (firefox) で確認してみます $ firefox fastqc_before/err038793_1_fastqc/fastqc_report.html fastqc_before/err038793_2_fastqc/fastqc_report.html ブラウザでタブが 2 つ開かれクオリティチェックの解析結果が確認できます Amelieff Corporation All Rights Reserved 36

クオリティコントロール FastQC の結果確認 1 問題なし注意 (warning) 問題あり (failure) Basic Statistics ファイルの基本的な情報

38 クオリティコントロール FastQC の結果確認 2 Per Base Sequence Quality 横軸はリード長縦軸は quality value を表すリードの位置における全体のクオリティの中央値や平均を確認できる赤線は中央値青線は平均値黄色のボックスは 25%~75% の領域を表す上下に伸びた黒いバーが 10%~90% の領域を意味する Amelieff Corporation All Rights Reserved 38

クオリティコントロール FastQC の結果確認 3 Per Sequence Quality Scores 縦軸がリード数横軸が

クオリティコントロール FastQC の結果確認 4 Per Base Sequence Content リードにおける位置での各塩基の割合を示すいずれかの位置で A と T の割合の差

クオリティコントロール FastQC の結果確認 5 Per Base GC Content リードにおける位置での GC 含量を表すいずれかの位置で全体での GC

42 クオリティコントロール FastQC の結果確認 6 Per sequence GC content リードの GC 含量の分布が示されているリードに含まれる GC 含量は一般に正規分布に従うとされている正規分布と比較しその残差が 15% 以上ならば Warning とされるまた 30% 以上ならば Failure とされる Amelieff Corporation All Rights Reserved 42

クオリティコントロール FastQC の結果確認 7 Per Base N Content N はシーケンサーの問題で ATGC いずれの塩基にも決定出来なかった場合に記述される

クオリティコントロール FastQC の結果確認 8 Sequence Length Distribution リード長の全体の分布

45 クオリティコントロール FastQC の結果確認 9 Sequence Duplication Levels リードの重複レベルを見ている 1~10 はそれぞれ重複のレベルで全体の 20% 以上がユニークでないものだと warning, 50% 以上がユニークでないと failure となる Overrepresented Sequences 重複している配列とその割合を表す特定の配列が全リードの 0.1% を超えると warning 1% を超えると failure となる Amelieff Corporation All Rights Reserved 45

クオリティコントロール FastQC の結果確認 10 K-mer Content K bp の任意の配列 (K-mer) を考えた時ライブラリに含まれる ATGC の割合を元に実際に観測された値 / 理論的に観測される期待値を計算している ( デフォルトは K=7) それぞれの任意の配列について実測が期待値を大きく上回っている時

46 クオリティコントロール FastQC の結果確認 10 K-mer Content K bp の任意の配列 (K-mer) を考えた時ライブラリに含まれる ATGC の割合を元に実際に観測された値 / 理論的に観測される期待値を計算している ( デフォルトは K=7) それぞれの任意の配列について実測が期待値を大きく上回っている時それはライブラリに配列的な偏りがあると解釈される実測値 / 期待値はリード長全体における計算とリードのある位置での計算を行い全体における値が 3 倍リードのある位置における値が 5 倍になると warning リードのある位置における値が 10 倍になると failure となる Amelieff Corporation All Rights Reserved 46

シーケンス技術が向上しクオリティの高いデータを目にする機会が増えましたが試料シーケンス

47 クオリティコントロール補足 ) クオリティの悪いデータ最初の 1 塩基の割合が不自然リード末端でクオリティが低下マッピング率の低下や変異の偽陽性が増加するなどの問題を引き起こすシーケンス技術が向上しクオリティの高いデータを目にする機会が増えましたが試料シーケンストリミングなどに問題がないか確認することをおすすめします Amelieff Corporation All Rights Reserved 47

48 クオリティコントロールクオリティを向上させるために (amelieff の場合 ) FASTQ 形式にマッチするかチェックデータクオリティチェック (FastQC) Illumina CASAVA filter [Y] を除去クオリティ20 未満が80% 以上のリードを除去クオリティ20 未満の末端をトリム未知の塩基 (N) が多いリード除去配列長が短いリード除去片側のみのリードを除外データクオリティチェック (FastQC) 様々な流儀が存在するが大切なのは処理の内容と処理の順番例えばロングリードの場合リードの大半が除外されてしまう可能性例えばペアエンドリードの場合ペアが揃っていないとマッピングソフトが停止する可能性 Amelieff Corporation All Rights Reserved 48

クオリティコントロール今回のデータに対する処理 (Trimmomatic を用いた一括処理 ) $ mkdir trimmed_data $ java -jar /usr/local/bin/trimmomatic-0.36.jar PE -threads 2 -phred33 -trimlog trimmed_data/log.txt data/err038793_1.

49 クオリティコントロール今回のデータに対する処理 (Trimmomatic を用いた一括処理 ) $ mkdir trimmed_data $ java -jar /usr/local/bin/trimmomatic-0.36.jar PE -threads 2 -phred33 -trimlog trimmed_data/log.txt data/err038793_1.fastq data/err038793_2.fastq trimmed_data/paired_output_err038793_1.fastq trimmed_data/unpaired_output_err038793_1.fastq trimmed_data/paired_output_err038793_2.fastq trimmed_data/unpaired_output_err038793_2.fastq SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:80 今回解析するデータは FastQC によるクオリティチェックの結果問題ありと判断された項目はありませんでしたそのため今回はリード末端のクオリティが悪い部分をトリムすることでクオリティの底上げを図ります Amelieff Corporation All Rights Reserved 49

50 クオリティコントロール Trimmomatic のオプション PE: ペアエンドの指定 -threds: 使用するスレッド数 -phred33: クオリティスコアの計算方法 -trimlog: logファイルの名前指定 SLIDINGWINDOW: ウィンドウサイズと平均クオリティの設定 LEADING: リードの先頭からトリム位置を探した時の下限クオリティ値 TRAILING: リードの末端からトリム位置を探した時の下限クオリティ値 MINLEN: リードそのものを除去する設定値 Amelieff Corporation All Rights Reserved 50

51 クオリティコントロール QC 後の結果確認 $ mkdir fastqc_after $ fastqc -o fastqc_after -f fastq trimmed_data/paired_output_err038793_1.fastq trimmed_data/paired_output_err038793_2.fastq $ firefox fastqc_after/paired_output_err038793_1_fastqc/ fastqc_report.html fastqc_after/paired_output_err038793_2_fastqc/ fastqc_report.html Amelieff Corporation All Rights Reserved 51

52 クオリティコントロール Trimmomatic によるクオリティコントロールの結果データクオリティチェック (FastQC) クオリティ20 未満のリード末端をトリム配列長が短いリード除去片側のみのリードを除外データクオリティチェック (FastQC) 解析するデータにアダプター配列が含まれている場合 Trimmomatic を用いてアダプター配列を除去することも出来るリード末端のクオリティが悪かった部分がトリムされました Amelieff Corporation All Rights Reserved 52

53 マッピング BWA mem によるマッピング準備 $ bwa mem -help Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq] : インデックスファイルと対象の fastq ファイルが要求されている例 ) リファレンスゲノムの FASTA ファイルに対するインデックスファイル search!! BWA インデックスファイル FASTA ファイル direct!! 参照したい配列 Amelieff Corporation All Rights Reserved 53

54 マッピング BWA memのインデックスファイル作成 $ mkdir -p Scerevisiae/BWAIndex $ cd Scerevisiae/BWAIndex $ ln -s /home/iu/genome/saccer3/genome.fa $ bwa index genome.fa $ ls genome.fa genome.fa.ann genome.fa.pac genome.fa.amb genome.fa.bwt genome.fa.sa インデックスファイルを格納するフォルダを作成し移動リファレンスゲノムのシンボリックリンクを作成しそれを用いてインデックスファイルを作成する Amelieff Corporation All Rights Reserved 54

55 マッピング BWA mem によるマッピング ( ここからは 1 ファイルで行います ) $ cd /home/iu/reseq $ mkdir mapping $ bwa mem -M -R '@RG tid:err038793_1 tsm:err tpl:illumina' Scerevisiae/BWAIndex/genome.fa trimmed_data/paired_output_err038793_1.fastq > mapping/err038793_mapped.sam -M: SAM/BAM ファイルの FLAG 列を他のソフトウェアに互換性のあるものに変更する -R: RG(read groups) の情報を付与する複数のサンプル情報が混在している場合に有用 GATK で BAM ファイルを扱うには platform (PL) および sample (SM) が必須 PL の例 :454, LS454, Illumina, Solid, ABI_Solid Amelieff Corporation All Rights Reserved 55

56 マッピング SAM ファイルを BAM ファイルに変換 $ samtools view -b mapping/err038793_mapped.sam > mapping/err038793_mapped.bam $ ls -lh mapping total 211M -rw-rw-r-- 1 iu iu 41M 7 月 13 17:14 ERR038793_mapped.bam -rw-rw-r-- 1 iu iu 171M 7 月 13 17:13 ERR038793_mapped.sam 171M の SAM ファイルが 41M のバイナリファイルに変換された Amelieff Corporation All Rights Reserved 56

57 マッピング BAM ファイルをソートインデックス作成 $ samtools sort mapping/err038793_mapped.bam -o mapping/err038793_sorted.bam $ samtools index mapping/err038793_sorted.bam $ ls mapping ERR038793_mapped.bam ERR038793_sorted.bam ERR038793_mapped.sam ERR038793_sorted.bam.bai BAM ファイルを高速に扱うためのインデックスファイルを作成 Amelieff Corporation All Rights Reserved 57

58 マッピングマッピングの結果を確認する $ samtools idxstats mapping/err038793_sorted.bam chri chrii chriii chriv chrix chrm : 1 列目 : コンティグ名 (fasta ファイルのヘッダ ) 2 列目 : コンティグの長さ 3 列目 : マッピングされたリード数 4 列目 : マッピングされなかったリード数 3 列目の総和 4 列目の総和マッピングされたリードの総数マッピングされなかったリードの総数 Amelieff Corporation All Rights Reserved 58

59 マッピングマッピングの結果を計算確認する 2 $ wc -l trimmed_data/paired_output_err038793_1.fastq awk '{print $1/4;}' fastq ファイルは 4 行 1 リードなので fastq ファイルの行数を 4 で割った数がリード数です $ samtools idxstats mapping/err038793_sorted.bam > multi.txt $ awk '{sum += $3} END {print sum}' multi.txt $ awk '{sum2 += $4} END {print sum2}' multi.txt Amelieff Corporation All Rights Reserved 59

61 マッピングマルチヒットしたリードを除きユニークリードのみにする $ samtools view -b -F 256 mapping/err038793_sorted.bam > mapping/err038793_unique.bam view : sam/bam を扱うサブコマンド -b : 出力を BAM ファイルにする -F : 指定されたフラグが付与されたリードを除外するマッピング状況を確認する $ samtools index mapping/err038793_unique.bam $ samtools idxstats mapping/err038793_unique.bam > unique.txt index : BAM ファイルのインデックスファイルを作成する idxstats : インデックスファイルのステータスを表示する Amelieff Corporation All Rights Reserved 61

62 マッピングマッピングの結果を再計算する $ awk '{sum += $3} END {print sum}' unique.txt $ awk '{sum2 += $4} END {print sum2}' unique.txt 全リード数 : マッピングされたリード数 : マッピングされなかったリード数 : 計 Amelieff Corporation All Rights Reserved 62

63 変異検出 GATK HaplotypeCaller による変異検出 $ mkdir variant_call $ java -jar /usr/local/bin/genomeanalysistk.jar -R /home/iu/genome/saccer3/genome.fa -T HaplotypeCaller -I mapping/err038793_sorted.bam -variant_index_type LINEAR -variant_index_parameter o variant_call/err vcf $ ls variant_call ERR vcf ERR vcf.idx VCF (Variant Call Format) が作成されました Amelieff Corporation All Rights Reserved 63

64 変異検出 GATK HaplotypeCaller で設定したオプション -R: リファレンスゲノムの場所 -T: 使用するアルゴリズム -I: 入力データ -variant_index_type: LINEAR で等間隔のインデックスを作成する -variant_index_parameter: インデックスの bin 幅 -o: 出力ファイル Amelieff Corporation All Rights Reserved 64

65 変異検出 VCF ファイルの確認 $ less variant_call/err vcf : #CHROM POS ID REF ALT QUAL ERR chri 111. C T /1:3,6:9:90:220,0,90 chri 136. G A /1:0,9:9:29:371,29,0 : CHROM : 染色体番号 POS : 変異箇所の1 塩基目の位置 ID : ID 情報 ( 情報がないので. と記載 ) REF : リファレンスゲノムの塩基配列 ALT : 変異のある塩基配列 QUAL : phred-scaleによるクオリティ値 FILTER : フィルタリング条件 ( 情報がないので. と記載 ) INFO : 変異情報 Amelieff Corporation All Rights Reserved 65

66 変異検出 VCFファイルの確認 FORMAT GT AD DP GQ PL genotype allelic depth ERR ( ファイル名 ) C/T のヘテロ REF のカバレージは 3 ALT のカバレージは 6 read depth 0/1:3,6:9:90:220,0,90 リード数は計 9 つ genotype quality phred-scaled genotype likelihoods C/T の最尤推定値が最も高い 10^(-PL/10) % の信頼性 #CHROM POS ID REF ALT chri 111. C T ひとつ目の SNP を例に Amelieff Corporation All Rights Reserved 66

変異検出検出した SNV INDEL の数を確認 $ grep -c -v '^#'

vcf 57869 57869 個の変異が検出されましたがこの中にはカバレージが低く

68 変異検出 Low coverage な SNV のカウント $ awk '{print $10;}' variant_call/err vcf grep '0/1' cut -d ':' -f 3 awk '{if($1 < 10){print $1;}}' wc -l 667 カバレージが 10 未満の信頼性の低い変異が 667 個存在しています 1: SAMPLE 列の抜きだし 2: SNV のみにフォーカス : SAMPLE 列の : 区切り 3 つめの要素の DP(coverage) を抜き出す 4: DP が 10 未満のもののみ出力する 5: 出力された行数を数える Amelieff Corporation All Rights Reserved 68

69 変異検出変異のフィルタリング (FILTER 列の活用 ) $ java -jar /usr/local/bin/genomeanalysistk.jar -R /home/iu/genome/saccer3/genome.fa -T VariantFiltration -V variant_call/err vcf -o variant_call/err038793_fil.vcf --clusterwindowsize 10 --filterexpression 'DP < 10' --filtername 'LowCoverage' -R: リファレンスゲノムの場所 -V: 入力 VCFファイル -o: 出力ファイル --filterexpression : フィルタリング条件 --filtername : フィルター名 Amelieff Corporation All Rights Reserved 69

70 変異検出変異のフィルタリング (FILTER 列の活用 ) $ less variant_call/err038793_fil.vcf : #CHROM POS ID REF ALT QUAL FILTER chri 111. C T LowCoverage chri 136. G A LowCoverage : カバレージが 10 以下の SNP を消すわけでなく LowCoverage というダグを付くことで後ほどフィルタリングできるようなっています Amelieff Corporation All Rights Reserved 70

71 アノテーション snpeff を用いたアノテーション方法 $ mkdir annotation $ cd annotation $ java -jar /usr/local/bin/snpeff.jar eff -c /usr/local/bin/snpeff.config -i vcf -o vcf R /variant_call/ERR038793_fil.vcf 1> ERR038793_eff.vcf $ less ERR038793_eff.vcf eff: 出力フォーマットの指定 -c: コンフィグファイルの場所様々なデータベースの情報が記載されている -i, -v: 入出力ファイルフォーマット R : Scerevisiae のデータベース SacCer3 に対応 Amelieff Corporation All Rights Reserved 71

アノテーション snpeff を用いたアノテーション方法 : #CHROM POS ID REF ALT QUAL FILTER INFO chri 111. C T 191.77 LowCoverage chri 136. G A 342.

72 アノテーション snpeff を用いたアノテーション方法 : #CHROM POS ID REF ALT QUAL FILTER INFO chri 111. C T LowCoverage chri 136. G A LowCoverage : 遺伝子名やコーディング情報翻訳後のタンパク質に与えるインパクト等の情報が付与される低 LOW MODIFIER MODERATE HIGH 高 Amelieff Corporation All Rights Reserved 72

73 アノテーション snpeff を用いたアノテーションエラーの回避 $ grep 'chrm' ERR038793_eff.vcf : ERROR_CHROMOSOME_NOT_FOUND : 今回作成した VCF ファイルではミトコンドリア DNA を chrm と記述していますしかしながら今回用いた snpeff のデータベース R ではミトコンドリアの DNA 情報が chrmito と記載されているために正しくマッチングが行われずエラーが起きています Amelieff Corporation All Rights Reserved 73

74 アノテーション snpeff を用いたアノテーションエラーの回避 $ sed -e 's/chrm/chrmito/g'../variant_call/err038793_fil.vcf >../variant_call/err038793_fil2.vcf $ java -jar /usr/local/bin/snpeff.jar eff -c /usr/local/bin/snpeff.config -i vcf R o vcf../variant_call/err038793_fil2.vcf 1> ERR038793_eff2.vcf $ grep 'chrm' ERR038793_eff2.vcf ミトコンドリア DNA もアノテーションされました sed コマンド : 文字列の全置換から行単位の抽出削除まで行えるテキスト加工コマンド Amelieff Corporation All Rights Reserved 74

75 解析結果の可視化 Integrative Genomics Viewer (IGV) を用いた解析結果の確認 1 $ igv.sh IGV を起動し Genomes タブから Load Genomes from File... を選択 /home/iu/ genome/ saccer3 の下にある genome.fa を選択し開く Amelieff Corporation All Rights Reserved 75

解析結果の可視化 Integrative Genomics Viewer (IGV) を用いた解析結果の確認 2 File タブから Load from File.

解析結果の可視化 Integrative Genomics Viewer (IGV) を用いた解析結果の確認 3 サーチウィンドウに

解析結果の可視化 Integrative Genomics Viewer (IGV) を用いた解析結果の確認 4 chrii:804,080-804,116 と入力し

79 まとめ本日行った解析のおさらい FastQCによるクオリティチェック Trimmomaticによるクオリティコントロール BWA memによるマッピング GATK HaplotypeCallerによる変異検出 GATK VariantFiltrationによる変異のフィルタリング snpeffによるアノテーション IGVによる解析結果の可視化 Amelieff Corporation All Rights Reserved 79

最後により高精度な解析を行うには本日行ってきたのはあくまで解析方法の一例ですツールの選択に正解はありません自身のデータに適したツールを選択しより良い解析手順を確立していってください

PowerPoint プレゼンテーション

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