Package 'piglet'

Description

A data.table of the allele similarity cluster table based on the HVGERM and hv_functionality germlie reference set. This is not the latest version of the allele similarity cluster table. For the latest version please refer either to the zenodo doi or you can use the recentAlleleClusters

Usage

allele_cluster_table

Format

An object of class data.table (inherits from data.frame) with 286 rows and 5 columns.

References

Peres, et al (2022) https://doi.org/10.1101/2022.12.26.521922

Description

Compare the sequences of two alleles (reference and sample alleles) and returns the differential nucleotide positions of the sample allele.

Usage

allele_diff(
  reference_allele,
  sample_allele,
  position_threshold = 0,
  snps = TRUE
)

Arguments

Details

The function utilizes c++ script to optimize the run time for large comparisons.

Value

A character vector of the differential nucleotide positions of the sample allele.

Examples

{
reference_allele = "AAGG"
sample_allele = "ATGA"

# setting position_threshold = 0 will return all differences
diff <- allele_diff(reference_allele, sample_allele)
# "A2T", "G4A"
print(diff)

# setting position_threshold = 3 will return the differences from position three onward
diff <- allele_diff(reference_allele, sample_allele, position_threshold = 3)
# "G4A"
print(diff)

# setting snps = FALSE will return the differences as indices
diff <- allele_diff(reference_allele, sample_allele, snps = FALSE)
# 2, 4
print(diff)

}

Description

Calculate differences between characters in columns of germs and return their indices as an int vector.

Usage

allele_diff_indices(germs, X = 0L, non_mismatch_chars_nullable = NULL)

Arguments

Value

A vector of integers containing indices of differing columns.

Examples

germs = c("ATCG", "ATCC") 
X = 3 
result = allele_diff_indices(germs, X)
# 1, 2, 3

Deprecated: use allele_diff_paired

Description

Renamed, and the counts change. The old implementation ignored gaps and ambiguous bases only on the germline side, so an aligned position where the germline carried a base and the input carried a gap was counted as a mismatch. It also looped to the length of the germline while indexing the input, reading past the end of the input whenever the input was shorter – undefined behaviour, and on real IGHV data it returned counts derived partly from adjacent memory.

Usage

allele_diff_indices_parallel(
  germs,
  inputs,
  X = 0L,
  parallel = FALSE,
  return_count = FALSE
)

Arguments

Details

Calls now route to allele_diff_paired, which ignores such positions on either side and pads the shorter sequence. Against the old behaviour, 194 of 400 real germline pairs differ, by a mean of 8 mismatches. Anything that depended on the old numbers needs re-checking rather than re-running.

Value

As allele_diff_paired.

See Also

allele_diff_paired

Deprecated: use allele_diff_paired

Description

Renamed. The old name said parallel, which was an argument rather than a behaviour and never had any effect (the package is not built with OpenMP), and indices, which is wrong whenever return_count = TRUE.

Usage

allele_diff_indices_parallel2(
  germs,
  inputs,
  X = 0L,
  parallel = FALSE,
  return_count = FALSE,
  non_mismatch_chars_nullable = NULL
)

Arguments

Details

Results are unchanged: this is the same implementation under a clearer name.

Value

As allele_diff_paired.

See Also

allele_diff_paired

Description

Compares each germline sequence in germs against the input sequence at the same position in inputs, and returns either the positions of the mismatches or how many there are. The two vectors are paired element by element, which is what distinguishes this from allele_diff_indices, where every sequence is compared against the first.

Usage

allele_diff_paired(
  germs,
  inputs,
  X = 0L,
  return_count = FALSE,
  non_mismatch_chars_nullable = NULL
)

Arguments

Details

A position is ignored when either sequence carries a non-mismatch character, by default a gap or an ambiguous base (N, ., -). Sequences of unequal length are padded to the longer with N, so the padded positions never count as mismatches.

Value

A list of integer vectors of one-based positions when return_count = FALSE, or an integer vector of counts when TRUE.

Examples

germs  <- c("ACGTACGT", "ACGTACGT")
inputs <- c("ACGTTCGT", "ACGTTCGA")

allele_diff_paired(germs, inputs, return_count = TRUE)   # 1, 2
allele_diff_paired(germs, inputs)                        # list(5), c(5, 8)

Description

Calculate differences between characters in columns of germs and return them as a string vector.

Usage

allele_diff_strings(germs, X = 0L, non_mismatch_chars_nullable = NULL)

Arguments

Value

A vector of strings containing differences between characters in columns.

Examples

germs = c("ATCG", "ATCC") 
X = 3 
result = allele_diff_strings(germs, X) 
# "A2T", "T3C", "C2G"

Description

A data.table of the allele thresholds table. The V alleles are based on the HVGERM and hv_functionality germline reference set. The D, and the J are based on the AIRR-C reference set (https://zenodo.org/records/10489725). The table contains these columns: allele - the IUIS allele name, asc_allele - the allele name based on allele similarity clusters (only for V), threshold = the genotype threshold for the alleles.

Usage

allele_threshold_table

Format

An object of class data.table (inherits from data.frame) with 262 rows and 4 columns.

References

Peres, et al (2022) https://doi.org/10.1101/2022.12.26.521922

Description

For a given cluster the function collapse similar sequences and renames the sequences based on the ASC name scheme

Usage

alleleClusterNames(
  cluster,
  allele.cluster.table,
  germ.dist,
  chain,
  segment,
  family_prefix = TRUE,
  retain_subgroup = FALSE
)

Arguments

Value

A data.frame with the clusters renamed alleles based on the ASC scheme.

Description

allelePalette returns a named vector of colors for a set of alleles, reusing the VDJbase allele color scheme (see the vdjbasevis package) for the curated allele numbers and generating additional, visually distinct colors when there are more alleles than the curated palette covers. This lets genotype plots scale to an arbitrary number of alleles while staying consistent with VDJbase for the common ones.

Usage

allelePalette(alleles)

Arguments

Details

The curated colors are taken from vdjbasevis::allelePalette (alleles 01-19). Alleles beyond that first draw from a pool of additional distinct colors and then, if still short, from evenly spaced hues generated in HCL space.

Value

a named character vector mapping each unique allele to a hex color.

Examples

allelePalette(c("01", "02", "06"))
# also covers many alleles without running out of colors
length(allelePalette(sprintf("%02d", 1:40)))

Description

Translates the ASC gene calls on a repertoire into the IUIS labels built by ascIUISVocabulary, adding v_gene_iuis, j_gene_iuis and, where the chain has one, d_gene_iuis.

Usage

annotateRepertoireIUIS(
  data,
  vocabulary,
  chain,
  v_gene = "v_gene",
  j_gene = "j_gene",
  d_gene = "d_gene",
  d_call = "d_call_new"
)

Arguments

Details

For IGH the D label is per allele, so it is mapped from the allele-level call rather than from the collapsed gene, and a call spanning more than one ASC group names no single gene and stays NA. V and J map straight from their ASC gene.

Rows in and rows out are reported: a call that maps to nothing is left NA rather than dropped, and the count is messaged so a labelling gap cannot pass unnoticed.

Value

data as a data.table with the *_gene_iuis columns added.

See Also

ascIUISVocabulary

Examples

husa <- data.frame(
  allele      = c("IGHV1-2*01", "IGHV1-2*02", "IGHD1-CO5H*01"),
  chain       = "IGH",
  gene_type   = c("IGHV", "IGHV", "IGHD"),
  seq         = c("ACGTACGT", "ACGTACGA", "GGTTGGTT"),
  iuis_allele = c("IGHV1-2*01", "IGHV1-2*02", "IGHD1-20*01"))
voc <- ascIUISVocabulary(husa)

rep_dt <- data.frame(
  v_gene     = c("IGHV1-2", "IGHV1-2"),
  d_gene     = c("IGHD1-CO5H", "IGHD1-CO5H"),
  j_gene     = c("IGHJ4", "IGHJ4"),
  d_call_new = c("IGHD1-CO5H*01", "IGHD1-CO5H*01"))

annotateRepertoireIUIS(rep_dt, voc, chain = "IGH")

Description

A function to artificially create an IGHV reference set with framework1 (FWR1) primers (see Details).

Usage

artificialFRW1Germline(
  germline_set,
  mask_primer = TRUE,
  trimm_primer = FALSE,
  quite = FALSE
)

Arguments

Details

The FRW1 primers used in this function were taken from the BIOMED-2 protocol. For more information on the protocol and primer design go to: van Dongen, J., Langerak, A., Brüggemann, M. et al. Design and standardization of PCR primers and protocols for detection of clonal immunoglobulin and T-cell receptor gene recombinations in suspect lymphoproliferations: Report of the BIOMED-2 Concerted Action BMH4-CT98-3936. Leukemia 17, 2257–2317 (2003). https://doi.org/10.1038/sj.leu.2403202Van Dongen, J. J. M., et al. "Design and standardization of PCR primers and protocols for detection of clonal immunoglobulin and T-cell receptor gene recombinations in suspect lymphoproliferations: report of the BIOMED-2 Concerted Action BMH4-CT98-3936." Leukemia 17.12 (2003): 2257-2317.

Value

A list with the input germline set allele and the trimmed/masked sequences.

Description

Derives the display labels that translate an Allele Similarity Cluster assignment into IUIS gene names, and returns the allele set behind each label.

Usage

ascIUISVocabulary(husa)

Arguments

Details

Two ASC subgroups that share any IUIS gene must carry the same label, so subgroups and IUIS parts are treated as a bipartite graph and every connected component gets one merged label. Where a component holds more than one subgroup the label alone does not identify the subgroup, so asc_and_iuis_group suffixes them A, B, C for display.

D is handled differently from V and J, and that asymmetry is the reason ASC labels look inconsistent downstream. D is relabelled at IUIS gene resolution by regroupASCByLabel: an ASC group holding alleles of more than one IUIS gene is broken apart and groups naming the same gene are joined, so the D label is a property of the allele rather than of the ASC gene. V and J keep the ASC grouping. This is why the resulting labels carry the locus for D (IGHD5-12) but not for V and J (V1-18, J4).

The label is a label. It is a display string, and it moves when the reference set behind it moves. What identifies a group is the set of member allele sequences, returned in identity; join on that, or compare it to detect that an upstream reference changed underneath a downstream table.

Value

A list of class ascIUISVocabulary:

labels

one row per ASC subgroup: chain, subgroup, iuis_group, asc_and_iuis_group, iuis_group_size.

asc_to_iuis

named character vector, ASC subgroup to IUIS label, for V and J.

d_allele_to_iuis

named character vector, ASC allele to IUIS label, for D.

identity

one row per chain and iuis_label: asc_genes, n_alleles, member_alleles and member_seq_md5. This is what the label means; the label itself is not.

See Also

annotateRepertoireIUIS, regroupASCByLabel

Examples

husa <- data.frame(
  allele      = c("IGHV1-2*01", "IGHV1-2*02", "IGHD1-CO5H*01"),
  chain       = "IGH",
  gene_type   = c("IGHV", "IGHV", "IGHD"),
  seq         = c("ACGTACGT", "ACGTACGA", "GGTTGGTT"),
  iuis_allele = c("IGHV1-2*01", "IGHV1-2*02", "IGHD1-20*01"))

voc <- ascIUISVocabulary(husa)
voc$labels

# What each label stands for. The label is a display string that moves when the
# reference set moves; these member alleles are what identify the group.
voc$identity

Description

A group's share of its own segment's repertoire in each subject, with a pseudocount so the logit is finite. This is the phenotype the usage scan explains.

Usage

ascUsagePhenotype(
  data,
  segments,
  subject = "subject",
  min_subject_fraction = 0.1,
  pseudocount = 0.5
)

Arguments

Details

Zero counts are kept deliberately. A variant that deletes a gene drives its usage to zero, and that is the canonical signal in this family, so dropping zeros would drop the effect. Every subject gets a row for every retained group, whether or not they used it.

Groups seen in fewer than min_subject_fraction of subjects are dropped before the shares are taken, and the count kept is reported. n_asc is how many groups the share was taken over and is part of the pseudocount, so usage cannot be recomputed from count and total without it.

Value

A data.table with subject, asc, count, total, segment, n_asc, usage and logit_usage.

See Also

runGeneUsageQTL

Examples

set.seed(1)
rep_dt <- data.frame(
  subject = rep(sprintf("s%02d", 1:70), each = 20),
  v_gene  = sample(c("V1-2", "V3-23", "V4-34"), 1400, replace = TRUE),
  j_gene  = sample(c("J4", "J6"), 1400, replace = TRUE))

usage <- ascUsagePhenotype(rep_dt, segments = c(V = "v_gene", J = "j_gene"))
head(usage)

# Zero counts are kept on purpose: a variant that deletes a gene drives its
# usage to zero, and that is the signal rather than missing data.
sum(usage$count == 0)

Description

assignAlleleClusters uses the allele clusters annotation to change the preliminary allele assignments to the new annotations before inferring a genotype.

Usage

assignAlleleClusters(
  data,
  alleleClusterTable,
  v_call = "v_call",
  from_col = "iuis_allele",
  to_col = "new_allele"
)

Arguments

Value

A modified input data.frame with the new assigned

Examples

# preferably obtain the latest ASC cluster table
# asc_archive <- recentAlleleClusters(doi="10.5281/zenodo.7429773", get_file = TRUE)

# allele_cluster_table <- extractASCTable(archive_file = asc_archive)

# example allele similarity cluster table
data(allele_cluster_table)

# loading TIgGER AIRR-seq b cell data
data <- tigger::AIRRDb

asc_data <- assignAlleleClusters(data, allele_cluster_table)

Description

Performs community detection on a weighted graph using the Leiden algorithm with CPM (Constant Potts Model) objective function.

Usage

detect_communities_leiden(g, resolution = 1)

Arguments

Details

The Leiden algorithm is a community detection method that optimizes a quality function (here CPM). It guarantees connected communities and is generally faster than Louvain while producing better quality partitions.

Value

An igraph communities object

See Also

distance_to_graph, optimize_resolution

Examples

data(HVGERM)
d <- igDistance(HVGERM[1:10], method = "hamming")
g <- distance_to_graph(d)
comm <- detect_communities_leiden(g, resolution = 0.5)

Description

Converts a distance matrix to a weighted igraph object using a log transform that spreads small distances and produces weights in [0,1].

Usage

distance_to_graph(distance_matrix)

Arguments

Details

The transformation uses a log-based similarity measure:

1. Normalize distances by max distance

2. Apply -log transform to convert to similarity

3. Normalize similarities to [0,1] range

4. Create weighted undirected graph

Value

An igraph object with weighted edges

See Also

detect_communities_leiden, igClust

Examples

data(HVGERM)
d <- igDistance(HVGERM[1:10], method = "hamming")
g <- distance_to_graph(d)

Description

Extracts the allele cluster table from the archive file.

Usage

extractASCTable(archive_file = NULL)

Arguments

Details

For downloading the latest archive file with the updated allele cluster table, use the function recentAlleleClusters.

Value

Returns the allele cluster table.

The table columns: new_allele - the ASC given allele name func_group - the ASC cluster number iuis_allele - the original IUIS allele name thresh - the allele threshold for ASC-based genotype inference amplicon_length - is the original length of the reference set.

Examples

## Not run: 
asc_archive <- recentAlleleClusters(doi="10.5281/zenodo.7429773", get_file = TRUE)

allele_cluster_table <- extractASCTable(archive_file = asc_archive)

## End(Not run)

Description

Generates the allele clusters reference set based on the clustering from ighvClust. The function collapse similar alleles and assign them into their respective allele clusters and family clusters. See details for naming scheme

Usage

generateReferenceSet(
  germline_distance,
  germline_set,
  alleleClusterTable,
  trim_3prime_side = NULL,
  family_prefix = TRUE,
  retain_subgroup = FALSE
)

Arguments

Details

Each allele is named by this scheme: IGHVF1-G1*01 - IGH = chain, V = region, F1 = family cluster numbering (the "F" prefix can be omitted by setting family_prefix = FALSE), G1 - allele cluster numbering, and 01 = allele numbering (given by clustering order, no connection to the expression)

In case there are alleles that are differentiated in a nucleotide position past the trimming position used for the clustering, then the alleles are separated and are annotated with the differentiating position as so: Say A101 and A102 are similar up to position 318, and thus collapsed in the clusters to G101. Upon checking the sequences past the trim position (318), a differentiating nucleotide was seen in position 319, A101 has a G, and A102 has a T. Then the alleles will be separated, and the new annotation will be as so: A101 = G101, and A102 = G1*01_G319T. Where the first nucleotide indicate the base, the following number the position, and the last nucleotide the one the base changed into.

Value

A list with the re-named germline set, and a table of the allele clusters and thresholds.

Description

genotypeToTigger reshapes the output of inferGenotypeAllele or inferGenotypeAllele_asc into the per-gene table layout used by TIgGER and VDJbase. The set of genotyped alleles is determined from the allele z-score and a z-score threshold.

Usage

genotypeToTigger(
  genotype,
  level = c("gene", "asc"),
  z_threshold = 0,
  file = NULL
)

Arguments

Value

A data.table with one row per gene (or ASC cluster), and columns:

• gene - the V gene or ASC cluster.

• alleles - the candidate allele numbers, comma-separated and ordered by descending count.

• counts - the read counts, comma-separated, matching alleles.

• total - the total read count for the gene.

• depth - the per-locus repertoire depth used as the denominator of the z-score.

• threshold - the per-allele presence thresholds, comma-separated, matching alleles.

• z_score - the per-allele z-scores, comma-separated, matching alleles.

• genotyped_alleles - the allele numbers with z_score >= z_threshold, comma-separated.

See Also

inferGenotypeAllele and inferGenotypeAllele_asc for producing the input genotype.

Examples

# loading TIgGER AIRR-seq b cell data
data <- tigger::AIRRDb

data(allele_threshold_table)
data(HVGERM)

genotype <- inferGenotypeAllele(
  data = data,
  allele_threshold_table = allele_threshold_table,
  germline_db = HVGERM, find_unmutated = TRUE)

# convert to the TIgGER/VDJbase table layout
geno_table <- genotypeToTigger(genotype)
head(geno_table)

Description

Converts IGHV germline set to ASC germline set.

Usage

germlineASC(allele_cluster_table, germline)

Arguments

Value

Returns the IGHV germline set with the ASC allele names.

Examples

# preferably obtain the latest ASC cluster table
# asc_archive <- recentAlleleClusters(doi="10.5281/zenodo.7429773", get_file = TRUE)

# allele_cluster_table <- extractASCTable(archive_file = asc_archive)

data(HVGERM)

# example allele similarity cluster table
data(allele_cluster_table)

asc_germline <- germlineASC(allele_cluster_table, germline = HVGERM)

Description

A data.table of all 498 human IGHV germline gene segment alleles in IMGT Gene-db release July 2022, with an additional 25 undocumented alleles from VDJbase. The first column is the allele name, the second column is the functionality annotation, the third column is the nt sequence and the last column is the aa sequence.

Usage

hv_functionality

Format

An object of class data.table (inherits from data.frame) with 521 rows and 4 columns.

References

Xochelli et al. (2014) Immunoglobulin heavy variable (IGHV) genes and alleles: new entities, new names and implications for research and prognostication in chronic lymphocytic leukaemia. Immunogenetics. 67(1):61-6.

Description

A character vector of all 498 human IGHV germline gene segment alleles in IMGT Gene-db release July 2022, with an additional 25 undocumented alleles from VDJbase.

Usage

HVGERM

Format

Values correspond to IMGT-gaped nuceltoide sequences (with nucleotides capitalized and gaps represented by '.').

References

Xochelli et al. (2014) Immunoglobulin heavy variable (IGHV) genes and alleles: new entities, new names and implications for research and prognostication in chronic lymphocytic leukaemia. Immunogenetics. 67(1):61-6.

Description

Cluster the distance matrix to create allele clusters. Supports both hierarchical clustering (default) and Leiden community detection.

Usage

igClust(
  germline_distance,
  method = c("hierarchical", "leiden"),
  distance_method = "decipher",
  family_threshold = 75,
  allele_cluster_threshold = 95,
  cluster_method = "complete",
  resolution = NULL,
  target_clusters = NULL,
  optimize_silhouette = TRUE,
  ncores = 1,
  quiet = FALSE
)

Arguments

Value

A named list that includes:

• alleleClusterTable: data.frame of allele clusters

• threshold: list of threshold parameters

• hclustAlleleCluster: hierarchical clustering object (both methods)

• communityObject: community detection result (Leiden method)

• graphObject: igraph object (Leiden method)

• silhouetteScore: silhouette score (Leiden method)

• resolutionParameter: resolution used (Leiden method)

See Also

igDistance, inferAlleleClusters

Description

Calculates the distance between pairs of alleles based on their aligned germline sequences. Supports multiple distance methods for different segment types.

Usage

igDistance(
  germline_set,
  AA = FALSE,
  method = c("decipher", "hamming", "lv"),
  trim_3prime = NULL,
  return_type = c("matrix", "dist"),
  quiet = TRUE
)

Arguments

Details

The aligned IMGT IGHV allele germline set can be downloaded from the IMGT site https://www.imgt.org/ under the section genedb.

For V segments, the "decipher" method is recommended as it handles alignment gaps properly. For D and J segments which may have variable lengths, the "lv" (Levenshtein) method is appropriate.

Value

A matrix or dist object of the computed distances between allele pairs.

See Also

ighvDistance for backward compatibility wrapper

Examples

data(HVGERM)
# Using DECIPHER method (default, for V segments)
d1 <- igDistance(HVGERM[1:10], method = "decipher")

# Using Hamming distance
d2 <- igDistance(HVGERM[1:10], method = "hamming")

# Using Levenshtein distance (good for D/J segments)
d3 <- igDistance(HVGERM[1:10], method = "lv")

Description

This function is deprecated. Use igClust instead.

Usage

ighvClust(
  germline_distance,
  family_threshold = 75,
  allele_cluster_threshold = 95,
  cluster_method = "complete"
)

Arguments

Value

A named list with clustering results.

See Also

igClust for the current implementation

Description

This function is deprecated. Use igDistance instead.

Usage

ighvDistance(germline_set, AA = FALSE)

Arguments

Value

A matrix of computed distances between allele pairs.

See Also

igDistance for the current implementation

Description

A wrapper function to infer the allele clusters. Supports both hierarchical clustering (default) and Leiden community detection.

Usage

inferAlleleClusters(
  germline_set,
  locus = NULL,
  clustering_method = c("hierarchical", "leiden"),
  distance_method = c("decipher", "hamming", "lv"),
  trim_3prime_side = 318,
  mask_5prime_side = 0,
  family_threshold = 75,
  allele_cluster_threshold = 95,
  cluster_method = "complete",
  resolution = NULL,
  target_clusters = NULL,
  optimize_silhouette = TRUE,
  ncores = 1,
  aa_set = FALSE,
  quiet = FALSE,
  family_prefix = TRUE,
  retain_subgroup = FALSE
)

Arguments

Details

The distance between pairs of allele sequences is calculated, then the alleles are clustered. For hierarchical clustering, two similarity thresholds define family and allele clusters. For Leiden clustering, community detection identifies clusters at a specified resolution.

The allele cluster names follow this scheme: IGHVF1-G1*01 - IGH = chain, V = region, F1 = family cluster numbering, G1 = allele cluster numbering, 01 = allele numbering (by clustering order)

For V segments, the "decipher" distance method is recommended. For D and J segments with variable lengths, "lv" (Levenshtein) is more appropriate.

Value

An object of class GermlineCluster containing:

• germlineSet: Modified germline set (3' trimming and 5' masking)

• alleleClusterSet: Renamed germline set with ASC names

• alleleClusterTable: data.frame of allele similarity clusters

• threshold: List of threshold parameters

• hclustAlleleCluster: hclust object (both methods)

• clusteringMethod: Method used ("hierarchical" or "leiden")

• communityObject: Community object (Leiden method)

• graphObject: igraph object (Leiden method)

• silhouetteScore: Silhouette score (Leiden method)

• resolutionParameter: Resolution used (Leiden method)

• locus: Locus identifier

See Also

igDistance, igClust, plot.GermlineCluster

Examples

# load the initial germline set

data(HVGERM)

germline <- HVGERM[!grepl("^[.]", HVGERM)]

# Hierarchical clustering (default)
asc <- inferAlleleClusters(germline)

# Leiden community detection
asc_leiden <- inferAlleleClusters(germline[1:50],
                                  clustering_method = "leiden",
                                  target_clusters = 10)

## plotting the clusters
plot(asc)

Description

inferGenotypeAllele infer an individual's genotype based on the allele-base method. The method utilize the allele specific threshold to determine the presence of an allele in the genotype. More specifically, based on the allele frequency, repertoire depth, and the specific allele threshold, a confidence level (Z score) is calculated for the presence of the allele in the genotype. The user can select the confidence level for the genotype inference.

Usage

inferGenotypeAllele(
  data,
  allele_threshold_table = NULL,
  call = "v_call",
  asc_annotation = FALSE,
  single_assignment = FALSE,
  translate_to_asc = FALSE,
  germline_db = NA,
  novel = NA,
  find_unmutated = FALSE,
  seq = "sequence_alignment",
  default_allele_threshold = 1e-04,
  depth_adjusted_threshold = FALSE,
  z_score_threshold = 0,
  quiet = TRUE
)

Arguments

Details

In naive repertoires, allele calls where more than one assignment is assigned is rare. Hence, in case the data represents the naive repertoire of a subject it is recommended to use the find_unmutated=TRUE option, to remove mutated sequences. For non-naive population, the allele calls in cases of multiple assignment are treated as belonging to all groups.

Value

A a data.frame with the inferred V genotype. The table contains the following columns:

• allele: The alleles in the allele_threshold_table.

• counts: The number of reads for each alleles.

• depth: The total number of reads in the genotype (Sum of counts).

• threshold: The population driven allele thresholds for genotype presence.

• z_score: The confidence level for the presence of the allele in the genotype.

• observed: Logical, whether the allele was seen in the data (in single or multiple assignment). Alleles that were never seen carry a smoothing pseudo-count and are marked FALSE.

• in_genotype: Logical, whether the allele passes z_score_threshold (z_score >= z_score_threshold).

• asc_allele: If translate_to_asc is true, the asc allele value from allele_threshold_table.

See Also

inferAlleleClusters will infer the allele clusters based on a supplied V reference set and set the default allele threshold of 1e-04. See recentAlleleClusters to obtain the latest version of the IGHV allele clusters and the naive population based allele threshold.

Examples

# loading TIgGER AIRR-seq b cell data
data <- tigger::AIRRDb

# allele threshold table
data(allele_threshold_table)

data(HVGERM)

# inferring the genotype
genotype <- inferGenotypeAllele(
data = data,
allele_threshold_table = allele_threshold_table,
germline_db = HVGERM, find_unmutated=TRUE)

# filter alleles with z_score >= 0 

head(genotype[genotype$z_score >= 0,])

Description

inferGenotypeAllele_asc infer an individual's genotype based on the allele-base method. The method utilize the allele specific threshold to determine the presence of an allele in the genotype. More specifically, the absolute frequency of each allele is calculated and checked against the threshold.

Usage

inferGenotypeAllele_asc(
  data,
  alleleClusterTable,
  v_call = "v_call",
  single_assignment = FALSE,
  germline_db = NA,
  novel = NA,
  find_unmutated = FALSE,
  seq = "sequence_alignment",
  confidence_level = NULL,
  default_allele_threshold = 1e-04
)

Arguments

Details

In naive repertoires, allele calls where more than one assignment is assigned is rare. Hence, in case the data represents the naive repertoire of a subject it is recommended to use the find_unmutated=TRUE option, to remove mutated sequences. For non-naive population, the allele calls in cases of multiple assignment are treated as belonging to all groups.

Value

A a data.frame with the inferred V genotype. The table contains the following columns:

See Also

inferAlleleClusters will infer the allele clusters based on a supplied V reference set and set the default allele threshold of 1e-04. See recentAlleleClusters to obtain the latest version of the IGHV allele clusters and the naive population based allele threshold.

Examples

# loading TIgGER AIRR-seq b cell data
data <- tigger::AIRRDb

# preferably obtain the latest ASC cluster table
# asc_archive <- recentAlleleClusters(doi="10.5281/zenodo.7429773", get_file = TRUE)

# allele_cluster_table <- extractASCTable(archive_file = asc_archive)

# example allele similarity cluster table
data(allele_cluster_table)

data(HVGERM)

# reforming the germline set
asc_germline <- germlineASC(allele_cluster_table, germline = HVGERM)

# assigning the ASC alleles
asc_data <- assignAlleleClusters(data, allele_cluster_table)

# inferring the genotype
asc_genotype <- inferGenotypeAllele_asc(
data = asc_data,
alleleClusterTable = allele_cluster_table,
germline_db = asc_germline, find_unmutated=TRUE)

Description

This function inserts gaps (e.g., . or -) into an ungapped sequence (ungapped) to match the positions of gaps in a reference sequence (gapped). It ensures that the aligned sequence has the same gap structure as the reference. Vectorised over both arguments, which are paired element by element.

Usage

insert_gaps(gapped, ungapped)

Arguments

Value

A vector of strings with gaps inserted to match the gapped reference.

Examples

gapped <- c("caggtc..aact", "caggtc---aact")
ungapped <- c("caggtcaact", "caggtcaact")

insert_gaps(gapped, ungapped)   # "caggtc..aact", "caggtc---aact"

Deprecated: use insert_gaps

Description

Renamed. The 2 marked a second attempt rather than anything about behaviour, and _vec restated that it is vectorised. Results are unchanged.

Usage

insert_gaps2_vec(gapped, ungapped, parallel = FALSE)

Arguments

Value

As insert_gaps.

See Also

insert_gaps

Description

GermlineCluster is an S3 class that stores the output of inferAlleleClusters. It contains the allele cluster table, clustering objects, and threshold parameters used for inference.

Usage

new_germline_cluster(
  germlineSet,
  alleleClusterSet,
  alleleClusterTable,
  threshold,
  hclustAlleleCluster = NULL,
  clusteringMethod = "hierarchical",
  communityObject = NULL,
  graphObject = NULL,
  distanceMatrix = NULL,
  silhouetteScore = NA_real_,
  resolutionParameter = NA_real_,
  locus = "IGHV",
  .familiesCut = NULL
)

Arguments

Value

An object of class "GermlineCluster".

See Also

inferAlleleClusters

GermlineCluster

Description

Performs a grid search over resolution parameters and selects the one that maximizes the silhouette score.

Usage

optimize_resolution(
  g,
  distance_matrix,
  target_clusters = 80,
  resolution_range_low = 0.1,
  resolution_range_high = 0.5,
  max_steps = 20,
  ncores = 1
)

Arguments

Value

A list containing:

• results: data.frame with Resolution, ClusterCount, Silhouette

• partitions: list of membership vectors for each resolution

• best_resolution: optimal resolution parameter

• best_partition: membership vector at optimal resolution

• best_clusters: number of clusters at optimal resolution

See Also

detect_communities_leiden, igClust

Description

Regresses each single anchor-partner cell on the variant, so a significant omnibus can be read as which partners moved and in which direction. Returns the fitted slope alongside the raw phenotype means in the lowest and highest observed genotype class, and how many subjects sat either side, which is what says whether a difference in means is worth anything.

Usage

pairingCellTests(
  pairs,
  dosage,
  omnibus,
  conditional,
  alpha = 0.05,
  min_subjects = 60L,
  min_complete_fraction = 0.9
)

Arguments

Details

A cell is only meaningful inside a row the omnibus already called significant, so the scan is restricted to those variants and omnibus_significant travels with every row. marked applies nominal alpha inside a cleared row; marked_strict divides alpha by the row width, the family being that row's partner genes. Both are returned so a downstream table and a figure make the same call rather than each re-deriving it.

Value

A data.table with conditional, variant, anchor_gene, partner_gene, n, beta, p_value, n_low, n_high, mean_low, mean_high, delta_mean, omnibus_p_value, omnibus_significant, min_genotype_group, marked and marked_strict.

See Also

pairingScan

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:70)
rep_dt <- data.frame(
  subject = rep(subjects, each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))
pairs <- pairingTable(rep_dt, anchor = "j_gene", partner = "d_gene")
dosage <- matrix(rep(c(0, 1, 2), length.out = 70), nrow = 1,
                 dimnames = list("v1", subjects))

omnibus <- pairingScan(pairs, dosage, conditional = "P(J|D)", min_subjects = 60)
omnibus$significant <- omnibus$p_value < 0.05

pairingCellTests(pairs, dosage, omnibus, conditional = "P(J|D)", min_subjects = 60)

Description

A significant omnibus says a variant is associated with an anchor gene's partner distribution. It does not say the partners were reallocated: the same p-value is produced by a rigid shift of the whole profile, by one unusual subject, and by a variant that simply changes how much the anchor or a partner is used. This runs the follow-ups that separate those, per lead, and returns a verdict.

Usage

pairingLeadCharacter(
  leads,
  pairs,
  dosage,
  positions = NULL,
  ancestry = NULL,
  replication_ancestry = "EUR",
  min_replication_n = 40L,
  min_genotype_group = 5L,
  alpha = 0.05,
  min_sep_d = 0.8,
  max_overlap = 0.7,
  min_mahalanobis = 1,
  min_subjects = 60L,
  min_complete_fraction = 0.9
)

Arguments

Value

A list of two data.tables: leads, the input with the summary columns and verdict merged on, and per_partner, one row per lead and partner with slopes and effect sizes.

What is tested

Rigid shift against reallocation

The mean across partners regressed on dosage is the shift every partner shares (p_common). Removing that mean and testing what is left asks whether partners moved by different amounts (p_heterogeneous). The second is the reallocation claim; the first is not. One deviation column is dropped before the test because deviations from a row mean sum to zero and the covariance would be singular.

Separation

mahalanobis, the distance between the highest and lowest genotype class centroids in units of the pooled within-genotype covariance, after residualising on the per-partner slopes. With overlap_joint beside it.

One subject

p_worst_loo refits the omnibus dropping each subject in turn and keeps the worst p-value. A signal that rests on one person shows up here and nowhere else, and a small genotype class is exactly when it happens.

Both margins

The enrichment divides both marginals out, so a lead is not pairing-specific until neither margin moved. Logit marginal usage is regressed on dosage for the anchor (anchor_usage_p) and for every partner (partner_usage_p_min, Bonferroni across partners).

Ancestry

ancestry_p, whether the dosage itself is stratified, and p_replication, whether the omnibus holds inside one ancestry group. A cis variant is expected to track ancestry, so neither is a rejection on its own.

Position

rho_position, Spearman of per-partner slope against partner position. A monotone trend along the locus is a different claim from a scattered one. Needs positions.

The verdict

verdict is assigned in order: single_subject_driven, overlapping_small_effect, pairing_plus_both_margins, pairing_plus_anchor_usage, pairing_plus_partner_usage, and otherwise pairing_only. Its cutoffs are arguments rather than constants, and their defaults are the ones the heavy-chain D/J analysis used. They are that analysis's policy, not facts about the method: a different cohort or a different pair of segments should set them deliberately rather than inherit them.

See Also

runPairingQTL, pairingCellTests

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:70)
rep_dt <- data.frame(
  subject = rep(subjects, each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))
pairs <- pairingTable(rep_dt, anchor = "j_gene", partner = "d_gene")
dosage <- matrix(rep(c(0, 1, 2), length.out = 70), nrow = 1,
                 dimnames = list("v1", subjects))

leads <- pairingScan(pairs, dosage, conditional = "P(J|D)", min_subjects = 60)
leads$threshold <- 0.05

out <- pairingLeadCharacter(leads, pairs, dosage, min_subjects = 60)
out$leads$verdict
out$per_partner

Description

Subjects by partner gene, of enrichment: the multivariate response pairingScan tests each variant against. Exported because the follow-up analyses that characterise a hit need the same matrix the scan saw, and rebuilding it slightly differently is how a follow-up stops describing the thing it followed up.

Usage

pairingPhenotype(
  pairs,
  anchor,
  subjects,
  min_complete_fraction = 0.9,
  min_subjects = 60L
)

Arguments

Details

Partner genes observed in too few subjects are dropped before complete cases are taken, so one rare partner cannot cost the anchor its whole subject set. Returns NULL when fewer than two partners or fewer than min_subjects subjects survive, which is the signal to skip this anchor rather than an error.

Value

A numeric matrix, subjects in rows (named) and partner genes in columns, or NULL.

See Also

pairingScan

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:70)
rep_dt <- data.frame(
  subject = rep(subjects, each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))
pairs <- pairingTable(rep_dt, anchor = "j_gene", partner = "d_gene")

# Subjects by partner gene, of enrichment: what one anchor is scanned on.
pheno <- pairingPhenotype(pairs, "J1", subjects, min_subjects = 60)
dim(pheno)

Description

One anchor gene at a time, tests every variant against the whole vector of partner enrichments by MANOVA (Pillai's trace, via qtlScanMultivariate). A variant that merely raises or lowers how much the anchor is used shifts every partner in the same direction and produces no signal; one that repartitions the anchor's partners does.

Usage

pairingScan(
  pairs,
  dosage,
  conditional,
  anchors = NULL,
  min_subjects = 60L,
  min_complete_fraction = 0.9,
  min_genotype_group = 5L
)

Arguments

Details

conditional has no default and must be stated. The two anchorings are separate scans over the same variants, so a result that does not say which one it is cannot be pooled, counted or plotted without double-reporting. The direction is easy to invert by mistake and has been inverted once: at fixed partner,

$enrichment = \log_2 P(\mathrm{partner} \mid \mathrm{anchor}) - \log_2 P(\mathrm{partner})$

and the second term is constant across the anchors, so the anchored response vector is $\log_2 P(\mathrm{partner} \mid \mathrm{anchor})$ shifted by a constant. Anchoring on J therefore gives the P(J|D) scan, and anchoring on D gives P(D|J).

min_genotype_group is counted on each anchor's own subject set, since complete cases differ between anchors and a variant can be well powered for one anchor and not another. It is returned on every row, not only on the leads: a p-value from this scan cannot be read without it.

Value

A data.table with conditional, variant, anchor_gene, n, n_response, pillai, f_stat, p_value, min_genotype_group and well_powered. No significance threshold is applied: set it from the number of independent variants (qtlLDGroups), not the number tested.

See Also

pairingTable, pairingCellTests

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:70)
rep_dt <- data.frame(
  subject = rep(subjects, each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))
pairs <- pairingTable(rep_dt, anchor = "j_gene", partner = "d_gene")

dosage <- matrix(rep(c(0, 1, 2), length.out = 70), nrow = 1,
                 dimnames = list("v1", subjects))

# `conditional` has no default: anchoring on J tests P(J|D), anchoring on D
# tests P(D|J), and the two are separate scans over the same variants.
pairingScan(pairs, dosage, conditional = "P(J|D)", min_subjects = 60)

Description

Builds the anchor-by-partner contingency of one subject's repertoire and compares it to what independence of the two marginals would give. This is the phenotype the conditional pairing scan explains.

Usage

pairingTable(data, anchor, partner, subject = "subject", pseudocount = 0.5)

Arguments

Details

Both marginals are divided out. expected is the product of the two margins over the subject's depth, so enrichment measures pairing preference rather than how much either gene is used. That distinction is the whole analysis: in the heavy chain the individual-level signal in the raw conditional p_partner_given_anchor was found to be almost entirely marginal usage, and it went to roughly zero once both marginals were removed. The raw conditionals are returned for inspection, but enrichment is what should be scanned.

Zero cells are created only where both genes occur somewhere in that subject, so a gene the subject does not carry stays a structural absence rather than a measured zero. This matters: a deletion is not the same observation as a pairing the subject never makes.

The enrichment cell is symmetric – swapping which margin conditions leaves the number unchanged – so to scan the other side, call this again with anchor and partner swapped. The omnibus test is not symmetric, because grouping the cells by one segment gives each anchor a response vector of a different width, which is a different multivariate test on a different covariance structure. Both sides are therefore separate scans and their results must not be pooled.

Value

A data.table with one row per subject, anchor gene and partner gene: subject, anchor_gene, partner_gene, count, depth, anchor_total, partner_total, expected, enrichment, p_anchor, p_partner, p_partner_given_anchor, p_anchor_given_partner. The anchor and partner column names are kept on the result as attributes.

Which combinations are established

IGH D-anchored and J-anchored pairing is the analysis this was developed for and the one that has been checked. IGK/IGL V–J, and IGH V–D and V–J, run through the same code and are exploratory: nobody has shown that pairing preference in those segments survives removing the marginals, and the heavy-chain result above is a reason to expect trouble rather than reassurance. Treat their output as a hypothesis.

See Also

pairingScan, pairingCellTests

Examples

set.seed(1)
rep_dt <- data.frame(
  subject = rep(sprintf("s%02d", 1:70), each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))

pairs <- pairingTable(rep_dt, anchor = "j_gene", partner = "d_gene")
head(pairs)

# The enrichment is symmetric, so swapping the roles gives the mirrored table.
# The omnibus that follows is not, which is why the two are separate scans.
mirrored <- pairingTable(rep_dt, anchor = "d_gene", partner = "j_gene")

Description

PIgLET is a suite of computational tools that improves genotype inference and downstream AIRR-seq data analysis. The package as two main tools. The first is Allele Clusters, this tool is designed to reduce the ambiguity within the IGHV alleles. The ambiguity is caused by duplicated or similar alleles which are shared among different genes. The second tool is an allele based genotype, that determined the presence of an allele based on a threshold derived from a naive population.

Allele Similarity Cluster

This section provides the functions that support the main tool of creating the allele similarity cluster form an IGHV germline set.

• inferAlleleClusters: The main function of the section to create the allele clusters based on a germline set.

• ighvDistance: Calculate the distance between IGHV aligned germline sequences.

• ighvClust: Hierarchical clustering of the distance matrix from ighvDistance.

• generateReferenceSet: Generate the allele clusters reference set.

• plotAlleleCluster: Plots the Hierarchical clustering.

• artificialFRW1Germline: Artificially create an IGHV reference set with framework1 (FWR1) primers.

Allele based genotype

This section provides the functions to infer the IGHV genotype using the allele based method and the allele clusters thresholds

• inferGenotypeAllele: Infer the IGHV genotype using the allele based method.

• assignAlleleClusters: Renames the v allele calls based on the new allele clusters.

• germlineASC: Converts IGHV germline set to ASC germline set.

• recentAlleleClusters: Download the most recent version of the allele clusters table archive from zenodo.

• extractASCTable: Extracts the allele cluster table from the zenodo archive file.

• zenodoArchive: An R6 object to query the zenodo api.

References

1. ##

Description

Plot method for GermlineCluster

Usage

## S3 method for class 'GermlineCluster'
plot(x, y = NULL, cex = 1, seed = 9999, ...)

Arguments

Value

A plot of the allele clusters dendrogram

Description

Plotting the dendrogram of the clusters

Usage

plotAlleleCluster(x, y = NULL, cex = 1, seed = 9999)

Arguments

Value

A plot of the allele clusters dendrogram

Description

Creates a comparison visualization showing cluster assignments from both methods.

Usage

plotClusterComparison(hierarchical_result, leiden_result, ...)

Arguments

Value

A ggplot object showing cluster agreement

See Also

inferAlleleClusters

Description

Creates a network visualization of allele clusters from community detection.

Usage

plotCommunityNetwork(
  x,
  layout = c("fr", "kk", "circle"),
  node_color = "cluster",
  node_size = "degree",
  edge_alpha = 0.3,
  show_labels = TRUE,
  label_size = 3,
  ...
)

Arguments

Details

This function creates a network visualization showing:

• Nodes representing alleles, colored by cluster

• Edges weighted by sequence similarity

• Layout optimized by specified algorithm

Value

A ggplot object

See Also

inferAlleleClusters, detect_communities_leiden

Examples

data(HVGERM)
asc <- inferAlleleClusters(HVGERM[1:30],
                           clustering_method = "leiden",
                           target_clusters = 5)
plotCommunityNetwork(asc)

Description

plotGenotypeAllele draws a genotype bar plot in the style of tigger::plotGenotype, with an aligned confidence panel. Unlike TIgGER, which carries a single confidence value per gene, PIgLET assigns a z-score to each allele; the confidence panel therefore mirrors the allele bars and colors each allele segment by its own z-score.

Usage

plotGenotypeAllele(
  genotype,
  level = c("gene", "asc"),
  z_threshold = 0,
  gene_sort = c("name", "position"),
  text_size = 12,
  silent = FALSE,
  ...
)

Arguments

Details

Only the genotyped alleles (those with z_score >= z_threshold) are shown. Each gene/ASC row is split into equal segments, one per genotyped allele; the left panel colors the segments by allele identity and the right panel colors the same segments by their z-score on a continuous blue scale.

Value

A gridExtra grob combining the allele panel and the per-allele z-score confidence panel (returned invisibly).

See Also

genotypeToTigger, inferGenotypeAllele, inferGenotypeAllele_asc

Examples

## Not run: 
data <- tigger::AIRRDb
data(allele_threshold_table)
data(HVGERM)
genotype <- inferGenotypeAllele(data, allele_threshold_table = allele_threshold_table,
                                germline_db = HVGERM, find_unmutated = TRUE)
plotGenotypeAllele(genotype)

## End(Not run)

Description

Creates a plot showing silhouette score and cluster count across resolution values.

Usage

plotSilhouetteOptimization(optimization_result, highlight_best = TRUE, ...)

Arguments

Value

A ggplot object

See Also

optimize_resolution, igClust

Examples

data(HVGERM)
d <- igDistance(HVGERM[1:30], method = "hamming")
g <- distance_to_graph(d)
opt <- optimize_resolution(g, d, target_clusters = 5)
plotSilhouetteOptimization(opt)

Description

Creates a circular or dendrogram tree visualization collapsed to ASC subgroup level, with optional heatmap annotations showing family assignments.

Usage

plotTruncatedTree(
  x,
  layout = c("circular", "dendrogram"),
  collapse_to = c("asc_subgroup", "iuis_subgroup", "family"),
  label_style = c("asc", "iuis", "both"),
  show_threshold_line = TRUE,
  threshold = 0.25,
  tip_size_by = "n_alleles",
  tip_color_by = "present",
  show_heatmap = TRUE,
  label_size = 7,
  ...
)

Arguments

Details

This function creates a publication-quality tree visualization that:

• Renames tree tips from original allele names to ASC names (new_allele)

• Collapses alleles to ASC subgroup level (single representative per ASC group)

• Shows tip point size by number of alleles in cluster

• Adds optional heatmap track showing IUIS vs ASC family assignments

• Draws threshold line at specified height

When using label_style = "iuis", if multiple ASC groups split a single IUIS subgroup, the labels are marked with superscript letters (e.g., IGHV1-2^A, IGHV1-2^B) to distinguish them.

Requires the ggtree package to be installed.

Value

A ggplot/ggtree object

See Also

inferAlleleClusters, plot.GermlineCluster

Examples

data(HVGERM)
asc <- inferAlleleClusters(HVGERM[1:50])

# Basic truncated tree with ASC labels
if (requireNamespace("ggtree", quietly = TRUE)) {
  plotTruncatedTree(asc, show_heatmap = FALSE)

  # With IUIS labels (marked if ASC splits IUIS group)
  plotTruncatedTree(asc, label_style = "iuis", show_heatmap = FALSE)
}

Description

Print method for GermlineCluster

Usage

## S3 method for class 'GermlineCluster'
print(x, ...)

Arguments

Value

Invisibly returns x

Description

Takes the strongest signal, absorbs everything correlated with it above r2, and repeats. What comes back is one row per independent signal rather than one row per significant variant, which is the number to quote.

Usage

qtlClump(assoc, dosage, r2 = 0.8)

Arguments

Details

Ties are broken on the variant name, not on the p-value alone. Variants in perfect linkage disequilibrium give the same fit to the last bit that the arithmetic carries, and their p-values then differ by around 1e-13 purely from the order of operations in the matrix product. Sorting on the p-value alone lets that noise decide which variant is named as the lead, and the lead variant name is what reaches a figure and a database label. So the sort is on the p-value and then the name, which is still an arbitrary choice between variants no data distinguishes, but the same arbitrary choice every run. A lead should be read as one representative of its linkage group.

Value

The subset of assoc that are leads, in increasing p-value order.

Examples

dosage <- rbind(lead   = c(0, 1, 2, 0, 1, 2, 1, 0),
                linked = c(0, 1, 2, 0, 1, 2, 1, 0),
                other  = c(0, 0, 1, 2, 2, 1, 0, 1))
colnames(dosage) <- sprintf("s%d", 1:8)

assoc <- data.table::data.table(
  variant = c("lead", "linked", "other"),
  p_value = c(1e-8, 2e-8, 1e-4))

# `linked` is absorbed into `lead`; two independent signals remain.
qtlClump(assoc, dosage, r2 = 0.8)

Description

Two variants carry the same information when their standardised dosage vectors match up to sign, so the mirrored case is folded in. This is exact linkage disequilibrium rather than a correlation threshold; use qtlClump for the latter.

Usage

qtlLDGroups(dosage, contig)

Arguments

Details

The number of distinct groups is the number of independent tests, which is what a Bonferroni threshold should be set from – the number of variants tested overstates it, often by half. The contig prefix stops identical patterns at different loci from merging, so positions from different loci are never pooled.

Value

A character vector named by variant, one group key per row of dosage. The key is opaque: use it to group and count, not to read. Variants with no variance each get their own group rather than collapsing together.

Examples

dosage <- rbind(a      = c(0, 1, 2, 0, 1, 2),
                copy   = c(0, 1, 2, 0, 1, 2),   # identical
                mirror = c(2, 1, 0, 2, 1, 0),   # identical up to sign
                other  = c(0, 0, 1, 2, 2, 1))
colnames(dosage) <- sprintf("s%d", 1:6)

groups <- qtlLDGroups(dosage, rep("chr2", 4))
length(unique(groups))   # three independent tests, not four

Description

Scores every variant against a multivariate response, one MANOVA per variant. For a one degree of freedom predictor Pillai's trace has a closed form,

$V = (u' S^{-1} u) / (g'g), \quad u = Y_c' g_c, \quad S = Y_c' Y_c$

and the F approximation is exact, so a whole missingness group is scored by one matrix product instead of a stats::manova() call per variant. The package tests check this against stats::manova().

Usage

qtlScanMultivariate(pheno, dosage, min_n = 60L)

Arguments

Details

Response columns with zero variance are dropped within each missingness group, and a group is skipped when fewer than two response columns survive or when n is too small for the covariance to be estimable (n < p + 3).

Value

A data.table with variant, n, n_response (response columns retained), pillai, f_stat and p_value.

See Also

pairingScan, which uses this for the conditional pairing scan

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:80)
dosage <- rbind(v_hit  = rep(c(0, 1, 2), length.out = 80),
                v_null = rep(c(0, 0, 1, 2), length.out = 80))
colnames(dosage) <- subjects
pheno <- matrix(rnorm(80 * 4), 80, 4,
                dimnames = list(subjects, paste0("partner", 1:4)))
pheno[, 1] <- pheno[, 1] + 0.8 * dosage["v_hit", ]

qtlScanMultivariate(pheno, dosage, min_n = 60)

Description

Ordinary least squares of each phenotype column on each variant, for the whole matrix at once. The slope of a simple regression is a ratio of cross products, so the entire scan is two matrix products per missingness group rather than one lm() call per variant-phenotype pair.

Usage

qtlScanUnivariate(pheno, dosage, min_n = 60L)

Arguments

Details

Variants are grouped by which subjects have no call, so every fit inside a group uses the same complete-case subject set and n is exact for each row rather than an upper bound. Groups smaller than min_n are skipped entirely.

A p-value from this scan must not be read without the smallest genotype class behind the fit: the extreme tail is anti-conservative, and a non-significant result means either no effect or no power. Use qtlSmallestGenotypeClass on the same dosage matrix, as runGeneUsageQTL does. Filtering associations on the p-value alone preferentially removes the underpowered rows, which leaves what survives looking better powered than it is.

Value

A data.table with one row per variant-phenotype pair that produced a finite p-value: variant, phenotype, n (subjects in the fit), beta, se, t_stat, p_value. Effect direction is relative to the coding of dosage, which is not necessarily the minor allele, so beta must not be described in terms of a named allele without resolving polarity first.

See Also

qtlScanMultivariate for a multivariate response, qtlSmallestGenotypeClass, runGeneUsageQTL

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:80)
dosage <- rbind(v_hit  = rep(c(0, 1, 2), length.out = 80),
                v_null = rep(c(0, 0, 1, 2), length.out = 80))
colnames(dosage) <- subjects
pheno <- cbind(gene_a = 0.7 * dosage["v_hit", ] + rnorm(80),
               gene_b = rnorm(80))
rownames(pheno) <- subjects

qtlScanUnivariate(pheno, dosage, min_n = 60)

# Read a p-value with the smallest genotype class beside it, never alone.
qtlSmallestGenotypeClass(dosage)

Description

How many subjects sit in the least populated genotype class at each variant. This is what says whether a fit is worth reading: an association driven by two homozygotes is not the same claim as one driven by forty, and the scan's extreme tail is anti-conservative precisely where this number is small.

Usage

qtlSmallestGenotypeClass(dosage)

Arguments

Details

Classes with no subject are not counted, so a variant seen only as 0 and 1 reports the smaller of those two rather than zero. Computed for the whole matrix at once; index the result by variant name rather than calling this per variant.

Value

A named numeric vector, one entry per row of dosage, NA where no class had a subject.

Examples

dosage <- rbind(all_three = c(0, 0, 0, 1, 1, 2),
                two_only  = c(0, 0, 0, 0, 1, 1))
colnames(dosage) <- sprintf("s%d", 1:6)

# The empty class is not counted as zero: two_only reports 2, not 0.
qtlSmallestGenotypeClass(dosage)

Description

A wrapper function for zenodoArchive, download the most recent allele similarity clusters and thresholds from the zenodo archive. The clusters and thresholds are based on https://yaarilab.github.io/IGHV_reference_book/ At the moment only available for human IGHV reference set.

Usage

recentAlleleClusters(
  doi = "10.5281/zenodo.7401189",
  path,
  get_file = FALSE,
  quite = FALSE
)

Arguments

Value

If get_file is TRUE, the function returns the path to the archive file

Examples

## Not run: 
recentAlleleClusters(doi="10.5281/zenodo.7401189")

## End(Not run)

Description

Post-processes an ASC assignment (e.g. the alleleClusterTable produced by inferAlleleClusters) by re-grouping alleles according to a per-allele label – typically the IMGT gene held in imgt_allele. It corrects the two ways an ASC grouping and a gene-level view disagree: an ASC group that merges several distinct genes, and a gene whose alleles are split across ASC groups. The clustering itself is not re-run; only the output is relabelled, so this drops onto the end of the pipeline without changing anything upstream.

Usage

regroupASCByLabel(
  x,
  group_col = "allele_cluster",
  label_col = "imgt_allele",
  duplicated_col = "duplicated",
  action = c("both", "split", "merge")
)

Arguments

Details

The label is not produced by the clustering; supply it on the table (the IMGT gene in imgt_allele, and any other genes that share the allele in duplicated). A gene is taken as the part of a label before *, so allele suffixes and mutation tags (IGHD1-20*01_a7g) are ignored.

Genes are grouped within a locus. A gene name is unique only inside its locus, so a label stripped of its prefix ("J1" rather than "IGHJ1") would otherwise put IGHJ1, IGKJ1 and IGLJ1 in one bin. The locus is read from the label where it still carries one and from group_col otherwise; if neither supplies it, the grouping falls back to plain gene names and warns. action = "split" is unaffected, since its output is prefixed by the ASC group.

Value

x with an added character column regrouped. A group made of several genes (a merge, e.g. a shared allele) is tagged with all of them joined by /, e.g. "IGHD2-2/IGHD2-8", so no gene name is dropped.

See Also

inferAlleleClusters

Examples

tb <- data.frame(
  allele         = paste0("a", 1:5),
  allele_cluster = c("G1", "G1", "G2", "G3", "G3"),
  imgt_allele    = c("IGHD1-1*01", "IGHD1-2*01", "IGHD1-1*02",
                     "IGHD2-2*01", "IGHD2-8*01"),
  duplicated     = c("", "", "", "IGHD2-8*01", ""),   # a4 shared: 2-2 and 2-8
  stringsAsFactors = FALSE)
regroupASCByLabel(tb, action = "both")$regrouped

Description

Runs the whole usage analysis in memory: builds the phenotype, scans every variant against every gene group, sets the significance threshold from the number of independent variants, and clumps the significant hits into leads.

Usage

runGeneUsageQTL(
  data,
  dosage,
  variants,
  segments,
  positions = NULL,
  locus = NA_character_,
  subject = "subject",
  min_subject_fraction = 0.1,
  pseudocount = 0.5,
  min_subjects = 60L,
  alpha = 0.05,
  r2_clump = 0.8,
  min_genotype_group = 5L,
  cis_window = 50000
)

Arguments

Details

The threshold is alpha over the number of independent variants after exact-LD collapse (qtlLDGroups), not over the number tested. Those differ by about half in this data, and using the number tested is the more conservative but wrong denominator.

Power columns are attached to the leads, where igqtl.R put them, and the counts behind every step are messaged: a filter that drops variants without saying so is a bug even when it returns cleanly.

Value

A list of data.tables: phenotype, associations, leads, per_asc and thresholds.

Ceiling

The scan is unweighted OLS on usage fractions, so its standard errors assume constant variance while binomial sampling variance across subjects spans roughly twentyfold. Any R-squared read off it is diluted and is an underestimate. Do not report an effect size from this as exact.

See Also

qtlScanUnivariate, runPairingQTL

Examples

set.seed(1)
n <- 70L
subjects <- sprintf("s%02d", seq_len(n))
dosage <- rbind(v_hit  = rep(c(0, 1, 2), length.out = n),
                v_null = rep(c(0, 0, 1, 2), length.out = n))
colnames(dosage) <- subjects

genes <- c("V1-2", "V3-23", "V4-34")
rep_dt <- do.call(rbind, lapply(seq_len(n), function(i) {
  w <- c(1 + dosage["v_hit", i], 1, 1)   # dosage raises V1-2 usage
  data.frame(subject = subjects[i],
             v_gene = sample(genes, 40, replace = TRUE, prob = w / sum(w)))
}))

variants <- data.frame(variant = rownames(dosage), contig = "igh",
                       pos = c(1000L, 50000L), maf = rowMeans(dosage) / 2)

res <- runGeneUsageQTL(rep_dt, dosage, variants,
                       segments = c(V = "v_gene"), min_subjects = 60)
res$associations
res$thresholds

Description

Runs one side of the pairing analysis: builds the enrichment table, scans every variant against each anchor's partner distribution, sets the threshold from the number of independent variants, and clumps the hits into leads.

Usage

runPairingQTL(
  data,
  dosage,
  variants,
  anchor,
  partner,
  conditional,
  subject = "subject",
  locus = NA_character_,
  pseudocount = 0.5,
  min_subjects = 60L,
  min_complete_fraction = 0.9,
  alpha = 0.05,
  r2_clump = 0.8,
  min_genotype_group = 5L
)

Arguments

Details

One anchoring per call, on purpose. Anchoring on each of the two segments gives two different multivariate tests over the same variants, and a function that ran both and returned them stacked would invite exactly the pooling that double-reports. Call it twice, with the conditional each side estimates.

Value

A list of data.tables: pairs, associations, leads and thresholds.

See Also

pairingTable, pairingScan, pairingCellTests

Examples

set.seed(1)
subjects <- sprintf("s%02d", 1:70)
rep_dt <- data.frame(
  subject = rep(subjects, each = 20),
  d_gene  = sample(c("D1", "D2", "D3"), 1400, replace = TRUE),
  j_gene  = sample(c("J1", "J2"), 1400, replace = TRUE))

dosage <- matrix(rep(c(0, 1, 2), length.out = 70), nrow = 1,
                 dimnames = list("v1", subjects))
variants <- data.frame(variant = "v1", contig = "igh", pos = 1000L, maf = 0.33)

# One anchoring per call. Running both and stacking them would double-report,
# so the caller states which conditional each side estimates.
jd <- runPairingQTL(rep_dt, dosage, variants,
                    anchor = "j_gene", partner = "d_gene",
                    conditional = "P(J|D)", min_subjects = 60)
jd$associations

Description

Summary method for GermlineCluster

Usage

## S3 method for class 'GermlineCluster'
summary(object, ...)

Arguments

Value

A list with summary statistics

Description

zenodoArchive

zenodoArchive

Format

R6Class object.

Value

Object of R6Class for modelling an zenodoArchive for ASC cluster files

Public fields

Methods

Public methods

• zenodoArchive$new()

• zenodoArchive$clean_doi()

• zenodoArchive$zenodo_query()

• zenodoArchive$get_versions()

• zenodoArchive$get_version_files()

• zenodoArchive$download_zenodo_files()

• zenodoArchive$clone()

Method new()

initializes the zenodoArchive

Usage

zenodoArchive$new(
  doi,
  page = 1,
  size = 20,
  all_versions = "true",
  sort = "mostrecent"
)

Arguments

Method clean_doi()

cleans the doi record for query

Usage

zenodoArchive$clean_doi(doi = self$doi)

Arguments

Returns

the clean doi

Method zenodo_query()

Query the zenodo archive according to the initial parameters.

Usage

zenodoArchive$zenodo_query(...)

Arguments

Returns

a list with the query values.

Method get_versions()

Extract all concept doi available versions.

Usage

zenodoArchive$get_versions(...)

Arguments

Returns

a data.frame of the available versions.

Method get_version_files()

get the chosen doi archive version available files

Usage

zenodoArchive$get_version_files(version = "latest")

Arguments

Returns

a list of the available files in the archive version.

Method download_zenodo_files()

get the chosen doi archive version available files

Usage

zenodoArchive$download_zenodo_files(
  file = NULL,
  path = tempdir(),
  version = "latest",
  all_files = F,
  get_file_path = F,
  quite = F
)

Arguments

Returns

If get_file_path is TRUE, the function returns the path to the archive file

Method clone()

The objects of this class are cloneable with this method.

Usage

zenodoArchive$clone(deep = FALSE)

Arguments

Examples

## Not run: 
  zenodo_archive <- zenodoArchive$new(
     doi = "10.5281/zenodo.7401189"
  )

  # view available version ins the archive
  archive_versions <- zenodo_archive$get_versions()

  # Getting the available files in the latest zenodo archive version
  files <- zenodo_archive$get_version_files()

  # downloading the first file from the latest archive version
  zenodo_archive$download_zenodo_files()

## End(Not run)