Change-O - Repertoire clonal assignment toolkit¶

changeo.Alignment.decodeBTOP(btop)¶

Parse a BTOP string into a list of tuples in CIGAR annotation.

Parameters: btop – BTOP string.
Returns: tuples of (operation, length) for each operation in the BTOP string using CIGAR annotation.
Return type: list

changeo.Alignment.decodeCIGAR(cigar)¶

Parse a CIGAR string into a list of tuples.

Parameters: cigar – CIGAR string.
Returns: tuples of (operation, length) for each operation in the CIGAR string.
Return type: list

changeo.Alignment.encodeCIGAR(alignment)¶

Encodes a list of tuple with alignment information into a CIGAR string.

Parameters: tuple – tuples of (type, length) for each alignment operation.
Returns: CIGAR string.
Return type: str

changeo.Alignment.gapV(seq, v_germ_start, v_germ_length, v_call, references, asis_calls=False)¶

Construction IMGT-gapped V segment sequences.

Parameters

seq (str) – V(D)J sequence alignment (SEQUENCE_VDJ).
v_germ_start (int) – start position V segment alignment in the germline (V_GERM_START_VDJ, 1-based).
v_germ_length (int) – length of the V segment alignment against the germline (V_GERM_LENGTH_VDJ, 1-based).
v_call (str) – V segment allele assignment (V_CALL).
references (dict) – dictionary of IMGT-gapped reference sequences.
asis_calls (bool) – if True do not parse v_call for allele names and just split by comma.

Returns

dictionary containing IMGT-gapped query sequences and germline positions.

Return type

Raises

KeyError – raised if the v_call is not found in the reference dictionary.

changeo.Alignment.getRegions(seq, junction_length)¶

Identify FWR and CDR regions by IMGT definition.

Parameters

seq – IMGT-gapped sequence.
junction_length – length of the junction region in nucleotides.

Returns

dictionary of FWR and CDR sequences.

Return type

changeo.Alignment.inferJunction(seq, j_germ_start, j_germ_length, j_call, references, asis_calls=False, regions='default')¶

Identify junction region by IMGT definition.

Parameters

seq (str) – IMGT-gapped V(D)J sequence alignment (SEQUENCE_IMGT).
j_germ_start (int) – start position J segment alignment in the germline (J_GERM_START, 1-based).
j_germ_length (int) – length of the J segment alignment against the germline (J_GERM_LENGTH).
j_call (str) – J segment allele assignment (J_CALL).
references (dict) – dictionary of IMGT-gapped reference sequences.
asis_calls (bool) – if True do not parse V_CALL for allele names and just split by comma.
regions (str) – name of the IMGT FWR/CDR region definitions to use.

Returns

dictionary containing junction sequence, translation and length.

Return type

changeo.Alignment.padAlignment(alignment, q_start, r_start)¶

Pads the start of an alignment based on query and reference positions.

Parameters

alignment – tuples of (operation, length) for each alignment operation.
q_start – query (input) start position (0-based)
r_start – reference (subject) start position (0-based)

Returns

updated list of tuples of (operation, length) for the alignment.

Return type

changeo.Applications¶

Application wrappers

changeo.Applications.getIgBLASTVersion(exec='igblastn')¶

Gets the version of the IgBLAST executable

Parameters: exec (str) – the name or path to the igblastn executable.
Returns: version number.
Return type: str

changeo.Applications.runASN(fasta, template=None, exec='tbl2asn')¶

Executes tbl2asn to generate Sequin files

Parameters

fasta (str) – fsa file name.
template (str) – sbt file name.
exec (str) – the name or path to the tbl2asn executable.

Returns

tbl2asn console output.

Return type

changeo.Applications.runIgBLASTN(fasta, igdata, loci='ig', organism='human', vdb=None, ddb=None, jdb=None, cdb=None, output=None, format='legacy', threads=1, exec='igblastn')¶

Runs igblastn on a sequence file

Parameters

fasta (str) – fasta file containing sequences.
igdata (str) – path to the IgBLAST database directory (IGDATA environment).
loci (str) – receptor type; one of ‘ig’ or ‘tr’.
organism (str) – species name.
vdb (str) – name of a custom V reference in the database folder to use.
ddb (str) – name of a custom D reference in the database folder to use.
jdb (str) – name of a custom J reference in the database folder to use.
cdb (str) – name of a custom C reference in the database folder to use.
output (str) – output file name. If None, automatically generate from the fasta file name.
format (str) – output format. One of ‘blast’ or ‘airr’.
threads (int) – number of threads for igblastn.
exec (str) – the name or path to the igblastn executable.

Returns

IgBLAST console output.

Return type

changeo.Applications.runIgBLASTP(fasta, igdata, loci='ig', organism='human', vdb=None, output=None, threads=1, exec='igblastp')¶

Runs igblastp on a sequence file

Parameters

fasta (str) – fasta file containing sequences.
igdata (str) – path to the IgBLAST database directory (IGDATA environment).
loci (str) – receptor type; one of ‘ig’ or ‘tr’.
organism (str) – species name.
vdb (str) – name of a custom V reference in the database folder to use.
output (str) – output file name. If None, automatically generate from the fasta file name.
threads (int) – number of threads for igblastp.
exec (str) – the name or path to the igblastp executable.

Returns

IgBLAST console output.

Return type

changeo.Applications.runIgPhyML(rep_file, rep_dir, model='HLP17', motifs='FCH', threads=1, exec='igphyml')¶

Run IgPhyML

Parameters

rep_file (str) – repertoire tsv file.
rep_dir (str) – directory containing input fasta files.
model (str) – model to use.
motif (str) – motifs argument.
threads – number of threads.
exec – the path to the IgPhyMl executable.

Returns

name of the output tree file.

Return type

changeo.Commandline¶

Commandline interface

class changeo.Commandline.CommonHelpFormatter(prog, indent_increment=2, max_help_position=24, width=None)¶

Bases: argparse.RawDescriptionHelpFormatter, argparse.ArgumentDefaultsHelpFormatter

Custom argparse.HelpFormatter

changeo.Commandline.checkArgs(parser)¶

Checks that arguments have been provided and prints help if they have not.

Parameters: parser – An argparse.ArgumentParser defining the commandline arguments.
Returns: True if arguments are present. Prints help and exits if not.
Return type: boolean

changeo.Commandline.getCommonArgParser(db_in=True, db_out=True, out_file=True, failed=True, log=True, format=True, multiproc=False, add_help=True)¶

Defines an ArgumentParser object with common pRESTO arguments

Parameters

db_in (bool) – if True include tab delimited database input arguments.
db_out (bool) – if True include explicit output file name argument.
out_file (bool) – if True add explicit output file name arguments.
failed (bool) – if True include arguments for output of failed results.
log (bool) – if True include log arguments.
format (bool) – input and output type arguments.
multiproc (bool) – if True include multiprocessing arguments.

Returns

an argument parser.

Return type

argparse.ArgumentParser

changeo.Commandline.parseCommonArgs(args, in_arg=None, in_types=None, in_list=False)¶

Checks common arguments from getCommonArgParser and transforms output options to a dictionary

Parameters

args – Argument Namespace defined by ArgumentParser.parse_args.
in_arg – String defining a non-standard input file argument to verify; by default ‘db_files’ and ‘seq_files’ are supported in that order.
in_types – List of types (file extensions as strings) to allow for files in file_arg; if None do not check type.
in_list – if True allow multiple input files with the out_name and log arguments.

Returns

Dictionary copy of args with output arguments embedded in the dictionary out_args

Return type

changeo.Commandline.setDefaultFields(args, defaults, format='airr')¶

Sets default field arguments by format

Parameters

args (dict) – parsed argument dictionary.
defaults (dict) – default variables to set with with keys as argument variables and values as AIRR field names.
format (str) – one of ‘changeo’ or ‘airr’ which defines the file format.

Returns

modified input args.

Return type

changeo.Distance¶

Distance calculations

changeo.Distance.calcDistances(sequences, n, dist_mat, sym='avg', norm=None)¶

Calculate pairwise distances between input sequences

Parameters

sequences – List of sequences for which to calculate pairwise distances
n – Length of n-mers to be used in calculating distance
dist_mat – pandas.DataFrame of mutation distances
norm – Normalization method. One of None, ‘len’, or ‘mut’.
sym – Symmetry method; one of ‘avg’ of ‘min.

Returns

numpy matrix of pairwise distances between input sequences

Return type

ndarray

changeo.Distance.formClusters(dists, link, distance)¶

Form clusters based on hierarchical clustering of input distance matrix with linkage type and cutoff distance

Parameters

dists – numpy matrix of distances
link – Linkage type for hierarchical clustering
distance – Distance at which to cut into clusters

Returns

List of cluster assignments

Return type

changeo.Distance.getAADistMatrix(mat=None, mask_dist=0, gap_dist=0)¶

Generates an amino acid distance matrix

Parameters

mat – Input distance matrix to extend to full alphabet; if unspecified, creates Hamming distance matrix that incorporates IUPAC equivalencies
mask_dict – Score for all matches against an X character
gap_dist – Score for all matches against a gap (-, .) character

Returns

pandas.DataFrame of distances

Return type

DataFrame

changeo.Distance.getDNADistMatrix(mat=None, mask_dist=0, gap_dist=0)¶

Generates a DNA distance matrix

Parameters

mat – Input distance matrix to extend to full alphabet; if unspecified, creates Hamming distance matrix that incorporates IUPAC equivalencies
mask_dist – Distance for all matches against an N character
gap_dist – Distance for all matches against a gap (-, .) character

Returns

pandas.DataFrame of distances

Return type

DataFrame

changeo.Distance.getNmers(sequences, n)¶

Breaks input sequences down into n-mers

Parameters

sequences – List of sequences to be broken into n-mers
n – Length of n-mers to return

Returns

Dictionary mapping sequence to a list of n-mers

Return type

changeo.Distance.zip_equal(*iterables)¶

Zips iterables and raises exception if different lengths

Parameters: iterables – pointer to iterables to zip together
Returns: A generator of tuples with combined elements from the iterables
Return type: iter

changeo.Gene¶

Gene annotations

changeo.Gene.buildClonalGermline(receptors, references, seq_field='sequence_imgt', v_field='v_call', d_field='d_call', j_field='j_call', amino_acid=False)¶

Determine consensus clone sequence and create germline for clone

Parameters

receptors (changeo.Receptor.Receptor) – list of Receptor objects
references (dict) – dictionary of IMGT gapped germline sequences
seq_field (str) – Receptor attribute in which to look for sequence
v_field (str) – Receptor attributein which to look for V call
d_field (str) – Receptor attributein which to look for D call
j_field (str) – Receptor attributein which to look for J call
amino_acid (bool) – if True then use the amino acid positional fields, otherwise use the nucleotide fields.

Returns

log dictionary, dictionary of {germline_type: germline_sequence},: dictionary of consensus {segment: gene call}

Return type

changeo.Gene.buildGermline(receptor, references, seq_field='sequence_imgt', v_field='v_call', d_field='d_call', j_field='j_call', amino_acid=False)¶

Join gapped germline sequences aligned with sample sequences

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object.
references (dict) – dictionary of IMGT gapped germline sequences.
seq_field (str) – Receptor attribute in which to look for sequence.
v_field (str) – Receptor attribute in which to look for V call.
d_field (str) – Receptor attribute in which to look for V call.
j_field (str) – Receptor attribute in which to look for V call.
amino_acid (bool) – if True then use the amino acid positional fields, otherwise use the nucleotide fields.

Returns

log dictionary, dictionary of {germline_type: germline_sequence}, dictionary of {segment: gene call}

Return type

changeo.Gene.getAllele(gene, action='first')¶

Extract allele from gene call string

Parameters

gene (str) – string with gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first allele calls when action is ‘first’. tuple: Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getAlleleNumber(gene, action='first')¶

Extract allele number from gene call string

Parameters

gene (str) – string with gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first allele number call when action is ‘first’. tuple: Tuple of allele numbers for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getCAllele(gene, action='first')¶

Extract C-region allele gene call string

Parameters

gene (str) – string with C-region gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first C-region allele call when action is ‘first’. tuple: Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getCGene(gene, action='first')¶

Extract C-region gene from gene call string

Parameters

gene (str) – string with C-region gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first C-region gene call when action is ‘first’. tuple: Tuple of gene calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getDAllele(gene, action='first')¶

Extract D allele gene from gene call string

Parameters

gene (str) – string with D gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first D allele call when action is ‘first’. tuple: Tuple of D allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getDGermline(receptor, references, d_field='d_call', amino_acid=False)¶

Extract D allele and germline sequence

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object
references (dict) – dictionary of germline sequences
d_field (str) – Receptor attribute containing the D allele assignment
amino_acid (bool) – if True then use the amino acid positional fields, otherwise use the nucleotide fields.

Returns

D allele name, D segment germline sequence

Return type

changeo.Gene.getFamily(gene, action='first')¶

Extract family from gene call string

Parameters

gene (str) – string with gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first family call when action is ‘first’. tuple: Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getGene(gene, action='first')¶

Extract gene from gene call string

Parameters

gene (str) – string with gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first gene call when action is ‘first’. tuple: Tuple of gene calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getJAllele(gene, action='first')¶

Extract J allele gene from gene call string

Parameters

gene (str) – string with J gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first J allele call when action is ‘first’. tuple: Tuple of J allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getJGermline(receptor, references, j_field='j_call', amino_acid=False)¶

Extract J allele and germline sequence

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object
references (dict) – dictionary of germline sequences
j_field (str) – Receptor attribute containing the J allele assignment
amino_acid (bool) – if True then use the amino acid positional fields, otherwise use the nucleotide fields.

Returns

J allele name, J segment germline sequence

Return type

changeo.Gene.getLocus(gene, action='first')¶

Extract locus from gene call string

Parameters

gene (str) – string with gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first locus call when action is ‘first’. tuple: Tuple of locus calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getVAllele(gene, action='first')¶

Extract V allele gene from gene call string

Parameters

gene (str) – string with V gene calls
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the first V allele call when action is ‘first’. tuple: Tuple of V allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.getVGermline(receptor, references, v_field='v_call', amino_acid=False)¶

Extract V allele and germline sequence

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object
references (dict) – dictionary of germline sequences
v_field (str) – Receptor attribute containing the V allele assignment
amino_acid (bool) – if True then use the amino acid positional fields, otherwise use the nucleotide fields.

Returns

V allele name, V segment germline sequence

Return type

changeo.Gene.parseGeneCall(gene, regex, action='first')¶

Extract alleles from strings

Parameters

gene (str) – string with gene calls
regex (re.Pattern) – compiled regular expression for allele match
action (str) – action to perform for multiple alleles; one of (‘first’, ‘set’, ‘list’).

Returns

String of the allele when action is ‘first’; tuple: Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

changeo.Gene.stitchRegions(receptor, v_seq, d_seq, j_seq, amino_acid=False)¶

Assemble full length region encoding

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object
v_seq (str) – V segment germline sequence as a string
d_seq (str) – D segment germline sequence as a string
j_seq (str) – J segment germline sequence as a string
amino_acid (bool) – if True use amino acid positional fields, otherwise use nucleotide fields.

Returns

string defining germline regions

Return type

changeo.Gene.stitchVDJ(receptor, v_seq, d_seq, j_seq, amino_acid=False)¶

Assemble full length germline sequence

Parameters

receptor (changeo.Receptor.Receptor) – Receptor object
v_seq (str) – V segment sequence as a string
d_seq (str) – D segment sequence as a string
j_seq (str) – J segment sequence as a string
amino_acid (bool) – if True use X for N/P regions and amino acid positional fields, otherwise use N and nucleotide fields.

Returns

full germline sequence

Return type

changeo.IO¶

File I/O and parsers

class changeo.IO.AIRRReader(handle)¶

Bases: changeo.IO.TSVReader

An iterator to read and parse AIRR formatted data.

class changeo.IO.AIRRWriter(handle, fields=['sequence_id', 'sequence', 'sequence_alignment', 'germline_alignment', 'rev_comp', 'productive', 'stop_codon', 'vj_in_frame', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction', 'junction_length', 'junction_aa', 'np1_length', 'np2_length', 'v_sequence_start', 'v_sequence_end', 'v_germline_start', 'v_germline_end', 'd_sequence_start', 'd_sequence_end', 'd_germline_start', 'd_germline_end', 'j_sequence_start', 'j_sequence_end', 'j_germline_start', 'j_germline_end'])¶

Bases: changeo.IO.TSVWriter

Writes AIRR formatted data.

writeReceptor(records)¶

Writes a row from a Receptor object

Parameters: records – a changeo.Receptor object to write or iterable of such objects.
Returns: None

class changeo.IO.ChangeoReader(handle)¶

Bases: changeo.IO.TSVReader

An iterator to read and parse Change-O formatted data.

class changeo.IO.ChangeoWriter(handle, fields=['SEQUENCE_ID', 'SEQUENCE_INPUT', 'FUNCTIONAL', 'IN_FRAME', 'STOP', 'MUTATED_INVARIANT', 'INDELS', 'LOCUS', 'V_CALL', 'D_CALL', 'J_CALL', 'SEQUENCE_VDJ', 'SEQUENCE_IMGT', 'V_SEQ_START', 'V_SEQ_LENGTH', 'V_GERM_START_VDJ', 'V_GERM_LENGTH_VDJ', 'V_GERM_START_IMGT', 'V_GERM_LENGTH_IMGT', 'NP1_LENGTH', 'D_SEQ_START', 'D_SEQ_LENGTH', 'D_GERM_START', 'D_GERM_LENGTH', 'NP2_LENGTH', 'J_SEQ_START', 'J_SEQ_LENGTH', 'J_GERM_START', 'J_GERM_LENGTH', 'JUNCTION', 'JUNCTION_LENGTH', 'GERMLINE_IMGT'], header=True)¶

Bases: changeo.IO.TSVWriter

Writes Change-O formatted data.

writeReceptor(records)¶

Writes a row from a Receptor object

Parameters: records – a changeo.Receptor.Receptor object to write or an iterable of such objects.
Returns: None

class changeo.IO.IHMMuneReader(ihmmune, sequences, references, receptor=True)¶

Bases: object

An iterator to read and parse iHMMune-Align output files.

__iter__()¶

Iterator initializer.

Returns: changeo.IO.IHMMuneReader

__next__()¶

Next method.

Returns: parsed IMGT/HighV-QUEST result as an Receptor (receptor=True) or dictionary (receptor=False).
Return type: changeo.Receptor.Receptor

static customFields(scores=False, regions=False, cell=False, schema=None)¶

Returns non-standard Receptor attributes defined by the parser

Parameters

scores – if True include alignment scoring fields.
regions – if True include IMGT-gapped CDR and FWR region fields.
schema – schema class to pass field through for conversion. If None, return changeo.Receptor.Receptor attribute names.

Returns

list of field names.

Return type

ihmmune_fields = ['SEQUENCE_ID', 'V_CALL', 'D_CALL', 'J_CALL', 'V_SEQ', 'NP1_SEQ', 'D_SEQ', 'NP2_SEQ', 'J_SEQ', 'V_MUT', 'D_MUT', 'J_MUT', 'NX_COUNT', 'J_INFRAME', 'V_SEQ_START', 'STOP_COUNT', 'D_PROB', 'HMM_SCORE', 'RC', 'COMMON_MUT', 'COMMON_NX_COUNT', 'V_SEQ_START', 'V_SEQ_LENGTH', 'A_SCORE']¶

parseRecord(record)¶

Parses a single row from each IMTG file.

Parameters: record – dictionary containing one row of iHMMune-Align file.
Returns: database entry for the row.
Return type: dict

class changeo.IO.IMGTReader(summary, gapped, ntseq, junction, receptor=True)¶

Bases: object

An iterator to read and parse IMGT output files.

__iter__()¶

Iterator initializer.

Returns: changeo.IO.IMGTReader

__next__()¶

Next method.

Returns: parsed IMGT/HighV-QUEST result as an Receptor (receptor=True) or dictionary (receptor=False).
Return type: changeo.Receptor.Receptor

static customFields(scores=False, regions=False, junction=False, schema=None)¶

Returns non-standard fields defined by the parser

Parameters

scores – if True include alignment scoring fields.
regions – if True include IMGT-gapped CDR and FWR region fields.
junction – if True include detailed junction annotation fields.
schema – schema class to pass field through for conversion. If None, return changeo.Receptor.Receptor attribute names.

Returns

list of field names.

Return type

parseRecord(summary, gapped, ntseq, junction)¶

Parses a single row from each IMTG file.

Parameters

summary – dictionary containing one row of the ‘1_Summary’ file.
gapped – dictionary containing one row of the ‘2_IMGT-gapped-nt-sequences’ file.
ntseq – dictionary containing one row of the ‘3_Nt-sequences’ file.
junction – dictionary containing one row of the ‘6_Junction’ file.

Returns

database entry for the row.

Return type

class changeo.IO.IgBLASTReader(igblast, sequences, references, asis_calls=False, regions='default', receptor=True, infer_junction=False)¶

Bases: object

An iterator to read and parse IgBLAST output files

__iter__()¶

Iterator initializer.

Returns: changeo.IO.IgBLASTReader

__next__()¶

Next method.

Returns: parsed IMGT/HighV-QUEST result as an Receptor (receptor=True) or dictionary (receptor=False).
Return type: changeo.Receptor.Receptor

static customFields(schema=None)¶

Returns non-standard fields defined by the parser

Parameters: schema – schema class to pass field through for conversion. If None, return changeo.Receptor.Receptor attribute names.
Returns: list of field names.
Return type: list

parseBlock(block)¶

Parses an IgBLAST result into separate sections

Parameters

block (iter) – an iterator from itertools.groupby containing a single IgBLAST result.

Returns

a parsed results block;: with the keys ‘query’ (sequence identifier as a string), ‘summary’ (dictionary of the alignment summary), ‘subregion’ (dictionary of IgBLAST CDR3 sequences), and ‘hits’ (VDJ hit table as a list of dictionaries). Returns None if the block has no data that can be parsed.

Return type

parseSections(sections)¶

Parses an IgBLAST sections into a db dictionary

Parameters: sections – dictionary of parsed sections from parseBlock.
Returns: db entries.
Return type: dict

class changeo.IO.IgBLASTReaderAA(igblast, sequences, references, asis_calls=False, regions='default', receptor=True, infer_junction=False)¶

Bases: changeo.IO.IgBLASTReader

An iterator to read and parse IgBLAST amino acid alignment output files

static customFields(schema=None)¶

Returns non-standard fields defined by the parser

Parameters: schema – schema class to pass field through for conversion. If None, return changeo.Receptor.Receptor attribute names.
Returns: list of field names.
Return type: list

parseSections(sections)¶

Parses an IgBLAST sections into a db dictionary

Parameters: sections – dictionary of parsed sections from parseBlock.
Returns: db entries.
Return type: dict

class changeo.IO.TSVReader(handle)¶

Bases: object

Simple csv.DictReader wrapper to read format agnostic TSV files.

reader¶

reader object.

Type: iter

fields¶

field names.

Type: list

__iter__()¶

Iterator initializer

Returns: changeo.IO.TSVReader

__next__()¶

Next method

Returns: row as a dictionary of field:value pairs.
Return type: dist

class changeo.IO.TSVWriter(handle, fields, header=True)¶

Bases: object

Simple csv.DictWriter wrapper to write format agnostic TSV files.

writeDict(records)¶

Writes a row from a dictionary

Parameters: records – dictionary of row data or an iterable of such objects.
Returns: None

writeHeader()¶

Writes the header

Returns: None

changeo.IO.checkFields(attributes, header, schema=<class 'changeo.Receptor.AIRRSchema'>)¶

Checks that a file header contains a required set of Receptor attributes

Parameters

attributes (list) – list of Receptor attributes to check for.
header (list) – list of fields names in the file header.
schema (object) – schema object to convert field names to Receptor attributes.

Returns

True if all attributes mapping fields are found.

Return type

bool

Raises

LookupError –

changeo.IO.countDbFile(file)¶

Counts the records in database files

Parameters: file – tab-delimited database file.
Returns: count of records in the database file.
Return type: int

changeo.IO.extractIMGT(imgt_output)¶

Extract necessary files from IMGT/HighV-QUEST results.

Parameters: imgt_output – zipped file or unzipped folder output by IMGT/HighV-QUEST.
Returns: (temporary directory handle, dictionary with names of extracted IMGT files).
Return type: tuple

changeo.IO.getDbFields(file, add=None, exclude=None, reader=<class 'changeo.IO.TSVReader'>)¶

Get field names from a db file

Parameters

file – db file to pull base fields from.
add – fields to append to the field set.
exclude – fields to exclude from the field set.
reader – reader class.

Returns

list of field names

Return type

changeo.IO.getFormatOperators(format)¶

Simple wrapper for fetching the set of operator classes for a data format

Parameters: format (str) – name of the data format.
Returns: a tuple with the reader class, writer class, and schema definition class.
Return type: tuple

changeo.IO.getOutputHandle(file, out_label=None, out_dir=None, out_name=None, out_type=None)¶

Opens an output file handle

Parameters

file – filename to base output file name on.
out_label – text to be inserted before the file extension; if None do not add a label.
out_type – the file extension of the output file; if None use input file extension.
out_dir – the output directory; if None use directory of input file
out_name – the short filename to use for the output file; if None use input file short name.

Returns

File handle

Return type

file

changeo.IO.getOutputName(file, out_label=None, out_dir=None, out_name=None, out_type=None)¶

Creates and output filename from an existing filename

Parameters

file – filename to base output file name on.
out_label – text to be inserted before the file extension; if None do not add a label.
out_type – the file extension of the output file; if None use input file extension.
out_dir – the output directory; if None use directory of input file
out_name – the short filename to use for the output file; if None use input file short name.

Returns

file name.

Return type

changeo.IO.readGermlines(references, asis=False, warn=False)¶

Parses germline repositories

Parameters

references (list) – list of strings specifying directories and/or files from which to read germline records.
asis (bool) – if True use sequence ID as record name and do not parse headers for allele names.
warn (bool) – print warning messages to standard error if True.

Returns

Dictionary of germlines in the form {allele: sequence}.

Return type

changeo.IO.splitName(file)¶

Extract the extension from a file name

Parameters: file (str) – file name.
Returns: tuple of the file directory, basename and extension.
Return type: tuple

changeo.IO.yamlDict(file)¶

Returns a dictionary from a yaml file

Parameters: file (str) – simple yaml file with rows in the form ‘argument: value’.
Returns: dictionary of key:value pairs in the file.
Return type: dict

changeo.Multiprocessing¶

Multiprocessing

class changeo.Multiprocessing.DbData(key, records)¶

Bases: object

A class defining data objects for worker processes

id¶: result identifier

data¶: list of data records

valid¶: True if preprocessing was successfull and data should be processed

class changeo.Multiprocessing.DbResult(key, records)¶

Bases: object

A class defining result objects for collector processes

id¶: result identifier

data¶: list of original data records

results¶: list of processed records

data_pass¶: list of records that pass filtering for workers that split data before processing

data_fail¶: list of records that failed filtering for workers that split data before processing

valid¶: True if processing was successful and results should be written

log¶: OrderedDict of log items

property data_count¶

changeo.Multiprocessing.collectDbQueue(alive, result_queue, collect_queue, db_file, label, fields, writer=<class 'changeo.IO.AIRRWriter'>, out_file=None, out_args={'failed': False, 'log_file': None, 'out_dir': None, 'out_name': None, 'out_type': 'tsv'})¶

Pulls from results queue, assembles results and manages log and file IO

Parameters

alive – multiprocessing.Value boolean controlling whether processing continues; when False function returns.
result_queue – multiprocessing.Queue holding worker results.
collect_queue – multiprocessing.Queue to store collector return values.
db_file – database file name.
label – task label used to tag the output files.
fields – list of output fields.
writer – writer class.
out_file – output file name. Automatically generated from the input file if None.
out_args – common output argument dictionary from parseCommonArgs.

Returns

Adds a dictionary with key value pairs to collect_queue containing: ’log’ defining a log object along with the ‘pass’ and ‘fail’ output file names.

Return type

None

changeo.Multiprocessing.feedDbQueue(alive, data_queue, db_file, reader=<class 'changeo.IO.AIRRReader'>, group_func=None, group_args={})¶

Feeds the data queue with Ig records

Parameters

alive – multiprocessing.Value boolean controlling whether processing continues if False exit process
data_queue – multiprocessing.Queue to hold data for processing
db_file – database file
reader – database reader class
group_func – function to use for grouping records
group_args – dictionary of arguments to pass to group_func

Returns

None

changeo.Multiprocessing.processDbQueue(alive, data_queue, result_queue, process_func, process_args={}, filter_func=None, filter_args={})¶

Pulls from data queue, performs calculations, and feeds results queue

Parameters

alive – multiprocessing.Value boolean controlling whether processing continues; when False function returns
data_queue – multiprocessing.Queue holding data to process
result_queue – multiprocessing.Queue to hold processed results
process_func – function to use for processing sequences
process_args – dictionary of arguments to pass to process_func
filter_func – function to use for filtering sequences before processing
filter_args – dictionary of arguments to pass to filter_func

Returns

None

changeo.Receptor¶

Receptor data structure

class changeo.Receptor.AIRRSchema¶

Bases: object

AIRR format to Receptor mappings

fields = ['sequence_id', 'sequence', 'sequence_alignment', 'germline_alignment', 'sequence_aa', 'sequence_aa_alignment', 'germline_aa_alignment', 'rev_comp', 'productive', 'stop_codon', 'vj_in_frame', 'v_frameshift', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction', 'junction_start', 'junction_end', 'junction_length', 'junction_aa', 'junction_aa_length', 'np1_length', 'np2_length', 'np1_aa_length', 'np2_aa_length', 'v_sequence_start', 'v_sequence_end', 'v_sequence_length', 'v_germline_start', 'v_germline_end', 'v_germline_length', 'v_sequence_aa_start', 'v_sequence_aa_end', 'v_sequence_aa_length', 'v_germline_aa_start', 'v_germline_aa_end', 'v_germline_aa_length', 'd_sequence_start', 'd_sequence_end', 'd_sequence_length', 'd_germline_start', 'd_germline_end', 'd_germline_length', 'd_sequence_aa_start', 'd_sequence_aa_end', 'd_sequence_aa_length', 'd_germline_aa_start', 'd_germline_aa_end', 'd_germline_aa_length', 'j_sequence_start', 'j_sequence_end', 'j_sequence_length', 'j_germline_start', 'j_germline_end', 'j_germline_length', 'j_sequence_aa_start', 'j_sequence_aa_end', 'j_sequence_aa_length', 'j_germline_aa_start', 'j_germline_aa_end', 'j_germline_aa_length', 'germline_alignment_d_mask', 'v_score', 'v_identity', 'v_support', 'v_cigar', 'd_score', 'd_identity', 'd_support', 'd_cigar', 'j_score', 'j_identity', 'j_support', 'j_cigar', 'vdj_score', 'cdr1', 'cdr2', 'cdr3', 'fwr1', 'fwr2', 'fwr3', 'fwr4', 'cdr1_aa', 'cdr2_aa', 'cdr3_aa', 'fwr1_aa', 'fwr2_aa', 'fwr3_aa', 'fwr4_aa', 'cdr1_start', 'cdr1_end', 'cdr2_start', 'cdr2_end', 'cdr3_start', 'cdr3_end', 'fwr1_start', 'fwr1_end', 'fwr2_start', 'fwr2_end', 'fwr3_start', 'fwr3_end', 'fwr4_start', 'fwr4_end', 'n1_length', 'n2_length', 'p3v_length', 'p5d_length', 'p3d_length', 'p5j_length', 'd_frame', 'cdr3_igblast', 'cdr3_igblast_aa', 'duplicate_count', 'consensus_count', 'umi_count', 'clone_id', 'cell_id']¶

static fromReceptor(field)¶

Returns an AIRR column name from a Receptor attribute name

Parameters: field – Receptor attribute name.
Returns: AIRR column name.
Return type: str

out_type = 'tsv'¶

required = ['sequence_id', 'sequence', 'sequence_alignment', 'germline_alignment', 'rev_comp', 'productive', 'stop_codon', 'vj_in_frame', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction', 'junction_length', 'junction_aa', 'np1_length', 'np2_length', 'v_sequence_start', 'v_sequence_end', 'v_germline_start', 'v_germline_end', 'd_sequence_start', 'd_sequence_end', 'd_germline_start', 'd_germline_end', 'j_sequence_start', 'j_sequence_end', 'j_germline_start', 'j_germline_end']¶

static toReceptor(field)¶

Returns a Receptor attribute name from an AIRR column name

Parameters: field – AIRR column name.
Returns: Receptor attribute name.
Return type: str

class changeo.Receptor.AIRRSchemaAA¶

Bases: changeo.Receptor.AIRRSchema

AIRR format to Receptor amino acid mappings

required = ['sequence_id', 'sequence', 'sequence_alignment', 'germline_alignment', 'sequence_aa', 'sequence_aa_alignment', 'germline_aa_alignment', 'rev_comp', 'productive', 'stop_codon', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'junction', 'junction_length', 'junction_aa', 'v_sequence_aa_start', 'v_sequence_aa_end', 'v_germline_aa_start', 'v_germline_aa_end']¶

class changeo.Receptor.ChangeoSchema¶

Bases: object

Change-O to Receptor mappings

fields = ['SEQUENCE_ID', 'SEQUENCE_INPUT', 'SEQUENCE_AA_INPUT', 'FUNCTIONAL', 'IN_FRAME', 'STOP', 'MUTATED_INVARIANT', 'INDELS', 'V_FRAMESHIFT', 'LOCUS', 'V_CALL', 'D_CALL', 'J_CALL', 'C_CALL', 'SEQUENCE_VDJ', 'SEQUENCE_IMGT', 'SEQUENCE_AA_VDJ', 'SEQUENCE_AA_IMGT', 'V_SEQ_START', 'V_SEQ_LENGTH', 'V_GERM_START_VDJ', 'V_GERM_LENGTH_VDJ', 'V_GERM_START_IMGT', 'V_GERM_LENGTH_IMGT', 'V_SEQ_AA_START', 'V_SEQ_AA_LENGTH', 'V_GERM_AA_START_VDJ', 'V_GERM_AA_LENGTH_VDJ', 'V_GERM_AA_START_IMGT', 'V_GERM_AA_LENGTH_IMGT', 'NP1_LENGTH', 'NP1_AA_LENGTH', 'D_SEQ_START', 'D_SEQ_LENGTH', 'D_GERM_START', 'D_GERM_LENGTH', 'D_SEQ_AA_START', 'D_SEQ_AA_LENGTH', 'D_GERM_AA_START', 'D_GERM_AA_LENGTH', 'NP2_LENGTH', 'NP2_AA_LENGTH', 'J_SEQ_START', 'J_SEQ_LENGTH', 'J_GERM_START', 'J_GERM_LENGTH', 'J_SEQ_AA_START', 'J_SEQ_AA_LENGTH', 'J_GERM_AA_START', 'J_GERM_AA_LENGTH', 'JUNCTION', 'JUNCTION_LENGTH', 'GERMLINE_IMGT', 'GERMLINE_AA_IMGT', 'JUNCTION_START', 'V_SCORE', 'V_IDENTITY', 'V_EVALUE', 'V_BTOP', 'V_CIGAR', 'D_SCORE', 'D_IDENTITY', 'D_EVALUE', 'D_BTOP', 'D_CIGAR', 'J_SCORE', 'J_IDENTITY', 'J_EVALUE', 'J_BTOP', 'J_CIGAR', 'VDJ_SCORE', 'FWR1_IMGT', 'FWR2_IMGT', 'FWR3_IMGT', 'FWR4_IMGT', 'CDR1_IMGT', 'CDR2_IMGT', 'CDR3_IMGT', 'FWR1_AA_IMGT', 'FWR2_AA_IMGT', 'FWR3_AA_IMGT', 'FWR4_AA_IMGT', 'CDR1_AA_IMGT', 'CDR2_AA_IMGT', 'CDR3_AA_IMGT', 'N1_LENGTH', 'N2_LENGTH', 'P3V_LENGTH', 'P5D_LENGTH', 'P3D_LENGTH', 'P5J_LENGTH', 'D_FRAME', 'CDR3_IGBLAST', 'CDR3_IGBLAST_AA', 'CONSCOUNT', 'DUPCOUNT', 'UMICOUNT', 'CLONE', 'CELL']¶

static fromReceptor(field)¶

Returns a Change-O column name from a Receptor attribute name

Parameters: field – Receptor attribute name.
Returns: Change-O column name.
Return type: str

out_type = 'tab'¶

required = ['SEQUENCE_ID', 'SEQUENCE_INPUT', 'FUNCTIONAL', 'IN_FRAME', 'STOP', 'MUTATED_INVARIANT', 'INDELS', 'LOCUS', 'V_CALL', 'D_CALL', 'J_CALL', 'SEQUENCE_VDJ', 'SEQUENCE_IMGT', 'V_SEQ_START', 'V_SEQ_LENGTH', 'V_GERM_START_VDJ', 'V_GERM_LENGTH_VDJ', 'V_GERM_START_IMGT', 'V_GERM_LENGTH_IMGT', 'NP1_LENGTH', 'D_SEQ_START', 'D_SEQ_LENGTH', 'D_GERM_START', 'D_GERM_LENGTH', 'NP2_LENGTH', 'J_SEQ_START', 'J_SEQ_LENGTH', 'J_GERM_START', 'J_GERM_LENGTH', 'JUNCTION', 'JUNCTION_LENGTH', 'GERMLINE_IMGT']¶

static toReceptor(field)¶

Returns a Receptor attribute name from a Change-O column name

Parameters: field – Change-O column name.
Returns: Receptor attribute name.
Return type: str

class changeo.Receptor.ChangeoSchemaAA¶

Bases: changeo.Receptor.ChangeoSchema

Change-O to Receptor amino acid mappings

required = ['SEQUENCE_ID', 'SEQUENCE_AA_INPUT', 'STOP', 'INDELS', 'LOCUS', 'V_CALL', 'SEQUENCE_AA_VDJ', 'SEQUENCE_AA_IMGT', 'V_SEQ_AA_START', 'V_SEQ_AA_LENGTH', 'V_GERM_AA_START_VDJ', 'V_GERM_AA_LENGTH_VDJ', 'V_GERM_AA_START_IMGT', 'V_GERM_AA_LENGTH_IMGT', 'GERMLINE_AA_IMGT']¶

class changeo.Receptor.Receptor(data)¶

Bases: object

A class defining a V(D)J sequence and its annotations

property d_germ_aa_end¶: Position of the last amino acid in the D germline amino acid alignment

property d_germ_end¶: Position of the last nucleotide in the D germline sequence alignment

property d_seq_aa_end¶: Position of the last D amino acid in the input amino acid sequence

property d_seq_end¶: Position of the last D nucleotide in the input sequence

getAIRR(field, seq=False)¶

Get an attribute from an AIRR field name

Parameters

field – AIRR column name as a string
seq – if True return the attribute as a Seq object

Returns

Value in the AIRR field. Returns None if the field cannot be found.

getAlleleCalls(calls, action='first')¶

Get multiple allele calls

Parameters

calls – iterable of calls to get; one or more of (‘v’,’d’,’j’)
actions – One of (‘first’,’set’)

Returns

List of requested calls in order

Return type

getAlleleNumbers(calls, action='first')¶

Get multiple allele numeric identifiers

Parameters

calls – iterable of calls to get; one or more of (‘v’,’d’,’j’)
actions – One of (‘first’,’set’)

Returns

List of requested calls in order

Return type

getChangeo(field, seq=False)¶

Get an attribute from a Change-O field name

Parameters

field – Change-O column name as a string
seq – if True return the attribute as a Seq object

Returns

Value in the Change-O field. Returns None if the field cannot be found.

getDAllele(action='first', field=None)¶

D segment allele getter

Parameters

actions – One of ‘first’, ‘set’ or ‘list’
field – attribute or annotation name containing the D call. Use d_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getDAlleleNumber(action='first', field=None)¶

D segment allele number getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the D call. Use d_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele numbers for ‘set’ or ‘list’ actions.

Return type

getDFamily(action='first', field=None)¶

D segment family getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the D call. Use d_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getDGene(action='first', field=None)¶

D segment gene getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the D call. Use d_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getFamilyCalls(calls, action='first')¶

Get multiple family calls

Parameters

calls – iterable of calls to get; one or more of (‘v’,’d’,’j’)
actions – One of (‘first’,’set’)

Returns

List of requested calls in order

Return type

getField(field)¶

Get an attribute or annotation value

Parameters: field – attribute name as a string
Returns: Value in the attribute. Returns None if the attribute cannot be found.

getGeneCalls(calls, action='first')¶

Get multiple gene calls

Parameters

calls – iterable of calls to get; one or more of (‘v’,’d’,’j’)
actions – One of (‘first’,’set’)

Returns

List of requested calls in order

Return type

getJAllele(action='first', field=None)¶

J segment allele getter

Parameters

actions – One of ‘first’, ‘set’ or ‘list’
field – attribute or annotation name containing the J call. Use j_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getJAlleleNumber(action='first', field=None)¶

J segment allele number getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the J call. Use j_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele numbers for ‘set’ or ‘list’ actions.

Return type

getJFamily(action='first', field=None)¶

J segment family getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the J call. Use j_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getJGene(action='first', field=None)¶

J segment gene getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the J call. Use j_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getSeq(field)¶

Get an attribute value converted to a Seq object

Parameters: field – variable name as a string
Returns: Value in the field as a Seq object
Return type: Bio.Seq.Seq

getVAllele(action='first', field=None)¶

V segment allele getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the V call. Use v_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getVAlleleNumber(action='first', field=None)¶

V segment allele number getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the V call. Use v_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele numbers for ‘set’ or ‘list’ actions.

Return type

getVFamily(action='first', field=None)¶

V segment family getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the V call. Use v_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

getVGene(action='first', field=None)¶

V segment gene getter

Parameters

actions – One of ‘first’, ‘set’ or list’
field – attribute or annotation name containing the V call. Use v_call attribute if None.

Returns

String of the allele when action is ‘first’; tuple : Tuple of allele calls for ‘set’ or ‘list’ actions.

Return type

property j_germ_aa_end¶: Position of the last amino acid in the J germline amino acid alignment

property j_germ_end¶: Position of the last nucleotide in the J germline sequence alignment

property j_seq_aa_end¶: Position of the last J amino acid in the input amino sequence

property j_seq_end¶: Position of the last J nucleotide in the input sequence

property junction_end¶: Position of the last junction nucleotide in the input sequence

setDict(data, parse=False)¶

Adds or updates multiple attributes and annotations

Parameters

data – a dictionary of annotations to add or update.
parse – if True pass values through string parsing functions for known fields.

Returns

updates attribute values and the annotations attribute.

Return type

None

setField(field, value, parse=False)¶

Set an attribute or annotation value

Parameters

field – attribute name as a string
value – value to assign
parse – if True pass values through string parsing functions for known fields.

Returns

None. Updates attribute or annotation.

toDict()¶

Convert the namespace to a dictionary

Returns: member fields with values converted to appropriate strings
Return type: dict

property v_germ_aa_end_imgt¶: Position of the last nucleotide in the IMGT-gapped V germline sequence alignment

property v_germ_aa_end_vdj¶: Position of the last nucleotide in the ungapped V germline sequence alignment

property v_germ_end_imgt¶: Position of the last nucleotide in the IMGT-gapped V germline sequence alignment

property v_germ_end_vdj¶: Position of the last nucleotide in the ungapped V germline sequence alignment

property v_seq_aa_end¶: Position of the last V nucleotide in the input sequence

property v_seq_end¶: Position of the last V nucleotide in the input sequence

class changeo.Receptor.ReceptorData¶

Bases: object

A class containing type conversion methods for Receptor data attributes

sequence_id¶

unique sequence identifier.

Type: str

rev_comp¶

whether the alignment is relative to the reverse compliment of the input sequence.

Type: bool

functional¶

whether sample V(D)J sequence is predicted to be functional.

Type: bool

in_frame¶

whether junction region is in-frame.

Type: bool

stop¶

whether a stop codon is present in the V(D)J sequence.

Type: bool

mutated_invariant¶

whether the conserved amino acids are mutated in the V(D)J sequence.

Type: bool

indels¶

whether the V(D)J nucleotide sequence contains insertions and/or deletions.

Type: bool

v_frameshift¶

whether the V segment contains a frameshift

Type: bool

sequence_input¶

input nucleotide sequence.

Type: Bio.Seq.Seq

sequence_vdj¶

Aligned V(D)J nucleotide sequence without IMGT-gaps.

Type: Bio.Seq.Seq

sequence_imgt¶

IMGT-gapped V(D)J nucleotide sequence.

Type: Bio.Seq.Seq

sequence_aa_input¶

input amino acid sequence.

Type: Bio.Seq.Seq

sequence_aa_vdj¶

Aligned V(D)J nucleotide sequence without IMGT-gaps.

Type: Bio.Seq.Seq

sequence_aa_imgt¶

IMGT-gapped V(D)J amino sequence.

Type: Bio.Seq.Seq

junction¶

ungapped junction region nucletide sequence.

Type: Bio.Seq.Seq

junction_aa¶

ungapped junction region amino acid sequence.

Type: Bio.Seq.Seq

junction_start¶

start positions of the junction in the input nucleotide sequence.

Type: int

junction_length¶

length of the junction in nucleotides.

Type: int

germline_vdj¶

full ungapped germline V(D)J nucleotide sequence.

Type: Bio.Seq.Seq

germline_vdj_d_mask¶

ungapped germline V(D)J nucleotides sequence with Ns masking the NP1-D-NP2 regions.

Type: Bio.Seq.Seq

germline_imgt¶

full IMGT-gapped germline V(D)J nucleotide sequence.

Type: Bio.Seq.Seq

germline_imgt_d_mask¶

IMGT-gapped germline V(D)J nucleotide sequence with ns masking the NP1-D-NP2 regions.

Type: Bio.Seq.Seq

germline_aa_vdj¶

full ungapped germline V(D)J amino acid sequence.

Type: Bio.Seq.Seq

germline_aa_imgt¶

full IMGT-gapped germline V(D)J amino acid sequence.

Type: Bio.Seq.Seq

v_call¶

V allele assignment(s).

Type: str

d_call¶

D allele assignment(s).

Type: str

j_call¶

J allele assignment(s).

Type: str

c_call¶

C region assignment.

Type: str

v_seq_start¶

position of the first V nucleotide in the input sequence (1-based).

Type: int

v_seq_length¶

number of V nucleotides in the input sequence.

Type: int

v_germ_start_imgt¶

position of the first V nucleotide in IMGT-gapped V germline sequence alignment (1-based).

Type: int

v_germ_length_imgt¶

length of the IMGT numbered germline V alignment.

Type: int

v_germ_start_vdj¶

position of the first nucleotide in ungapped V germline sequence alignment (1-based).

Type: int

v_germ_length_vdj¶

length of the ungapped germline V alignment.

Type: int

v_seq_aa_start¶

position of the first V amino acid in the amino acid input sequence (1-based).

Type: int

v_seq_aa_length¶

number of V amino acid in the amino acid input sequence.

Type: int

v_germ_aa_start_imgt¶

position of the first V amino acid in IMGT-gapped V germline amino acid alignment (1-based).

Type: int

v_germ_aa_length_imgt¶

length of the IMGT numbered germline V amino acid alignment.

Type: int

v_germ_aa_start_vdj¶

position of the first amino acid in ungapped V germline amino acid alignment (1-based).

Type: int

v_germ_aa_length_vdj¶

length of the ungapped germline V amino acid alignment.

Type: int

np1_start¶

position of the first untemplated nucleotide between the V and D segments in the input sequence (1-based).

Type: int

np1_length¶

number of untemplated nucleotides between the V and D segments.

Type: int

np1_aa_start¶

position of the first untemplated amino acid between the V and D segments in the input amino acid sequence (1-based).

Type: int

np1_aa_length¶

number of untemplated amino acids between the V and D segments.

Type: int

d_seq_start¶

position of the first D nucleotide in the input sequence (1-based).

Type: int

d_seq_length¶

number of D nucleotides in the input sequence.

Type: int

d_germ_start¶

position of the first nucleotide in D germline sequence alignment (1-based).

Type: int

d_germ_length¶

length of the germline D alignment.

Type: int

d_seq_aa_start¶

position of the first D amino acid in the input amino acidsequence (1-based).

Type: int

d_seq_aa_length¶

number of D amino acids in the input amino acid sequence.

Type: int

d_germ_aa_start¶

position of the first amino acid in D germline amino acid alignment (1-based).

Type: int

d_germ_aa_length¶

length of the germline D amino acid alignment.

Type: int

np2_start¶

position of the first untemplated nucleotide between the D and J segments in the input sequence (1-based).

Type: int

np2_length¶

number of untemplated nucleotides between the D and J segments.

Type: int

np2_aa_start¶

position of the first untemplated amino acid between the D and J segments in the input amino acid sequence (1-based).

Type: int

np2_aa_length¶

number of untemplated amino acid between the D and J segments.

Type: int

j_seq_start¶

position of the first J nucleotide in the input sequence (1-based).

Type: int

j_seq_length¶

number of J nucleotides in the input sequence.

Type: int

j_germ_start¶

position of the first nucleotide in J germline sequence alignment (1-based).

Type: int

j_germ_length¶

length of the germline J alignment.

Type: int

j_seq_aa_start¶

position of the first J amino acid in the input amino acidsequence (1-based).

Type: int

j_seq_aa_length¶

number of J amino acid in the input amino acidsequence.

Type: int

j_germ_aa_start¶

position of the first amino acid in J germline amino acid alignment (1-based).

Type: int

j_germ_aa_length¶

length of the germline J amino acid alignment.

Type: int

v_score¶

alignment score for the V.

Type: float

v_identity¶

alignment identity for the V.

Type: float

v_evalue¶

E-value for the alignment of the V.

Type: float

v_btop¶

BTOP for the alignment of the V.

Type: str

v_cigar¶

CIGAR for the alignment of the V.

Type: str

d_score¶

alignment score for the D.

Type: float

d_identity¶

alignment identity for the D.

Type: float

d_evalue¶

E-value for the alignment of the D.

Type: float

d_btop¶

BTOP for the alignment of the D.

Type: str

d_cigar¶

CIGAR for the alignment of the D.

Type: str

j_score¶

alignment score for the J.

Type: float

j_identity¶

alignment identity for the J.

Type: float

j_evalue¶

E-value for the alignment of the J.

Type: float

j_btop¶

BTOP for the alignment of the J.

Type: str

j_cigar¶

CIGAR for the alignment of the J.

Type: str

vdj_score¶

alignment score for the V(D)J.

Type: float

fwr1_imgt¶

IMGT-gapped FWR1 nucleotide sequence.

Type: Bio.Seq.Seq

fwr2_imgt¶

IMGT-gapped FWR2 nucleotide sequence.

Type: Bio.Seq.Seq

fwr3_imgt¶

IMGT-gapped FWR3 nucleotide sequence.

Type: Bio.Seq.Seq

fwr4_imgt¶

IMGT-gapped FWR4 nucleotide sequence.

Type: Bio.Seq.Seq

cdr1_imgt¶

IMGT-gapped CDR1 nucleotide sequence.

Type: Bio.Seq.Seq

cdr2_imgt¶

IMGT-gapped CDR2 nucleotide sequence.

Type: Bio.Seq.Seq

cdr3_imgt¶

IMGT-gapped CDR3 nucleotide sequence.

Type: Bio.Seq.Seq

cdr3_igblast¶

CDR3 nucleotide sequence assigned by IgBLAST.

Type: Bio.Seq.Seq

fwr1_aa_imgt¶

IMGT-gapped FWR1 amino acid sequence.

Type: Bio.Seq.Seq

fwr2_aa_imgt¶

IMGT-gapped FWR2 amino acid sequence.

Type: Bio.Seq.Seq

fwr3_aa_imgt¶

IMGT-gapped FWR3 amino acid sequence.

Type: Bio.Seq.Seq

fwr4_aa_imgt¶

IMGT-gapped FWR4 amino acid sequence.

Type: Bio.Seq.Seq

cdr1_aa_imgt¶

IMGT-gapped CDR1 amino acid sequence.

Type: Bio.Seq.Seq

cdr2_aa_imgt¶

IMGT-gapped CDR2 amino acid sequence.

Type: Bio.Seq.Seq

cdr3_aa_imgt¶

IMGT-gapped CDR3 amino acid sequence.

Type: Bio.Seq.Seq

cdr3_igblast_aa¶

CDR3 amino acid sequence assigned by IgBLAST.

Type: Bio.Seq.Seq

n1_length¶

M nucleotides 5’ of the D segment.

Type: int

n2_length¶

nucleotides 3’ of the D segment.

Type: int

p3v_length¶

palindromic nucleotides 3’ of the V segment.

Type: int

p5d_length¶

palindromic nucleotides 5’ of the D segment.

Type: int

p3d_length¶

palindromic nucleotides 3’ of the D segment.

Type: int

p5j_length¶

palindromic nucleotides 5’ of the J segment.

Type: int

d_frame¶

D segment reading frame.

Type: int

conscount¶

number of reads contributing to the UMI consensus sequence.

Type: int

dupcount¶

copy number of the sequence.

Type: int

umicount¶

number of UMIs representing the sequence.

Type: int

clone¶

clonal cluster identifier.

Type: str

cell¶

origin cell identifier.

Type: str

annotations¶

dictionary containing all unknown fields.

Type: dict

static aminoacid(v, deparse=False)¶

static double(v, deparse=False)¶

end_fields = {'cdr1_end': ('cdr1_start', 'cdr1_length'), 'cdr2_end': ('cdr2_start', 'cdr2_length'), 'cdr3_end': ('cdr3_start', 'cdr3_length'), 'd_germ_aa_end': ('d_germ_aa_start', 'd_germ_aa_length'), 'd_germ_end': ('d_germ_start', 'd_germ_length'), 'd_seq_aa_end': ('d_seq_aa_start', 'd_seq_aa_length'), 'd_seq_end': ('d_seq_start', 'd_seq_length'), 'fwr1_end': ('fwr1_start', 'fwr1_length'), 'fwr2_end': ('fwr2_start', 'fwr2_length'), 'fwr3_end': ('fwr3_start', 'fwr3_length'), 'fwr4_end': ('fwr4_start', 'fwr4_length'), 'j_germ_aa_end': ('j_germ_aa_start', 'j_germ_aa_length'), 'j_germ_end': ('j_germ_start', 'j_germ_length'), 'j_seq_aa_end': ('j_seq_aa_start', 'j_seq_aa_length'), 'j_seq_end': ('j_seq_start', 'j_seq_length'), 'junction_end': ('junction_start', 'junction_length'), 'v_alignment_aa_end': ('v_alignment_aa_start', 'v_alignment_aa_length'), 'v_alignment_end': ('v_alignment_start', 'v_alignment_length'), 'v_germ_aa_end_imgt': ('v_germ_aa_start_imgt', 'v_germ_aa_length_imgt'), 'v_germ_aa_end_vdj': ('v_germ_aa_start_vdj', 'v_germ_aa_length_vdj'), 'v_germ_end_imgt': ('v_germ_start_imgt', 'v_germ_length_imgt'), 'v_germ_end_vdj': ('v_germ_start_vdj', 'v_germ_length_vdj'), 'v_seq_aa_end': ('v_seq_aa_start', 'v_seq_aa_length'), 'v_seq_end': ('v_seq_start', 'v_seq_length')}¶

static identity(v, deparse=False)¶

static integer(v, deparse=False)¶

length_fields = {'cdr1_length': ('cdr1_start', 'cdr1_end'), 'cdr2_length': ('cdr2_start', 'cdr2_end'), 'cdr3_length': ('cdr3_start', 'cdr3_end'), 'd_germ_aa_length': ('d_germ_aa_start', 'd_germ_aa_end'), 'd_germ_length': ('d_germ_start', 'd_germ_end'), 'd_seq_aa_length': ('d_seq_aa_start', 'd_seq_aa_end'), 'd_seq_length': ('d_seq_start', 'd_seq_end'), 'fwr1_length': ('fwr1_start', 'fwr1_end'), 'fwr2_length': ('fwr2_start', 'fwr2_end'), 'fwr3_length': ('fwr3_start', 'fwr3_end'), 'fwr4_length': ('fwr4_start', 'fwr4_end'), 'j_germ_aa_length': ('j_germ_aa_start', 'j_germ_aa_end'), 'j_germ_length': ('j_germ_start', 'j_germ_end'), 'j_seq_aa_length': ('j_seq_aa_start', 'j_seq_aa_end'), 'j_seq_length': ('j_seq_start', 'j_seq_end'), 'junction_length': ('junction_start', 'junction_end'), 'v_alignment_aa_length': ('v_alignment_aa_start', 'v_alignment_aa_end'), 'v_alignment_length': ('v_alignment_start', 'v_alignment_end'), 'v_germ_aa_length_imgt': ('v_germ_aa_start_imgt', 'v_germ_aa_end_imgt'), 'v_germ_aa_length_vdj': ('v_germ_aa_start_vdj', 'v_germ_aa_end_vdj'), 'v_germ_length_imgt': ('v_germ_start_imgt', 'v_germ_end_imgt'), 'v_germ_length_vdj': ('v_germ_start_vdj', 'v_germ_end_vdj'), 'v_seq_aa_length': ('v_seq_aa_start', 'v_seq_aa_end'), 'v_seq_length': ('v_seq_start', 'v_seq_end')}¶

static logical(v, deparse=False)¶

static nucleotide(v, deparse=False)¶

parsers = {'c_call': 'identity', 'cdr1_aa_imgt': 'aminoacid', 'cdr1_imgt': 'nucleotide', 'cdr2_aa_imgt': 'aminoacid', 'cdr2_imgt': 'nucleotide', 'cdr3_aa_imgt': 'aminoacid', 'cdr3_igblast': 'nucleotide', 'cdr3_igblast_aa': 'aminoacid', 'cdr3_imgt': 'nucleotide', 'cell': 'identity', 'clone': 'identity', 'conscount': 'integer', 'd_btop': 'identity', 'd_call': 'identity', 'd_cigar': 'identity', 'd_evalue': 'double', 'd_frame': 'integer', 'd_germ_aa_length': 'integer', 'd_germ_aa_start': 'integer', 'd_germ_length': 'integer', 'd_germ_start': 'integer', 'd_identity': 'double', 'd_score': 'double', 'd_seq_aa_length': 'integer', 'd_seq_aa_start': 'integer', 'd_seq_length': 'integer', 'd_seq_start': 'integer', 'dupcount': 'integer', 'functional': 'logical', 'fwr1_aa_imgt': 'aminoacid', 'fwr1_imgt': 'nucleotide', 'fwr2_aa_imgt': 'aminoacid', 'fwr2_imgt': 'nucleotide', 'fwr3_aa_imgt': 'aminoacid', 'fwr3_imgt': 'nucleotide', 'fwr4_aa_imgt': 'aminoacid', 'fwr4_imgt': 'nucleotide', 'germline_aa_imgt': 'aminoacid', 'germline_aa_vdj': 'aminoacid', 'germline_imgt': 'nucleotide', 'germline_imgt_d_mask': 'nucleotide', 'germline_vdj': 'nucleotide', 'germline_vdj_d_mask': 'nucleotide', 'in_frame': 'logical', 'indels': 'logical', 'j_btop': 'identity', 'j_call': 'identity', 'j_cigar': 'identity', 'j_evalue': 'double', 'j_germ_aa_length': 'integer', 'j_germ_aa_start': 'integer', 'j_germ_length': 'integer', 'j_germ_start': 'integer', 'j_identity': 'double', 'j_score': 'double', 'j_seq_aa_length': 'integer', 'j_seq_aa_start': 'integer', 'j_seq_length': 'integer', 'j_seq_start': 'integer', 'junction': 'nucleotide', 'junction_aa': 'aminoacid', 'junction_length': 'integer', 'junction_start': 'integer', 'locus': 'identity', 'mutated_invariant': 'logical', 'n1_length': 'integer', 'n2_length': 'integer', 'np1_aa_length': 'integer', 'np1_aa_start': 'integer', 'np1_length': 'integer', 'np1_start': 'integer', 'np2_aa_length': 'integer', 'np2_aa_start': 'integer', 'np2_length': 'integer', 'np2_start': 'integer', 'p3d_length': 'integer', 'p3v_length': 'integer', 'p5d_length': 'integer', 'p5j_length': 'integer', 'rev_comp': 'logical', 'sequence_aa_imgt': 'aminoacid', 'sequence_aa_input': 'aminoacid', 'sequence_aa_vdj': 'aminoacid', 'sequence_id': 'identity', 'sequence_imgt': 'nucleotide', 'sequence_input': 'nucleotide', 'sequence_vdj': 'nucleotide', 'stop': 'logical', 'umicount': 'integer', 'v_btop': 'identity', 'v_call': 'identity', 'v_cigar': 'identity', 'v_evalue': 'double', 'v_frameshift': 'logical', 'v_germ_aa_length_imgt': 'integer', 'v_germ_aa_length_vdj': 'integer', 'v_germ_aa_start_imgt': 'integer', 'v_germ_aa_start_vdj': 'integer', 'v_germ_length_imgt': 'integer', 'v_germ_length_vdj': 'integer', 'v_germ_start_imgt': 'integer', 'v_germ_start_vdj': 'integer', 'v_identity': 'double', 'v_score': 'double', 'v_seq_aa_length': 'integer', 'v_seq_aa_start': 'integer', 'v_seq_length': 'integer', 'v_seq_start': 'integer', 'vdj_score': 'double'}¶

start_fields = {'cdr1_start': ('cdr1_length', 'cdr1_end'), 'cdr2_start': ('cdr2_length', 'cdr2_end'), 'cdr3_start': ('cdr3_length', 'cdr3_end'), 'd_germ_aa_start': ('d_germ_aa_length', 'd_germ_aa_end'), 'd_germ_start': ('d_germ_length', 'd_germ_end'), 'd_seq_aa_start': ('d_seq_aa_length', 'd_seq_aa_end'), 'd_seq_start': ('d_seq_length', 'd_seq_end'), 'fwr1_start': ('fwr1_length', 'fwr1_end'), 'fwr2_start': ('fwr2_length', 'fwr2_end'), 'fwr3_start': ('fwr3_length', 'fwr3_end'), 'fwr4_start': ('fwr4_length', 'fwr4_end'), 'j_germ_aa_start': ('j_germ_aa_length', 'j_germ_aa_end'), 'j_germ_start': ('j_germ_length', 'j_germ_end'), 'j_seq_aa_start': ('j_seq_aa_length', 'j_seq_aa_end'), 'j_seq_start': ('j_seq_length', 'j_seq_end'), 'junction_start': ('junction_length', 'junction_end'), 'v_alignment_aa_start': ('v_alignment_aa_length', 'v_alignment_aa_end'), 'v_alignment_start': ('v_alignment_length', 'v_alignment_end'), 'v_germ_aa_start_imgt': ('v_germ_aa_length_imgt', 'v_germ_aa_end_imgt'), 'v_germ_aa_start_vdj': ('v_germ_aa_length_vdj', 'v_germ_aa_end_vdj'), 'v_germ_start_imgt': ('v_germ_length_imgt', 'v_germ_end_imgt'), 'v_germ_start_vdj': ('v_germ_length_vdj', 'v_germ_end_vdj'), 'v_seq_aa_start': ('v_seq_aa_length', 'v_seq_aa_end'), 'v_seq_start': ('v_seq_length', 'v_seq_end')}¶

Using IgBLAST¶

Example data¶

We have hosted a small example data set resulting from the UMI barcoded MiSeq workflow described in the pRESTO documentation. In addition to the example FASTA files, we have included the standalone IgBLAST results. The files can be downloded from here:

Configuring IgBLAST¶

A collection of scripts for setting up the standalone IgBLAST database from the IMGT reference sequences are available on the Immcantation repository. To use these scripts, copy all the tools in the /scripts folder to a location in your PATH. At a minimum, you’ll need the following scripts:

fetch_igblastdb.sh
fetch_imgtdb.sh
clean_imgtdb.py
imgt2igblast.sh

Download and configure the IgBLAST and IMGT reference databases as follows, adjusting the version number to taste:

# Download and extract IgBLAST
VERSION="1.17.0"
wget ftp://ftp.ncbi.nih.gov/blast/executables/igblast/release/${VERSION}/ncbi-igblast-${VERSION}-x64-linux.tar.gz
tar -zxf ncbi-igblast-${VERSION}-x64-linux.tar.gz
cp ncbi-igblast-${VERSION}/bin/* ~/bin
# Download reference databases and setup IGDATA directory
fetch_igblastdb.sh -o ~/share/igblast
cp -r ncbi-igblast-${VERSION}/internal_data ~/share/igblast
cp -r ncbi-igblast-${VERSION}/optional_file ~/share/igblast
# Build IgBLAST database from IMGT reference sequences
fetch_imgtdb.sh -o ~/share/germlines/imgt
imgt2igblast.sh -i ~/share/germlines/imgt -o ~/share/igblast

Note

Several Immcantation tools require the observed V(D)J sequence (sequence_alignment) and associated germline fields (germline_alignment or germline_alignment_d_mask) to have gaps inserted to conform to the IMGT numbering scheme. Thus, when a tool such as MakeDb.py or CreateGermlines.py requires a reference sequence set as input, it will required the IMGT-gapped reference set. Meaning, the reference sequences that were downloaded using the fetch_imgtdb.sh script, or downloaded manually from the IMGT reference directory, rather than the final upgapped reference set required by IgBLAST.

See also

The provided scripts download only the mouse and human IMGT reference databases. See the IgBLAST documentation for instructions on how to build the database in a more general case. Shown below is an example of how to performed the same steps as the Immcantation scripts using a separately downloaded IMGT reference set and the scripts provided by IgBLAST. You must have all of the associated commands in your PATH and the appropriate directories created:

# V segment database
edit_imgt_file.pl IMGT_Human_IGHV.fasta > ~/share/igblast/fasta/imgt_human_ig_v.fasta
makeblastdb -parse_seqids -dbtype nucl -in ~/share/igblast/fasta/imgt_human_ig_v.fasta \
    -out ~/share/igblast/database/imgt_human_ig_v
# D segment database
edit_imgt_file.pl IMGT_Human_IGHD.fasta > ~/share/igblast/fasta/imgt_human_ig_d.fasta
makeblastdb -parse_seqids -dbtype nucl -in ~/share/igblast/fasta/imgt_human_ig_d.fasta \
    -out ~/share/igblast/database/imgt_human_ig_d
# J segment database
edit_imgt_file.pl IMGT_Human_IGHJ.fasta > ~/share/igblast/fasta/imgt_human_ig_j.fasta
makeblastdb -parse_seqids -dbtype nucl -in ~/share/igblast/fasta/imgt_human_ig_j.fasta \
    -out ~/share/igblast/database/imgt_human_ig_j

Once these databases are built for each segment they can be referenced when running IgBLAST.

Running IgBLAST¶

Change-O provides a simple wrapper script to run IgBLAST with the required options as the igblast subcommand of AssignGenes.py. This wrapper can be run as follows using the database built using the Immcantation scripts:

AssignGenes.py igblast -s HD13M.fasta -b ~/share/igblast \
    --organism human --loci ig --format blast

The optional --format blast argument defines the output format of IgBLAST. The default, blast, is the blocked tabular output provided by specifying the -outfmt '7 std qseq sseq btop' argument to IgBLAST. Specifying --format airr will output a tab-delimited file compliant with the AIRR Rearrangement schema defined by the AIRR Community. AIRR format support requires IgBLAST v1.9.0 or higher.

The -b ~/share/igblast argument specifies the path containing the database, internal_data, and optional_file directories required by IgBLAST. This option sets the IGDATA environment variable that controls where IgBLAST looks for internal database files. See the IgBLAST documentation for more details regarding the IGDATA environment variable.

See also

The AssignGenes.py IgBLAST wrapper provides limited functionality. For more control, IgBLAST should be run directly. The only strict requirement for compatibility with Changeo-O is that the output must either be an AIRR tab-delimited file (--outfmt 19) or a blast-style tabular output with the optional query sequence, subject sequence and BTOP fields (-outfmt '7 std qseq sseq btop'). An example of how to run IgBLAST directly is shown below:

export IGDATA=~/share/igblast
igblastn \
    -germline_db_V ~/share/igblast/database/imgt_human_ig_v\
    -germline_db_D ~/share/igblast/database/imgt_human_ig_d \
    -germline_db_J ~/share/igblast/database/imgt_human_ig_j \
    -auxiliary_data ~/share/igblast/optional_file/human_gl.aux \
    -domain_system imgt -ig_seqtype Ig -organism human \
    -outfmt '7 std qseq sseq btop' \
    -query HD13M.fasta \
    -out HD13M.fmt7

Processing the output of IgBLAST¶

Standalone IgBLAST blast-style tabular output is parsed by the igblast subcommand of MakeDb.py to generate the standardized tab-delimited database file on which all subsequent Change-O modules operate. In addition to the IgBLAST output (-i HD13M.fmt7), both the FASTA files input to IgBLAST (-s HD13M.fasta) and the IMGT-gapped reference sequences (-r IMGT_Human_IGHV.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta) must be provided to MakeDb.py:

MakeDb.py igblast -i HD13M.fmt7 -s HD13M.fasta \
    -r IMGT_Human_IGHV.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta \
    --extended

The optional --extended argument adds extra columns to the output database containing IMGT-gapped CDR/FWR regions and alignment metrics.

Warning

The references sequences you provide to MakeDb.py must contain IMGT-gapped V segment references, and these reference must be the same sequences used to build the IgBLAST reference database. If your IgBLAST germlines are not IMGT-gapped and/or they are not identical to those provided to MakeDb.py, then sequences which were assigned missing germlines will fail the parsing operation and the junction (CDR3) sequences will not be correct.

Parsing IMGT output¶

Example data¶

We have hosted a small example data set resulting from the UMI barcoded MiSeq workflow described in the pRESTO documentation. In addition to the example FASTA files, we have included the IMGT/HighV-QUEST results. The files can be downloded from here:

Reducing file size for submission to IMGT/HighV-QUEST¶

IMGT/HighV-QUEST currently limits the size of uploaded files to 500,000 sequences. To accomodate this limit, you can use the count subcommand of the pRESTO tool SplitSeq to divide your files into small pieces:

SplitSeq.py count -s file.fastq -n 500000 --fasta

The -n 500000 argument sets the maximum number of sequences in each file and the --fasta argument tells the tool to output a FASTA, rather than FASTQ, formatted file suitable for upload to IMGT/HighV-QUEST.

See also

For additional details see the corresponding example in the pRESTO documentation

Processing the output of IMGT/HighV-QUEST¶

The output from IMGT/HighV-QUEST may be parsed via the imgt subcommand of MakeDb.py to generate the standardized tab-delimited database file on which all subsequent Change-O modules operate. Processing the IMGT output requires either the compressed output file (.zip or .txz) or an uncompressed folder containing the 1_Summary, 2_IMGT-gapped, 3_Nt-sequences and 6_Junction files (-i HD13M.txz). Additionally, it is recommended that you provide the FASTA file that was submitted to HighV-QUEST (-s HD13M.fasta), as this will allow MakeDb.py to correct the changes HighV-QUEST makes to the sequence identifier and add additional columns corresponding any annotations generated by pRESTO:

MakeDb.py imgt -i HD13M.txz -s HD13M.fasta --extended

The optional --extended argument add extra columns to the output database containing IMGT-gapped CDR/FWR regions and alignment metrics.

Merging processed IMGT/HighV-QUEST output¶

If you previously split files for submission to IMGT/HighV-QUEST, you can run each partition through MakeDb.py individually and merge the resulting output files using the merge subcommand of ParseDb.py:

MakeDb.py imgt -i part1.txz -s part1.fasta -o part1.tsv
MakeDb.py imgt -i part2.txz -s part2.fasta -o part2.tsv
ParseDb.py merge -d part1.tsv part2.tsv -o merged.tsv

Parsing 10X Genomics V(D)J data¶

New

We have an updated tutorial covering the processing of 10x Genomics VDJ data with Change-O and SCOPer. You can also follow the steps below to process 10x VDJ data using methods available in Change-O.

Example data¶

10X Genomics provides an example data set of Ig V(D)J processed by the Cell Ranger pipeline, which is available for download from their Single Cell Immune Profiling support site.

Converting 10X V(D)J data into the AIRR Community standardized format¶

To process 10X V(D)J data, a combination of AssignGenes.py and MakeDb.py can be used to generate a TSV file compliant with the AIRR Community Rearrangement schema that incorporates annotation information provided by the Cell Ranger pipeline. The --10x filtered_contig_annotations.csv specifies the path of the contig annotations file generated by cellranger vdj, which can be found in the outs directory.

Generate AIRR Rearrangement data from the 10X V(D)J FASTA files using the steps below:

AssignGenes.py igblast -s filtered_contig.fasta -b ~/share/igblast \
   --organism human --loci ig --format blast
MakeDb.py igblast -i filtered_contig_igblast.fmt7 -s filtered_contig.fasta \
   -r IMGT_Human_*.fasta --10x filtered_contig_annotations.csv --extended

all_contig.fasta can be exchanged for filtered_contig.fasta, and all_contig_annotations.csv can be exchanged for filtered_contig_annotations.csv.

Warning

The resulting table overwrites the V, D and J gene assignments generated by Cell Ranger and uses those generated by IgBLAST or IMGT/HighV-QUEST instead.

See also

To process mouse data and/or TCR data alter the --organism and --loci arguments to AssignGenes.py accordingly (e.g., --organism mouse, --loci tcr) and use the appropriate V, D and J IMGT reference databases (e.g., IMGT_Mouse_TR*.fasta)

See the IgBLAST usage guide for further details regarding the setup and use of IgBLAST with Change-O.

Identifying clones from B cells in AIRR formatted 10X V(D)J data¶

Splitting into separate light and heavy chain files¶

To group B cells into clones from AIRR Rearrangement data, the output from MakeDb.py must be parsed into a light chain file and a heavy chain file:

ParseDb.py select -d 10x_igblast_db-pass.tsv -f locus -u "IGH" \
        --logic all --regex --outname heavy
ParseDb.py select -d 10x_igblast_db-pass.tsv -f locus -u "IG[LK]" \
        --logic all --regex --outname light

Assign clonal groups to the heavy chain data¶

The heavy chain file must then be clonally clustered separately. See Clustering sequences into clonal groups for how to use DefineClones.py to assign clonal cluster annotations to the IGH file.

Correct clonal groups based on light chain data¶

DefineClones.py currently does not support light chain cloning. However, cloning can be performed after heavy chain cloning using light_cluster.py provided on the Immcantation Bitbucket repository in the scripts directory:

light_cluster.py -d heavy_select-pass_clone-pass.tsv -e light_select-pass.tsv \
        -o 10X_clone-pass.tsv

Here, heavy_select-pass_clone-pass.tsv refers to the cloned heavy chain AIRR Rearrangement file, light_select-pass.tsv refers to the light chain file, and 10X_clone-pass.tsv is the resulting output file.

The algorithm will (1) remove cells associated with more than one heavy chain and (2) correct heavy chain clone definitions based on an analysis of the light chain partners associated with the heavy chain clone.

Note

By default, light_chain.py expects the AIRR Rearrangement columns:

v_call
j_call
junction_length
umi_count
cell_id
clone_id

To process legacy Change-O formatted data add the --format changeo argument:

light_cluster.py -d heavy_select-pass_clone-pass.tab -e light_select-pass.tab \
    -o 10X_clone-pass.tab --format changeo

Which expects the following Change-O columns:

V_CALL
J_CALL
JUNCTION_LENGTH
UMICOUNT
CELL
CLONE

Filtering records¶

The ParseDb.py tool provides a basic set of operations for manipulating Change-O database files from the commandline, including removing or updating rows and columns.

Removing non-productive sequences¶

After building a Change-O database from either IMGT/HighV-QUEST or IgBLAST output, you may wish to subset your data to only productive sequences. This can be done in one of two roughly equivalent ways using the ParseDb.py tool:

ParseDb.py select -d HD13M_db-pass.tsv -f productive -u T
ParseDb.py split -d HD13M_db-pass.tsv -f productive

The first line above uses the select subcommand to output a single file labeled parse-select containing only records with the value of T (-u T) in the productive column (-f productive).

Alternatively, the second line above uses the split subcommand to output multiple files with each file containing records with one of the values found in the productive column (-f productive). This will generate two files labeled productive-T and productive-F.

Removing disagreements between the C-region primers and the reference alignment¶

If you have data that includes both heavy and light chains in the same library, the V-segment and J-segment alignments from IMGT/HighV-QUEST or IgBLAST may not always agree with the isotype assignments from the C-region primers. In these cases, you can filter out such reads with the select subcommand of ParseDb.py. An example function call using an imaginary file db.tsv is provided below:

ParseDb.py select -d db.tsv -f v_call j_call c_call -u "IGH" \
    --logic all --regex --outname heavy
ParseDb.py select -d db.tsv -f v_call j_call c_call -u "IG[LK]" \
    --logic all --regex --outname light

These commands will require that all of the v_call, j_call and c_call fields (-f v_call j_call c_call and --logic all) contain the string IGH (lines 1-2) or one of IGK or IGL (lines 3-4). The --regex argument allows for partial matching and interpretation of regular expressions. The output from these two commands are two files, one containing only heavy chains (heavy_parse-select.tsv) and one containg only light chains (light_parse-select.tsv).

Exporting records to FASTA files¶

You may want to use external tools, or tools from pRESTO, on your Change-O result files. The ConvertDb.py tool provides two options for exporting data from tab-delimited files to FASTA format.

Standard FASTA¶

The fasta subcommand allows you to export sequences and annotations to FASTA formatted files in the pRESTO annototation scheme:

ConvertDb.py fasta -d HD13M_db-pass.tsv --if sequence_id \
    --sf sequence_alignment --mf v_call duplicate_count

Where the column containing the sequence identifier is specified by --if sequence_id, the nucleotide sequence column is specified by --sf sequence_id, and additional annotations to be added to the sequence header are specified by --mf v_call duplicate_count.

BASELINe FASTA¶

The baseline subcommand generates a FASTA derivative format required by the BASELINe web tool. Generating these files is similar to building standard FASTA files, but requires a few more options. An example function call using an imaginary file db.tsv is provided below:

ConvertDb.py baseline -d db.tsv --if sequence_id \
    --sf sequence_alignment --mf v_call duplicate_count \
    --cf clone_id --gf germline_alignment_d_mask

The additional arguments required by the baseline subcommand include the clonal grouping (--cf clone_id) and germline sequence (--gf germline_alignment_d_mask) columns added by the DefineClones and CreateGermlines tasks, respectively.

Note

The baseline subcommand requires the CLONE column to be sorted. DefineClones.py generates a sorted CLONE column by default. However, you needed to alter the order of the CLONE column at some point, then you can re-sort the clonal assignments using the sort subcommand of ParseDb.py. An example function call using an imaginary file db.tsv is provided below:

ParseDb.py sort -d db.tsv -f clone_id

Which will sort records by the value in the clone_id column (-f clone_id).

Clustering sequences into clonal groups¶

Example data¶

We have hosted a small example data set resulting from the UMI barcoded MiSeq workflow described in the pRESTO documentation. The files can be downloded from here:

The following examples use the HD13M_db-pass.tsv database file provided in the example bundle, which has already undergone the IMGT/IgBLAST parsing and filtering operations.

Determining a clustering threshold¶

Before running DefineClones.py, it is important to determine an appropriate threshold for trimming the hierarchical clustering into B cell clones. The distToNearest function in the SHazaM R package calculates the distance between each sequence in the data and its nearest-neighbor. The resulting distribution should be bimodal, with the first mode representing sequences with clonal relatives in the dataset and the second mode representing singletons. The ideal threshold for separating clonal groups is the value that separates the two modes of this distribution and can be found using the findThreshold function in the SHazaM R package. The distToNearest function allows selection of all parameters that are available in DefineClones.py. Using the length normalization parameter ensures that mutations are weighted equally regardless of junction sequence length. The distance to nearest-neighbor distribution for the example data is shown below. The threshold is approximately 0.16 - indicated by the red dotted line.

See also

For additional details see the vignette on tuning clonal assignment thresholds.

Assigning clones¶

There are several parameter choices when grouping Ig sequences into B cell clones. The argument --act set accounts for ambiguous V gene and J gene calls when grouping similar sequences. The distance metric --model ham is nucleotide Hamming distance. Because the threshold was generated using length normalized distances, the --norm len argument is selected with the previously determined threshold --dist 0.16:

DefineClones.py -d HD13M_db-pass.tsv --act set --model ham \
    --norm len --dist 0.16

Note

Because T cells don’t undergo SHM, non-zero nucleotide distances suggest sequences orginate from a different ancestor. To identify TCR clones, use –dist 0 or a very low distance value to allow for sequencing error.

Reconstructing germline sequences¶

Example data¶

We have hosted a small example data set resulting from the UMI barcoded MiSeq workflow described in the pRESTO documentation. The files can be downloaded from here:

The following examples use the HD13M_db-pass.tsv AIRR Rearrangement file provided in the example bundle, which has already undergone the IMGT/IgBLAST parsing and filtering operations.

Adding germline sequences to the database¶

The CreateGermlines.py tool is used to reconstruct the germline V(D)J sequence, from which the Ig lineage and mutations can be inferred. In addition to the alignment information parsed by MakeDb.py to generate the initial database, CreateGermlines.py also requires the set of germline sequences that were used for the alignment passed to the -r argument. In the case of V-segment germlines, the reference sequences must be IMGT-gapped. Because the D-segment call for B cell receptor alignments is often low confidence, the default germline format (-g dmask) places Ns in the N/P and D-segments of the junction region rather than using the D-segment assigned during reference alignment; this can be modified to generate a complete germline (-g full) or a V-segment only germline (-g vonly) if you wish. The command below adds the germline sequence to the germline_alignment_d_mask column of the output database:

CreateGermlines.py -d HD13M_db-pass.tsv -g dmask \
    -r IMGT_Human_IGHV.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta

Alternatively, if you have run the clonal assignment task prior to invoking CreateGermlines.py, then adding the --cloned argument is recommended, as this will generate a single germline of consensus length for each clone:

CreateGermlines.py -d HD13M_db-pass_clone-pass.tsv -g dmask --cloned \
    -r IMGT_Human_IGHV.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta

Important

The germline set passed to -r must contain the complete set of germlines used by the reference alignment software (IMGT/HighV-QUEST or IgBLAST). If alleles called by the aligner are missing from the reference set, they will not be successfully processed. Additionally, the V-segment reference set must contain IMGT-gapped sequences to properly reconstruct germlines, even if the reference alignment was performed on ungapped sequences.

Note

While MakeDb.py provides the ihmm subcommand to parse alignment output generated by iHMMuneAlign, there is insufficient information to successfully reconstruct germline sequences for all cases using CreateGermlines.py.

See also

The TIgGER R package provided tools for identifying novel polymorphisms and building a personalized germline database. To use the germline corrections provided by TIgGER you would replace the V-segment germline file with the one generated by genotypeFasta (-r IGHV_genotype.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta) and specify the genotyped V-segment column (--vf v_call_genotyped):

CreateGermlines.py -d genotyped.tsv -g dmask --vf v_call_genotyped \
    -r IGHV_genotype.fasta IMGT_Human_IGHD.fasta IMGT_Human_IGHJ.fasta

IgPhyML lineage tree analysis¶

IgPhyML is a program designed to build phylogenetic trees and test evolutionary hypotheses regarding B cell affinity maturation.

The biology of B cell somatic hypermutation (SHM) violates important assumptions in most standard phylogenetic substitution models; further, while most phylogenetics programs are designed to analyze single lineages, B cell repertoires typically contain thousands of lineages. IgPhyML addresses both of these issues by implementing substitution models that correct for the context-sensitive nature of SHM, and combines information from multiple lineages to give more precisely estimated repertoire-wide model parameter estimates.

An in-depth description of IgPhyML installation and usage can be found at the IgPhyML website.

Quick start¶

Once installed, IgPhyML can be run through BuildTrees by specifying the --igphyml option. IgPhyML is easiest to run through the Immcantation Docker image. If this is not possible, these instructions require Change-O 0.4.6 or higher, Alakazam 0.3.0 or higher, and IgPhyML to be installed, with the executable in your PATH variable.

The following commands should work as a first pass on many reasonably sized datasets, but if you really want to understand what’s going on or make sure what you’re doing makes sense, please check out the IgPhyML website.

Build trees and estimate model parameters¶

Download the IgPhyML repository, move to the examples folder, and run BuildTrees:

# Clone IgPhyML repository to get example files
git clone https://bitbucket.org/kleinstein/igphyml

# Move to examples directory
cd igphyml/examples

# Run BuildTrees
BuildTrees.py -d example.tsv --outname ex --log ex.log --collapse \
    --sample 3000 --igphyml --clean all --nproc 1

This command processes an AIRR-formatted dataset of BCR sequences that have been clonally clustered with germlines reconstructed. It then quickly builds trees using the GY94 model and, using these fixed topologies, estimates HLP19 model parameters. This can be sped up by increasing the --nproc option. Subsampling using the --sample option in isn’t strictly necessary, but IgPhyML will run slowly when applied to large datasets. Here, the --collapse flag is used to collapse identical sequences. This is highly recommended because identical sequences slow down calculations without affecting likelihood values in IgPhyML.

Visualize results¶

The output file of the above command can be read using the readIgphyml function of Alakazam. After opening an R session in the examples subfolder, enter the following commands. Note that when using the Docker container, you’ll need to run dev.off() after plotting the tree to create a pdf plot in the examples directory:

library(alakazam)
library(igraph)

db = readIgphyml("ex_igphyml-pass.tab")

# Plot largest lineage tree
plot(db$trees[[1]],layout=layout_as_tree)

# Show HLP10 parameters
print(t(db$param[1,]))
CLONE         "REPERTOIRE"
NSEQ          "4"
NSITE         "107"
TREE_LENGTH   "0.286"
LHOOD         "-290.7928"
KAPPA_MLE     "2.266"
OMEGA_FWR_MLE "0.5284"
OMEGA_CDR_MLE "2.3324"
WRC_2_MLE     "4.8019"
GYW_0_MLE     "3.4464"
WA_1_MLE      "5.972"
TW_0_MLE      "0.8131"
SYC_2_MLE     "-0.99"
GRS_0_MLE     "0.2583"

Lineage tree of example clone.¶

To visualize a larger dataset with bigger trees, and bifurcating tree topologies, again open an R session in the examples directory:

library(alakazam)
library(ape)

db = readIgphyml("sample1_igphyml-pass.tab",format="phylo")

# Plot largest lineage tree
plot(ladderize(db$trees[[1]]),cex=0.7,no.margin=TRUE)

Phylo-formatted lineage tree of a larger B cell clone.¶

Generating MiAIRR compliant GenBank/TLS submissions¶

MiAIRR¶

The MiAIRR standard (minimal information about adaptive immune receptor repertoires) is a minimal reporting standard for experiments using sequencing-based technologies to study adaptive immune receptors (T and B cell receptors). The current version (1.0) of the standard was published in Rubelt et al, 2017 and accepted by the general assembly at the annual AIRR Community meeting in December 2017.

MiAIRR recommends submission of raw read data to the Sequence Read Archive (SRA) and submission of processed and annotated data to the Targeted Locus Study (TLS) section of GenBank.

This example will cover generation of files for submission to TLS starting from Change-O formatted data. For complete details of the required and optional elements of the TLS submission see the AIRR Standards documentation site.

Special attention should be paid to the REQUIRED elements. Note that GenBank expects there to be a CDS element that corresponds to the JUNCTION. If submitting single-cell heavy:light paired BCR data, GenBank expects separate files for the heavy, the kappa, and the lambda chains. Note that even though the kappa and the lambda chain sequences should be in separate files, their misc_feature comments should both read immunoglobulin light chain variable region, per AIRR standard requirements. In addition, every effort should be made to make sure that the values of the attributes for GenBank submission match those of the BioSample attributes. In particular, if the BioSample specifies a strain value (e.g. for mouse data), then a strain attribute MUST be included when preparing GenBank submission, and that value MUST match the BioSample value.

Example data¶

We have hosted a small example data set resulting from the UMI barcoded MiSeq workflow described in the pRESTO documentation. The files can be downloded from here: