2.1.1 Loading in the data

This dataset was generated by Grun et al. (2016) using the CEL-seq protocol with unique molecular identifiers (UMIs) and ERCC spike-ins. Count tables were obtained from the NCBI Gene Expression Omnibus using the accession number above.

library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
grun.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE81nnn/GSE81076/suppl/GSE81076%5FD2%5F3%5F7%5F10%5F17%2Etxt%2Egz"))

We first read the table into memory.

gse81076.df <- read.table(grun.fname, sep='\t', 
    header=TRUE, stringsAsFactors=FALSE, row.names=1)
dim(gse81076.df)

## [1] 20148  1728

head(rownames(gse81076.df))

## [1] "A1BG-AS1__chr19" "A1BG__chr19"     "A1CF__chr10"     "A2M-AS1__chr12" 
## [5] "A2ML1__chr12"    "A2MP1__chr12"

head(colnames(gse81076.df))

## [1] "D2ex_1" "D2ex_2" "D2ex_3" "D2ex_4" "D2ex_5" "D2ex_6"

Unfortunately, the data and metadata are all mixed together in this file. As a result, we need to manually extract the metadata from the column names.

donor.names <- sub("^(D[0-9]+).*", "\\1", colnames(gse81076.df))
table(donor.names)

## donor.names
##    D101    D102  D10631     D17   D1713 D172444      D2      D3     D71 
##      96      96      96     288      96      96      96     480      96 
##     D72     D73     D74 
##      96      96      96

plate.id <- sub("^D[0-9]+(.*)_.*", "\\1", colnames(gse81076.df))
table(plate.id)

## plate.id
##      All1 All2 TGFB  en1  en2  en3  en4   ex 
##  864   96   96   96   96   96   96   96  192

Another irritating feature of this dataset is that gene symbols were supplied, rather than stable identifiers such as Ensembl. Stable identifiers are desirable as they facilitate reliable cross-referencing of gene identities between data sets11 The two data sets used in this workflow come from the same provider, who at least uses gene symbols consistently. So, technically, we could have skipped the conversion here. Nonetheless, we have chosen to perform it to demonstrate how one would deal with more heterogeneous data sources (e.g., mixtures of Ensembl identifiers and gene symbols, or multiple synonymous gene symbols).. We convert all row names to Ensembl identifiers, removing NA or duplicated entries (with the exception of spike-in transcripts).

gene.symb <- gsub("__chr.*$", "", rownames(gse81076.df))
is.spike <- grepl("^ERCC-", gene.symb)
table(is.spike)

## is.spike
## FALSE  TRUE 
## 20064    84

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.symb, keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.symb[is.spike]

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
gse81076.df <- gse81076.df[keep,]
rownames(gse81076.df) <- gene.ids[keep]
summary(keep)

##    Mode   FALSE    TRUE 
## logical    2372   17776

We create a SingleCellExperiment object to store the counts and metadata together. This reduces the risk of book-keeping errors in later steps of the analysis. Note that we re-identify the spike-in rows, as the previous indices would have changed after the subsetting.

library(SingleCellExperiment)
sce.gse81076 <- SingleCellExperiment(list(counts=as.matrix(gse81076.df)),
    colData=DataFrame(Donor=donor.names, Plate=plate.id),
    rowData=DataFrame(Symbol=gene.symb[keep]))
isSpike(sce.gse81076, "ERCC") <- grepl("^ERCC-", rownames(gse81076.df)) 
sce.gse81076  

## class: SingleCellExperiment 
## dim: 17776 1728 
## metadata(0):
## assays(1): counts
## rownames(17776): ENSG00000268895 ENSG00000121410 ... ENSG00000074755
##   ENSG00000036549
## rowData names(1): Symbol
## colnames(1728): D2ex_1 D2ex_2 ... D17TGFB_95 D17TGFB_96
## colData names(2): Donor Plate
## reducedDimNames(0):
## spikeNames(1): ERCC

2.1.2 Quality control and normalization

We compute quality control (QC) metrics for each cell (McCarthy et al. 2017) and identify cells with low library sizes, low numbers of expressed genes, or high ERCC content.

library(scater)
sce.gse81076 <- calculateQCMetrics(sce.gse81076, compact=TRUE)
QC <- sce.gse81076$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3)
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
    HighSpike=sum(high.spike, na.rm=TRUE))

##   LowLib LowNgenes HighSpike
## 1     55       130       388

Cells with extreme values for these QC metrics are presumed to be of low quality and are removed. A more thorough analysis would examine the distributions of these QC metrics beforehand, but we will skip that step for brevity here.

discard <- low.lib | low.genes | high.spike
sce.gse81076 <- sce.gse81076[,!discard]
summary(discard)

##    Mode   FALSE    TRUE 
## logical    1292     436

We compute size factors for the endogenous genes using the deconvolution method (Lun, Bach, and Marioni 2016). This is done with pre-clustering by quickCluster() to avoid pooling together very different cells.

library(scran)
library(BiocSingular)
set.seed(1000) # for irlba. 
clusters <- quickCluster(sce.gse81076, use.ranks=FALSE, BSPARAM=IrlbaParam())
table(clusters)

## clusters
##   1   2   3   4   5   6   7 
## 106 196 273 299 153 163 102

sce.gse81076 <- computeSumFactors(sce.gse81076, min.mean=0.1, clusters=clusters)
summary(sizeFactors(sce.gse81076))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.003541 0.438547 0.791872 1.000000 1.290387 8.275084

We also compute size factors for the spike-in transcripts (Lun et al. 2017). Recall that we set general.use=FALSE to ensure that the spike-in size factors are only applied to the spike-in transcripts.

sce.gse81076 <- computeSpikeFactors(sce.gse81076, general.use=FALSE)
summary(sizeFactors(sce.gse81076, "ERCC"))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01042 0.57782 0.88699 1.00000 1.27765 7.43998

We then compute normalized log-expression values for use in downstream analyses.

sce.gse81076 <- normalize(sce.gse81076)

2.1.3 Identifying highly variable genes

We identify highly variable genes (HVGs) using trendVar() and decomposeVar(), using the variances of spike-in transcripts to model technical noise. We set block= to ensure that uninteresting differences between plates or donors do not inflate the variance. The small discrepancy in the fitted trend in Figure 1 is caused by the fact that the trend is fitted robustly to the block-wise variances of the spike-ins, while the variances shown are averaged across blocks and not robust to outliers.

block <- paste0(sce.gse81076$Plate, "_", sce.gse81076$Donor)
fit <- trendVar(sce.gse81076, block=block, parametric=TRUE) 
dec <- decomposeVar(sce.gse81076, fit)

plot(dec$mean, dec$total, xlab="Mean log-expression", 
    ylab="Variance of log-expression", pch=16)
is.spike <- isSpike(sce.gse81076)
points(dec$mean[is.spike], dec$total[is.spike], col="red", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)

We order genes by decreasing biological component, revealing some usual suspects such as insulin and glucagon. We will be using this information later when performing feature selection prior to running mnnCorrect().

dec.gse81076 <- dec
dec.gse81076$Symbol <- rowData(sce.gse81076)$Symbol
dec.gse81076 <- dec.gse81076[order(dec.gse81076$bio, decreasing=TRUE),]
head(dec.gse81076)

## DataFrame with 6 rows and 7 columns
##                             mean            total              bio
##                        <numeric>        <numeric>        <numeric>
## ENSG00000254647 2.84438696097646 6.37771955862257 5.93492163812922
## ENSG00000129965 1.89026042119725   6.029667887451 5.58238903990959
## ENSG00000115263 4.01989463994372 5.74235119909272 5.50960551193726
## ENSG00000118271 3.67399583064614 5.77781565535565 5.47705987083213
## ENSG00000115386 4.24266345753652 5.39869761580333 5.16308284971761
## ENSG00000164266 3.03518659963559  5.4170378405534 5.00991839667255
##                              tech   p.value       FDR      Symbol
##                         <numeric> <numeric> <numeric> <character>
## ENSG00000254647 0.442797920493351         0         0         INS
## ENSG00000129965 0.447278847541415         0         0    INS-IGF2
## ENSG00000115263 0.232745687155467         0         0         GCG
## ENSG00000118271 0.300755784523523         0         0         TTR
## ENSG00000115386 0.235614766085716         0         0       REG1A
## ENSG00000164266 0.407119443880852         0         0      SPINK1

2.2.1 Loading in the data

This dataset was generated by Muraro et al. (2016) using the CEL-seq2 protocol with unique molecular identifiers (UMIs) and ERCC spike-ins. Count tables were obtained from the NCBI Gene Expression Omnibus using the accession number above.

muraro.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE85nnn/GSE85241/suppl",
    "GSE85241%5Fcellsystems%5Fdataset%5F4donors%5Fupdated%2Ecsv%2Egz"))

We first read the table into memory.

gse85241.df <- read.table(muraro.fname, sep='\t', 
    header=TRUE, row.names=1, stringsAsFactors=FALSE)
dim(gse85241.df)

## [1] 19140  3072

head(rownames(gse85241.df))

## [1] "A1BG-AS1__chr19" "A1BG__chr19"     "A1CF__chr10"     "A2M-AS1__chr12" 
## [5] "A2ML1__chr12"    "A2M__chr12"

head(colnames(gse85241.df))

## [1] "D28.1_1" "D28.1_2" "D28.1_3" "D28.1_4" "D28.1_5" "D28.1_6"

We extract the metadata from the column names.

donor.names <- sub("^(D[0-9]+).*", "\\1", colnames(gse85241.df))
table(donor.names)

## donor.names
## D28 D29 D30 D31 
## 768 768 768 768

plate.id <- sub("^D[0-9]+\\.([0-9]+)_.*", "\\1", colnames(gse85241.df))
table(plate.id)

## plate.id
##   1   2   3   4   5   6   7   8 
## 384 384 384 384 384 384 384 384

Yet again, gene symbols were supplied instead of Ensembl or Entrez identifiers. We convert all row names to Ensembl identifiers, removing NA or duplicated entries (with the exception of spike-in transcripts).

gene.symb <- gsub("__chr.*$", "", rownames(gse85241.df))
is.spike <- grepl("^ERCC-", gene.symb)
table(is.spike)

## is.spike
## FALSE  TRUE 
## 19059    81

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.symb, keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.symb[is.spike]

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
gse85241.df <- gse85241.df[keep,]
rownames(gse85241.df) <- gene.ids[keep]
summary(keep)

##    Mode   FALSE    TRUE 
## logical    2223   16917

We create a SingleCellExperiment object to store the counts and metadata together.

sce.gse85241 <- SingleCellExperiment(list(counts=as.matrix(gse85241.df)),
    colData=DataFrame(Donor=donor.names, Plate=plate.id),
    rowData=DataFrame(Symbol=gene.symb[keep]))
isSpike(sce.gse85241, "ERCC") <- grepl("^ERCC-", rownames(gse85241.df)) 
sce.gse85241  

## class: SingleCellExperiment 
## dim: 16917 3072 
## metadata(0):
## assays(1): counts
## rownames(16917): ENSG00000268895 ENSG00000121410 ... ENSG00000074755
##   ENSG00000036549
## rowData names(1): Symbol
## colnames(3072): D28.1_1 D28.1_2 ... D30.8_95 D30.8_96
## colData names(2): Donor Plate
## reducedDimNames(0):
## spikeNames(1): ERCC

2.2.2 Quality control and normalization

We compute QC metrics for each cell and identify cells with low library sizes, low numbers of expressed genes, or high ERCC content.

sce.gse85241 <- calculateQCMetrics(sce.gse85241, compact=TRUE)
QC <- sce.gse85241$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3)
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
    HighSpike=sum(high.spike, na.rm=TRUE))

##   LowLib LowNgenes HighSpike
## 1    577       669       695

Low-quality cells are defined as those with extreme values for these QC metrics and are removed.

discard <- low.lib | low.genes | high.spike
sce.gse85241 <- sce.gse85241[,!discard]
summary(discard)

##    Mode   FALSE    TRUE 
## logical    2346     726

We compute size factors for the endogenous genes and spike-in transcripts, and use them to compute log-normalized expression values.

set.seed(1000)
clusters <- quickCluster(sce.gse85241, use.ranks=FALSE, BSPARAM=IrlbaParam())
table(clusters)

## clusters
##   1   2   3   4   5   6   7   8   9 
## 257 393 224 316 226 426 129 203 172

sce.gse85241 <- computeSumFactors(sce.gse85241, min.mean=0.1, clusters=clusters)
summary(sizeFactors(sce.gse85241))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0845  0.5328  0.8151  1.0000  1.2090 15.3318

sce.gse85241 <- computeSpikeFactors(sce.gse85241, general.use=FALSE)
summary(sizeFactors(sce.gse85241, "ERCC"))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.09295 0.61309 0.88902 1.00000 1.27519 4.04643

sce.gse85241 <- normalize(sce.gse85241)

2.2.3 Identifying highly variable genes

We fit a trend to the spike-in variances as previously described, allowing us to model the technical noise for each gene (Figure 2). Again, we set block= to ensure that uninteresting differences between plates or donors do not inflate the variance.

block <- paste0(sce.gse85241$Plate, "_", sce.gse85241$Donor)
fit <- trendVar(sce.gse85241, block=block, parametric=TRUE) 
dec <- decomposeVar(sce.gse85241, fit)
plot(dec$mean, dec$total, xlab="Mean log-expression", 
    ylab="Variance of log-expression", pch=16)
is.spike <- isSpike(sce.gse85241)
points(dec$mean[is.spike], dec$total[is.spike], col="red", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)

We order genes by decreasing biological component, as described above.

dec.gse85241 <- dec
dec.gse85241$Symbol <- rowData(sce.gse85241)$Symbol
dec.gse85241 <- dec.gse85241[order(dec.gse85241$bio, decreasing=TRUE),]
head(dec.gse85241)

## DataFrame with 6 rows and 7 columns
##                             mean            total              bio
##                        <numeric>        <numeric>        <numeric>
## ENSG00000115263 7.67110494913732 6.64869412141223 6.61993009624089
## ENSG00000089199 4.63141777116535 6.49101561206233 6.37005001267999
## ENSG00000169903 3.02363719480826 6.60374469250074 6.25171241896403
## ENSG00000254647 2.00812251287305 6.46803316161483 5.91300126053973
## ENSG00000118271 7.34417012391726 5.82255630100473  5.7911444440125
## ENSG00000171951  4.2002996520299 5.63389093078673 5.47030577732233
##                               tech   p.value       FDR      Symbol
##                          <numeric> <numeric> <numeric> <character>
## ENSG00000115263 0.0287640251713411         0         0         GCG
## ENSG00000089199  0.120965599382341         0         0        CHGB
## ENSG00000169903  0.352032273536714         0         0      TM4SF4
## ENSG00000254647  0.555031901075099         0         0         INS
## ENSG00000118271 0.0314118569922332         0         0         TTR
## ENSG00000171951   0.16358515346441         0         0        SCG2

6.2.1 Motivation

In more complex experiments, we may know beforehand that certain sets of batches are more similar to each other. We might then want to merge those similar batches before attempting the more difficult merges involving batches with different cell type composition and/or gene expression. Examples include:

• Merging batches that represent replicate experiments from the same condition, prior to merging across conditions.
• Merging batches generated with the same scRNA-seq technology prior to merging across technologies.

This strategy limits encourages detection of correct MNN pairs as similar batches should have more shared populations, By comparison, performing the more difficult merges first is more likely to introduce errors whereby distinct subpopulations are incorrectly placed together. This unnecessarily propagates the error to later steps as the initial merge is used as a reference for subsequent merges.

6.2.2 Setting up the inputs

Hierarchical merging can be achieved through multiple calls to fastMNN() with progressively merged batches. To illustrate, assume that we want to remove Donor effects within each batch prior to merging across batches. We split up the cells in GSE85241 according to the donor of origin, creating one SingleCellExperiment object for each donor.

all.donors <- unique(rescaled.gse85241$Donor)
table(rescaled.gse85241$Donor)

## 
## D28 D29 D30 D31 
## 340 604 689 713

by.donor.85241 <- vector("list", length(all.donors))
names(by.donor.85241) <- sort(all.donors)
for (x in all.donors) {
    by.donor.85241[[x]] <- rescaled.gse85241[,rescaled.gse85241$Donor==x]
}

We repeat this process for GSE81076. For demonstration purposes, we will aggregate some of the donors together to ensure that there are enough cells in each level for MNN detection.

adj.donors <- c(D101="A", D102="A", D10631="A",
    D17="B", D1713="B", D172444="B",
    D2="C", D3="C",
    D71="D", D72="D", D73="D", D74="D")[rescaled.gse81076$Donor]
table(adj.donors)

## adj.donors
##   A   B   C   D 
## 162 464 320 346

all.donors <- unique(adj.donors)
by.donor.81076 <- vector("list", length(all.donors))
names(by.donor.81076) <- sort(all.donors)
for (x in all.donors) {
    by.donor.81076[[x]] <- rescaled.gse81076[,adj.donors==x]
}

We use the multiBatchPCA() function to perform a PCA across all batches to be merged. This ensures that all cells are placed onto the same coordinate space, which would obviously not be possible if a PCA was performed for each batch separately. Specifically, multiBatchPCA() performs a modified PCA to ensure that each supplied matrix contributes equally to the definition of the PC space. This avoids problems with imbalances in the number of cells across batches, meanining that smaller batches (possibly with unique cell types) are not ignored.

all.batches <- c(by.donor.85241, by.donor.81076)

# Cosine normalizing for consistency with fastMNN() defaults.
all.logcounts <- lapply(all.batches, logcounts)
scaled <- lapply(all.logcounts, batchelor::cosineNorm) 

set.seed(1000) # for irlba.
pc.all <- do.call(batchelor::multiBatchPCA, c(scaled, 
    list(d=50, BSPARAM=IrlbaParam(deferred=TRUE))
))
names(pc.all)

## [1] "D28" "D29" "D30" "D31" "A"   "B"   "C"   "D"

Comments from Aaron:

• Here, we have applied multiBatchPCA() to the batch-level inputs for convenience. It is also possible to supply donor-level matrices to equalize contributions across donors.
• Recall that we ran multiBatchNorm() earlier to generate the rescaled.gse* objects. As a result, cosine normalization is not technically necessary, as all batches should be on the same scale already (see ?cosineNorm for a discussion of this). Nonetheless, we run it here for consistency with our previous fastMNN() call where cosine normalization is turned on by default.

6.2.3 Performing progressive merges

We pass the matrices corresponding to the donors in GSE85241 to fastMNN(), setting pc.input=TRUE to indicate that dimensionality reduction has already been performed. This uses the first donor to define the reference space22 Which can be changed with auto.order=, if so desired. and merges cells from all other donors to the first. In this manner, we remove donor effects within the GSE85241 batch.

pcs.85241 <- pc.all[seq_along(by.donor.85241)]
mnn.out.85241 <- do.call(batchelor::fastMNN, c(pcs.85241, list(pc.input=TRUE)))
Rle(mnn.out.85241$batch)

## character-Rle of length 2346 with 4 runs
##   Lengths:   340   604   689   713
##   Values : "D28" "D29" "D30" "D31"

Note that when pc.input=TRUE, a DataFrame is returned instead of a SingleCellExperiment. This reflects the fact that per-gene identities are lost when PCs are supplied instead of per-gene expression vectors33 Though it is entirely possible to reconstruct these using rotation vectors, see below.. We repeat this for the donors in the GSE81076 batch.

pcs.81076 <- tail(pc.all, length(by.donor.81076))
mnn.out.81076 <- do.call(batchelor::fastMNN, c(pcs.81076, list(pc.input=TRUE)))
Rle(mnn.out.81076$batch)

## character-Rle of length 1292 with 4 runs
##   Lengths: 162 464 320 346
##   Values : "A" "B" "C" "D"

The next step is to merge the two batches together. To do this, we simply repeat the fastMNN() call with the donor-corrected values for each batch. Again, we need to set pc.input=TRUE to prevent the function from unnecessarily (and incorrectly) repeating the cosine normalization and PCA steps on the corrected values. We use a larger k in the final fastMNN() call to improve the robustness of MNN detection to outlier cells on the edges of each cluster.

mnn.out3 <- batchelor::fastMNN(
    GSE81076=mnn.out.81076,
    GSE85241=mnn.out.85241, 
    pc.input=TRUE, k=100 # see comments below.
) 

Rle(mnn.out3$batch) # by dataset

## character-Rle of length 3638 with 2 runs
##   Lengths:       1292       2346
##   Values : "GSE81076" "GSE85241"

Rle(c(mnn.out.81076$batch, mnn.out.85241$batch)) # by donor

## character-Rle of length 3638 with 8 runs
##   Lengths:   162   464   320   346   340   604   689   713
##   Values :   "A"   "B"   "C"   "D" "D28" "D29" "D30" "D31"

This yields a final corrected expression matrix where both within-batch donor effects and batch effects have been corrected. We examine the quality of each of the merge steps with t-SNE plots (Figure 5). Within each batch, the donors are generally well-mixed, and the final merge is consistent with Figure 3.

set.seed(1000)
par(mfrow=c(1,3))

library(Rtsne)
tout.85241 <- Rtsne(mnn.out.85241$corrected, pca=FALSE)
plot(tout.85241$Y[,1], tout.85241$Y[,2], main="GSE85241 donors",
    col=as.factor(mnn.out.85241$batch), xlab="tSNE1", ylab="tSNE2")

tout.81076 <- Rtsne(mnn.out.81076$corrected, pca=FALSE)
plot(tout.81076$Y[,1], tout.81076$Y[,2], main="GSE81076 donors",
    col=as.factor(mnn.out.81076$batch), xlab="tSNE1", ylab="tSNE2")

tout.all <- Rtsne(mnn.out3$corrected, pca=FALSE)
plot(tout.all$Y[,1], tout.all$Y[,2], main="Final",
    col=as.factor(mnn.out3$batch), xlab="tSNE1", ylab="tSNE2")

Comments from Aaron:

• Here, we have applied multiBatchPCA() to the batch-level inputs for convenience. It is also possible to supply donor-level matrices to equalize contributions across donors.
• One might ask why a larger k was not needed in Figure 3 as well. This was probably because the outliers were masked by inter-donor heterogeneity when the t-SNE plots were generated without removing the donor effects.

• In this specific example, cells from the same donor will occupy contiguous rows in the mnn.out3$corrected matrix. However, this may not have been the case for the original ordering of cells in each SingleCellExperiment. This requires some extra account-keeping to match up the final corrected matrix to the original ordering, e.g., when cross-referencing to metadata.

  original.plate <- unlist(lapply(rescaled, "[[", i="Plate"))
  original.names <- unlist(lapply(rescaled, colnames))
  
  # Needs unique names: trigger error otherwise.
  stopifnot(anyDuplicated(original.names)==0L)
  
  m <- match(rownames(mnn.out3$corrected), original.names)
  new.plate <- original.plate[m]