Upgrade to Pro — share decks privately, control downloads, hide ads and more …

【第7回】ゼロから始めるゲノム解析.pdf

B57079e5c504d22b7fda28c5db13bcac?s=47 nkimoto
December 18, 2020

 【第7回】ゼロから始めるゲノム解析.pdf

2020/11/18 (金) 【第7回】ゼロから始めるゲノム解析(R編) 資料
https://bioalgorithms.connpass.com/event/198815/

B57079e5c504d22b7fda28c5db13bcac?s=128

nkimoto

December 18, 2020
Tweet

Transcript

  1. Quality Check, Processing and Alignment of High-throughput Sequencing Reads Naoki

    Kimoto (@_kimoton) 【第7回】ゼロから始めるゲノム解析
  2. 今回扱う範囲 Chapter 6では、アライメントしたリードの情報をどのように扱うかを見てきた。 今回はその前段階としてQCを行いアライメントを行うまでのステップを学習する Preprocessing Quality Check & Trimming Alignment

    (Mapping) Further Analysis Chapter 6 で 扱った範囲 Chapter 8, 9, 10 で 扱う範囲 Chapter 7(今回)で 扱う範囲!
  3. description シーケンス識別子と 配列に関するメタデータ FASTAフォーマット 塩基配列データをテキスト形式で格納するためのフォーマット。 元々FASTAというソフトウェアにおいて使用されていたことに由来している。 1行目:「>」から始まるシーケンス識別子と、その他配列のメタデータ。 2行目以降:塩基配列

  4. FASTQフォーマット FASTAフォーマットの拡張。シーケンサから出力された配列データを格納する際に一般 的に使用されている。 1行目:「@」から始まるシーケンス識別子と、その他配列のメタデータ。 2行目:塩基配列 3行目:「+」文字で始まる。オプションで1行目と同じシーケンス識別子が続く。 4行目:エンコードされたPhredクオリティスコア(2行目の配列に対応) ベースコールの誤り率pを変換したもの

  5. RでQuality Check Rqc パッケージを使用すると、様々なシーケンスクオリティを可視化することができる。ま ずは、fastqファイルをrqc() 関数に与え、シーケンスクオリティに関連する情報を持つオ ブジェクトを取得する。 library(Rqc) # ShortRead

    パッケージからテストデータを読み込み folder = system.file(package="ShortRead", "extdata/E-MTAB-1147") # 対象のfastq.qzファイルをrqc() 関数に与え、シーケンスクオリティに関連する情報を # 持つオブジェクトを取得 qcRes=rqc(path = folder, pattern = ".fastq.gz", openBrowser=FALSE)
  6. RでQuality Check - リードクオリティの確認 読み取りの終わりに向かって品質 が低下していることがわかる(ボック スは四分位範囲)。 良好なサンプルでは、 ベースあたり のQuality

    Scoreの中央値が28(約 1/631のエラー率)を超える。 低クオリティのサンプルでは20 (1/100のエラー率)を下回る。 28 20 # 各ポジションのクオリティスコアを描画 rqcCycleQualityBoxPlot(qcRes)
  7. RでQuality Check - 塩基割合の確認 ランダムな配列のライブラリでは、ヌ クレオチドの偏りがなく、線が互いに ほぼ平行となる。 RNA-Seqでは、ライブラリーの調製中 にプライマー配列が読み取りの開始 位置にアニーリングするため、バイア

    スがかかる。 バイサルファイトシークエンス法で は、メチル化されていないシトシン (C)はチミン(T)に変換されるため、 バイアスがかかる。 # 各ポジションの塩基割合を描画 rqcCycleBaseCallsLinePlot(qcRes)
  8. RでQuality Check - 重複リードの確認 配列フラグメントを増幅する際の PCR アーティファクトによってリードの重複 が起こりうる。 RNA-seqにおいて、リードの重複は、 特定の遺伝子が高レベルで発現して

    いることに起因している可能性があ る。 ⇒ RNA-seqではこの段階で重複した リードを削除すべきではない。 ChIP-seqでは、リードの重複はPCR バイアスが原因である可能性が高 い。 # リードの重複度を描画 rqcReadFrequencyPlot(qcRes)
  9. RでQuality Check - FastQCの利用 library(fastqcr) # FastQC をインストール fastqc_install() #

    FASTQC を実行し、fastqc_resultsフォルダに結果レポートを出力 fastqc(fq.dir = folder,qc.dir = "fastqc_results") fastqcr パッケージを使用すると、Java製のQCツール、FastQCをR経由でインストール ・実行することができる。 FastQCを使うと、ここまでに出力したプロットの他に、過剰に観測されたk-merと過剰なア ダプター配列に関するプロットを出力することができる。
  10. RでQuality Check - FastQCの利用

  11. RでQuality Check - FastQC結果からレポート作成 # Rmarkdown形式にレンダリングしたレポートを表示 qc_report(qc.path="fastqc_results", result.file="reportFile", preview =

    TRUE) gc_report() 関数を使うと、 FastQC結果からRmarkdownベースのレポートを作成可能。
  12. RでQuality Check - FastQC結果からレポート作成

  13. RでQuality Check - FastQC結果の一部を描画 # FastQCの結果を読み込み qc <- qc_read("fastqc_results/ERR127302_1_subset_fastqc.zip" )

    # "Per base sequence quality plot" についての結果をプロット qc_plot(qc, "Per base sequence quality" ) qc_read() 関数を使用してFastQC結果を読み取り、 qc_plot() 関数を使用して関心のある特定のプロットを作成することが可能。
  14. RでQuality Check - FastQC結果の一部を描画

  15. QCを行うことが可能なパッケージ • Rqc (de Souza, Carvalho, and Lopes-Cendes 2018) •

    QuasR (Gaidatzis, Lerch, Hahne, et al. 2015) • systemPipeR (Backman and Girke 2016) • ShortRead (Morgan, Anders, Lawrence, et al. 2009) QCのための様々なRパッケージがbioconductor上で公開されている。 これらはFastQCと同様のレポートを作成するが、その内容はパッケージによって 異なる。
  16. RでPreprocessing クオリティスコアが低いリードは、ベースコールの誤りのため、 正常にアライメントされない可能性が高い。 ⇒ QCの結果を参考に、リードをトリミング・フィルタリングする。 Rでは、QuasR パッケージやShortReadパッケージを使用することでリードを除去 することができる。 QuasR::preprocessReads() 関数では、単一のインターフェースから以下が可能

    • リード長、「N」が含まれている数 • アダプター配列の除去 • 読み取りの開始または終了を事前定義された長さでトリミング
  17. RでPreprocessing - QuasR パッケージ QuasR::preprocessReads() 関数では各フィルタリング条件について、下記のように引数 を指定して実行する。 library(QuasR) # Fastqファイルのパスリストを作成

    fastqFiles <- system.file(package="ShortRead" , "extdata/E-MTAB-1147" , c("ERR127302_1_subset.fastq.gz" , "ERR127302_2_subset.fastq.gz" )) # 前処理後のfastqファイル名を定義 outfiles <- paste(tempfile( pattern=c("processed_1_" , "processed_2_" )), ".fastq",sep="") # 前処理を実行 preprocessReads(fastqFiles, outfiles, nBases=1, 引数 意味 nBases 指定した数よりNが多い リードを除去 truncateEndBases 末端から除去する塩基数 Lpattern 指定した配列から開始 されるリードを除去 minLength 指定した長さより短いリード を除去
  18. RでPreprocessing - ShortRead パッケージ QuasR::preprocessReads() は ShortRead パッケージに依存している。 ShortRead パッケージではより低レベルな関数を使用してフィルタリングが可能

    library(ShortRead) # Fastqファイルのパスリストを作成 fastqFile <- system.file(package="ShortRead", "extdata/E-MTAB-1147" , "ERR127302_1_subset.fastq.gz" ) # fastqファイルを読み込み fq = readFastq(fastqFile) # 塩基ごとのクオリティスコアを matrix形式で取得 qPerBase = as(quality(fq), "matrix") # リードごとに、クオリティスコアが 20以下の塩基数をカウント qcount = rowSums( qPerBase <= 20) # 全ての塩基でクオリティスコアが 20以上のリード数を出力 fq[qcount == 0] ## class: ShortReadQ ## length: 10699 reads; width: 72 cycles
  19. RでPreprocessing - ShortRead パッケージ ShortRead::writeFastq() 関数を使用すると、フィルタリング後のfastqファイルを出 力することができる。 # 全ての塩基でクオリティスコアが 20以上のリード数を出力

    writeFastq(fq[qcount == 0], paste(fastqFile, "Qfiltered", sep="_"))
  20. RでPreprocessing - ShortRead パッケージ fastqファイルは一般的に巨大なファイルである場合が多い(~30GB) ShortRead::FastqStreamer() 関数を使用すると、fastqファイルを分割して「ストリー ミング」することができる。 # ブロックサイズ1000でストリーミングを設定

    f <- FastqStreamer(fastqFile,readerBlockSize=1000) while(length(fq <- yield(f))) { # 塩基ごとのクオリティスコアをmatrix形式で取得 qPerBase = as(quality(fq), "matrix") # リードごとに、クオリティスコアが20以下の塩基数をカウント qcount = rowSums(qPerBase <= 20) # 全ての塩基でクオリティスコアが20以上のリードについて、 # fastqファイルを mode="a" で書き込み writeFastq(fq[qcount == 0], paste(fastqFile, "Qfiltered", sep="_"), mode="a") } yield() 関数を呼び 出すたびに、ファ イルの新しいブ ロックがメモリに 読み込まれる
  21. アライメント QC、Preprocessiongを終えたら、リード配列をリファレンスゲノムにマッピングまたはアラ インメントする。 アライメントでは通常、一致するリードを効率的に探索するためのデータ構造(ゲノムイ ンデックス)を作成する。 ゲノムインデックスの作成方法には、以下の2種類が存在する。 1. ハッシュテーブルを作成 2. Burrows-Wheeler変換を作成

    BWA(H.Li and Durbin 2009)、Bowtie1 / 2(Langmead and Salzberg 2012)、 SOAP(R.Li、Yu、Li、etal. 2009)など
  22. アライメントを行うことが可能なパッケージ • gmapR (Barr, Wu, and Lawrence 2019), • QuasR

    (Gaidatzis, Lerch, Hahne, et al. 2015) • Rsubread (Liao, Smyth, and Shi 2013) • systemPipeR (Backman and Girke 2016) アライメントを行うための様々なRパッケージがbioconductor上で公開されている。
  23. Rでアライメント - QuasRパッケージ qAlign() 関数を使用すると、リファレンスゲノムのファイル、アライメントするリードに関す る情報を記載したテキストファイルの2つの引数を与えるだけでアライメントを実行でき る。 library(QuasR) # copy

    example data to current working directory file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) # fasta形式のゲノムファイル genomeFile <- "extdata/hg19sub.fa" # サンプル名とパスを記載したファイル sampleFile <- "extdata/samples_chip_single.txt" # アライメントを実行 proj <- qAlign(sampleFile, genomeFile) 1)fasta形式または配列情報を格納できる BSgenomeクラスのいずれかに格納したリ ファレンスゲノムのファイル 2)fastqファイルへのファイルパスとサンプ ル名を含むテキストファイル FileName SampleName chip_1_1.fq.bz2 Sample1 chip_2_1.fq.bz2 Sample2
  24. Rでアライメント - QuasRパッケージ # アライメントを実行 proj <- qAlign(sampleFile, genomeFile) Creating

    .fai file for: /home/kimoton/extdata/hg19sub.fa alignment files missing - need to: create alignment index for the genome create 2 genomic alignment(s) ・ ・ Genomic alignments have been created successfully # bamファイルが作成されていることを確認 > dir(pattern="\\.bam$", "extdata/") [1] "chip_1_1_3b0d53a66292.bam" "chip_2_1_3b0d3e161e85.bam" [3] "phiX_paired_withSecondary.bam" ゲノムインデックスが存在しな い場合は自動で生成。
  25. アライメントした後は?? バイサルファイトシーケンシングによって得られたメチル化データの場合、 ⇒ C-> Tミスマッチをカウント! RNA-seqによって得られた遺伝子発現データの場合、 ⇒ 転写産物と重複するリードをカウント! Preprocessing Quality

    Check & Trimming Alignment (Mapping) Further Analysis Chapter 6 で 扱った範囲 Chapter 8, 9, 10 で 扱う範囲 Chapter 7(今回)で 扱った範囲