Introduction
Installation
- External tools to be installed by the user
Input Files
- 1) Samples file
- 2) FASTQ files
Running the Pipeline
- 1) Default run
- 2) Custom run
Re-running a single Step
Output folder structure
Troubleshooting
Session Info

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 VH gene from within the FR1 to the 5’ part of the JH region (note: because the FR1 primer does not always bind at the very start of FR1, the full VH gene sequence is not always recovered), and is sequenced on Illumina platforms such as the NextSeq.

This makes vgene applicable where V gene usage, clonality, and somatic hypermutation (SHM) status are important, for example in the study of B-cell malignancies. vgene takes the raw paired-end FASTQ and performs: quality filtering, primer detection and JH trimming, VH-group clustering, BLAST-based germline assignment, somatic mutation frequency, clone clustering, and per-clone lineage tree reconstruction.

For each clone it reports an 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
Protein translation status	`stop_codon`	Whether translation is in or out of frame
Segment annotation	—	V, D, and CDR3 per-segment annotations

Additionally, phylogeny plots per clone are saved in the lineage_plots/ folder inside each patient directory.

Note on scope and customisation: Custom VH FR1 primers can be substituted via the primers argument in run_pipeline(). However, currently the JH trimming motif (TGGGGC) is fixed and cannot be changed. Additionally, vgene is designed for IGH heavy chain analysis only and is not currently compatible with kappa (IGK) or lambda (IGL) light chain rearrangements, as the germline reference, VH-group clustering, and JH trimming steps are all IGH-specific.

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	What it does	Installation	How to verify
NCBI BLAST+	Aligns reads against the IGHV germline reference to assign V-gene.	`sudo apt-get install ncbi-blast+` (Linux) or `brew install blast` (macOS) or download (Windows)	`blastn -version`
GNU parallel	Runs VH clustering job in parallel.	`sudo apt-get install parallel` (Linux) or `brew install parallel` (macOS)	`parallel --version`
PHYLIP `dnapars`	Reconstructs maximum-parsimony phylogenetic trees for each B-cell clone.	`sudo apt-get install phylip` (Linux) or `brew install phylip` (macOS) or download installer (Windows)	`dnapars`

Input Files

Two inputs are required.

1) Samples file

A .csv file named samples.csv placed in the working directory (the folder you set as the workdir argument in run_pipeline()) with one row per sample. Three columns are required:

Example samples.csv

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	Unique sample identifier	`Patient_77_PB-CSW`
`folder`	Character	Donor or group label	`Patient_77`

One donor (column folder) can have multiple samples(column name), and all the samples with the same donor are analysed together for BLAST, mutation analysis, and clone clustering.

Sample naming convention. The B-cell sample classification 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 stored in a single directory called R1_and_R2/ inside the working directory. If such a folder is not found, run_pipeline() creates an R1_and_R2/ folder inside the working directory and copies them there automatically. This folder serves as the input for quality filtering in step 1.

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 only the required parameters supplied by the user. The file samples.csv and the R1_and_R2 folder containing the FASTQ files must be inside the working directory.

Paths below are shown relative to the working directory unless noted as absolute (full) paths.

library(vgene)
run_pipeline(
  workdir      = "<YOUR VALUE>",  # e.g. "/home/aditee/v-gene"
  samples_file = "<YOUR VALUE>",  # e.g. "samples.csv"
  fastq_dir    = "<YOUR VALUE>",  # e.g. "R1_and_R2"
  ighv_path    = system.file("extdata", "makeblastdb",
                             "IGHV_reference.fasta", package = "vgene"),
                                    # do not edit
  blast_dir    = "<YOUR VALUE>",  # e.g. "ncbi-blast-2.16.0/bin"
  phylip_exec  = "<YOUR VALUE>"   # e.g. "/opt/phylip-3.697/exe/dnapars"
)

Arguments:

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

2) Custom run

Override only the parameters you need to change.

library(vgene)

run_pipeline(
  # --- Required parameters (always supply these) ---
  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 = TRUE, 

  # 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	Accepts	Description
`qual_mean_cutoff`	`20`	Positive integer	Minimum mean read quality (Phred)
`qual_nuc_cutoff`	`35`	Positive integer	Per-nucleotide quality threshold (Phred)
`min_length`	`200`	Positive integer	Minimum read length after trimming (bp)
`max_length`	`200`	Positive integer ≥ `min_length`	Maximum read length for quality control (bp)
`primers`	`c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA", "CTCACCTG", "CTGTAAGG", "ACTCACCT")`	Character vector	Primer sequences for detection and JH trimming
`multicores`	`TRUE`	`TRUE` / `FALSE` / integer ≥ 1	`TRUE` = auto-detect cores, `FALSE` = single core, integer = exact number of cores
`vh_diff_subdir`	`"VH_sorted/diff"`	Character	Path to intermediate vgene clones directory
`igblast_output_path`	`"airr_igblast_output.tsv"`	Character	Path to IgBLAST AIRR output file from step 14
`output_vgene_clone_airr_file`	`"output_clones.tsv"`	Character	Path to annotated clones for lineage tree building

Re-running a single Step

The pipeline is built from individual step functions, one R function per step (e.g. quality filtering using the process_fastq_pipeline() function, germline mapping using the get_antibody_family_info_via_blast(), etc, as shown below). The main run using run_pipeline() simply calls these step functions in order. So, if one step fails, you can run that step on its own to diagnose and fix the issue, then either restart the main run, or continue by running the remaining step functions manually.

Before re-running the step functions, set your working directory and read the folders variable:

setwd("v-gene")   # <YOUR VALUE>: replace with your workdir
samples <- read.csv("samples.csv", stringsAsFactors = FALSE)
folders <- unique(read.csv("samples.csv")$folder)

# Step 1: quality filter 
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 
find_primers(
  input_dir  = "R1_and_R2_corrected/",
  output_dir = "With_primer/",
  primers    = c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA",
                 "CTCACCTG", "CTGTAAGG", "ACTCACCT")
)

# Step 3: Use JH to remove sequences outside the primer pair
trim_vh_files_by_jh(
  primers = c("CTGCAAGG", "ACAGAGAC", "CCTGTGCA",
              "CTCACCTG", "CTGTAAGG", "ACTCACCT")
)

# Step 4: clusters sequences by VH group in parallel
if (!dir.exists("logs/cluster_logs")) {
  dir.create("logs/cluster_logs", recursive = TRUE)
}
system(paste0(
  "parallel --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: computes per-sample cluster size distributions
analyze_clustered_fasta(input_dir = "Collapsed")

# Step 6: merges per-VH clusters into per-sample files
combine_sample_sequences(
  input_dir = "Collapsed/new_more_one/",
  samples   = samples
)

# Step 7: assigns VH sequences using IGHV germline reference with BLAST
get_antibody_family_info_via_blast(
  blast_dir = "<YOUR VALUE>",     # e.g. "ncbi-blast-2.16.0/bin"
  path_IGHV = system.file("extdata", "makeblastdb",
                           "IGHV_reference.fasta", package = "vgene"),
  folders   = folders,
  verbose   = TRUE
)

# Step 8: allele collapsing 
process_folders_collapse_alleles(folders)

# Step 10: mutation data cleaning 
clean_mutation_rate_files(folders)

# Step 11: mutation frequency splitting, one file per sample
split_mutation_frequency_by_sample_type(folders)

# Step 12: use sequences between samples to find a clone
reduce_similar_sequences_in_no_diff(folders)

# Step 13: use sequences within samples to find a clone
reduce_similar_sequences_in_diff(folders)

# Step 14: IgBLAST produces AIRR-formatted output for lineage reconstruction
run_igblast_for_folders(folders)

# Step 15: lineage tree reconstruction builds maximum parsimony trees per clone
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 representative tree shown below. Results are written to the directory given in workdir. Inside it, the pipeline creates a directory for each unique Donor in your samples.csv file. In each Donor folder, clone information and annotations are written to output_clones.tsv, and the lineage plots for each clone are saved in the lineage_plots/ subfolder.

The plots/ folder contains per-patient summary plots of VH family usage, CDR3 length distribution, clone size distributions, and somatic hypermutation.

v-gene/
├── R1_and_R2/                          # copies of input FASTQ files
├── 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 Donor (folder column in CSV)
│   ├── VH_sorted/                      # sequences sorted by VH group
│   │   ├── 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 Donor
│   ├── VH_sorted/
│   │   ├── diff/
│   │   └── no_diff/
│   ├── output_clones.tsv
│   └── lineage_plots/
└── Patient_N/                          # repeated for all Donors

Troubleshooting

1) BLAST returns no hits

Confirm that ighv_path points to the correct FASTA file provided with the vgene package, and that the BLAST database index files (.nhr, .nin, .nsq) exist alongside it. The pipeline builds the database automatically on first run. If the index files are missing or corrupted, delete them and re-run step 7:

get_antibody_family_info_via_blast(
  blast_dir = "ncbi-blast-2.16.0/bin",   # replace with your blast_dir
  path_IGHV = system.file("extdata", "makeblastdb",
                           "IGHV_reference.fasta", package = "vgene"),
  folders   = unique(read.csv("samples.csv")$folder),
  verbose   = TRUE
)

2) Finding your BLAST installation

The helper function find_blast() from the vgene package searches for an existing BLAST+ installation automatically:

library(vgene)
blast_dir <- find_blast()
blast_dir

find_blast() searches these locations in order:

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

To check the path manually, run in a terminal:

which blastn             # Linux/macOS — shows blastn location
dirname $(which blastn)  # shows the bin directory to use as blast_dir

You can then use it as the blast_dir argument in run_pipeline():

blast_dir <- find_blast()
run_pipeline(
  ...,
  blast_dir = blast_dir
)

If your BLAST is in a non-standard location, set BLAST_DIR first:

Sys.setenv(BLAST_DIR = "/your/actual/blast/bin/path")  # << YOUR VALUE
blast_dir <- find_blast()

3) Mutation frequencies split incorrectly

The splitting regex is derived from the folder column of your samples CSV. If folder 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. Folder naming rules are also described in the Input Files section above. See ?vgene::split_mutation_frequency_by_sample_type for details.

4) `dnapars` not found

phylip_exec must be the full path to the dnapars binary from the PHYLIP package. To find the path on your system, run in a terminal:

which dnapars                          # Linux/macOS
find / -name "dnapars" 2>/dev/null    # searches the whole system

Then pass it explicitly:

run_pipeline(
  ...,
  phylip_exec = "<YOUR VALUE>"  # example "/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

Getting Started with vgene

Aditee Kadam

2026-09-11

Introduction

Installation

External tools to be installed by the user

Input Files

1) Samples file

2) FASTQ files

Running the Pipeline

1) Default run

Arguments:

2) Custom run

Arguments:

Re-running a single Step

Output folder structure

Troubleshooting

1) BLAST returns no hits

2) Finding your BLAST installation

3) Mutation frequencies split incorrectly

4) `dnapars` not found

Session Info

Getting Started with vgene

Aditee Kadam

2026-09-11

Introduction

Installation

External tools to be installed by the user

Input Files

1) Samples file

2) FASTQ files

Running the Pipeline

1) Default run

Arguments:

2) Custom run

Arguments:

Re-running a single Step

Output folder structure

Troubleshooting

1) BLAST returns no hits

2) Finding your BLAST installation

3) Mutation frequencies split incorrectly

4) dnapars not found

Session Info

4) `dnapars` not found