[DADA2] 16S rRNA V3V4에서 최적 primer trim 길이 구하기

2024.03.13. 수

 

개요

박테리아 엠플리콘 분석에서 마커 유전자로 16S rRNA의 V3 V4 영역이 가장 많이 사용된다. 대략적인 길이는 약 ~465 bp이다. 엠플리콘 데이터는 DADA2의 denoising과정을 통해 퀄리티를 필터링하고 paired-end 서열을 합쳐준다. 이때 merge를 위해서는 적어도 12(최적 20nt) nt 이상 겹치는 영역이 존재해야 한다.

 

그러므로 V3-V4 서열을 필터링할 때, 너무 짧지도 혹은 길게 자르지 않도록 주의해야 한다. 그러면 어떤 길이로 잘라야 최적의 퀄리티를 얻을 수 있을까? 이에 대한 실행 코드를 소개한다. 

(참고: "Many loss in V3-V4 filterAndTrim",  https://github.com/benjjneb/dada2/issues/1476)

 

참고글의 borexino는 본인의 데이터에서 최적의 퀄리티를 얻기 위해 각 길이로 잘라서 결과를 비교하는 코드를 제시했다. 만약 본인의 데이터셋이 30을 넘어간다면,  데이터가 많을 때, 랜덤 추출 후 돌리는 것을 추천한다. 

 

 

코드

library(dada2); packageVersion("dada2")

path <-("~/test") # CHANGE ME to the directory containing the fastq files after unzipping.
list.files(path)

# Forward and reverse fastq filenames have format: SAMPLENAME_R1_001.fastq and SAMPLENAME_R2_001.fastq
fnFs <- sort(list.files(path, pattern="_1.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_2.fastq", full.names = TRUE))

# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

# Set a little function
getN <- function(x) sum(getUniques(x))

# Set summary_tab dataset
summary_tab <- data.frame()

# Set my lengths
MyTruncF <- c(240, 260, 280, 290)
MyTruncR <- c(220, 230, 240, 260)

# Benchmarking loops
for (Fwr in MyTruncF) {
  for (Rev in MyTruncR) {
    print("==============================================")
    print(paste0("--->  Testing TruncF:", Fwr, " TruncR:", Rev))
    print("==============================================")
    
    # Filtering and Trimming
    out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
                         truncLen=c(Fwr, Rev), trimLeft=c(20,20),
                         #minLen=200, 
                         maxN=0, truncQ=2, maxEE=c(2,5),
                         rm.phix=TRUE, compress=TRUE, verbose=TRUE, multithread=TRUE)
    
    # Learn errors
    set.seed(101)
    print("Starting learnErrors")
    errF <- learnErrors(filtFs, multithread=TRUE, verbose = TRUE)
    errR <- learnErrors(filtRs, multithread=TRUE, verbose = TRUE)
    
    # Dereplecating
    print("Starting derepFastq")
    derepFs <- derepFastq(filtFs, verbose=TRUE)
    derepRs <- derepFastq(filtRs, verbose=TRUE)
    names(derepFs) <- sample.names    # Name the derep-class objects by the sample names
    names(derepRs) <- sample.names    # Name the derep-class objects by the sample names
    
    # DADA step
    print("Starting dada")
    dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
    dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
    
    # Merging F/R
    print("Starting mergePairs")
    mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE) 
    
    # Sequences 
    seqtab <- makeSequenceTable(mergers)
    
    # Removing kimeras
    print("Starting removeBimeraDenovo")
    seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
    
    # Making a little table
    foo <- data.frame(row.names=sample.names, dada2_input=out[,1],
                      filtered=out[,2], dada_f=sapply(dadaFs, getN),
                      dada_r=sapply(dadaRs, getN), merged=sapply(mergers, getN),
                      nonchim=rowSums(seqtab.nochim),
                      ratioNonChimericChimeric = sum(seqtab.nochim)/sum(seqtab),
                      final_perc_reads_retained=round(rowSums(seqtab.nochim)/out[,1]*100, 1),
                      TruncL = Fwr,
                      TruncR = Rev)
    
    # Storing results
    summary_tab <- rbind(summary_tab, foo) ## 전체 결과 합산
    
  }
}

# Results hard saving
saveRDS(summary_tab, file = "~/test/SummaryStat.rds")

 

 

 

실행 결과 

실제 데이터를 사용하여 수행 후 시각화 해보자. 

summary_tab$set <- paste0("F", summary_tab$TruncL, "_R", summary_tab$TruncR)
summary_tab2 <- aggregate(final_perc_reads_retained ~ set, data = summary_tab, FUN = mean)
order <- summary_tab2 %>% arrange(final_perc_reads_retained) %>% select(set) %>% .$set
summary_tab2$set <- factor(summary_tab2$set, levels = order)

ggplot(summary_tab2, aes(x = set, y = final_perc_reads_retained)) + 
  geom_point()+
  theme(axis.text.x=element_text(angle=45, hjust=1))

 

 

지금 분석하고 있는 데이터는 adapter trimming단계에서 QC30을 설정한 후 dada2 분석을 수행하였다. 이때 각 forward와 reverse의 길이가 240/220일 때 read의 수율이 가장 좋았다.

 

=> 240/220> 250/220> 260/220> 270/220 

=> reverse를 최대한 적게 잡아야 퀄리티가 올라갔다.  

 

 

그러나 240/220은 겹치는 부위가 약 20nt라고 했을 때, 460 -20 = 440 은 465에 비해 많이 모자라다. 그러므로 260/220 혹은 270/210(자료에는 없지만 260/220보다 퀄리티가 좋음)으로 설정하는 것이 좋을 듯하다.