processing_all.Rmd

---
title: "ProcessingPigBAL run"
author: "Arianne Richard"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Intro

This code will run an initial processing of the porcine BAL 10x Genomics scRNAseq data generated at the Pirbright Institute in July and October 2022. bcl files were demultiplexed using `cellranger mkfastq` and aligned/counted using `cellranger count` (cellranger-7.0.0) using default parameters and Sus scrofa genome (genome assembly 11.1, Ensembl release 107).

Note saturation estimated by cellranger as follows:
1. 60.8
2. 55.9
3. 40.3
4. 62.1
5. 69.6
6. 64.1
7. 77.2
8. 86.0
9. 61.4
10. 75.7
11. 78.5
12. 75.8

# Processing droplets

We will read in raw counts matrices generated by cellranger to allow us to filter out empty droplets, remove ambient RNA and removed barcode-swapped droplets using tools from DropletUtils. (cellranger filtered matrices are also probably fine, but it's nice to be able to do a bit of extra filtering ourselves.)

## Reading in the data

In this data set, we have 12 samples, run in pairs on individual lanes of the Chromium controller. The two libraries in each pair were then pooled for a run of sequencing. The first set (1-6) were run on 3 separate days in July 2022, and the second set (7-12) were run on 3 days in October 2022. 

```{r}

library(SingleCellExperiment)
library(DropletUtils)

samps <- paste0('sample', as.character(c(1:12)))

dirs <- paste0(samps, '/raw_feature_bc_matrix')

sce_list <- vector('list', length=length(samps))

for(i in 1:length(samps)){
  sce_list[[i]] <- read10xCounts(samples=dirs[i], sample.names=samps[i], col.names=TRUE)
}


```

## Removing swapped drops

Now we want to look for any droplets that have swapped cell barcodes due to use of patterned flow cell sequencing. Two samples were sequenced in each run, and so swapping could occur anywhere within the sequencing run but not outside of it.

```{r}
## run 1
fs <- paste0('sample', as.character(c(1,2)), '/molecule_info.h5')
run1 <- swappedDrops(samples=fs)
sum(colSums(run1$cleaned[[1]]))/sum(colSums(counts(sce_list[[1]])))
sum(colSums(run1$cleaned[[2]]))/sum(colSums(counts(sce_list[[2]])))

## run 2
fs <- paste0('sample', as.character(c(3,4)), '/molecule_info.h5')
run2 <- swappedDrops(fs)
sum(colSums(run2$cleaned[[1]]))/sum(colSums(counts(sce_list[[3]])))
sum(colSums(run2$cleaned[[2]]))/sum(colSums(counts(sce_list[[4]])))

## run 3
fs <- paste0('sample', as.character(c(5,6)), '/molecule_info.h5')
run3 <- swappedDrops(fs)
sum(colSums(run3$cleaned[[1]]))/sum(colSums(counts(sce_list[[5]])))
sum(colSums(run3$cleaned[[2]]))/sum(colSums(counts(sce_list[[6]])))

## run 4
fs <- paste0('sample', as.character(c(7,8)), '/molecule_info.h5')
run4 <- swappedDrops(fs)
sum(colSums(run4$cleaned[[1]]))/sum(colSums(counts(sce_list[[7]])))
sum(colSums(run4$cleaned[[2]]))/sum(colSums(counts(sce_list[[8]])))

## run 5
fs <- paste0('sample', as.character(c(9,10)), '/molecule_info.h5')
run5 <- swappedDrops(fs)
sum(colSums(run5$cleaned[[1]]))/sum(colSums(counts(sce_list[[9]])))
sum(colSums(run5$cleaned[[2]]))/sum(colSums(counts(sce_list[[10]])))

## run 6
fs <- paste0('sample', as.character(c(11,12)), '/molecule_info.h5')
run6 <- swappedDrops(fs)
sum(colSums(run6$cleaned[[1]]))/sum(colSums(counts(sce_list[[11]])))
sum(colSums(run6$cleaned[[2]]))/sum(colSums(counts(sce_list[[12]])))

cleaned <- list(run1[[1]][[1]], run1[[1]][[2]], run2[[1]][[1]], run2[[1]][[2]], run3[[1]][[1]], run3[[1]][[2]], run4[[1]][[1]], run4[[1]][[2]], run5[[1]][[1]], run5[[1]][[2]], run6[[1]][[1]], run6[[1]][[2]])

for(i in 1:length(cleaned)){
  colnames(cleaned[[i]]) <- paste0(colnames(cleaned[[i]]), '-1')
}

for(i in 1:length(sce_list)){
  counts(sce_list[[i]])[,colnames(cleaned[[i]])] <- cleaned[[i]]
}

rm(cleaned, run1, run2, run3, run4, run5, run6)

```

Swapping was extremely minimal (< 0.01). We've removed the swapped reads anyways, but it was hardly consequential.

## Filtering out empty drops

Now we filter each sample separately for empty droplets. First we'll look at the Count by rank plot for UMI counts within each droplet just to get a visual of whether anything is odd.

```{r}

brs <- lapply(sce_list, function(x) {
  temp <- counts(x)
  return(barcodeRanks(temp))})

for(i in 1:length(brs)){
  # returning only unique for plotting speed
  uniq <- !duplicated(brs[[i]]$rank)
  plot(brs[[i]]$rank[uniq], brs[[i]]$total[uniq], log="xy",
       xlab="Rank", ylab="Total UMI count", main=unique(sce_list[[i]]$Sample), cex.lab=1.2)
  abline(h=metadata(brs[[i]])$inflection, col="darkgreen", lty=2)
  abline(h=metadata(brs[[i]])$knee, col="dodgerblue", lty=2)
  legend("bottomleft", legend=c("Inflection", "Knee"), 
         col=c("darkgreen", "dodgerblue"), lty=2, cex=1.2)
}

rm(brs)

```

Sample 3 is a bit odd in shape - there is a strangely shallow slope after the knee that makes the inflection point lower. Let's see how this affects emptyDrops filtering. 

```{r}

e.out <- vector('list', length=length(sce_list))

for(i in 1:length(sce_list)){
  set.seed(100)
  e.out[[i]] <- emptyDrops(counts(sce_list[[i]]))
  print(summary(e.out[[i]]$FDR <= 0.001))
  print(table(Sig=e.out[[i]]$FDR <= 0.001, Limited=e.out[[i]]$Limited))
  
  set.seed(100)
  limit <- 100   ## set at 100 because this is the default for emptyDrops used above
  all.out <- emptyDrops(counts(sce_list[[i]]), lower=limit, test.ambient=TRUE)
  hist(all.out$PValue[all.out$Total <= limit & all.out$Total > 0],
       xlab="P-value", main=unique(sce_list[[i]]$Sample), col="grey80") 
}

```

Number of iterations looks appropriate as Limited=TRUE and Sig=FALSE is always 0.

P-values of assumed ambient reads sometimes dip near 1 but have no peaks near 0, so it is presumed that the default assumption that droplets containing fewer than 100 UMIs do not contain any real cells is accurate.

Thus we can filter the droplets accordingly.

```{r}
for(i in 1:length(sce_list)){
  sce_list[[i]] <- sce_list[[i]][,which(e.out[[i]]$FDR <= 0.001)]
}

lapply(sce_list, dim)
```

So the odd slope in Sample 3 corresponded to dramatically more called cells than the other samples. This might be an issue, and we will revisit this as we continue to process and QC.

## Removing ambient RNA

While removing ambient RNA is not always necessary, the data we are dealing with comes from very different immunological settings of control versus infected versus vaccinated. As such, there will be dramatically different cell states in each animal and these will contribute to background RNA. The worry is that it might appear that *all* cells are expressing e.g. specific inflammatory genes as opposed to the main cell type in that condition. Thus, we will remove ambient RNA from each sample. 

We start by running a quick clustering.

```{r}

library(scran)
library(scater)
library(uwot)

for(i in 1:length(sce_list)){
  set.seed(100) 
  clust <- quickCluster(sce_list[[i]])
  colLabels(sce_list[[i]]) <- clust
}

```

Next we will use this rough clustering to remove ambient RNA.

```{r}

for(i in 1:length(sce_list)){
  amb <- metadata(e.out[[i]])$ambient[,1]
  amb_full <- rep(0, nrow((sce_list[[i]])))
  names(amb_full) <- rownames(sce_list[[i]])
  amb_full[names(amb)] <- amb
  assays(sce_list[[i]], withDimnames=FALSE)$counts <-
    removeAmbience(counts(sce_list[[i]]), 
                   ambient=amb_full, groups=sce_list[[i]]$label)
}

rm(e.out, amb_full, amb)
```

## Cellranger filtering

For comparison, we will also read in the cellranger filtered output to see how results compare in terms of number of cells. 

```{r}
samps <- paste0('sample', as.character(c(1:12)))
dirs <- paste0(samps, '/filtered_feature_bc_matrix')
sce_list_cr <- vector('list', length=length(samps))

for(i in 1:length(samps)){
  sce_list_cr[[i]] <- read10xCounts(samples=dirs[i], sample.names=samps[i], col.names=TRUE)
}

lapply(sce_list, dim)

lapply(sce_list_cr, dim)

for(i in 1:length(sce_list)){
  print((dim(sce_list[[i]])[2] - dim(sce_list_cr[[i]])[2])/
          ((dim(sce_list[[i]])[2] + dim(sce_list_cr[[i]])[2])/2)*100)
}


rm(sce_list_cr)

save(list=ls(), file='ongoing.RData')
```

So this wasn't terribly different - cellranger pipelines treated the samples similarly within 5%. We'll stick with our pipeline for now, and may very well filter to a similar number as we go on with the QC. 

# QC of cells

## Removing residual low quality cells

Here we will look for outlying cells in terms of detected gene number and log-ratio of mitochondrial reads. 

In order to look at specific types of genes, we need to annotate the dataset. To make this easier, we'll first remove genes with no counts in any dataset.

```{r}
# Remove genes with zero counts
detected <- lapply(sce_list, function(x){rowSums(counts(x)) > 0})
names(detected) <- paste0('Sample', as.character(c(1:12)))
detected <- data.frame(detected)
undetected <- which(rowSums(detected) == 0)
sce_list <- lapply(sce_list, function(x) {
  return(x[-undetected,])
})
rm(detected, undetected)
# Annotate the remaining non-zero genes
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "sscrofa_gene_ensembl")

map <- getBM(attributes=c('ensembl_gene_id', 'hgnc_symbol', 'chromosome_name', 'start_position', 'end_position'),
             filters='ensembl_gene_id', 
             values=rownames(sce_list[[1]]),
             mart=ensembl)

rownames(map) <- as.character(map$ensembl_gene_id)

sce_list <- lapply(sce_list, function(x){
  rowData(x)$hgnc_symbol <- map[rownames(x), 'hgnc_symbol']
  rowData(x)$chr <- map[rownames(x), 'chromosome_name']
  rowData(x)$start <- map[rownames(x), 'start_position']
  rowData(x)$end <- map[rownames(x), 'end_position']
  return(x)
})
# Identify and plot incidence of mitochondrial genes
length(which(rowData(sce_list[[1]])$chr == 'MT'))
rowSums(counts(sce_list[[1]])[which(rowData(sce_list[[1]])$chr == 'MT'),])
hist(colSums(counts(sce_list[[1]])[which(rowData(sce_list[[1]])$chr == 'MT'),]), main='Mitochondrial UMI counts', xlab='UMI counts of mitochondrial transcripts per cell', ylab= 'Number of cells', breaks=40)

```

There are very few annotated mitochondrial genes detected here, but the UMI count for each is very high, so it might be all right. 

Now we can check various QC paramters.

```{r}
sce_list <- lapply(sce_list, function(x){
  stats <- perCellQCMetrics(x, subsets=list(Mito=which(rowData(x)$chr=="MT")))
  colData(x)$sum_counts <- stats$sum
  colData(x)$genes_detected <- stats$detected
  colData(x)$sum_mito <- stats$subsets_Mito_sum
  colData(x)$mito_ratio <- log10(stats$subsets_Mito_sum) - 
    log10(stats$sum - stats$subsets_Mito_sum)
  return(x)
}) 

lapply(sce_list, function(x){
  hist(log10(x$sum_counts), main=paste0(unique(x$Sample),', total reads per cell'), xlab='log total counts')
  hist(x$genes_detected, main=paste0(unique(x$Sample),', total genes detected per cell'), xlab='genes detected')
  plot(log10(x$sum_mito), log10(x$sum_counts - x$sum_mito), main=paste0(unique(x$Sample),', mitochondrial reads versus non-mitochondrial reads'), xlab='log mitochondrial counts', ylab='log non-mitochondrial counts')
  hist(x$mito_ratio, main=paste0(unique(x$Sample),', ratio of mitochondrial reads to non-mitochondrial reads'), xlab='log mito:non mito count ratio')
})

```

There appear to be a subset of cells in most samples that are very low in total counts but high in mitochondrial reads, suggesting poor quality cells. There are also a very large number of cells in Sample 3 with a low number of detected genes, suggesting that the funny slope of the UMI count versus rank plot represents retention of cells that are either lower quality of or simply lower sequencing depth. As everything is shifted in Sample 3 and the cell number is so high, it seems most likely that this is due to a sequencing depth issue.

Let's now visualize the data and colour the plot by various QC metrics and population markers to see if the suspect cells represent truly low quality cells or an odd cell type.

```{r}

for(i in 1:length(sce_list)){
  set.seed(100) 
  clust <- quickCluster(sce_list[[i]])
  sce_list[[i]] <- computeSumFactors(sce_list[[i]], cluster=clust)
  sce_list[[i]] <- logNormCounts(sce_list[[i]])
  sce_list[[i]] <- fixedPCA(sce_list[[i]], subset.row=NULL) ## note intentionally not using HVGs here as this is a rough visualization to see all effects on cells
  sce_list[[i]] <- runUMAP(sce_list[[i]], dimred='PCA')
}

```

Let's have a look at samples 1, 3 and 7 and check QC parameters and marker genes for sorted subsets to determine whether QC-outliers are of one specific cell type.

```{r}

x <- sce_list[[1]]
plotReducedDim(x, 'UMAP', colour_by='genes_detected') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by='mito_ratio') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'HNMT')]) + ggtitle(paste0(unique(x$Sample), ', highlighting myeloid cells')) ## myeloid
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'CD3E')]) + ggtitle(paste0(unique(x$Sample), ', highlighting T cells')) ## T cells
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'PAX5')]) + ggtitle(paste0(unique(x$Sample), ', highlighting B cells')) ## B cells

x <- sce_list[[3]]
plotReducedDim(x, 'UMAP', colour_by='genes_detected') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by='mito_ratio') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'HNMT')]) + ggtitle(paste0(unique(x$Sample), ', highlighting myeloid cells')) ## myeloid
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'CD3E')]) + ggtitle(paste0(unique(x$Sample), ', highlighting T cells')) ## T cells
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'PAX5')]) + ggtitle(paste0(unique(x$Sample), ', highlighting B cells')) ## B cells

x <- sce_list[[7]]
plotReducedDim(x, 'UMAP', colour_by='genes_detected') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by='mito_ratio') + ggtitle(unique(x$Sample))
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'HNMT')]) + ggtitle(paste0(unique(x$Sample), ', highlighting myeloid cells')) ## myeloid
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'CD3E')]) + ggtitle(paste0(unique(x$Sample), ', highlighting T cells')) ## T cells
plotReducedDim(x, 'UMAP', colour_by=rownames(sce_list[[1]])[which(rowData(sce_list[[1]])$Symbol == 'PAX5')]) + ggtitle(paste0(unique(x$Sample), ', highlighting B cells')) ## B cells

```

From this examination, it appears that none of the major cell types are specifically characterized by high mitochondrial ratios, so I think that we can safely filter these out. It would honestly be difficult to justify retaining these outliers. We can always revisit this later to look at marker genes for the clusters if desired.

Sample 3 remains tricky as there is a population of cells with myeloid markers that has a much lower rate of gene detection than the other myeloid cluster, which seems to match that found in Sample 1. This could be simply a cluster of poorly sequenced cells that clustered independently due to lower sequencing depth, but it also might represent a novel population of cells that appears in this biological condition (though it is unique among the 3 samples from the same biological condition, so this seems unlikely and definitely unprovable). We may want to consider additional sequencing of this library in the future. 

So now we filter out cells with outlying mitochondrial ratios within each dataset.

```{r, fig.width=10, fig.height=6} 

for(i in 1:length(sce_list)){
  sce_list[[i]]$outlier <- isOutlier(sce_list[[i]]$mito_ratio, type='both')
  print(table(sce_list[[i]]$outlier))
  cols <- rep('black', nrow(sce_list[[i]]))
  cols[sce_list[[i]]$outlier] <- 'red'
  plot(log10(sce_list[[i]]$sum_mito), log10(sce_list[[i]]$sum_counts - sce_list[[i]]$sum_mito), 
       col=cols, main=paste0(unique(sce_list[[i]]$Sample), ', highlighting cells with outlying mito to non-mito ratios'), xlab='log mito counts', ylab='log non-mito counts')
  plotReducedDim(sce_list[[i]], 'UMAP', colour_by='outlier')
  sce_list[[i]] <- sce_list[[i]][,!sce_list[[i]]$outlier]
}

```

## Checking for doublets

Here we annotate suspected doublets within individual samples. We will wait to remove these cells until after we have merged and clustered everything because it can be the case that a whole cluster is primarily doublets and we want to be sure we detect that.

```{r}

library(scds)
for(i in 1:length(sce_list)){
  sce_list[[i]] <- cxds(sce_list[[i]])
  sce_list[[i]]$query_doublet <- isOutlier(sce_list[[i]]$cxds_score, nmads=5, type='higher')
  plotReducedDim(sce_list[[i]], 'UMAP', colour_by='cxds_score')
}

```

These are occasionally clustered, but the scores are fairly incomparable across samples. For this reason we've marked those cells within each sample that are more than 5 MAD above the median. We will revisit whether these are doublets after merging, normalizing and clustering.


# Normalization and merging

Now we will properly normalize the data and merge all of the samples together. We will exclude Sample 3 due to its low sequencing saturation and excessive cell number as it is a very strong outlier. 

Note that we also considered using multiBatchNorm instead to account for differential sequencing depth between samples, but a quick visualization showed that it did not make particularly much difference. Hence, we don't do this.

```{r}
for(i in 1:length(sce_list)){
  colnames(sce_list[[i]]) <- paste0(colnames(sce_list[[i]]), '-', as.character(i))
}
names(sce_list) <- paste0('sample', c(1:12))
sce_all <- do.call('cbind', sce_list[c(1:2, 4:12)])

## add additional annotation to the cells for subsequent visualization

run <- as.character(sce_all$Sample)
run[run %in% c('sample1', 'sample2')] <- 'run1'
run[run %in% c('sample3', 'sample4')] <- 'run2'
run[run %in% c('sample5', 'sample6')] <- 'run3'
run[run %in% c('sample7', 'sample8')] <- 'run4'
run[run %in% c('sample9', 'sample10')] <- 'run5'
run[run %in% c('sample11', 'sample12')] <- 'run6'
sce_all$run <- factor(run)

set <- as.character(sce_all$run)
set[set %in% c('run1', 'run2', 'run3')] <- 'set1'
set[set %in% c('run4', 'run5', 'run6')] <- 'set2'
sce_all$set <- factor(set)

annot <- read.csv('sample_info.csv')
annot$sample <- paste0('sample', as.character(annot$AR_no.))
rownames(annot) <- annot$sample
sce_all$condition <- annot[as.character(sce_all$Sample), 'Condition']
sce_all$sex <- annot[as.character(sce_all$Sample), 'Sex']

```

Now we perform quick clustering and normalization by pooled factors. 

```{r}
clust <- quickCluster(sce_all)
size_factors <- calculateSumFactors(sce_all, cluster=clust)
sce_all <- logNormCounts(sce_all, size_factors=size_factors)

```

Now we look at the data to determine whether there are batch effects that ought to be removed between samples. Here we will look for the top 2500 genes for biological variability and use these for PCA and other dimensionality reduction as these are the genes we will ultimately be interested in.

Whether or not sample should be specified in the design for determining biological variability is an interesting question. In some ways, sample corresponds to a batch effect in that technical effects of e.g. sequencing depth will be captured here. In other ways, it is a biological effect as different samples underwent different experimental treatments and these would be expected to impact gene expression in quite a profound way. As a compromise, we will include sequencing run (i.e. sample pair) in the design as this should pick up batch effects but will hopefully avoid excluding biologically variable genes. Of note, using sequencing set (i.e. 6 samples per batch) left 2 samples (sample 1 and sample 6) as obvious outliers. Using Sample would obviously merge things even more. As the goal here is identifying cell types, we'll stick with the pairs as the correction that results in sufficient overlap with the lowest degree of forced merging possible.


```{r}

varmod <- modelGeneVar(sce_all, block=sce_all$run)

for(i in 1:length(levels(sce_all$run))){
  plot(metadata(varmod$per.block[[i]])$mean, metadata(varmod$per.block[[i]])$var,
       main=levels(sce_all$run)[i],
       xlab="Mean of log-expression", 
       ylab="Variance of log-expression")
  curve(metadata(varmod$per.block[[i]])$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}

vargenes <- getTopHVGs(varmod, n=5000)
length(vargenes)

#original
#vargenes <- getTopHVGs(varmod, n=2500)


sce_all <- fixedPCA(sce_all, subset.row=vargenes) 
set.seed(100)
sce_all <- runUMAP(sce_all, dimred='PCA')
set.seed(100)
sce_all <- runTSNE(sce_all, dimred='PCA')

col_list <- c(sample1='#A6CEE3', sample2='#1F78B4', 
                                  sample4='#B15928', 
                                  sample5='#B2DF8A', sample6='#33A02C',
                                  sample7='#FB9A99', sample8='#E31A1C',
                                  sample9='#FDBF6F', sample10='#FF7F00',
                                  sample11='#CAB2D6', sample12='#6A3D9A')

set.seed(100)
ord <- sample.int(ncol(sce_all), ncol(sce_all), replace=FALSE)
pl <- plotUMAP(sce_all[,ord], colour_by="Sample")
pl + scale_colour_manual(values=col_list) + ggtitle('UMAP using HVGs, uncorrected')
plotUMAP(sce_all[,ord], colour_by="condition") + ggtitle('UMAP using HVGs, uncorrected')
plotUMAP(sce_all[,ord], colour_by="sex") + ggtitle('UMAP using HVGs, uncorrected')


plotReducedDim(sce_all, 'UMAP', colour_by=rownames(sce_all)[which(rowData(sce_all)$Symbol == 'HNMT')]) + ggtitle('UMAP using HVGs, myeloid highlighted (HMNT), uncorrected') ## myeloid
plotReducedDim(sce_all, 'UMAP', colour_by=rownames(sce_all)[which(rowData(sce_all)$Symbol == 'CD3E')]) + ggtitle('UMAP using HVGs, T cells highlighted (CD3E), uncorrected') ## T cells
plotReducedDim(sce_all, 'UMAP', colour_by=rownames(sce_all)[which(rowData(sce_all)$Symbol == 'PAX5')]) + ggtitle('UMAP using HVGs, B cells highlighted (PAX5), uncorrected') ## B cells


```


Now we have to make a decision about batch correction. Will we batch correct between all samples, between run pairs, or just between run sets? As mentioned above, batch-correcting by set left 2 severe outliers. The best option appears to be to correct between run pairs.

```{r}
library(batchelor)
set.seed(100)
sce_all_corrected <- fastMNN(sce_all, batch=sce_all$run, subset.row = vargenes, get.variance = TRUE)

set.seed(100)
sce_all_corrected <- runUMAP(sce_all_corrected, dimred="corrected")
sce_all_corrected <- runTSNE(sce_all_corrected, dimred="corrected")
sce_all_corrected$Sample <- sce_all$Sample
sce_all_corrected$run <- sce_all$run
sce_all_corrected$condition <- sce_all$condition
sce_all_corrected$sex <- sce_all$sex
sce_all_corrected$cxds_score <- sce_all$cxds_score
sce_all_corrected$query_doublet <- sce_all$query_doublet
sce_all_corrected$sum_counts <- sce_all$sum_counts
sce_all_corrected$genes_detected <- sce_all$genes_detected
sce_all_corrected$mito_ratio <- sce_all$mito_ratio

rowdat <- rowData(sce_all_corrected)
rowdat <- cbind(rowdat, rowData(sce_all)[rownames(sce_all_corrected),])

rowData(sce_all_corrected) <- rowdat

sce_all_corrected$batch <- factor(sce_all_corrected$batch)

set.seed(100)
ord <- sample.int(ncol(sce_all_corrected), ncol(sce_all_corrected), replace=FALSE)
pl <- plotUMAP(sce_all_corrected[,ord], colour_by="Sample")
pl + scale_colour_manual(values=col_list) + ggtitle('UMAP using HVGs, corrected by run')

plotUMAP(sce_all_corrected[,ord], colour_by="run") + ggtitle('UMAP using HVGs, corrected by run')
plotUMAP(sce_all_corrected[,ord], colour_by="condition") + ggtitle('UMAP using HVGs, corrected by run')
plotUMAP(sce_all_corrected[,ord], colour_by="sex") + ggtitle('UMAP using HVGs, corrected by run')
plotUMAP(sce_all_corrected[,ord], colour_by="cxds_score") + ggtitle('UMAP using HVGs, corrected by run')
plotUMAP(sce_all_corrected[,ord], colour_by="query_doublet") + ggtitle('UMAP using HVGs, corrected by run')

```

Interesting to see groups of suspected doublets in subclusters of cluster 1 and cluster 6. Must decide if going to remove these or move on clustering as is. 

We will leave this analysis here and start a new script where we cluster, subcluster, filter doublets, etc.

Save everything

```{r}
sessionInfo()
save(list=ls(), file='ongoing_full.RData')
```