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

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

Transcript

1. Quality Check, Processing and Alignment of High-throughput Sequencing Reads Naoki Kimoto (@_kimoton) 【第７回】ゼロから始めるゲノム解析
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（今回）で 扱った範囲