Tutoriel DADA2

library(dada2); packageVersion("dada2")

[1] ‘1.28.0’

path <- "~/MiSeq_SOP" # CHANGE ME to the 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] "HMP_MOCK.v35.fasta"            "Mock_S280_L001_R1_001.fastq"  
[41] "Mock_S280_L001_R2_001.fastq"   "mouse.dpw.metadata"           
[43] "mouse.time.design"             "stability.batch"              
[45] "stability.files"              

# Forward and reverse fastq filenames have format: SAMPLENAME_R1_001.fastq and SAMPLENAME_R2_001.fastq
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

plotQualityProfile(fnFs[1:2])

plotQualityProfile(fnRs[1:2])

# Place filtered files in filtered/ subdirectory
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names

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

Creating output directory: /home/rstudio/MiSeq_SOP/filtered

head(out)

                              reads.in
F3D0_S188_L001_R1_001.fastq       7793
F3D1_S189_L001_R1_001.fastq       5869
F3D141_S207_L001_R1_001.fastq     5958
F3D142_S208_L001_R1_001.fastq     3183
F3D143_S209_L001_R1_001.fastq     3178
F3D144_S210_L001_R1_001.fastq     4827
                              reads.out
F3D0_S188_L001_R1_001.fastq        7113
F3D1_S189_L001_R1_001.fastq        5299
F3D141_S207_L001_R1_001.fastq      5463
F3D142_S208_L001_R1_001.fastq      2914
F3D143_S209_L001_R1_001.fastq      2941
F3D144_S210_L001_R1_001.fastq      4312

errF <- learnErrors(filtFs, multithread=TRUE)

33514080 total bases in 139642 reads from 20 samples will be used for learning the error rates.

errR <- learnErrors(filtRs, multithread=TRUE)

22342720 total bases in 139642 reads from 20 samples will be used for learning the error rates.

plotErrors(errF, nominalQ=TRUE)

dadaFs <- dada(filtFs, err=errF, multithread=TRUE)

Sample 1 - 7113 reads in 1979 unique sequences.
Sample 2 - 5299 reads in 1639 unique sequences.
Sample 3 - 5463 reads in 1477 unique sequences.
Sample 4 - 2914 reads in 904 unique sequences.
Sample 5 - 2941 reads in 939 unique sequences.
Sample 6 - 4312 reads in 1267 unique sequences.
Sample 7 - 6741 reads in 1756 unique sequences.
Sample 8 - 4560 reads in 1438 unique sequences.
Sample 9 - 15637 reads in 3590 unique sequences.
Sample 10 - 11413 reads in 2762 unique sequences.
Sample 11 - 12017 reads in 3021 unique sequences.
Sample 12 - 5032 reads in 1566 unique sequences.
Sample 13 - 18075 reads in 3707 unique sequences.
Sample 14 - 6250 reads in 1479 unique sequences.
Sample 15 - 4052 reads in 1195 unique sequences.
Sample 16 - 7369 reads in 1832 unique sequences.
Sample 17 - 4765 reads in 1183 unique sequences.
Sample 18 - 4871 reads in 1382 unique sequences.
Sample 19 - 6504 reads in 1709 unique sequences.
Sample 20 - 4314 reads in 897 unique sequences.

dadaRs <- dada(filtRs, err=errR, multithread=TRUE)

Sample 1 - 7113 reads in 1660 unique sequences.
Sample 2 - 5299 reads in 1349 unique sequences.
Sample 3 - 5463 reads in 1335 unique sequences.
Sample 4 - 2914 reads in 853 unique sequences.
Sample 5 - 2941 reads in 880 unique sequences.
Sample 6 - 4312 reads in 1286 unique sequences.
Sample 7 - 6741 reads in 1803 unique sequences.
Sample 8 - 4560 reads in 1265 unique sequences.
Sample 9 - 15637 reads in 3414 unique sequences.
Sample 10 - 11413 reads in 2522 unique sequences.
Sample 11 - 12017 reads in 2771 unique sequences.
Sample 12 - 5032 reads in 1415 unique sequences.
Sample 13 - 18075 reads in 3290 unique sequences.
Sample 14 - 6250 reads in 1390 unique sequences.
Sample 15 - 4052 reads in 1134 unique sequences.
Sample 16 - 7369 reads in 1635 unique sequences.
Sample 17 - 4765 reads in 1084 unique sequences.
Sample 18 - 4871 reads in 1161 unique sequences.
Sample 19 - 6504 reads in 1502 unique sequences.
Sample 20 - 4314 reads in 732 unique sequences.

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

mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)

6540 paired-reads (in 107 unique pairings) successfully merged out of 6891 (in 197 pairings) input.
5028 paired-reads (in 101 unique pairings) successfully merged out of 5190 (in 157 pairings) input.
4986 paired-reads (in 81 unique pairings) successfully merged out of 5267 (in 166 pairings) input.
2595 paired-reads (in 52 unique pairings) successfully merged out of 2754 (in 108 pairings) input.
2553 paired-reads (in 60 unique pairings) successfully merged out of 2785 (in 119 pairings) input.
3646 paired-reads (in 55 unique pairings) successfully merged out of 4109 (in 157 pairings) input.
6079 paired-reads (in 81 unique pairings) successfully merged out of 6514 (in 198 pairings) input.
3968 paired-reads (in 91 unique pairings) successfully merged out of 4388 (in 187 pairings) input.
14233 paired-reads (in 143 unique pairings) successfully merged out of 15355 (in 352 pairings) input.
10528 paired-reads (in 120 unique pairings) successfully merged out of 11165 (in 278 pairings) input.
11154 paired-reads (in 137 unique pairings) successfully merged out of 11797 (in 298 pairings) input.
4349 paired-reads (in 85 unique pairings) successfully merged out of 4802 (in 179 pairings) input.
17431 paired-reads (in 153 unique pairings) successfully merged out of 17812 (in 272 pairings) input.
5850 paired-reads (in 81 unique pairings) successfully merged out of 6095 (in 159 pairings) input.
3716 paired-reads (in 86 unique pairings) successfully merged out of 3894 (in 147 pairings) input.
6865 paired-reads (in 99 unique pairings) successfully merged out of 7191 (in 187 pairings) input.
4426 paired-reads (in 67 unique pairings) successfully merged out of 4603 (in 127 pairings) input.
4576 paired-reads (in 101 unique pairings) successfully merged out of 4739 (in 174 pairings) input.
6092 paired-reads (in 109 unique pairings) successfully merged out of 6315 (in 173 pairings) input.
4269 paired-reads (in 20 unique pairings) successfully merged out of 4281 (in 28 pairings) input.

# Inspect the merger data.frame from the first sample
head(mergers[[1]])

seqtab <- makeSequenceTable(mergers)
dim(seqtab)

[1]  20 293

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

251 252 253 254 255 
  1  88 196   6   2 

seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)

Identified 61 bimeras out of 293 input sequences.

dim(seqtab.nochim)

[1]  20 232

sum(seqtab.nochim)/sum(seqtab)

[1] 0.9640374

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)

       input filtered denoisedF
F3D0    7793     7113      6976
F3D1    5869     5299      5227
F3D141  5958     5463      5331
F3D142  3183     2914      2799
F3D143  3178     2941      2822
F3D144  4827     4312      4151
       denoisedR merged nonchim
F3D0        6979   6540    6528
F3D1        5239   5028    5017
F3D141      5357   4986    4863
F3D142      2830   2595    2521
F3D143      2868   2553    2519
F3D144      4228   3646    3507

taxa <- assignTaxonomy(seqtab.nochim, "~/silva_nr_v132_train_set.fa.gz", multithread=TRUE)

taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)

     Kingdom    Phylum         
[1,] "Bacteria" "Bacteroidetes"
[2,] "Bacteria" "Bacteroidetes"
[3,] "Bacteria" "Bacteroidetes"
[4,] "Bacteria" "Bacteroidetes"
[5,] "Bacteria" "Bacteroidetes"
[6,] "Bacteria" "Bacteroidetes"
     Class         Order          
[1,] "Bacteroidia" "Bacteroidales"
[2,] "Bacteroidia" "Bacteroidales"
[3,] "Bacteroidia" "Bacteroidales"
[4,] "Bacteroidia" "Bacteroidales"
[5,] "Bacteroidia" "Bacteroidales"
[6,] "Bacteroidia" "Bacteroidales"
     Family          
[1,] "Muribaculaceae"
[2,] "Muribaculaceae"
[3,] "Muribaculaceae"
[4,] "Muribaculaceae"
[5,] "Bacteroidaceae"
[6,] "Muribaculaceae"
     Genus        
[1,] NA           
[2,] NA           
[3,] NA           
[4,] NA           
[5,] "Bacteroides"
[6,] NA           

unqs.mock <- seqtab.nochim["Mock",]
unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE) # Drop ASVs absent in the Mock
cat("DADA2 inferred", length(unqs.mock), "sample sequences present in the Mock community.\n")

DADA2 inferred 20 sample sequences present in the Mock community.

mock.ref <- getSequences(file.path(path, "HMP_MOCK.v35.fasta"))
match.ref <- sum(sapply(names(unqs.mock), function(x) any(grepl(x, mock.ref))))
cat("Of those,", sum(match.ref), "were exact matches to the expected reference sequences.\n")

Of those, 20 were exact matches to the expected reference sequences.

library(phyloseq); packageVersion("phyloseq")

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

[1] ‘1.44.0’

library(Biostrings); packageVersion("Biostrings")

Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, aperm,
    append, as.data.frame,
    basename, cbind, colnames,
    dirname, do.call,
    duplicated, eval, evalq,
    Filter, Find, get, grep,
    grepl, intersect,
    is.unsorted, lapply, Map,
    mapply, match, mget,
    order, paste, pmax,
    pmax.int, pmin, pmin.int,
    Position, rank, rbind,
    Reduce, rownames, sapply,
    setdiff, sort, table,
    tapply, union, unique,
    unsplit, which.max,
    which.min

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: ‘S4Vectors’

The following object is masked from ‘package:utils’:

    findMatches

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:phyloseq’:

    distance

Loading required package: XVector
Loading required package: GenomeInfoDb

Attaching package: ‘Biostrings’

The following object is masked from ‘package:base’:

    strsplit

[1] ‘2.68.1’

library(ggplot2); packageVersion("ggplot2")

Need help? Try Stackoverflow:
https://stackoverflow.com/tags/ggplot2

[1] ‘3.4.3’

theme_set(theme_bw())

samples.out <- rownames(seqtab.nochim)
subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
gender <- substr(subject,1,1)
subject <- substr(subject,2,999)
day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
samdf$When <- "Early"
samdf$When[samdf$Day>100] <- "Late"
rownames(samdf) <- samples.out

ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE), 
               sample_data(samdf), 
               tax_table(taxa))
ps <- prune_samples(sample_names(ps) != "Mock", ps) # Remove mock sample

dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- taxa_names(ps)
ps <- merge_phyloseq(ps, dna)
taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 232 taxa and 19 samples ]
sample_data() Sample Data:       [ 19 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 232 taxa by 6 taxonomic ranks ]
refseq()      DNAStringSet:      [ 232 reference sequences ]

plot_richness(ps, x="Day", measures=c("Shannon", "Simpson"), color="When")

Warning: The data you have provided does not have
any singletons. This is highly suspicious. Results of richness
estimates (for example) are probably unreliable, or wrong, if you have already
trimmed low-abundance taxa from the data.

We recommended that you find the un-trimmed data and retry.

# Transform data to proportions as appropriate for Bray-Curtis distances
ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")

Run 0 stress 0.08043117 
Run 1 stress 0.08043117 
... New best solution
... Procrustes: rmse 5.720551e-07  max resid 1.172521e-06 
... Similar to previous best
Run 2 stress 0.1212044 
Run 3 stress 0.08616061 
Run 4 stress 0.09477104 
Run 5 stress 0.1327107 
Run 6 stress 0.08076337 
... Procrustes: rmse 0.01049156  max resid 0.03228709 
Run 7 stress 0.08076336 
... Procrustes: rmse 0.01047811  max resid 0.0322432 
Run 8 stress 0.09477199 
Run 9 stress 0.08616061 
Run 10 stress 0.1228545 
Run 11 stress 0.08616061 
Run 12 stress 0.08616061 
Run 13 stress 0.08043117 
... Procrustes: rmse 4.784089e-06  max resid 1.209747e-05 
... Similar to previous best
Run 14 stress 0.08616061 
Run 15 stress 0.1010631 
Run 16 stress 0.08076337 
... Procrustes: rmse 0.01050448  max resid 0.0323295 
Run 17 stress 0.08043116 
... New best solution
... Procrustes: rmse 7.89584e-07  max resid 1.37892e-06 
... Similar to previous best
Run 18 stress 0.08076338 
... Procrustes: rmse 0.01053278  max resid 0.03242151 
Run 19 stress 0.08043117 
... Procrustes: rmse 1.24138e-06  max resid 2.668817e-06 
... Similar to previous best
Run 20 stress 0.08076337 
... Procrustes: rmse 0.01049622  max resid 0.03230052 
*** Best solution repeated 2 times

plot_ordination(ps.prop, ord.nmds.bray, color="When", title="Bray NMDS")

top20 <- names(sort(taxa_sums(ps), decreasing=TRUE))[1:20]
ps.top20 <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
plot_bar(ps.top20, x="Day", fill="Family") + facet_wrap(~When, scales="free_x")