Genetic diversity

Genetic diversity is provided by two lightweight extensions, so you only load what your data needs:

  • sequences — load BioSequences (using Diversity, BioSequences) to build a GeneticType from a vector of aligned BioSequences;
  • VCF — load PopGen (using Diversity, PopGen) to build a GeneticType from a PopGen.PopData object read from a VCF file.

Similarity between types is derived from pairwise genetic distances, mirroring the gen2dist() / dist2sim() pipeline in the R package rdiversity.

Usage

Using the functionality in the package is simple:

  • Create genetic data, either a vector of aligned BioSequences or a PopGen.PopData object (e.g. read from a VCF file with PopGen.vcf)
  • Create a GeneticType (an AbstractTypes subtype) from it
  • Create a Metacommunity from that
  • Calculate diversity!

Sequences

julia> using Diversity, BioSequences
julia> seqs = [dna"ACGTACGT", dna"ACGAACGT", dna"TTTTTTTT"]3-element Vector{BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}: ACGTACGT ACGAACGT TTTTTTTT
julia> gt = GeneticType(seqs; names = ["a", "b", "c"])DiversityBioSequencesExt.GeneticFASTA{Vector{BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}}(BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}[ACGTACGT, ACGAACGT, TTTTTTTT], ["a", "b", "c"], 3, [1.0 0.8571428571428572 0.1428571428571429; 0.8571428571428572 1.0 0.0; 0.1428571428571429 0.0 1.0])
julia> calcsimilarity(gt, 1.0)3×3 Matrix{Float64}: 1.0 0.857143 0.142857 0.857143 1.0 0.0 0.142857 0.0 1.0
julia> metagen = Metacommunity([0.3, 0.3, 0.4], gt)Metacommunity{Float64, Vector{Float64}, Matrix{Float64}, DiversityBioSequencesExt.GeneticFASTA{Vector{BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}}, Onecommunity} with 3 species in 1 subcommunity measuring Genetic (sequence) diversity. Species names: a, b, c Subcommunity names: 1
julia> meta_gamma(metagen, 0)1×8 DataFrame Row div_type measure q type_level type_name partition_le String String Int64 String String String ⋯ ─────┼────────────────────────────────────────────────────────────────────────── 1 │ Genetic (sequence) Gamma 0 types metacommunit ⋯ 3 columns omitted

Sequences a and b differ at one site out of eight and are correspondingly similar; both are far from c, which shares no site with either. The metacommunity therefore holds rather less than three types' worth of diversity.

VCF

The example below reads the small biallelic VCF that ships with the package, and builds a similarity matrix using biallelic Manhattan distances — matching rdiversity's gen2dist(vcf, biallelic = TRUE):

julia> using Diversity, PopGen
julia> pd = PopGen.vcf(Diversity.path("data", "biallelic.vcf"); silent = true)PopData{Diploid, 3 SNP loci} Samples: 3 Populations: 1
julia> gt = GeneticType(pd; distance = :manhattan)DiversityPopGenExt.GeneticVCF{PopGenCore.PopData}(PopData{Diploid, 3 SNP loci} Samples: 3 Populations: 1, ["s1", "s2", "s3"], 3, [1.0 0.0 0.25; 0.0 1.0 0.25; 0.25 0.25 1.0])
julia> gettypenames(gt, true)3-element Vector{String}: "s1" "s2" "s3"
julia> calcsimilarity(gt, 1.0)3×3 Matrix{Float64}: 1.0 0.0 0.25 0.0 1.0 0.25 0.25 0.25 1.0
julia> metagen = Metacommunity([0.2 0.1; 0.1 0.3; 0.2 0.1], gt)Metacommunity{Float64, Matrix{Float64}, Matrix{Float64}, DiversityPopGenExt.GeneticVCF{PopGenCore.PopData}, Subcommunities} with 3 species in 2 subcommunities measuring Genetic (VCF) diversity. Species names: s1, s2, s3 Subcommunity names: 1, 2
julia> meta_gamma(metagen, 0)1×8 DataFrame Row div_type measure q type_level type_name partition_level String String Int64 String String String ⋯ ─────┼────────────────────────────────────────────────────────────────────────── 1 │ Genetic (VCF) Gamma 0 types metacommunity ⋯ 2 columns omitted

Here the samples are the types, so a metacommunity is a matrix of sample abundances with one column per subcommunity.

The distance (:manhattan or :hamming), transform (:linear or :exponential), k and normalise keyword arguments control how the pairwise distances and the resulting similarity matrix are calculated.

vcf_dataframe(pd) converts a PopData back into the VCF-body layout that rdiversity's gen2dist() consumes, so the same data can drive both the Julia and the R calculation — which is how the two are cross-validated against each other in test/run_rcall.jl.