Vgene: Usage guide

Introduction

Vgene is an R package for analyzing bulk targeted DNA sequencing data of B-cell receptors (BCR) generated from the Invivoscribe LymphoTrack IGH FR1 assay kit. The assay targets the conserved framework region 1 (FR1) of the immunoglobulin heavy chain (IGH) locus, spanning the full VH and part of the JH region, and is sequenced on Illumina platforms such as the MiSeq. Therefore, this makes vgene applicable where clonality and somatic hypermutation (SHM) status are important, for example in the study of B-cell malignancies. In particular, vgene takes the raw paired-end FASTQ and performs: quality filtering, primer detection and JH trimming, VH-family clustering, BLAST-based germline assignment, mutation rate calculation, clone clustering, and per clone lineage tree reconstruction.

For each clone it reports in output_clones.tsv:

Output	Column	Description
Clonality	clone_id, sample_type	Clone identity and sample of origin
Somatic hypermutation	v_identity	V-gene identity to germline
Frame status	stop_codon	Whether protein translation is in or out of frame
Segment annotation	—	VD, CDR3 and other per-segment annotations

And phylogeny plots per clone are present in in the lineage_plots folder

Installation

# install.packages("remotes")
remotes::install_gitlab(
  "nfdi4immuno/tools/v-gene",
  host = "https://codebase.helmholtz.cloud"
)

External tools to be installed by the user

The following command-line tools must be installed and on your PATH before running the pipeline:

Tool	How to verify
NCBI BLAST+	blastn -version
GNU parallel	parallel --version
PHYLIP dnapars	dnapars

---

Input Files

Two inputs are required.

1) Samples file

A .csv file named samples.csv in the working directory with one row per sample. Three columns are required:

number,name,folder
1,Patient_77_PB-CSW,Patient_77
2,Patient_77_PB-MD27,Patient_77
3,Patient_78_PB-CSW,Patient_78
4,Patient_78_PB-MD27,Patient_78

Column	Type	Description	Example
number	Integer	Index of the sample	1
name	Character	Sample name	Patient_77_PB-CSW
folder	Character	Donor or group label	Patient_77

Multiple samples can share the same folder value when they come from the same donor. All samples with the same folder are analysed together for BLAST, mutation analysis, and clone clustering.

Sample naming convention. The B-cell sample name for each sample is derived automatically by separating the folder prefix from the name. For example, a sample named Patient_77_PB-CSW in folder Patient_77 yields the sample name PB-CSW. Keep sample names consistent with this pattern.

2) FASTQ files

Paired-end FASTQ files (*.fastq.gz) should all be in a single directory. The main function run_pipeline() creates an R1_and_R2/ folder inside the working directory and copies them there automatically.

---

Running the Pipeline

Here is an example run. Replace each path with the one that matches your own files and installations.

1) Default run

The default run uses the required parameters supplied by the user. The rest of the parameters have a default value that can be customized by the user.

library(vgene)

run_pipeline(
  workdir       = "v-gene",
  samples_file  = "samples.csv",
  fastq_dir     = "R1_and_R2",
  ighv_path     = system.file("extdata", "makeblastdb",
                               "IGHV_reference.fasta", package = "vgene"),
  blast_dir     = "ncbi-blast-2.16.0/bin",
  phylip_exec   = "/opt/phylip-3.697/exe/dnapars"
)

Arguments

Parameter	Description	Example
workdir	Output directory (created if absent)	"v-gene"
samples_file	Path to samples CSV	"samples.csv"
fastq_dir	Directory containing input FASTQ files	"R1_and_R2"
blast_dir	Directory containing BLAST+ binaries	"ncbi-blast-2.16.0/bin"
phylip_exec	Full path to PHYLIP dnapars executable	"/opt/phylip-3.697/exe/dnapars"

2) Custom run

Override only the parameters you need to change.

library(vgene)

run_pipeline(
  # --- Required parameters ---
  workdir       = "v-gene",
  samples_file  = "samples.csv",
  fastq_dir     = "R1_and_R2",
  ighv_path     = system.file("extdata", "makeblastdb",
                               "IGHV_reference.fasta", package = "vgene"),
  blast_dir     = "ncbi-blast-2.16.0/bin",
  phylip_exec   = "/opt/phylip-3.697/exe/dnapars",
  
  # --- Custom parameters and their defaults ---
  
  # Relax quality thresholds, allow shorter reads ---
  qual_mean_cutoff = 25,
  qual_nuc_cutoff  = 30,
  min_length       = 150,
  max_length       = 250,

  # Use another primer list ---
  primers = c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA", "CTCACCTG", "CTGTAAGG", "ACTCACCT"),

  # Parallize the jobs
  multicores = NULL, 

  # Lineage tree paths ---
  vh_diff_subdir                = file.path("VH_sorted", "diff"), 
  igblast_output_path           = "airr_igblast_output.tsv",
  output_vgene_clone_airr_file  = "output_clones.tsv"
)

Arguments

Parameter	Default	Description
qual_mean_cutoff	20	Minimum mean read quality (Phred)
qual_nuc_cutoff	35	Per-nucleotide quality threshold (Phred)
min_length	200	Minimum read length after trimming (bp)
max_length	200	Maximum read length to be used for quality control (bp)
primers	c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA", "CTCACCTG", "CTGTAAGG", "ACTCACCT")	Invivoscribe primer list
multicores	NULL	Use parallel processing
vh_diff_subdir	"VH_sorted/diff"	Path to clones with shared sequences between the samples
igblast_output_path	airr_igblast_output.tsv	Path to clones with shared sequences between the samples
output_vgene_clone_airr_file	"output_clones.tsv"	Name of annotated output file

---

Re-running a Single Step

Every step function is exported, so you can run the functions individually as well. The step functions are:

# Step 1
process_fastq_pipeline(
  input_dir        = "R1_and_R2",
  output_dir       = "R1_and_R2_corrected",
  qual_mean_cutoff = 20,
  qual_nuc_cutoff  = 35,
  min_length       = 200,
  max_length       = 200,
  logfile          = "log.txt"
)

# Step 2
find_primers(
  input_dir  = "R1_and_R2_corrected/",
  output_dir = "With_primer/",
  primers    = c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA",
                 "CTCACCTG", "CTGTAAGG", "ACTCACCT")
)

# Step 3
trim_vh_files_by_jh(
  primers = c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA",
              "CTCACCTG", "CTGTAAGG", "ACTCACCT")
)

# Step 4 — create log directory first
if (!dir.exists("logs/cluster_logs")) {
  dir.create("logs/cluster_logs", recursive = TRUE)
}
system(paste0(
"parallel -j ", n_parallel_jobs,
  " --results logs/cluster_logs --quote Rscript -e",
  " \"renv::load(); vgene::run_cluster_all(",
  "VH='{}',",
  " out_dir='With_primer/',",
   " multicores=", multicores, ")\"",
  " ::: ", paste(vh_list, collapse = " ")
))

# Step 5
analyze_clustered_fasta(input_dir = "Collapsed", filter_number = filter_number)

# Step 6
samples <- read.csv("samples.csv", stringsAsFactors = FALSE)
combine_sample_sequences(
  input_dir = "Collapsed/new_more_one/",
  samples   = samples,
  verbose   = TRUE,
  recursive = TRUE
)

folders <- unique(read.csv("samples.csv")$folder)

# Step 7
get_antibody_family_info_via_blast(
  blast_dir = "ncbi-blast-2.16.0/bin",
  path_IGHV = system.file("extdata", "makeblastdb",
                           "IGHV_reference.fasta", package = "vgene"),
  folders   = folders,
  verbose   = TRUE
)

# Step 8
process_folders_collapse_alleles(folders)

# Step 10
clean_mutation_rate_files(folders)

# Step 11
split_mutation_rate_by_sample_type(folders)

# Step 12
reduce_similar_sequences_in_no_diff(folders)

# Step 13
reduce_similar_sequences_in_diff(folders)

# Step 14
run_igblast_for_folders(folders)

# Step 15
typesB  <- unique(unname(mapply(
  function(name, folder) sub(paste0(folder, "_"), "", name),
  samples$name, samples$folder
)))

run_lineage_tree_pipeline(
  folders,
  typesB,
  vh_diff_subdir               = file.path("VH_sorted", "diff"),
  phylip_exec                  = "/opt/phylip-3.697/exe/dnapars",
  igblast_output_path          =  "airr_igblast_output.tsv",
  output_vgene_clone_airr_file = "output_clones.tsv"
)

Output Folder Structure

After a successful run, the working directory contains the tree shown below.Results are written to the directory given in workdir. Inside it, the pipeline creates one folder for each unique folder value in your samples file. Clone information and annotations are written to a output_clones.tsv in each folder, and the lineage plots for each clone are saved in the lineage_plotssubfolder.

v-gene/
├── R1_and_R2/                          # copies of input FASTQ 
├── R1_and_R2_corrected/                # QC-filtered reads
├── With_primer/                        # primer-positive, JH-trimmed reads
├── Collapsed/                          # clustered sequences
│   └── new_more_one/                   # clusters with more than one member
├── plots/                              # per-patient summary plots
├── table/                              # clone summary statistics across all patients
├── logs/
│   └── cluster_logs/                   # GNU parallel logs from step 4
├── Patient_1/                          # one folder per patient (folder column in CSV)
│   ├── VH_sorted/                      # sequences sorted by VH family
│   │   ├── diff/                       # clones with sequences shared between samples
│   │   └── no_diff/                    # clones with no sequences shared between samples
│   ├── output_clones.tsv               # annotated clone table (step 15)
│   └── lineage_plots/                  # one phylogeny plot per clone (step 15)
│       ├── clone_001.pdf
│       └── clone_002.pdf
├── Patient_2/                          # same structure for every patient
│   ├── VH_sorted/
│   │   ├── diff/
│   │   └── no_diff/
│   ├── output_clones.tsv
│   └── lineage_plots/
└── Patient_N/                          # ... repeated for all patients

---

Troubleshooting

1) BLAST returns no hits

Confirm that ighv_path points to the correct FASTA file, which is provided in the package, and that the BLAST database index files (.nhr, .nin, .nsq) exist alongside it.

2) Finding your BLAST installation

The helper function find_blast() searches for an existing BLAST+ installation automatically. Use it to verify BLAST is accessible before running the pipeline, or to retrieve the correct blast_dir path to pass to run_pipeline().

library(vgene)

# Automatically searches common installation locations
blast_dir <- find_blast()
blast_dir

find_blast() searches the following locations in order:

1. A path you supply directly via the blast_dir argument
2. The BLAST_DIR environment variable
3. A ncbi-blast-*/bin folder in your current working directory
4. /usr/local/ncbi/blast/bin
5. /opt/ncbi/blast/bin

If BLAST is found, the path to the bin directory is returned and printed. You can then pass it directly to run_pipeline():

If your BLAST installation is in a non-standard location, set the BLAST_DIR environment variable before calling find_blast():

Sys.setenv(BLAST_DIR = "~/myblast/bin")
blast_dir <- find_blast()

---

3) Mutation rates split incorrectly

The splitting regex is derived from the folder column of your samples CSV. If sample names contain special regex characters such as ., +, or (, the regex will not match correctly. Simplify or rename the entries in the folder column before running. See ?split_mutation_rate_by_sample_type for details.

---

4)dnapars not found

phylip_exec must be the full path to the dnapars binary. There is no default for this parameter. Check the path and pass it explicitly:

run_pipeline(
  ...,
  phylip_exec = "/opt/phylip-3.697/exe/dnapars"
)

---

Session Info

sessionInfo()
#> R version 4.5.2 (2025-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: Europe/Berlin
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices datasets  utils     methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39       R6_2.6.1            fastmap_1.2.0      
#>  [4] xfun_0.57           cachem_1.1.0        knitr_1.51         
#>  [7] htmltools_0.5.9     rmarkdown_2.31      lifecycle_1.0.5    
#> [10] cli_3.6.6           sass_0.4.10         renv_1.2.2         
#> [13] jquerylib_0.1.4     compiler_4.5.2      rstudioapi_0.18.0  
#> [16] tools_4.5.2         evaluate_1.0.5      bslib_0.10.0       
#> [19] yaml_2.3.12         otel_0.2.0          BiocManager_1.30.27
#> [22] jsonlite_2.0.0      rlang_1.3.0