tidyGenR core functions cover all the steps necessary to determine sequence variants and genotypes from sequencing reads coming from multilocus amplicon libraries sequenced by high throughput sequencing technologies.
This vignette contains a standard workflow for paired-end data.
library(tidyGenR)
#> Warning: replacing previous import 'BiocGenerics::detail' by
#> 'Biostrings::detail' when loading 'GenomicAlignments'
#> Warning: replacing previous import 'BiocGenerics::detail' by
#> 'Biostrings::detail' when loading 'cigarillo'
#> Warning: replacing previous import 'BiocGenerics::detail' by
#> 'Biostrings::detail' when loading 'pwalign'The raw and processed data used in this vignette come from amplicon libraries of Rattus baluensis (Camacho-Sanchez et al. 2026).
Naming of FASTQ must meet a given format.
Reads from different loci are mixed together in the same sample
FASTQ. The known primer sequences at the ends of amplicons can be used
to demultiplex reads by locus using demultiplex(). This
function produces an executable for cutadapt. A
data.frame with primer sequences can be used to feed
demultiplex() with primer sequences.
# data.frame with primers
data("primers")
head(primers, 3)
#> locus fw rv
#> 1 nfkbia GCCTCCAAACACACAGTCAT TGAGGAGAGCTATGACACGG
#> 2 chrna9 TTATCTGGGAGAGCGTGACC TTGGGAAARGATGAACCGGC
#> 3 rogdi AGAARCCGGCTCACTACCC GAGGCACAGCTTGTTGAGGp_sh_out <- tempfile(fileext = ".sh")
# demultiplex reads by locus using 3 primer pairs
demultiplex(
freads = freads,
rreads = rreads,
primers = primers[1:3, ],
sh_out = p_sh_out,
run = FALSE
)
#> Warning in system2("which", "cutadapt", stdout = TRUE): running command
#> ''which' cutadapt' had status 1 and error message 'Function not implemented'demultiplex() creates a script as below:
readLines(p_sh_out)
#> [1] "#!/bin/bash"
#> [2] "cat /tmp/RtmpH20xBv/samples | while read sample"
#> [3] "do"
#> [4] " \\"
#> [5] "--discard-untrimmed --pair-adapters --no-indels \\"
#> [6] "--overlap 15 \\"
#> [7] "-e 0.15 \\"
#> [8] "-g file:/tmp/RtmpH20xBv/fw_primers \\"
#> [9] "-G file:/tmp/RtmpH20xBv/rv_primers \\"
#> [10] "-o /tmp/RtmpH20xBv/$sample.{name}.1.fastq.gz \\"
#> [11] "-p /tmp/RtmpH20xBv/$sample.{name}.2.fastq.gz \\"
#> [12] "/tmp/Rtmp0p1S7G/Rinst13e55c316639/tidyGenR/extdata/raw/\"$sample\".1.fastq.gz \\"
#> [13] "/tmp/Rtmp0p1S7G/Rinst13e55c316639/tidyGenR/extdata/raw/\"$sample\".2.fastq.gz"
#> [14] "done > /tmp/RtmpH20xBv/file149b446d2a0.log"If problems arise when running cutadapt
run = TRUE, the script written to sh_out can
be edited manually and run outside R. FASTQ with zero or few reads can
be removed using remove_poor_fastq()
Variants are determined on the basis of 1-nt resolution. Thus, reads
to be compared need to have the same length and minimum quality
stantards (eg. no ‘Ns’) (see DADA2 documentation, Callahan et al. 2016).
Sequence quality can be checked with FASTQC,
dada2::plotQualityProfile() or other QC checking software
to guide truncation lengths. Truncation and minimal quality filtering
are performed with trunc_amp(), a wrapper function of
dada2::filterAndTrim(). Global or locus-specific truncation
lengths can be supplied in a data.frame.
# Download per-locus demultiplexed FASTQ
dem <-
system.file("extdata", "demultiplexed", package = "tidyGenR")
# truncate reads for one locus
p_trun <- file.path(tempdir(), "truncated")
one_locus <-
trunc_amp(
loci = "chrna9",
mode_trun = "pe",
in_dir = dem,
fw_pattern = "1.fastq.gz",
rv_pattern = "2.fastq.gz",
trunc_fr = c(250, 180),
max_ee = c(3, 3),
outdir = p_trun
)
#>
#> Samples detected:
#> BOR1061, BOR1063, BOR1069
#>
#> Loci detected:
#> chrna9, nfkbia, rogdi
#>
#> Sequencing files matching patterns seem to conform the required naming format 'sample.locus.[1|2].fastq.gz'
#>
#> Truncating chrna9
#> Creating output directory: /tmp/RtmpH20xBv/truncated
#>
#> Forwards reads of chrna9 truncated at 250 nt.
#>
#> Reverse reads of chrna9 truncated at 180 nt.
#> Filtered/truncated reads written to/tmp/RtmpH20xBv/truncatedIt is possible to use locus-specific truncation lengths specified in
a data.frame.
# dataframe with locus-specific truncation lengths
data("trunc_fr")
head(trunc_fr, 3)
#> locus trunc_f trunc_r
#> 1 nfkbia 198 197
#> 2 chrna9 196 182
#> 3 rogdi 183 185# truncate reads for all loci using locus-specific truncation lengths
all_loci <-
trunc_amp(
mode_trun = "pe",
in_dir = dem,
fw_pattern = "1.fastq.gz",
rv_pattern = "2.fastq.gz",
trunc_fr = trunc_fr,
max_ee = c(3, 3),
outdir = p_trun
)trunc_amp() writes truncated FASTQ to a directory and
returns a list with in and out reads.
Truncated FASTQ:
head(list.files(p_trun))
#> [1] "BOR1061_chrna9_F_filt.fastq.gz" "BOR1061_chrna9_R_filt.fastq.gz"
#> [3] "BOR1061_nfkbia_F_filt.fastq.gz" "BOR1061_nfkbia_R_filt.fastq.gz"
#> [5] "BOR1061_rogdi_F_filt.fastq.gz" "BOR1061_rogdi_R_filt.fastq.gz"IN and OUT reads after trunc_amp():
Real sequence variants and their frequency are determined with
variant_call(). It utilizes the ‘DADA2’ ASV determination
pipeline, but for multiple loci. Determination of variants and
potentially, merging paired-end reads, and removal of chimeras are
performed altogether with variant_call(). Binned read
qualities (i.e. from NovaSeq) require the use of
error_function = loess_err_mod4 (details in
?loess_err_mod4()). Variants are returned in
tidy format.
head(variants)
#> # A tibble: 6 × 7
#> sample locus variant reads nt md5 sequence
#> <chr> <fct> <chr> <int> <dbl> <chr> <chr>
#> 1 BOR1061 chrna9 01 28 196 9598f07291c660d1a00c1497350381ad TGCAGTGTG…
#> 2 BOR1061 chrna9 02 34 196 a189edcac75e0853e8fac55df31acede TGCAGTGTG…
#> 3 BOR1061 nfkbia 1 46 198 cf753e1049c5f09f805df1918309a34b CGTAGGGCA…
#> 4 BOR1061 rogdi 1 50 183 6ddfd64d697d724f09d5f0ba1c909de6 CAGCCACCC…
#> 5 BOR1063 chrna9 01 52 196 9598f07291c660d1a00c1497350381ad TGCAGTGTG…
#> 6 BOR1063 nfkbia 1 44 198 cf753e1049c5f09f805df1918309a34b CGTAGGGCA…Variants can be further filtered with filter_variants()
based on frequency (maf) and depth (ad)
thresholds.
Genotypes can be determined from variants. So far, the implemented
method for genotyping is based on the expected maximum number of alleles
per sample and locus and the read threshold (ADt) for
discriminating homozygotes from hemizygotes.
genotypes <-
genotype(variants, ADt = 10, ploidy = 2)
head(genotypes)
#> # A tibble: 6 × 8
#> sample locus allele allele_no reads nt md5 sequence
#> <chr> <chr> <chr> <int> <dbl> <dbl> <chr> <chr>
#> 1 BOR1061 chrna9 01 1 28 196 9598f07291c660d1a00c1497… TGCAGTG…
#> 2 BOR1061 chrna9 02 2 34 196 a189edcac75e0853e8fac55d… TGCAGTG…
#> 3 BOR1061 nfkbia 1 1 23 198 cf753e1049c5f09f805df191… CGTAGGG…
#> 4 BOR1061 nfkbia 1 2 23 198 cf753e1049c5f09f805df191… CGTAGGG…
#> 5 BOR1061 rogdi 1 1 25 183 6ddfd64d697d724f09d5f0ba… CAGCCAC…
#> 6 BOR1061 rogdi 1 2 25 183 6ddfd64d697d724f09d5f0ba… CAGCCAC…