macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

DADA2

Preprintより

 微生物群集の人間および環境への健康への重要性は、それらの効率的な特徴付けのための方法に動機を与えている。最も一般的で費用効果の高い方法は、標的遺伝子エレメントの増幅および配列決定である。 16S rRNA [ref.1]、ITS領域[ref.2]または18S rRNA [ref.3]のような分類学的マーカー遺伝子のアンプリコンシーケンシングは、コミュニティの国勢調査を提供する。機能的多様性は、機能的遺伝子を標的とすることによって探ることができる[ref.4]。
 アンプリコンシークエンシングデータの生物学的変異による誤差の解消は、アンプリコン特異的エラー訂正法の開発を促す独自の課題を提示している[ref.5,6,7,8]。これらの方法のほとんどは、パイロシーケンシングされたアンプリコン用に設計されており、Illuminaシーケンシングには適用できない。
 現在、イルミナシーケンシングされたアンプリコンデータのエラーは、低品質のリードをフィルタリングし、Operational Taxonomic Units(OTU)を構築することで最も頻繁に対処されている。類似の配列を一括してまとめることで、エラーが生物学的な変化として誤って解釈される割合は減少するが、OTUはストレインレベルの変異を解決する可能性を排除しておりシーケンシングの質を完全には活用していない[ref.7,12,13,14 、15]。最近の研究では、微妙な変化は生態学的なニッチ[12,13]、時間的ダイナミクス[ref.15]、および人口構造[ref.4]について有益であることが示されている。微妙な変異は、いくつかの症例では共生株と病原性を区別し[ref.16、17]、より複雑な微生物関連疾患の臨床関連情報を含む可能性がある[ref.18,19,20]。
 DADA(Divisive Amplicon Denoising Algorithm)は、OTUを構築することなく、パイロシーケンシングアンプリコンエラーを修正するために導入された[ref.7]。 DADAは、454系列のアンプリコンデータの中で最も優れたスケールで真の変異を同定する一方で、偽陽性は少ない[ref.7、4]ことが示された。
 このPreprintでは、Illuminaシークエンシングでの使用に適したDADAの拡張と再実装であり、オープンソースRパッケージとしてhttps://github.com/benjjneb/dada2で入手可能なDADA2を紹介する。 DADA2は、クオリティ情報を組み込んだIlluminaシーケンシングアンプリコンエラーの新しいモデルを実装している。

 

DADA2 マニュアル

f:id:kazumaxneo:20180721120826j:plain

 

チュートリアル

https://benjjneb.github.io/dada2/tutorial.html

 

DADA2に関するツイート

 

インストール

本体 Github

各インストール方法の詳細はこちらを参照

DADA2 Installation

ここではDockerhubに公開されているdada2イメージを使った(Docker Hub)。

docker pull golob/dada2

 

解析の流れ

https://benjjneb.github.io/dada2/tutorial.htmlに従う。データはオーサーらが準備した mothur MiSeq SOP.を使う(リンク先からダウンロードして解凍)。19サンプルある。2x250 V4 シーケンシングデータと記載されている。

 

1、アンプリコンのシーケンシングデータが置いてあるディレクトリに移動し、dockerイメージをスタートさせる。

docker run -itv $PWD:/data/ -w /data golob/dada2

 docker run -itv $PWD:/data/ -w /data/ blekhmanlab/dada2

 

R version 3.5.1 (2018-07-02) -- "Feather Spray"

Copyright (C) 2018 The R Foundation for Statistical Computing

Platform: x86_64-pc-linux-gnu (64-bit)

 

R is free software and comes with ABSOLUTELY NO WARRANTY.

You are welcome to redistribute it under certain conditions.

Type 'license()' or 'licence()' for distribution details.

 

R is a collaborative project with many contributors.

Type 'contributors()' for more information and

'citation()' on how to cite R or R packages in publications.

 

Type 'demo()' for some demos, 'help()' for on-line help, or

'help.start()' for an HTML browser interface to help.

Type 'q()' to quit R.

 

Bioconductor version 3.7 (BiocInstaller 1.30.0), ?biocLite for help

A newer version of Bioconductor is available for this version of R,

  ?BiocUpgrade for help

2、ここからRの中で作業する。dada2とデータの読み込み

library(dada2); packageVersion("dada2")
path <- "/data/MiSeq_SOP/"
#ファイル確認
list.files(path)

> list.files(path)

 [1] "F3D0_S188_L001_R1_001.fastq"   "F3D0_S188_L001_R2_001.fastq"  

 [3] "F3D1_S189_L001_R1_001.fastq"   "F3D1_S189_L001_R2_001.fastq"  

 [5] "F3D141_S207_L001_R1_001.fastq" "F3D141_S207_L001_R2_001.fastq"

 [7] "F3D142_S208_L001_R1_001.fastq" "F3D142_S208_L001_R2_001.fastq"

 [9] "F3D143_S209_L001_R1_001.fastq" "F3D143_S209_L001_R2_001.fastq"

[11] "F3D144_S210_L001_R1_001.fastq" "F3D144_S210_L001_R2_001.fastq"

[13] "F3D145_S211_L001_R1_001.fastq" "F3D145_S211_L001_R2_001.fastq"

[15] "F3D146_S212_L001_R1_001.fastq" "F3D146_S212_L001_R2_001.fastq"

[17] "F3D147_S213_L001_R1_001.fastq" "F3D147_S213_L001_R2_001.fastq"

[19] "F3D148_S214_L001_R1_001.fastq" "F3D148_S214_L001_R2_001.fastq"

[21] "F3D149_S215_L001_R1_001.fastq" "F3D149_S215_L001_R2_001.fastq"

[23] "F3D150_S216_L001_R1_001.fastq" "F3D150_S216_L001_R2_001.fastq"

[25] "F3D2_S190_L001_R1_001.fastq"   "F3D2_S190_L001_R2_001.fastq"  

[27] "F3D3_S191_L001_R1_001.fastq"   "F3D3_S191_L001_R2_001.fastq"  

[29] "F3D5_S193_L001_R1_001.fastq"   "F3D5_S193_L001_R2_001.fastq"  

[31] "F3D6_S194_L001_R1_001.fastq"   "F3D6_S194_L001_R2_001.fastq"  

[33] "F3D7_S195_L001_R1_001.fastq"   "F3D7_S195_L001_R2_001.fastq"  

[35] "F3D8_S196_L001_R1_001.fastq"   "F3D8_S196_L001_R2_001.fastq"  

[37] "F3D9_S197_L001_R1_001.fastq"   "F3D9_S197_L001_R2_001.fastq"  

[39] "filtered"                      "HMP_MOCK.v35.fasta"           

[41] "Mock_S280_L001_R1_001.fastq"   "Mock_S280_L001_R2_001.fastq"  

[43] "mouse.dpw.metadata"            "mouse.time.design"            

[45] "stability.batch"               "stability.files"              

全てペアエンドfastq

 

3、 ファイル名を取り出し、 サンプル1のクオリティ分布を視覚化する。

fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

#R1
postscript("qualtiy1_R1.eps")
plotQualityProfile(fnFs[1:1])
dev.off()

#R2
postscript("qualtiy1_R2.eps")
plotQualityProfile(fnRs[1:1])
dev.off()

f:id:kazumaxneo:20190126231254j:plain
左がR1、右がR2。緑の線が各ポジションのクオリティスコア中央値を表す。R2はR1よりクオリティの落ち込みが激しい。

 

4、全サンプル(1-38)R1クオリティ分布

postscript("qualtiy1-38_R1.eps")
plotQualityProfile(fnFs[1:19]) #R2ならplotQualityProfile(fnRs[1:19])
dev.off()

f:id:kazumaxneo:20190126142406j:plain

 

5、トリミングの実行。オーサーらが提案している通り、R1は3'末端10bpのみトリミングし、R2は急激にクオリティが下がる160付近でトリミングする。注意点として、マージするため十分長くないといけない。本データはMiseqの250x2 (v4 kit)シーケンシングで、オーバーラップも多いため、大胆なトリミングを行っている。

#フィルタリング後のファイル名は~F_filt.fastq.gz、~R_filt.fastq.gzとする
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

#クオリティトリミング、R1は240、R2は160
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE

#リード数確認
head(out)

 filterAndTrim詳細リンク 

 

6、エラーモデル作成と視覚化

#エラー率を学習(本データでは数分かかる)
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

#視覚化
postscript("errorrates.eps")
plotErrors(errF, nominalQ=TRUE)
dev.off()

learnErrors詳細リンク

f:id:kazumaxneo:20190126235831j:plain

各塩基から次の塩基に変わる時のエラー率がクオリティを横軸にとってグラフにプロットされる。黒線は学習したエラー率(観測されたデータ)、赤線が定義上のエラー率変化。今回のデータについては、観測されたエラー率はエラー率として定義された分布(公式ではこれを名目上(または名ばかりの)のクオリティ分布、と記載している)とよくフィッティングしている。すなわち横軸のクオリティが増えるとエラー率も下がる。この傾向は全ての2塩基組み合わせで共通している。

 

7、duplication除去。同一の配列を1つにまとめる。

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
# Name the derep-class objects by the sample names
names(derepFs) <- sample.names
names(derepRs) <- sample.names

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

derepFastq詳細リンク 。454とIontorrentではパラメータ変更が推奨されている(リンク)。

 

8、ペアエンドのマージ。最低12-bpのオーバーラップが必要。

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)

#確認
head(mergers1)

 

9、amplicon sequence variant (ASV) テーブルの作成

seqtab <- makeSequenceTable(mergers)
#行列のサイズ(列数と行数)確認
dim(seqtab)

# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

 

10、キメラ除去。

seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
#行列のサイズ(列数と行数)確認
dim(seqtab.nochim)

sum(seqtab.nochim)/sum(seqtab)

removeBimeraDenovo詳細リンク 。

大半の配列は保持される。仮にこのステップで大半の配列が除去されてしまうようなら、プライマー配列除去が不十分なことに起因している可能性が高い。上流の行程に戻ってプロセスを見直す。

 

11、これまでの作業でのリード数の変化を調べる。

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)

> head(track)

       input filtered denoisedF denoisedR merged nonchim

F3D0    7793     7113      6976      6979   6540    6528

F3D1    5869     5299      5227      5239   5028    5017

F3D141  5958     5463      5331      5357   4986    4863

F3D142  3183     2914      2799      2830   2595    2521

F3D143  3178     2941      2822      2867   2552    2518

F3D144  4827     4312      4151      4228   3627    3488

 

12、系統アサイン。

実行前に管理されているトレーニングデータセットをダウンロードする。このリンクの"Silva version 132"リンク先にある"silva_nr_v132_train_set.fa.gz"をダウンロードする。そのままfastqがあるMiSeq_SOP/下に置く。fungiのシーケンシング解析の場合はトレーニングデータセットを変える(上にリンクを張ったチュートリアルを参照)。

#読み込み
taxa <- assignTaxonomy(seqtab.nochim, "/data/MiSeq_SOP/silva_nr_v132_train_set.fa.gz", multithread=TRUE)
#うまくアサインされないなら、アンチセンスが検索されていない可能性がある。tryRC=TRUEも加える。

taxa.print <- taxa
rownames(taxa.print) <- NULL
head(taxa.print)

##      Kingdom    Phylum          Class         Order          

## [1,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

## [2,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

## [3,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

## [4,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

## [5,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

## [6,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"

##      Family                      Genus         Species       

## [1,] "Bacteroidales_S24-7_group" NA            NA            

## [2,] "Bacteroidales_S24-7_group" NA            NA            

## [3,] "Bacteroidales_S24-7_group" NA            NA            

## [4,] "Bacteroidales_S24-7_group" NA            NA            

## [5,] "Bacteroidaceae"            "Bacteroides" "acidifaciens"

やはり種レベルでのアサインはうまくいかないものが多い。

 

チュートリアルでは、この後phyloseqを読み込んで、アルファ多様性の変化を調べる流れまで丁寧に説明されています。

引用

DADA2: High-resolution sample inference from Illumina amplicon data

Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP

Nat Methods. 2016 Jul;13(7):581-3

 

Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses

Ben J. Callahan, Kris Sankaran, Julia A. Fukuyama, Paul J. McMurdie

https://f1000research.com/articles/5-1492

 

参考 

R3.5のインストール