TCR beta repertoire exploratory analysis

In tcR: Advanced Data Analysis of Immune Receptor Repertoires

Before analysis

Before running this pipeline you should do next steps:

1. Save your parsed MiTCR data to the immdata variable (it could be a mitcr data frame or a list).

immdata <- parse.folder('/home/username/mitcrdata/')

1. Save the immdata variable to the some folder as the .rda file.

save(immdata, file = '/home/username/immdata.rda')

1. In the code block below change the path string to the path to yours immdata.rda file. After that click the Knit HTML button to start analysis and the script will make a html report.

load('../data/twb.rda')
immdata <- twb[1:2]
library(tcR)

1. Friendly advice: run the pipeline on the top N sequences first and configure sizes of figures.

N <- 50000
immdata <- head(immdata, N)

Data statistics

cloneset.stats(immdata)
repseq.stats(immdata)

Segments' statistics

V-segment usage

if (has.class(immdata, 'list')) {
  for (i in 1:length(immdata)) {
    plot(vis.gene.usage(immdata[[i]], HUMAN_TRBV, .main = paste0(names(immdata)[i], ' ', 'V-usage')))
  }
} else {
  vis.gene.usage(immdata, HUMAN_TRBV, .coord.flip=F)
}

J-segment usage

if (has.class(immdata, 'list')) {
  for (i in 1:length(immdata)) {
    plot(vis.gene.usage(immdata[[i]], HUMAN_TRBJ, .main = paste0(names(immdata)[i], ' ', 'J-usage')))
  }
} else {
  vis.gene.usage(immdata, HUMAN_TRBJ, .coord.flip=F)
}

CDR3 length and read distributions

if (has.class(immdata, 'list')) {
  for (i in 1:length(immdata)) {
    plot(vis.count.len(immdata[[i]], .name = paste0(names(immdata)[i], ' ', 'CDR3 length')))
  }
} else {
  vis.count.len(immdata)
}

if (has.class(immdata, 'list')) {
  for (i in 1:length(immdata)) {
    plot(vis.number.count(immdata[[i]], .name = paste0(names(immdata)[i], ' ', 'read histogram')))
  }
} else {
  vis.number.count(immdata)
}

Proportions of the most abundant clones

vis.top.proportions(immdata)

Most frequent kmers

kms <- get.kmers(immdata, .verbose = F)
vis.kmer.histogram(kms, .position = 'fill')

Rarefaction analysis

clmn <- 'Read.count'
if (has.class(immdata, 'list')) {
  if (!is.na(immdata[[1]]$Umi.count[1])) {
    clmn <- 'Umi.count'
  }
} else {
  if (!is.na(immdata$Umi.count[1])) {
    clmn <- 'Umi.count'
  }
}
vis.rarefaction(rarefaction(immdata, .col = clmn, .verbose = F), .log = T)

tcR documentation built on July 2, 2020, 3:18 a.m.