▶ Dada2란? R을 기반으로 qiime2(리눅스 기반)와 같이 미생물 분석에 사용되는 R 패키지이다.
▶ 이전 글 : [R] Dada2 설치
▶ 튜토리얼 목표 : Dada2로 미생물 데이터 전 처리와 phyloseq을 통한 시각화
▶ 본인의 Dada2 버전에 맞는 튜토리얼을 선택
- Dada2 ver 1.2 tutorial : https://benjjneb.github.io/dada2/tutorial_1_2.html
- Dada2 ver 1.8 tutorial : https://benjjneb.github.io/dada2/tutorial_1_8.html
▶ 참고할 강의 영상 : https://www.youtube.com/watch?v=wV5_z7rR6yw (브라운 대학 워크샵)
▶ 예제 데이터 source : Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. (2013): Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Applied and Environmental Microbiology. 79(17):5112-20.
- 시퀀싱 기계 : illumina MiSeq platform -> paired reads 생성
- 우리가 주목할것 : 박테리아 16s rRNA의 V4영역
🟦 1. 준비 과정
#install.packaged( 패키지_이름)
library(dada2) # sequence processing
library(ShortRead) # FASTQ input and manipulation
# if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#
# BiocManager::install("DECIPHER")
library(DECIPHER) # sequence alignment
#install.packages("phangorn")
library(phangorn) # phylogenetic tree generation
library(ggplot2) # data visualization and analysis
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("phyloseq")
library(phyloseq) # visualization튜토리얼에 사용하는 데이터는 http://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip 여기에서 다운 가능 하다
이를 압축 풀고 Dada2에 파일이 위치한 경로를 입력해 준다
path <- "~/MiSeq_SOP" #directory containing the fastq files after unzipping.
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"# 파일 경로와 파일의 이름을 합해준다
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
fnFs
# "D:~/MiSeq_SOP/F3D0_S188_L001_R1_001.fastq"
# ...
# 파일 이름에서 샘플 이름을 추출한다
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
sample.names
#[1] "F3D0" "F3D1" "F3D141" "F3D142" "F3D143" "F3D144" "F3D145" "F3D146" "F3D147" "F3D148" "F3D149"
#[12] "F3D150" "F3D2" "F3D3" "F3D5" "F3D6" "F3D7" "F3D8" "F3D9" "Mock"
🟦 2. Quality filtering and trimming
이제 시퀀싱 된 read의 퀄리티를 봐야 한다
plotQualityProfile(fnFs[[n]]) // n번째 샘플의 퀄리티 점수를 볼 수 있다
plotQualityProfile(c(fnFs[[1]], fnRs[[1]]))
확연하게 reverse 가 퀄리티가 좋지 않은것이 보인다.
이는 일루미나 시퀀싱 시계에서 양 말단을 읽어낼때 흔하게 발생하는 문제이다.
그렇기 때문에 우리는 forward와 reverse를 겹치면서(overlap) 나쁜 퀄리티의 염기를 조금이나마 보상하는 것이다.
merge전에 나쁜 퀄리티의 염기를 제거하기 위해 이 예제에서는 160번째 염기에서 자르도록 한다!
또한 앞에 10bp정도 퀄리티가 나쁘니 같이 제거해준다.
!! 이때 주의해야 할 점은 V4영역을 보기 때문에 두 염기가 overlap되어야 merge(join)이 가능함으로 실제 데이터는 이러한 부분을 고려해서 잘라야 한다.
# 일단 퀄리티 점수대로 잘라낼 output이 들어갈 경로를 만들어 준다
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
out <- filterAndTrim(fnFs, # forward 서열 경로
filtFs, # forward 서열 경로
fnRs,
filtRs,
truncLen=c(240,160), # 잘라냉 bp길이
maxN=0, # basepair number Cut
maxEE=c(2,2), # reads wirh more EE’s than this number will be discarded
truncQ=2, # 잘라낼 QC점수 기준 (2 = QC 20을 뜻함, 그 이하는 버림)
rm.phix=TRUE, # TRUE라면 phiX genome(박테리로 파지 라이브러리) 과 매치되는 서열은 버림
compress=TRUE, # 결과물이 gzip으로 압축해서 나옴( 속도 빨라지니 좋은것 )
multithread=TRUE) # On Windows set multithread=FALSE
# terminal 환경에서는 "FALSE"로 바꿔주어야 함
head(out)
# reads.in reads.out
#F3D0_S188_L001_R1_001.fastq 7793 7113
#F3D1_S189_L001_R1_001.fastq 5869 5299
#F3D141_S207_L001_R1_001.fastq 5958 5463
#F3D142_S208_L001_R1_001.fastq 3183 2914
#F3D143_S209_L001_R1_001.fastq 3178 2941
#F3D144_S210_L001_R1_001.fastq 4827 4312out파일로 각 샘플에서 몇개의 read가 필터링을 통과하고 남았는지 확인 가능하다
!! 주의할 점은 나중에 merge(dada2::mergePaires)할 때 filter된 forward와 reverse read의 겹치는 부분이 있어야 한다! 다 잘려나가면 큰일난다...그럴땐 maxEE를 조금 완화시켜야 한다.
🟦 3. Error 확률 추정
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)주어진 샘플데이터를 기반으로 에러 모델을 만든다.(설명 봐도 잘 모르겠다...)
결과를 그래프로 확인해보자.
plotErrors(errF, nominalQ=TRUE)
plotErrors(errR, nominalQ=TRUE)
그래프 해석하기
- 검정색 점 : 실제 데이터
- 빨간줄 : 예상되는 QC점수와 error
- 검정 점 : 머신러닝 알고리즘으로 실제 데이터를 기반으로 계산한 에러율
- 각 A2A나 A2C의 의미 : A염기가 C로 바뀔 확률을 말한 것
=> 그래프가 우하향인 것은 퀄리티 점수가 높을수록 그러한 error확률이 적다는 것을 의미함
- 그래서 이 그래프가 의미하는것이 뭐지?
=> 검정 선과 빨간선이 대략 일치한다 -> 예상된 에러율 정도네
-> 데이터가 나쁘지 않군! 하고 다음 단계로 데이터를 신뢰하고 넘어갈 수 있음
🟦 4. 중복 서열 제거 (Dereplicate Reads)
unique한 서열만 남기고 제거하는 단계이다.
필수 단계는 아니지만 이 단계를 시행하는 이유는
중복을 제거해 컴퓨터 계산 시간을 줄이기 위함이며 Dada2의 정확성을 높일 수 있다
derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
# 이름 바꿔주기
names(derepFs) <- sample.names
names(derepRs) <- sample.names
🟦 5. Sample Inference
Dada2에서 가장 주요한 단계로 Unique한 서열을 위해 abundance의 p-value를 계산한다
즉 통계적 분석으로 이 ASV가 시퀀싱 에러인지 아니면 진짜 unique해서 하나의 ASV로 봐도 되는지 계산하는 것이다.
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
dadaFs[[1]]
# dada-class: object describing DADA2 denoising results
# 128 sequence variants were inferred from 1979 input unique sequences.
# Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16
🟦 6. Paired read서열 Merge(=join) 하기
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)default로 20nt가 겹쳐야 merge가 가능하다.
defualt값으로 돌리면 forward와 reverse서열이 정확하게 겹치는 부분이 생겨야 merge가 가능하다.
아니라면 maxMismatch = 1 or 2로 바꾸면 1~2개 틀려도 merge가능 하다
🟦 7. chimera 서열 제거
chimera서열이란 PCR의 인공 생산물(artifact)로써 primer가 우리가 원하는 영역이 아닌 다른 영역에 붙어서 그 서열이 시퀀싱 되었을때 나 등등의 경우에 나타난다
서열을 제거 하기 전에 merge한 데이터를 otu table형식으로 바꾸어 주자
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# [1] 20 29320개의 샘플이 행, 293개의 ASV가 열로 표시된것을 볼 수 있다
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim) # [1] 20 232
sum(seqtab.nochim)/sum(seqtab) # 0.964064293개 중 61개의 ASV가 키메라 서열로 인식되어 제거됨
이제까지 단계를 거친 파일들을 요약해서 보고자 한다면 아래 script를 참고
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
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 2868 2553 2519
#F3D144 4827 4312 4151 4228 3646 3507
#F3D145 7377 6741 6592 6627 6079 5820
#F3D146 5021 4560 4450 4470 3968 3879
#F3D147 17070 15637 15433 15506 14233 13006
#F3D148 12405 11413 11250 11268 10528 9935
#F3D149 13083 12017 11857 11898 11154 10653
#F3D150 5509 5032 4880 4925 4349 4240
#F3D2 19620 18075 17907 17939 17431 16835
#F3D3 6758 6250 6145 6176 5850 5486
#F3D5 4448 4052 3930 3994 3716 3716
#F3D6 7989 7369 7231 7294 6865 6678
#F3D7 5129 4765 4646 4671 4426 4215
#F3D8 5294 4871 4786 4802 4576 4547
#F3D9 7070 6504 6341 6442 6092 6015
#Mock 4779 4314 4287 4286 4269 4269
만약 대량의 서열이 키메라로 인식되어 삭제 되었을때 주의해야 할 점은 프라이머가 잘 제거되었는지이다.
아니라면 다시 filterAndTrim()으로 제거해 주어야 하고
만약 그래도 고쳐지지 않으면
1. 더 낮은 퀄리티에서 서열 trimming하기
2. 자신이 보고자 하는 영역이 어디인지 보기
- 만약 V4-V5사이를 보고자 한다면 각 영역 사이에 constant한 부분이 chimera형성의 주요 타겟이다.
- 그러므로 V4만 보는 것이 V4V5보다 키메라 형성이 더 적다
3. 연구하고자 하는 환경이 복잡한지
- 다양한 생물이 살 고 있다면 여러 다양한 시퀀스가 에러라고 인식될 것이다
- 이는 불가항력이다ㅠ
다음글에선 silva database를 이용해 taxonomy classification과 phyloseq을 통한 시각화 튜토리얼을 정리해보았다.
