Beállítás és bevitel

  • Kezdjük párosított végű FASTQ fájlokkal (pl. sample1_R1.fastq, sample1_R2.fastq).

  • A fájloknak a következőknek kell lenniük:

    • Demultiplexelt (minta szerint felosztva)
    • Primer-trimmelt

  • A DADA2 az előremenő (R1) és a fordított (R2) fájlokat a minta neve alapján automatikusan illeszti egymáshoz.

Miért az előre (R1) és a fordított (R2) olvasatok?

A párosított végű szekvenálás (pl. Illumina MiSeq) a DNS-fragmentum két végét olvassa, és mintánként két FASTQ-fájlt generál:

  • R1.fastq → Előre olvasás

  • R2.fastq → Visszaolvasás

  • Az előre olvasás a felső szál 5′ végéről származik

  • A fordított olvasás az alsó szál 5′ végéről származik (fordított komplementer)

  • A két leolvasás átfedhet, ami lehetővé teszi a teljes DNS-darab jobb rekonstrukcióját.

A DADA2 páros végű leolvasásokat használ:

  • Függetlenül kiszűri és modellezi az egyes irányok hibáit.

  • Elvegyíti a leolvasásokat az átfedésük segítségével:

    • A szekvenálási hibák kijavítása
    • A teljes vonalkód (amplikon) rekonstruálása
    • A pontosság és a felbontás javítása
Fájl Leírás
R1.fastq Előre olvasott adatok (5′→3′ egy szálból)
R2.fastq Fordított olvasatok (5′→3′ a másik szálból)
  • Mindkettőre szükség van a teljes hosszúságú, jó minőségű szekvenciák létrehozásához.

A minőség vizsgálata

  • A DADA2 a minőségprofilok vizsgálatával kezdi, hogy meghatározza a csonkítási hosszokat és a szűrési paramétereket.
plotQualityProfile(fwd)
plotQualityProfile(rev)
  • Előre olvasott adatok: általában jó minőség
  • Fordított olvasatok: a minőség gyakran csökken a vége felé.
  • Használja ezt a csonkítási hossz meghatározásához (pl. 240 bp előre, 160 bp hátra).

Szűrés és levágás

  • A DADA2 a filterAndTrim() függvényt használja a minőségellenőrzéshez és a csonkításhoz.
filterAndTrim(fwd, filtF,
              rev, filtR,
              truncLen = c(240,160),
              maxN = 0,
              maxEE = c(2,2),
              truncQ = 2,
              rm.phix = TRUE,
              compress = TRUE,
              multithread = TRUE)
  • Eltávolítja a rossz minőségű és kétértelmű olvasásokat.
  • Eltávolítja a PhiX szennyeződéseket
  • Mintánként visszatartott olvasásszámokat ad ki
Mi az a PhiX és miért távolítják el?

A PhiX (vagy PhiX174) egy kis vírusgenom, amelyet gyakran adnak hozzá kontrollként az Illumina szekvenálási futtatások során.

  • Segít a szekvenáló kalibrálásában és a teljesítmény javításában alacsony diverzitású minták (például amplikon alapú eDNS) esetén.
  • Ez nem része a tényleges biológiai mintának, ezért ki kell szűrni.

A DADA2-ben ez automatikusan történik:

rm.phix = TRUE

Hibaarányok tanulása

  • A DADA2 a learnErrors() függvényt használja a szekvenálási hibák modellezésére.
errF <- learnErrors(filtFs, multithread = TRUE)
errR <- learnErrors(filtRs, multithread = TRUE)
  • Felépít egy parametrikus hibamodellt a megfigyelt adatokból.
  • Külön-külön az előre és a hátrafelé történő olvasásra

Zajszűrés / Mintavételi következtetés

  • A DADA2 a dada() függvényt használja a zajszűréshez és a mintavételi következtetéshez.
dadaFs <- dada(filtFs, err = errF, multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, multithread = TRUE)
  • Az ASV-k származtatása.
  • Korrigálja a szekvenálási hibákat a valódi biológiai szekvenciák azonosításához.

Páros olvasatok egyesítése

  • A DADA2 a mergePairs() függvényt használja az előre és a hátrafelé olvasott adatok egyesítésére.
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs)
  • Összevonja az előre és a hátrafelé olvasott adatokat
  • Átfedés (pl. ≥12 bp) és egyezés szükséges az átfedési régióban.

Szekvencia táblázat készítése

  • A DADA2 a makeSequenceTable() függvényt használja az ASV-k táblázatának létrehozásához.
seqtab <- makeSequenceTable(mergers)
  • Sorok = minták
  • Oszlopok = ASV-k
  • Cellák = olvasott adatok száma

Kimérák eltávolítása

  • A DADA2 a removeBimeraDenovo() függvényt használja a kimérák eltávolítására.
seqtab.nochim <- removeBimeraDenovo(seqtab, method = "consensus", multithread = TRUE)
  • Eltávolítja a PCR leletek (kimérák)
  • Egy megtisztított ASV táblázatot eredményez.
Mik azok a kimérák és miért kell eltávolítani őket?

A kimerák az amplifikáció során létrehozott mesterséges DNS-szekvenciák.

  • Akkor jönnek létre, amikor az egyik templátból származó, részben meghosszabbított DNS-szál véletlenül egy másik, de hasonló templáthoz kötődik a következő PCR-ciklusban.

  • Ezáltal egy olyan hibrid szekvencia jön létre, amely nem létezik egyetlen valódi szervezetben sem.

  • Gyakoriak az amplikon-alapú munkafolyamatokban.

  • Helhúzhatják a biológiai sokféleség becslését vagy téves fajmeghatározáshoz vezethetnek.

Olvasásvesztés nyomon követése

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")
  • Figyeli, hogy hány olvasás halad át az egyes lépéseken
  • Segít észlelni a nagyobb adatvesztést vagy adatcsatorna problémákat

Taxonómia hozzárendelése

  • A DADA2 a assignTaxonomy() függvényt használja a taxonómiai osztályozáshoz.
taxa <- assignTaxonomy(seqtab.nochim, "silva_nr_v138_train_set.fa.gz", multithread = TRUE)
taxa <- addSpecies(taxa, "silva_species_assignment_v138.fa.gz")
  • Az ASV-k osztályozása egy referenciaadatbázissal (pl. SILVA) való összehasonlítással.
  • Hozzáadja a fajszintű azonosítókat a pontos egyezésekhez

Választható: fajszintű megfeleltetés

  • Használjon kurátori referenciákat, mint például a MiFish a hal-specifikus azonosításhoz.

  • Alternatívák:

    • BLAST a fajok kézi ellenőrzéséhez
    • DECIPHER IdTaxa() függvénye.

Validálás ismert közösséggel

  • Az adatcsatorna futtatása egy ismert mintán
  • A várt fajokat nagy pontossággal vissza kell állítania

Elmentes a phyloseq-be

  • A DADA2 eredményeit a phyloseq csomagba importálhatja, amely egy ökoszisztéma az ökológiai adatok kezelésére.
library(phyloseq)
ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows = FALSE),
               tax_table(taxa),
               sample_data(metadata))

Példa kód

library(dada2)
library(phyloseq)
library(Biostrings)
library(ggplot2)

path <- "MiSeq_SOP"
list.files(path)

fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))

sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

plotQualityProfile(fnFs[1:2])
plotQualityProfile(fnRs[1:2])

filtFs <- file.path(path, "filtered", basename(fnFs))
filtRs <- file.path(path, "filtered", basename(fnRs))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
                     truncLen=c(240,200), maxN=0, maxEE=c(2,2),
                     truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=FALSE)
head(out)

errF <- learnErrors(filtFs, multithread=FALSE)
errR <- learnErrors(filtRs, multithread=FALSE)
plotErrors(errF, nominalQ=TRUE)

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

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

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

seqtab <- makeSequenceTable(mergers)
dim(seqtab)
table(nchar(getSequences(seqtab)))

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

getN <- function(x) sum(getUniques(x))
track <- data.frame(
  input=out[,1],
  filtered=out[,2],
  denoisedF=sapply(dadaFs, getN),
  denoisedR=sapply(dadaRs, getN),
  merged=sapply(mergers, getN),
  nonchim=rowSums(seqtab.nochim)
)
rownames(track) <- sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
print(track)

taxa <- assignTaxonomy(seqtab.nochim, "silva_nr_v138_train_set.fa.gz", multithread=FALSE)
print(head(taxa))

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

rownames(seqtab.nochim)
sample.names <- sapply(strsplit(rownames(seqtab.nochim), "_"), `[`, 1)
rownames(seqtab.nochim) <- sample.names

unqs.mock <- seqtab.nochim["Mock",]
unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE) 
cat("DADA2 inferred", length(unqs.mock), "minta szekvenciák jelen vannak a Mock közösségben.\n")

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("Ezek közül", sum(match.ref), "pontosan egyeznek a várt referencia szekvenciákkal.\n")

######

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)

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

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

ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")

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")

A relatív gyakoriság megértése az DNS vizsgálatokban

Relatív gyakoriság

A DNS használatakor a teljes mintának csak egy részét használjuk fel, és az amplifikáció pedig torzítást eredményez azáltal, hogy több millió DNS másolatot hoz létre. Ez felveti a kérdést:

Hogyan bízhatunk abban, hogy a DNS-ek száma (vagy a “leolvasások száma”) tükrözi a fajok valódi gyakoriságát a környezetben?

Az amplifikáció exponenciálisan sokszorozza a DNS-t - ideális esetben minden egyes ciklusban megduplázza a DNS másolatok számát. Ez a folyamat azonban az amplifikációs torzítást eredményez, mint például:

  • Primerhatékonysági torzítás: egyes szekvenciák jobban kötődnek a primerekhez és könnyebben amplifikálódnak.
  • Sablonpreferencia: bizonyos DNS-szekvenciák előnyben részesülnek a másolás során.
  • Plateau-effektus: az erőforrások kimerülésével az amplifikáció lelassul és leáll.
  • Részmintavételezés: a teljes DNS-nek csak egy kis részét használják fel a PCR-hez és a szekvenáláshoz.

Ez azt jelenti, hogy az amplifikáció nem őrzi meg pontosan a DNS eredeti arányait - a DNS jelenlétét igazolja, de nem a tökéletes gyakoriságot.

A leolvasásszámlálás azonban még így is hasznos, de óvatosan kell értelmezni. Míg az amplifikáció torzítja az abszolút mennyiségeket, a leolvasott számok még mindig ökológiai jeleket hordoznak, különösen a domináns fajok esetében.

A kutatók jellemzően a következőket használják:

  • Relatív leolvasásszám (Relative Read Abundance, RRA):

    • A leolvasások 40%-a = európai tengeri sügér
    • A leolvasások 10%-a = háromágú pikó

Az RRA lehetővé teszi a mintázatok összehasonlítását a különböző helyszíneken vagy körülmények között, nem pedig a pontos organizmusszámok összehasonlítását.

A szekvenálás után egy olvasatszámtáblázat mutatja, hogy az egyes fajok esetében mintánként hány olvasatot detektáltunk.

Ahhoz, hogy ezek az adatok összehasonlíthatóak legyenek, a következő technikákat alkalmazzuk:

Normalizálás

Cél: A különböző szekvenálási mélységek (azaz a mintánkénti különböző teljes olvasatszámok) kiigazítása.

  • Arányos skálázás: a nyers számokat relatív arányokká (az összes olvasat %-ában) alakítja át.
  • Log-transzformáció: tompítja a nagyon gyakori fajok hatását.
  • Modellalapú megközelítések: pl. DESeq2, edgeR - a számlálási variabilitást statisztikailag figyelembe veszi.

Ritkítás

Cél: A szekvenálási ráfordítások standardizálása a minták között.

  • Minden adatkészletet közös olvasatszámra (pl. 10 000 olvasat mintánként) vesz részmintát.
  • Biztosítja a diverzitás tisztességes összehasonlítását.
  • Megakadályozza a gazdagság túlbecslését a mélyebben szekvenált mintákban.

Kompromisszum: egyes adatok (pl. ritka olvasatok) elvesznek, de az eredmények kevésbé torzítottak.

Statisztikai eszközök

Miután az adatokat normalizáltuk és ritkítottuk, különböző ökológiai és statisztikai elemzéseket alkalmazhatunk:

  • Diverzitási indexek: Shannon, Simpson - a fajgazdagság és az egyenletesség mérésére.
  • Rendezés: PCA, NMDS, PCoA - a helyszínek közötti mintázatok megjelenítésére.
  • Hipotézisvizsgálat: PERMANOVA, ANOSIM - csoportok vagy kezelések összehasonlítására
  • Ökológiai modellezés: A környezeti tényezők és a közösségek közötti kapcsolatok feltárása