## ----setup, include=FALSE----------------------------------------------------- library(BiocStyle) library(igblastr) reset_germline_dbs() reset_c_region_dbs() ## ----eval=FALSE--------------------------------------------------------------- # if (!require("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # BiocManager::install("igblastr") ## ----message=FALSE------------------------------------------------------------ library(igblastr) ## ----------------------------------------------------------------------------- if (!has_igblast()) install_igblast() ## ----------------------------------------------------------------------------- igblast_info() ## ----------------------------------------------------------------------------- list_germline_dbs(preinstalled.only=TRUE) ## ----------------------------------------------------------------------------- head(list_IMGT_releases()) ## ----------------------------------------------------------------------------- list_IMGT_organisms("202631-1") ## ----message=FALSE------------------------------------------------------------ install_IMGT_germline_db("202631-1", organism="Homo sapiens") ## ----------------------------------------------------------------------------- use_germline_db("IMGT-202631-1.Homo_sapiens.IGH+IGK+IGL") ## ----------------------------------------------------------------------------- list_germline_dbs() ## ----------------------------------------------------------------------------- head(load_intdata("human")) head(load_auxdata("human")) ## ----------------------------------------------------------------------------- head(load_intdata("IMGT-202631-1.Homo_sapiens.IGH+IGK+IGL")) head(load_auxdata("IMGT-202631-1.Homo_sapiens.IGH+IGK+IGL")) ## ----------------------------------------------------------------------------- list_c_region_dbs() ## ----------------------------------------------------------------------------- use_c_region_db("_IMGT.human.IGH+IGK+IGL.202605") ## ----------------------------------------------------------------------------- query <- system.file(package="igblastr", "extdata", "BCR", "1279067_1_Paired_sequences.fasta.gz") ## ----------------------------------------------------------------------------- json <- system.file(package="igblastr", "extdata", "BCR", "1279067_1_Paired_All.json") query_metadata <- jsonlite::fromJSON(json) query_metadata ## ----------------------------------------------------------------------------- use_germline_db() use_c_region_db() ## ----------------------------------------------------------------------------- AIRR_df <- igblastn(query, num_alignments_V=1) ## ----------------------------------------------------------------------------- AIRR_df ## ----message=FALSE------------------------------------------------------------ library(ggplot2) ## ----------------------------------------------------------------------------- AIRR_df |> ggplot(aes(locus, 100 - v_identity)) + theme_bw(base_size=14) + geom_point(position = position_jitter(width = 0.3), alpha = 0.1) + geom_boxplot(color = "blue", fill = NA, outliers = FALSE, alpha = 0.3) + ggtitle("Distribution of V percent mutation by locus at the nucleotide level") + xlab(NULL) ## ----------------------------------------------------------------------------- perc_mut_aa <- percent_mutation(AIRR_df, for.aa=TRUE) head(perc_mut_aa) ## ----------------------------------------------------------------------------- perc_mut_aa |> ggplot(aes(locus, v_perc_mut_aa)) + theme_bw(base_size=14) + geom_point(position = position_jitter(width = 0.3), alpha = 0.1) + geom_boxplot(color = "blue", fill = NA, outliers = FALSE, alpha = 0.3) + ggtitle("Distribution of V percent mutation by locus at the amino acid level") + xlab(NULL) ## ----message=FALSE------------------------------------------------------------ library(dplyr) library(scales) ## ----------------------------------------------------------------------------- plot_gene_dist <- function(AIRR_df, loc) { df_v_gene <- AIRR_df |> filter(locus == loc) |> mutate(v_gene = allele2gene(v_call)) |> # drop allele info group_by(v_gene) |> summarize(n = n(), .groups = "drop") |> mutate(frac = n / sum(n)) df_v_gene |> ggplot(aes(frac, v_gene)) + theme_bw(base_size=13) + geom_col() + scale_x_continuous('Percent of sequences', labels = scales::percent) + ylab("Germline gene") + ggtitle(paste0(loc, "V gene prevalence")) } ## ----fig.height=8------------------------------------------------------------- plot_gene_dist(AIRR_df, "IGH") ## ----fig.height=5.9----------------------------------------------------------- plot_gene_dist(AIRR_df, "IGK") ## ----fig.height=5.3----------------------------------------------------------- plot_gene_dist(AIRR_df, "IGL") ## ----message=FALSE------------------------------------------------------------ library(ggseqlogo) ## ----fig.height=4.5----------------------------------------------------------- AIRR_df$cdr3_aa_length <- nchar(AIRR_df$cdr3_aa) AIRR_df |> group_by(locus, cdr3_aa_length) |> summarize(n = n(), .groups = "drop") |> ggplot(aes(cdr3_aa_length, n)) + theme_bw(base_size=14) + facet_wrap(~locus) + geom_col() + ggtitle("Histograms of CDR3 length by locus") ## ----------------------------------------------------------------------------- AIRR_df |> filter(locus == "IGK", cdr3_aa_length == 9) |> pull(cdr3_aa) |> ggseqlogo(method = "probability") + theme_bw(base_size=14) + ggtitle("Logo plot of kappa chain CDR3 sequences that are 9 AA long") ## ----------------------------------------------------------------------------- igblastn_help() ## ----eval=FALSE--------------------------------------------------------------- # V_alleles <- c("IGHV3-23*01", "IGHV3-23*04") # igblastn(query, germline_db_V_seqidlist=V_alleles) ## ----eval=FALSE--------------------------------------------------------------- # igblastn(query, germline_db_V_seqidlist=file("path/to/my_V_gene_alleles.txt")) ## ----------------------------------------------------------------------------- fasta_file <- "Influenza-KV7015_heavy.fasta" query <- system.file(package="igblastr", "extdata", "BCR", "KyDab", fasta_file) ## ----------------------------------------------------------------------------- db_name1 <- "IMGT-202631-1.Homo_sapiens.IGH+IGK+IGL" db_name2 <- install_IMGT_germline_db("202631-1", "Mus_musculus") db_name <- "comb.IMGT-202631-1.human+mouse.IGH+IGK+IGL" combine_germline_dbs(db_name, db_name1, db_name2, suffix1="_Hs", suffix2="_Mm") ## ----------------------------------------------------------------------------- use_germline_db(db_name) list_germline_dbs() ## ----------------------------------------------------------------------------- AIRR_df <- igblastn(query, num_alignments_V=1, num_alignments_D=1, num_alignments_J=1) ## ----------------------------------------------------------------------------- # extract_suffix() assumes suffixes of length 3, e.g. "_Hs"/"_Mm": extract_suffix <- function(x) {nc <- nchar(x); factor(substr(x, nc-2, nc), levels=c("_Hs", "_Mm"))} table(extract_suffix(AIRR_df$v_call)) table(extract_suffix(AIRR_df$d_call)) table(extract_suffix(AIRR_df$j_call)) ## ----echo=FALSE, fig.cap="Number of human vs. mouse V/D/J allele assignments
for the 3,950 Kymouse heavy-chain sequences
in Influenza-KV7015_heavy.fasta", fig.align="center"---- library(ggplot2) plot_human_vs_mouse_assignments <- function(AIRR_df) { is_human <- function(allele_name) igblastr:::has_suffix(allele_name, "_Hs") is_mouse <- function(allele_name) igblastr:::has_suffix(allele_name, "_Mm") gene_type <- factor(rep(c("V", "D", "J"), each=2), levels=c("V", "D", "J")) allele_organism <- rep(c("human (_Hs)", "mouse (_Mm)"), 3) counts <- c(sum(is_human(AIRR_df$v_call)), sum(is_mouse(AIRR_df$v_call)), sum(is_human(AIRR_df$d_call), na.rm=TRUE), sum(is_mouse(AIRR_df$d_call), na.rm=TRUE), sum(is_human(AIRR_df$j_call)), sum(is_mouse(AIRR_df$j_call))) labs <- counts counts[counts == 0L] <- max(counts) * 0.003 data <- data.frame(`gene type`=gene_type, `allele organism`=allele_organism, `allele assignments`=counts, lab=labs, check.names=FALSE) ggplot(data, aes(x=`gene type`, y=`allele assignments`, fill=`allele organism`)) + geom_bar(position="dodge", stat="identity") + coord_cartesian(ylim=c(0, max(counts) * 1.15)) + theme(axis.text=element_text(face="bold")) + geom_text(data=data, aes(x=`gene type`, y=`allele assignments`, label=lab), position=position_dodge(width=0.9), vjust=-1) } plot_human_vs_mouse_assignments(AIRR_df) ## ----------------------------------------------------------------------------- filename <- "SRR11341217.fasta.gz" query <- system.file(package="igblastr", "extdata", "TCR", filename) ## ----warning=FALSE, message=FALSE--------------------------------------------- db_name <- install_IMGT_germline_db("202631-1", organism="Homo sapiens", tcr.db=TRUE) ## ----------------------------------------------------------------------------- list_germline_dbs() ## ----message=FALSE------------------------------------------------------------ use_germline_db(db_name) ## ----------------------------------------------------------------------------- list_c_region_dbs() ## ----message=FALSE------------------------------------------------------------ use_c_region_db("_IMGT.human.TRA+TRB+TRG+TRD.202605") ## ----------------------------------------------------------------------------- use_germline_db() use_c_region_db() ## ----------------------------------------------------------------------------- AIRR_df <- igblastn(query) AIRR_df ## ----------------------------------------------------------------------------- igblast_info() ## ----------------------------------------------------------------------------- sessionInfo()