Data loading

AnalyzAIRR proposes two different data loading processes depending on the alignment files format.

l <- list.files(system.file(file.path('extdata/MiAIRR'),
                            package='AnalyzAIRR'),
                            full.names = T)

metadata <- read.table(system.file(file.path('extdata/sampledata.txt'),
                            package='AnalyzAIRR'), 
                            header=T,
                            row.names = 1)
metadata$cell_subset <- factor(metadata$cell_subset)
metadata$sex <- factor(metadata$sex)

RepSeqData <- readAIRRSet(fileList = l,
                          fileFormat = "MiAIRR",
                          chain = "TRA",
                          sampleinfo = metadata,
                          keep.ambiguous = FALSE,
                          keep.unproductive=FALSE,
                          filter.singletons=FALSE,
                          aa.th=8,
                          outFiltered=TRUE,
                          raretab=FALSE,
                          cores=1L)
#> Running on 1 cores...
#> Loading and filtering sequences...
#> Assembling sequences...
#> Creating a RepSeqExperiment object...
#> Done.

The RepSeqExperiment object

assay(RepSeqData)

mData(RepSeqData)

oData(RepSeqData)$filtered

History(RepSeqData)

Filtering

• Filter out sequences based on their count in each sample

  filterCount() filters out sequences having a count below the chosen threshold. The function returns a RepSeqExperiment object that can be used to perform in-depth analyses on the remaining dataset

• Extract sequences that are either shared across samples, or private to a single sample

  getPublic() allows to subset a RepSeqExperiment object in order to keep sequences that are shared by:

  • samples belonging to a specified group
  • samples within the whole dataset if the group parameter is not specified.
    

The function returns a RepSeqExperiment object with the selected clones.

# Get clones present in at least 50% of the samples belonging to the amTreg group
getPublic(x=RepSeqData, level="aaClone", group=c("cell_subset", "amTreg"))

Similarly, it is possible to extract private sequences with getPrivate(). If the parameter singletons is set to TRUE, only private sequences with a count of 1 will be returned.

Note: All the above-mentioned functions can be applied at any of the following repertoire level: aaClone, ntClone, aaCDR3 and ntCDR3.

• Filter out one or multiple sequence(s)

  filterSequence() filters out specified sequences in all samples or a specified group of samples.

• Extract most frequent sequences

  getTopSequences() allows the extraction of the top most expressed repertoire levels. Users can choose the proportion of the top sequences to extract using the prop parameter as well as the group of samples to analyze.

• Filter productive or unproductive sequences

  getProductive() and getUnproductive() allow the extraction of productive and unproductive sequences respectively in case no filters were applied during the building of the RepSeqExperiment object.

• Exclude one or multiple samples from a RepSeqExperiment object

  dropSamples() offers the possibility to filter out one or multiple repertoires by specifying their corresponding sample_id. These repertoires can be ones identified as outliers during the exploratory analysis routine described further below.

Normalization

AnalyzAIRR proposes two normalization strategies:

• Down-sampling

This strategy can be applied when the studied samples largely differ in their repertoire sizes.

By using the sampleRepSeqExp() function, users can choose the value to which all the samples are downsampled. If not specified, the lowest number of sequences across the dataset will be used. This function returns a new RepSeqExperiment object with the downsized data.

RepSeqData_ds<- sampleRepSeqExp(x=RepSeqData, sample.size=50000)

• Shannon-based normalization

This strategy adapted in (Chaara et al., 2018) can be used to eliminate “uninformative” sequences resulting from experimental noise. It is particularly efficient when applied on small samples as it corrects altered count distributions caused by a high-sequencing depth.

The function ShannonNorm() allows the application of this strategy without the need to specify any parameter and returns a new RepSeqExperiment object with the corrected data.

RepSeqData_sh<- ShannonNorm(x=RepSeqData)

Eliminated sequences in each sample are stored in the otherData slot.

Descriptive statistics

• Visualize the metaData statistics

Statistics in the metaData slot, calculated during the building of the RepSeqExperiment object, can be visualized for each sample using the plotStatistics() function with the colorBy parameter set to the sample id column. Groups can be separated using the facetBy parameter.

plotStatistics(x = RepSeqData,
               stat = "nSequences",
               colorBy = "sample_id")

Figure 1: A barplot showing the number of TCR sequences identified in each sample.

plotStatistics(x = RepSeqData,
               stat = "ntClone",
               colorBy = "sample_id")

Figure 2: A barplot showing the number of ntClones identified in each sample.

Comparison of descriptive statistics

The majority of the functions described in the Exploratory analysis section can be applied to compare groups of samples and perform statistical tests.

• Compare metaData statistics

  The calculated statistics in the metaData slot can be compared between groups of samples. This is possible using the plotStatistics() function by setting the grouped parameter to TRUE and specifying one group column from the metadata slot to compare in the colorBy parameter, and up to two columns in the facetBy parameter.

  In case multiple groups are specified, the first group will be compared and the others will be used as facets.

plotStatistics(x = RepSeqData,
               stat = "ntClone",
               colorBy = "cell_subset",
               grouped=TRUE)

Figure 14: Comparison of the the number of ntClones between amTreg and nTreg samples. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

plotStatistics(x = RepSeqData,
               stat = "ntClone",
               colorBy = "cell_subset", facetBy="sex",
               grouped=TRUE)

Figure 15: Comparison of the the number of ntClones between amTreg and nTreg samples within each sex. Boxplots represent the median across all samples belonging to the same specified group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

• Compare gene usage

  plotGeneUsage() allows the comparison of V or J gene usages between groups of samples using the colorBy parameter

plotGeneUsage(x=RepSeqData, level="J", scale="frequency", colorBy="cell_subset")

• Compare repertoire diversities

  For The plotDiversity() function, the grouped parameter should be set to TRUE, and the groups to be analysed should be specified in the colorBy parameter.

plotDiversity(x = RepSeqData,
              level = "aaClone",
              index = "gini",
              colorBy = "cell_subset",
              grouped=TRUE)

Figure 17: The Gini index is calculated at the aaClone level and compared between cell subsets. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

plotDiversity(x = RepSeqData,
              level = "aaClone",
              index = "gini",
              colorBy = "cell_subset", facetBy = "sex",
              grouped = TRUE)

Figure 18: The Gini index is calculated at the aaClone level and compared between cell subsets within each sex. Boxplots represent the median across all samples belonging to the same group. A Wilcoxon test is applied and adjusted p-values using the Holm method are shown.

The Renyi index also supports the grouped parameter.

plotRenyiIndex(x = RepSeqData,
                  level = "aaClone",
                  colorBy = "cell_subset",
                  shapeBy = "sex",
                  grouped = TRUE)

• Compare clonal distributions

  Similarly to the analysis applied at the single-sample level, plotIntervals() allows to plot and compare clonal distributions within defined count intervals between experimental groups.

plotIntervals(x=RepSeqData, level="aaClone", colorBy = "cell_subset", facetBy="sex", fractions="count", grouped=TRUE)

Similarity analysis

The repertoire sharing at any level evaluates the degree of convergence between repertoires and experimental conditions. Different statistical methods are proposed herein to evaluate this convergence.

• Compute the number of shared sequences

  The number of shared sequences, at any level, between samples belonging for instance to the same experimental group can be obtained using the plotVenn() function. If sampleNames is not specified, the first 3 samples in the datasets will be analyzed.

ctrnames <- rownames(mData(RepSeqData))[which(mData(RepSeqData)[,"cell_subset" ] %in% "nTreg")]
plotVenn(x = RepSeqData,
           level = "aaClone",
           sampleNames = ctrnames)

• Assess the correlation between a pair of samples

  The correlation between a pair of repertoires can be calculated using the plotScatter() function by simply specifying a pair of sampleNames to compare.

plotScatter(x = RepSeqData,
            level = "V",
            scale = "frequency",
            sampleNames = c("tripod-30-813","tripod-31-846"))

Figure 22: A biplot comparing the V gene usage between two defined samples. A linear model is fitted on the data. The V usage is significantly similar between the two samples.

plotScatter(x = RepSeqData,
            level = "aaClone",
            scale = "frequency",
            sampleNames = c("tripod-30-813","tripod-31-846"))

Figure 23: A biplot comparing the clonal expression between two defined samples. A linear model is fitted on the data. Clonal overlap and expression are significantly low between the compared samples.

• Compute dissimilarity distances

  AnalyzAIRR proposes a list of dissimilarity indices, each taking into account different parameters. The proposed methods include:

  • The Jaccard similarity: a measure of similarity between sample sets defined as the size of the intersection divided by the size of the union of the sample sets.

  • The Morisita-Horn similarity: a measure of similarity that tends to be over-sensitive to abundant species.

  Details on these indices and others can be found here

  These distances can be calculated at any level of the repertoire and can be:

  • Visualized As a dissimilarity heatmap by setting plot to Heatmap. This function performs a hierarchical clustering on the calculated distance scores using the method specified in the clustering parameter.
  • Used to perform a multidimensional scaling analysis by setting plot to MDS. In this case, no clustering method is needed.

plotDissimilarity(x = RepSeqData,
                  level = "aaClone",
                  method = "morisita",
                  clustering = "ward.D",
                  colorBy = c("sex","cell_subset"),
                  plot="Heatmap")

Figure 24: A heatmap showing the Morisita distances calculated at the aaClone level between all pairwise samples. The ward D method was used to perform the hierarchical clustering.

plotDissimilarity(x = RepSeqData,
                 level = "aaClone",
                 method = "morisita",
                 colorBy = "cell_subset",
                 plot="MDS"
                 )

Figure 25: Morisita distances were calculated at the aaClone level between all pairwise samples within the dataset and used to perform a multidimensional scaling analysis.

Differential analysis

Differentially expressed genes or sequences can be identified using diffExpGroup. The experimental groups to be compared can be specified with the group parameter, and the function outputs the statistics calculated for each gene/sequence.

In the following example, we identify over-expressed VJ combinations within the amTreg samples compared to the nTregs.

DS <- diffExpGroup(x = RepSeqData,
                   colGrp = "cell_subset" ,
                   level = "VJ",
                   group = c("cell_subset", "amTreg", "nTreg"))