Shazam: Quantification of selection pressure

Constructing clonal consensus sequences

Individual sequences within clonal groups are not, strictly speaking, independent events and it is generally appropriate to only analyze selection pressures on an effective sequence for each clonal group. The collapseClones function provides one strategy for generating an effective sequences for each clone. It reduces the input database to one row per clone and appends clonal_sequence and clonal_germline columns which contain the consensus sequences for each clone.

# Collapse clonal groups into single sequences
clones <- collapseClones(ExampleDb, cloneColumn="clone_id", 
                         sequenceColumn="sequence_alignment", 
                         germlineColumn="germline_alignment_d_mask", 
                         regionDefinition=IMGT_V, 
                         method="thresholdedFreq", minimumFrequency=0.6,
                         includeAmbiguous=FALSE, breakTiesStochastic=FALSE, 
                         nproc=1)

Incorporating lineage information

For each clone, lineage information can be incorporated following these steps:

# Subset to sequences with clone_id=3170
db_3170 <- subset(ExampleDb, clone_id == 3170)
dim(db_3170)
colnames(db_3170)

# Generate a ChangeoClone object for lineage construction 
clone_3170 <- makeChangeoClone(db_3170, seq="sequence_alignment", germ="germline_alignment")

# Run the lineage reconstruction
dnapars_exec <- "/usr/local/bin/dnapars"
graph_3170 <- buildPhylipLineage(clone_3170, dnapars_exec, rm_temp=TRUE)

# Generating a data.frame from the lineage tree graph object, 
# and merge it with clone data.frame
graph_3170_df <- makeGraphDf(graph_3170, clone_3170)
dim(graph_3170_df)
colnames(graph_3170_df)

makeGraphDf creates a data.frame with the column parent_sequence, which can be used to analyze mutations for each sequence relative to their parent_sequence.

Calculating selection in multiple steps

Following construction of an effective sequence for each clone, the observed and expected mutation counts are calculated for each sequence in the clonal_sequence column relative to the clonal_germline. observedMutations is used to calculate the number of observed mutations and expectedMutations calculates the expected frequency of mutations. The underlying targeting model for calculating expectations can be specified using the targetingModel parameter. In the example below, the default HH_S5F is used. Column names for sequence and germline sequence may also be passed in as parameters if they differ from the Change-O defaults.

Mutations are counted by these functions separately for complementarity determining (CDR) and framework (FWR) regions. The regionDefinition argument defines whether these regions are handled separately, and where the boundaries lie. There are several built-in region definitions in the shazam package, both dependent upon the V segment being IMGT-gapped:

• IMGT_V: All regions in the V segment, excluding CDR3, grouped as either CDR or FWR.
• IMGT_V_BY_REGIONS: The CDR1, CDR2, FWR1, FWR and FWR3 regions in the V segment (no CDR3) treated as individual regions.
• IMGT_VDJ: All regions, including CDR3 and FWR4, grouped as either CDR or FWR. This RegionDefinition is initially empty, and one is created on the fly for each set of clonally related sequences.
• IMGT_VDJ_BY_REGIONS: CDR1, CDR2, CDR3, FWR1, FWR, FWR3 and FWR4 regions treated as individual regions. This RegionDefinition is initially empty, and one is created on the fly for each set of clonally related sequences.

Users may define other region sets and boundaries by creating a custom RegionDefinition object.

# Count observed mutations and append mu_count columns to the output
observed <- observedMutations(clones, 
                              sequenceColumn="clonal_sequence",
                              germlineColumn="clonal_germline",
                              regionDefinition=IMGT_V, nproc=1)
# Count expected mutations and append mu_expected columns to the output
expected <- expectedMutations(observed, 
                              sequenceColumn="clonal_sequence",
                              germlineColumn="clonal_germline",
                              targetingModel=HH_S5F,
                              regionDefinition=IMGT_V, nproc=1)

The counts of observed and expected mutations can be combined to test for selection using calcBaseline. The statistical framework used to test for selection based on mutation counts can be specified using the testStatistic parameter.

# Calculate selection scores using the output from expectedMutations
baseline <- calcBaseline(expected, testStatistic="focused", 
                         regionDefinition=IMGT_V, nproc=1)

## calcBaseline will use existing observed and expected mutations, in the fields: mu_count_cdr_r, mu_count_cdr_s, mu_count_fwr_r, mu_count_fwr_s and mu_expected_cdr_r, mu_expected_cdr_s, mu_expected_fwr_r, mu_expected_fwr_s

Calculating selection in one step

It is not required for observedMutation and expectedMutations to be run prior to calcBaseline. If the output of these two steps does not appear in the input data.frame, then calcBaseline will call the appropriate functions prior to calculating selection scores.

# Calculate selection scores from scratch
baseline <- calcBaseline(clones, testStatistic="focused", 
                         regionDefinition=IMGT_V, nproc=1)

Using alternative mutation definitions and models

The default behavior of observedMutations and expectedMutations, and by extension calcBaseline, is to define a replacement mutation in the usual way - any change in the amino acid of a codon is considered a replacement mutation. However, these functions have a mutationDefinition argument which allows these definitions to be changed by providing a MutationDefinition object that contains alternative replacement and silent criteria. shazam provides the following built-in MutationDefinition objects:

• CHARGE_MUTATIONS: Amino acid mutations are defined by changes in side chain charge class.
• HYDROPATHY_MUTATIONS: Amino acid mutations are defined by changes in side chain hydrophobicity class.
• POLARITY_MUTATIONS: Amino acid mutations are defined by changes in side chain polarity class.
• VOLUME_MUTATIONS: Amino acid mutations are defined by changes in side chain volume class.

The default behavior of expectedMutations is to use the human 5-mer mutation model, HH_S5F. Alternative SHM targeting models can be provided using the targetingModel argument.

# Calculate selection on charge class with the mouse 5-mer model
baseline_mk_rs5nf <- calcBaseline(clones, testStatistic="focused", 
                         regionDefinition=IMGT_V, 
                         targetingModel=MK_RS5NF,
                         mutationDefinition=CHARGE_MUTATIONS,
                         nproc=1)

Grouping by a single annotation

The following example generates a single selection PDF for each unique annotation in the sample_id column.

# Combine selection scores by time-point
grouped_1 <- groupBaseline(baseline, groupBy="sample_id")

Subsetting and grouping by multiple annotations

Grouping by multiple annotations follows the sample procedure as a single annotation by simply adding columns to the groupBy argument. Subsetting the data can be performed before or after generating selection PDFs via calcBaseline. However, note that subsetting may impact the clonal representative sequences generated by collapseClones. In the following example, subsetting precedes the collapsing of clonal groups.

# Subset the original data to switched isotypes
db_sub <- subset(ExampleDb, c_call %in% c("IGHM", "IGHG"))

# Collapse clonal groups into single sequence
clones_sub <- collapseClones(db_sub, cloneColumn="clone_id",
                             sequenceColumn="sequence_alignment",
                             germlineColumn="germline_alignment_d_mask",
                             regionDefinition=IMGT_V, 
                             method="thresholdedFreq", minimumFrequency=0.6,
                             includeAmbiguous=FALSE, breakTiesStochastic=FALSE, 
                             nproc=1)

# Calculate selection scores from scratch
baseline_sub <- calcBaseline(clones_sub, testStatistic="focused", 
                             regionDefinition=IMGT_V, nproc=1)

## calcBaseline will calculate observed and expected mutations for clonal_sequence using clonal_germline as a reference.

# Combine selection scores by time-point and isotype
grouped_2 <- groupBaseline(baseline_sub, groupBy=c("sample_id", "c_call"))

Convolving variables at multiple levels

To make selection comparisons using two levels of variables, you would need two iterations of groupings, where the first iteration of groupBaseline groups on both variables, and the second iteration groups on the “outer” variable. For example, if a data set has both case and control subjects annotated in status and subject columns, then generating convolved PDFs for each status would be performed as:

# First group by subject and status
subject_grouped <- groupBaseline(baseline, groupBy=c("status", "subject"))

# Then group the output by status
status_grouped <- groupBaseline(subject_grouped, groupBy="status")

Testing the difference in selection PDFs between groups

The testBaseline function will perform significance testing between two grouped BASELINe PDFs, by region, and return a data.frame with the following information:

• region: The sequence region, such as cdr and fwr.
• test: The name of the two groups compared.
• pvalue: Two-sided p-value for the comparison.
• fdr: FDR corrected p-value.

testBaseline(grouped_1, groupBy="sample_id")

##   region       test     pvalue        fdr
## 1    cdr -1h != +7d 0.04494357 0.08988715
## 2    fwr -1h != +7d 0.49922881 0.49922881