#'
#'
plot.scRNA.features.complexHeatmap=function(s.anno, features, name,width=300, cores){
log=mclapply(seq(dim(s.anno)[1]), plot.scRNA.features, s.anno,features, name, mc.cores = cores, width=width)
}
plot.scRNA.features=function(i, s.anno, features, name, add.annotation=NULL, width=300){
# eccite PBMC dataset:
load(s.anno[i,2])
#HSC, myeloid, B-solu, T-solu, erythroi
# STEP1 find the cluster order, and set them based on lineage:
clusters=Idents(scmat)
if(s.anno[i,3]!=""){
order.clu=as.numeric(unlist(strsplit(s.anno[i,3], ",|, | ,"))) # MODIFY
ord.names=names(clusters[order(match(clusters,order.clu))])
}else{
order.clu=levels(clusters)
# order samples within cluster by umap component with more spread
# this way there is still information within cluster, for heatmaps
coords=Embeddings(scmat[["umap"]])
fit.cluster=levels(Idents(scmat))
ord.names=unlist(lapply(fit.cluster, function(clu){
d=coords[Idents(scmat)%in%clu,]
ind=which.max(c(max(d[,1])-min(d[,1]), max(d[,2])-min(d[,2]))) # take the axis with most spread
names(sort(d[,ind]))
}))
}
features=data.frame(features, stringsAsFactors = F)
if(!grepl("^B:|^N:", features[1,1])){
features[,1]=paste0("N:GEXP:", features[,1])
allfeats=rownames(fm)[!rowSums(fm>0)==0]
}else{
allfeats=gsub("N:GEXP:|N:RPPA:", "", rownames(fm)[!rowSums(fm>0)==0])
}
# text annot created if more columns than 1.
if(dim(features)[2]>1){
features=features[features[,1]%in%allfeats,]
feats=features[,1]
t.annot=apply(features[,2:dim(features)[2], drop=F], 1, paste, collapse=" ")
names(t.annot)=feats
}else{
t.annot=NULL
feats=features[features[,1]%in%allfeats,1]
}
# make a data frame for annotations on top of the heatmap
annotdf=data.frame("cluster"=clusters,"HLAIScore"=as.numeric(fm["N:SAMP:HLAIScore",]), "HLAIIScore"=as.numeric(fm["N:SAMP:HLAIIScore",]), "CytolyticScore"=as.numeric(fm["N:SAMP:CytolyticScore",]))
if(!is.null(add.annotation))annotdf$annotation=add.annotation
name=gsub("__", "_", gsub("[[:punct:]]", "_", paste(name, s.anno[i,1], sep="_")))
# plot using a complexheatmap package using a featurematrix format and custom plotting function
r1=plot.complexHM.fm(feats = feats, text_annot = t.annot, fm.f = fm, annotdf = annotdf, feats.barplot = c("HLAIScore", "HLAIIScore", "CytolyticScore"), order_columns=ord.names, order_rows=F, split.columns = T, plotting.param="/research/groups/sysgen/PROJECTS/students/HEMAP_IMMUNOLOGY/minna_work/scRNA_genelist_plotting/plotting_param_fm.txt", NAME=name, WIDTH = width)
}
plot.gene.seurat=function(s.anno, feature, name,cols=c("grey75", "red"), max=1, min=0, cores=1){
plot.feature=function(i, feature){
# plot proportions as a barplot next to 2D points
load(s.anno[i,2])
if(!feature%in%rownames(scmat))return(NULL)
p <- FeaturePlot(scmat, features = feature, cols = cols, max.cutoff = max, min.cutoff = min, order = T, combine = FALSE)
p=lapply(p, `+`, labs(title = s.anno[i,1]))
if(i>1){
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend() + NoAxes()
}
}else{
p[[1]] <- p[[1]] + NoAxes()
}
cowplot::plot_grid(plotlist = p)
}
# plot here lineage proportion per patient
plot.list=mclapply(seq(dim(s.anno)[1]), plot.feature, feature, mc.cores=cores)
plot.list=plot.list[!sapply(plot.list, is.null)]
if(!length(plot.list))return(NULL)
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170, height = ceiling(length(plot.list)/2)*75+75, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
plot.fm.feature=function(s.anno, feature, name, max=3, min=0, cores=1){
plot.feature=function(i, feature){
# plot proportions as a barplot next to 2D points
load(s.anno[i,2])
scmat[[feature]]=fm[rownames(fm)%in%feature,]
p <- FeaturePlot(scmat, features = feature, cols = c("grey75", "red"), max.cutoff = max, min.cutoff = min, order = T, combine = FALSE)
p=lapply(p, `+`, labs(title = s.anno[i,1]))
if(i>1){
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend() + NoAxes()
}
}else{
p[[1]] <- p[[1]] + NoAxes()
}
return(p)
}
# plot here lineage proportion per patient
plot.list=unlist(mclapply(seq(dim(s.anno)[1]), plot.feature, feature, mc.cores=cores), recursive = F)
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170, height = ceiling(length(plot.list)/2)*75+75, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170, height = ceiling(length(plot.list)/2)*75+75, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
plot.several.fm.feature=function(s.anno.object, features, name, min=0, max=5, cores=5){
load(s.anno.object)
features=features[features%in%rownames(fm)]
plot.feature=function(i, features){
# plot proportions as a barplot next to 2D points
scmat[[features[i]]]=fm[rownames(fm)%in%features[i],]
p <- FeaturePlot(scmat, features = features[i], cols = c("grey75", "red"), max.cutoff = max, min.cutoff = min, order = T, combine = FALSE)
p=lapply(p, `+`, labs(title = features[i]))
if(i>1){
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend() + NoAxes()
}
}else{
p[[1]] <- p[[1]] + NoAxes()
}
return(p)
}
# plot here lineage proportion per patient
plot.list=unlist(mclapply(seq(features), plot.feature, features, mc.cores=cores), recursive = F)
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170, height = ceiling(length(plot.list)/2)*75+75, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
plot.seurat.ident=function(s.anno, colorv, name, cores=1){
plot.feature=function(i){
# plot proportions as a barplot next to 2D points
load(s.anno[i,2])
cells=gsub(" \\d+$", "", levels(Idents(scmat)))
p=DimPlot(scmat, group.by = c("ident"), cols = colorv[match(cells, colorv[,1]),2], ncol = 1, label = T, repel = T) + NoLegend() + NoAxes()
}
# plot here lineage proportion per patient
plot.list=mclapply(seq(dim(s.anno)[1]), plot.feature, mc.cores=cores)
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170, height = ceiling(length(plot.list)/2)*75+75, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
#'
plot.multi.scatter.scRNA=function(s.anno, scmat=NULL, colors.group, identityVector.samples.list=NULL, seurat.feature=NULL, name="scRNA", cores=7, SIZE=0.1, rasterize=T, width=77, height = 74, text.size=10){
# tools:
library(Matrix)
library(Seurat)
source("/research/users/ppolonen/git_home/common_scripts/visualisation/plotting_functions.R")
library(data.table)
library(ComplexHeatmap)
library(circlize)
library(parallel)
library(ggplot2)
# plot here lineage proportion per patient
plot.list=mclapply(seq(dim(s.anno)[1]), function(i){
# plot proportions as a barplot next to 2D points
load(s.anno[i,2])
coords=Embeddings(scmat[["umap"]])
identityVector=unlist(strsplit(s.anno$Type[i], ", |,"))
clusters=unlist(strsplit(s.anno$Clusters[i], ", |,"))
clusters.samples=as.character(Idents(scmat))
if(is.null(identityVector.samples.list)){
identityVector.samples=clusters.samples
for(j in seq(clusters)){
identityVector.samples[clusters.samples%in%clusters[j]]=identityVector[j]
}
}else{
identityVector.samples=as.character(identityVector.samples.list[[i]]) # must be in a list
}
if(!is.null(seurat.feature)){
identityVector.samples=as.character(scmat[[seurat.feature]][,1]) # must be in a list
}
name.sample=s.anno$Sample[i]
# take the colors:
colorv=colors.group[match(unique(identityVector.samples), colors.group[,1]),2]
a=plot.scatter(x = coords[,1], y = coords[,2], group = identityVector.samples, colorv = colorv, namev=name, main = name, add.proportions = T, SIZE = SIZE, add.legend=F, rasterize=rasterize, width=width, height = height, text.size=text.size)
return(a)
}, mc.cores=cores)
# add legend last
add.legend=get.only.legend(group = colors.group[,1], colorv = colors.group[,2])
nrows=ceiling(length(plot.list)/2)+1
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 2*width+width*0.2, height = ceiling(length(plot.list)/2)*height+120, rows = nrows, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
rowi=1
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
figure <- multipanelfigure::fill_panel(figure,row = nrows,column = 1:2, add.legend)
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
plot.scatter.seurat=function(scmat, colors.group,clusters="", identityVector="", identityVector.samples.list=NULL, seurat.feature=NULL, lv=NULL, name="scRNA", cores=7, SIZE=0.1, rasterize=T, width=77, height=74, text.size=10, add.proportions=T, add.density=F, add.labels=F){
# tools:
library(Matrix)
library(Seurat)
source("/research/users/ppolonen/git_home/common_scripts/visualisation/plotting_functions.R")
library(data.table)
library(ComplexHeatmap)
library(circlize)
library(parallel)
library(ggplot2)
coords=Embeddings(scmat[["umap"]])
clusters.samples=as.character(Idents(scmat))
if(is.null(identityVector.samples.list)){
identityVector.samples=clusters.samples
for(j in seq(clusters)){
identityVector.samples[clusters.samples%in%clusters[j]]=identityVector[j]
}
}else{
identityVector.samples=as.character(identityVector.samples.list[[i]]) # must be in a list
}
if(!is.null(seurat.feature)){
identityVector.samples=as.character(scmat[[seurat.feature]][,1]) # must be in a list
ind=order(match(identityVector.samples, colors.group[,1]))
# order everything:
identityVector.samples=identityVector.samples[ind]
coords=coords[ind,]
lv=lv[ind]
}
# take the colors:
colorv=colors.group[match(unique(identityVector.samples), colors.group[,1]),2]
plot.list=list(plot.scatter(x = coords[,1], y = coords[,2], group = identityVector.samples, colorv = colorv, lv = lv, namev=name, main = name, add.proportions = add.proportions, add.density = add.density, add.labels = add.labels, SIZE = SIZE, add.legend=T, rasterize=rasterize, width=width, height = height, text.size=text.size))
figure <-multipanelfigure:: multi_panel_figure(width = width+50, height = height, rows = 1, columns = 1, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))
for(i in seq(plot.list)){
figure <- multipanelfigure::fill_panel(figure, plot.list[[i]])
}
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_fm_scatterplot.pdf"), limitsize = FALSE)
}
#'
#'
compute.FDR.scRNA=function(i, s.anno, nperm=1000, geneset, nCores=6){
library(AUCell)
library(Matrix)
library(Seurat)
load(s.anno[i,2])
# take only genes that are in the matrix
geneset=geneset[geneset%in%rownames(scmat)]
# Random
geneset.random=lapply(seq(nperm), function(i, genes){
sample(rownames(scmat), length(geneset))
})
names(geneset.random)=paste0("Random",seq(nperm))
geneset.observed=list("MDS_signature"=geneset)
data_n <- data.matrix(GetAssayData(scmat, assay = "RNA"))
cells_rankings <- AUCell_buildRankings(data_n, nCores=nCores, plotStats=TRUE)
AUC.random <- AUCell_calcAUC(geneset.random, cells_rankings, nCores = nCores, nrow(cells_rankings)*0.05)
AUC.observed <- AUCell_calcAUC(geneset.observed, cells_rankings, nCores = nCores, nrow(cells_rankings)*0.05)
# Alternative: assign cells according to the 'automatic' threshold
cells_assignment <- AUCell_exploreThresholds(AUC.observed, plotHist=T, assignCells=TRUE)
# how many times random higher than observed, divide by the number of permutations == eFDR
y=as.numeric(getAUC(AUC.observed))
x=getAUC(AUC.random)
FDR=sapply(seq(y), function(i){
sum(x[,i]>y[i])/nperm
})
empirical.FDR=data.frame(t(FDR), check.names = F)
rownames(empirical.FDR)=paste0(names(geneset.observed), "_eFDR_", nperm, "_", length(geneset))
colnames(empirical.FDR)=colnames(scmat)
a=scmat[["SingleR.label"]]
table(a[rownames(a)%in%cells_assignment$MDS_signature$assignment,])
table(a[empirical.FDR<0.05,])
scmat[["empirical.FDR"]]=as.numeric(empirical.FDR)
DimPlot(scmat, group.by = "empirical.FDR")
scmat[["y"]]=as.numeric(y)>0.04
DimPlot(scmat, group.by = "y")
return(list(empirical.FDR, y))
}
get.only.legend=function(colorv, group, position="right", direction="vertical"){
library(ggplot2)
dat=data.frame(group, colorv)
levels(dat$group)=group
myColors <- colorv
names(myColors) <- levels(group)
colScale <- scale_colour_manual(name = "group",values = myColors)
dat$x=dim(dat)[1]
dat$y=dim(dat)[1]
legend <- cowplot::get_legend(ggplot(dat, aes(x, y, colour = group)) + theme_bw() +
geom_point(size=0.6) + colScale +theme(legend.position = position, legend.direction = direction, legend.title = element_blank(),
legend.key=element_blank(), legend.box.background=element_blank(),
legend.text=element_text(size = 10, face = "bold", family="Helvetica"))+
theme(plot.margin=unit(c(1,1,1,1), "mm"))+
guides(colour = guide_legend(override.aes = list(size=3))))
return(legend)
}
statistical.comparisons.scRNA=function(i, s.anno, genelist.statistical.analysis=NULL, test.use="wilcox"){
# plot proportions as a barplot next to 2D points
load(s.anno[i,2])
coords=Embeddings(scmat[["umap"]])
identityVector=gsub(" ", "", unlist(strsplit(s.anno$Type[i], ", |,")))
clusters=gsub(" ", "", unlist(strsplit(s.anno$Clusters[i], ", |,")))
clusters.samples=as.character(Idents(scmat))
identityVector.samples=clusters.samples
for(j in seq(clusters)){
identityVector.samples[clusters.samples%in%clusters[j]]=identityVector[j]
}
name.sample=s.anno$Sample[i]
if(is.null(genelist.statistical.analysis))genelist.statistical.analysis=rownames(scmat)
# Find the DE genes per factor and within factor:
statistical.analysis=do.call(rbind, mclapply(unique(identityVector.samples), function(id){
clusters2test=unique(clusters.samples[identityVector.samples%in%id])
# clusters2test vs Rest:
markers=FindMarkers(scmat, ident.1 = clusters2test, group.by = "seurat_clusters", test.use = test.use, logfc.threshold =0.01, only.pos = T, features = genelist.statistical.analysis)
#DE in clusters within category
cluster.markers=FindAllMarkers(scmat[,clusters.samples%in%clusters2test], test.use = test.use, logfc.threshold = 0.01, only.pos = T, features = genelist.statistical.analysis)
# make genelists with annotation:
markers$test=paste0(id, " vs. Rest")
markers$gene=rownames(markers)
cluster.markers$test=paste0(id, " ", cluster.markers$cluster, " vs. Rest")
cluster.markers=cluster.markers[,!colnames(cluster.markers)%in%"cluster"]
df.statistics=plyr::rbind.fill(list(markers, cluster.markers))
df.statistics$test.used=test.use
return(df.statistics)
}))
return(statistical.analysis)
}
#' @param seurat.object Seurat class object or data matrix that can be converted to Seurat object. Can also be a list of data matrices (RNA+citeseq etc.), Must be named by assay type with four uppercase letters (PROT, etc) (also found from seurat object with this name).
#' @param regress.cell.label Character vector with as many columns as seurat.object Contains batch information for each cell. Uses fastMNN/SCTtransform to batch correct the data (usually different samples or experiments).
#' @param check.pcs Plot jackstraw and elbowplot to optimize the number of PCA componens
#' @param plot.umap Umap plot for clusters and for sample ID, if batch corrected also batch clustering
#' @param resolution Resolution parameter for clustering
#' @param nFeature.min Filter cells with < nFeature.min
#' @param nFeature.max Filter cells with > nFeature.max
#' @param percent.mitoDNA Filter cells with mitochondrial genes > percent.mitoDNA
#' @param singleR Annotate each cell using singleR built in annotations. Mainly uses Encode+blueprint annotations.
#' @param cores Number of cores to use for Seurat/singleR/MNNcorrect
#' @param add.gm.score Adds geneset score (using geometric mean) to fm, must be a named list of genes. Geneset name is added to fm and seurat object.
#' @param ... Pass parameters to Seurat tools
sc.data.analysis=function(scmat, name="scRNA_object", nr.pcs=30, regress.cell.label=NULL, batch.correction.method=c("SCTtransform","MNNcorrect"), auto.order=F, normalize.input=T, check.pcs=T, plot.umap=T, resolution=0.5, nFeature.min=0, nFeature.max=10000, percent.mitoDNA=25, singleR=T, singleR.reference=c("blueprint_encode", "HPCA"), cores=6, add.gm.score=NULL, ...){
# determine computational resources:
future::plan("multiprocess", workers = cores)
options(future.globals.maxSize= 5e+09)
# set method for batch correction
batch.correction.method=match.arg(batch.correction.method)
# list of data matrix, can be citeseq data:
if(class(scmat)=="list"){
cat("Input is list, splitting to data matrix by type and adding to Seurat", sep="\n\n")
list.additional=lapply(2:length(scmat), function(i)CreateAssayObject(scmat[[i]]))
names(list.additional)=names(scmat)[2:length(scmat)]
scmat <- CreateSeuratObject(scmat[[1]])
}else{
list.additional=NULL
# if class is not seurat object
if(!class(scmat)=="Seurat"){
cat("Input is data matrix, converting to Seurat object", sep="\n\n")
scmat <- CreateSeuratObject(scmat)
}else{
cat("Input is Seurat object", sep="\n\n")
}
}
DefaultAssay(object = scmat) <- "RNA"
# is it a vector of pcs or single value?
if(length(nr.pcs)==1)nr.pcs=seq(nr.pcs)
# QC
scmat=PercentageFeatureSet(scmat, pattern = "^MT-", col.name = "percent.mt")
r1=VlnPlot(scmat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggplot2::ggsave(r1, filename = paste0(name, "_QC.pdf"), width = 12)
# Filtering
nFeature_RNA <- FetchData(object = scmat, vars = "nFeature_RNA")
pct.mt <- FetchData(object = scmat, vars = "percent.mt")
filt=which(x = nFeature_RNA > nFeature.min & nFeature_RNA < nFeature.max & pct.mt < percent.mitoDNA)
cat(paste("Dimension before/after filtering:", paste(paste(dim(scmat), collapse="x"), paste(dim(scmat[, filt]), collapse="x"), collapse = "/")), sep="\n\n")
scmat=scmat[, filt]
if(normalize.input){
assay="SCT"
scmat=SCTransform(scmat, variable.features.n = 3000) #,...)
# seurat standard processing
scmat=RunPCA(scmat, npcs = 100)
reduction="pca"
}else{
assay="RNA"
# assume that the data is normalized
scmat <- FindVariableFeatures(object = scmat)
scmat <- ScaleData(object = scmat)
scmat=RunPCA(scmat, npcs = 100)
reduction="pca"
}
# set default assay to use:
DefaultAssay(object = scmat) <- assay
# use SCTransform to remove batch effect, used in umap and clustering
# else use SCT to normalize data
if(!is.null(regress.cell.label)&batch.correction.method=="SCTtransform"){
batch.df=data.frame("batch"=regress.cell.label[filt])
rownames(batch.df)=colnames(scmat)
scmat[["batch"]] = batch.df
scmat=SCTransform(scmat,vars.to.regress="batch", variable.features.n = 3000)
# seurat standard processing
scmat=RunPCA(scmat, npcs = 100)
name=paste0(name, "_", batch.correction.method)
reduction="pca"
}
# optimize pcs number:
if(check.pcs){
# setting here original RNA to define PCAs, SCT not working yet, wait for a fix
# also the issue here how to deal with the mnnCorrect?
# DefaultAssay(object = scmat) <- "RNA"
# scmat <- NormalizeData(scmat,display.progress = FALSE)
# scmat <- FindVariableFeatures(scmat,do.plot = F,display.progress = FALSE)
# scmat <- ScaleData(scmat)
# scmat=RunPCA(scmat, npcs = 100)
# determine how many pcs should be used:
r3=ElbowPlot(scmat, ndims = 100, reduction = "pca")
ggplot2::ggsave(r3, filename = paste0(name, "_elbow.pdf"))
# find optimal number of PCs to use:
scmat <- JackStraw(scmat, num.replicate = 100, dims = 100, reduction = "pca")
scmat <- ScoreJackStraw(scmat, dims = 1:100, reduction = "pca")
r4=JackStrawPlot(scmat, dims = 1:100)
ggplot2::ggsave(r4, filename = paste0(name, "_Jackstraw.pdf"), width = 12)
cat("PCs checking done", sep="\n\n")
}
# use MNN to remove batch effect, used in umap and clustering:
if(!is.null(regress.cell.label)&batch.correction.method=="MNNcorrect"){
library(BiocParallel)
library(batchelor)
cat("FastMNN batch correction...", sep="\n\n")
if(is.factor(regress.cell.label)&!auto.order){
order=levels(regress.cell.label)
cat("FastMNN batch correction order:", sep="\n\n")
cat(order, sep="-->")
cat("", sep="\n\n")
}
if(!is.factor(regress.cell.label)&!auto.order){
order=unique(regress.cell.label)
cat("FastMNN batch correction order:", sep="\n\n")
cat(order, sep="-->")
cat("", sep="\n\n")
}
batch.df=data.frame("batch"=regress.cell.label[filt])
rownames(batch.df)=colnames(scmat)
scmat[["batch"]] = batch.df
# http://htmlpreview.github.io/?https://github.com/satijalab/seurat-wrappers/blob/master/docs/fast_mnn.html
scmat <- NormalizeData(scmat)
scmat <- FindVariableFeatures(scmat)
object.list = SplitObject(scmat, split.by = "batch")[order(order)]
scmat <- SeuratWrappers::RunFastMNN(object.list, reduction.name = "mnnCorrect", auto.order=auto.order, BPPARAM=MulticoreParam(cores))
# set defaults
reduction="mnnCorrect"
cat("FastMNN batch correction... done", sep="\n\n")
}
# Basic seurat processing with umap:
scmat=RunUMAP(scmat, dims = nr.pcs, reduction = reduction)
scmat=FindNeighbors(scmat, dims = nr.pcs, reduction = reduction)
scmat=FindClusters(scmat, resolution = resolution)
# add celltype automatically
if(singleR){
cat("Annotating with singleR", sep="\n\n")
if(!file.exists(paste0(name, "_singler_object.Rdata"))){
library("SingleR")
counts <- GetAssayData(scmat, assay = assay, slot = "counts")
if(singleR.reference=="blueprint_encode"){
if(dim(counts)[2]<100000){
singler = CreateSinglerObject(counts, annot = NULL, clusters=Idents(scmat), ref.list = list(blueprint_encode), project.name=name, fine.tune = T,do.signatures = T, numCores=cores)
}else{
cat("SingleR in batches", sep="\n\n")
singler = CreateBigSingleRObject.custom(counts, annot = NULL, clusters=Idents(scmat), N=30000, ref.list = list(blueprint_encode), xy = Embeddings(scmat[["umap"]]), project.name=name, fine.tune = T, do.signatures = T, numCores=cores)
}
}
if(singleR.reference=="HPCA"){
if(dim(counts)[2]<100000){
singler = CreateSinglerObject(counts, annot = NULL, clusters=Idents(scmat), ref.list = list(hpca), project.name=name, fine.tune = T,do.signatures = T, numCores=cores)
}else{
cat("SingleR in batches", sep="\n\n")
singler = CreateBigSingleRObject.custom(counts, annot = NULL, clusters=Idents(scmat), N=30000, ref.list = list(hpca), project.name=name, xy = Embeddings(scmat[["umap"]]), fine.tune = T,do.signatures = T, numCores=cores)
}
}
# add original identifiers
singler$meta.data$orig.ident = scmat@meta.data$orig.ident # the original identities, if not supplied in 'annot'
## if using Seurat v3.0 and over use:
singler$meta.data$xy = scmat@reductions$umap@cell.embeddings # the tSNE coordinates
singler$meta.data$clusters = scmat@active.ident # the Seurat clusters (if 'clusters' not provided)
# add annot from hpca
# singler2 = CreateSinglerObject(counts, annot = NULL, clusters=Idents(scmat), ref.list = list(hpca), project.name=name, fine.tune = T,do.signatures = T, numCores=cores)
# singler$singler[[1]]$SingleR.single$labels[singler2$singler[[1]]$SingleR.single$labels%in%"Pre-B_cell_CD34-"]="pre-B-cells"
# singler$singler[[1]]$SingleR.single$labels[singler2$singler[[1]]$SingleR.single$labels%in%"Pro-B_cell_CD34+"]="pro-B-cells"
# singler$singler[[1]]$SingleR.clusters$labels[singler2$singler[[1]]$SingleR.clusters$labels%in%"Pre-B_cell_CD34-"]="pre-B-cells"
# singler$singler[[1]]$SingleR.clusters$labels[singler2$singler[[1]]$SingleR.clusters$labels%in%"Pro-B_cell_CD34+"]="pro-B-cells"
# these names are not intuitive, replace them, mentioned that these are actually naive http://comphealth.ucsf.edu/SingleR/SupplementaryInformation2.html:
singler$singler[[1]]$SingleR.single$labels[singler$singler[[1]]$SingleR.single$labels%in%"CD4+ T-cells"]="CD4+ Tn"
singler$singler[[1]]$SingleR.single$labels[singler$singler[[1]]$SingleR.single$labels%in%"CD8+ T-cells"]="CD8+ Tn"
singler$singler[[1]]$SingleR.clusters$labels[singler$singler[[1]]$SingleR.clusters$labels%in%"CD4+ T-cells"]="CD4+ Tn"
singler$singler[[1]]$SingleR.clusters$labels[singler$singler[[1]]$SingleR.clusters$labels%in%"CD8+ T-cells"]="CD8+ Tn"
save(singler, file=paste0(name, "_singler_object.Rdata"))
}else{
library("SingleR")
load(paste0(name, "_singler_object.Rdata"))
cat(paste0("Loaded existing SingleR object (", paste0(name, "_singler_object.Rdata"), "), remove/rename it if you want to re-compute."), sep="\n\n")
}
# can be added as cluster id
cluster=singler$singler[[1]]$SingleR.clusters$labels
cluster.main=singler$singler[[1]]$SingleR.clusters.main$labels
# order based on cell type and add celltype to cluster:
lineage=c("HSC","MPP","CLP","CMP","GMP","MEP","Monocytes","DC","Macrophages","Macrophages M1","Macrophages M2","CD4+ Tn","CD4+ T-cells", "CD4+ Tcm", "CD4+ Tem","CD8+ Tn","CD8+ T-cells","CD8+ Tcm","CD8+ Tem","NK cells","Tregs","naive B-cells","Memory B-cells","Class-switched memory B-cells","Plasma cells","Endothelial cells","Neutrophils","Eosinophils","Fibroblasts","Smooth muscle","Erythrocytes","Megakaryocytes")
lineage=unique(c(lineage,blueprint_encode$types, hpca$types))
lineage=lineage[lineage%in%singler$singler[[1]]$SingleR.single$labels]
scmat[["SingleR.label.main"]]=singler$singler[[1]]$SingleR.single.main$labels
scmat[["SingleR.label.cluster"]]=paste(singler$singler[[1]]$SingleR.single$labels, Idents(scmat))
scmat[["SingleR.label"]]=factor(singler$singler[[1]]$SingleR.single$labels, levels=lineage[lineage%in%unique(singler$singler[[1]]$SingleR.single$labels)])
identityVector.samples=as.character(Idents(scmat))
clusters.samples=Idents(scmat)
for(j in seq(cluster)){
identityVector.samples[clusters.samples%in%rownames(cluster)[j]]=paste0(cluster[j], ":", rownames(cluster)[j])
}
# order and cluster identity;
cluster=cluster[order(match(cluster[,1],lineage)),,drop=F]
cat(paste(levels(Idents(scmat)), collapse=","), paste(cluster, collapse=","), sep="\t\t")
scmat[["SingleR.cluster"]]=gsub(":.*.", "", Idents(scmat)) # just the "cluster label" per cell
# order seurat clusters too:
Idents(scmat)=factor(identityVector.samples, levels=paste0(cluster, ":", rownames(cluster)))
r4=DimPlot(object = scmat, group.by = "SingleR.label", reduction = 'umap' , label = TRUE, pt.size = 0.5, repel = T) + NoLegend() + NoAxes()
ggplot2::ggsave(r4, filename = paste0(name, "_UMAP_singleR.pdf"), width = 6, height = 6)
cat("\nSingleR done", sep="\n\n")
# new singleR:
# counts <- GetAssayData(scmat, assay = "RNA", slot = "counts")
#
# library(SingleR)
# library(scater)
#
# ref=NovershternHematopoieticData()
#
# common <- intersect(rownames(counts), rownames(ref))
# ref <- ref[common,]
# counts <- counts[common,]
# counts.log <- logNormCounts(counts)
#
# singler <- SingleR(test = counts.log, ref = ref, labels = ref$label.main, numCores=cores)
#
# # order based on cell type and add celltype to cluster:
# lineage=c("HSC","MPP","CLP","CMP","GMP","MEP","Monocytes","DC","Macrophages","Macrophages M1","Macrophages M2","CD4+ Tn","CD4+ T-cells", "CD4+ Tcm", "CD4+ Tem","CD8+ Tn","CD8+ T-cells","CD8+ Tcm","CD8+ Tem","NK cells","Tregs","naive B-cells","Memory B-cells","Class-switched memory B-cells","Plasma cells","Endothelial cells","Neutrophils","Eosinophils","Fibroblasts","Smooth muscle","Erythrocytes","Megakaryocytes")
# lineage=unique(c(lineage,ref$label.main))
#
# scmat[["SingleR.label.main"]]=singler$first.labels
# scmat[["SingleR.label"]]=factor(singler$labels, levels=lineage[lineage%in%unique(singler$labels)])
#
# r4=DimPlot(object = scmat, group.by = "SingleR.label", reduction = 'umap' , label = TRUE, pt.size = 0.5, repel = T) + NoLegend() + NoAxes()
# ggplot2::ggsave(r4, filename = paste0(name, "_UMAP_singleR.pdf"), width = 6, height = 6)
#
# cat("\nSingleR done", sep="\n\n")
}
# plot clusters and samples
if(plot.umap){
if(!is.null(regress.cell.label)){
r6=DimPlot(scmat, group.by = c("batch"), ncol = 2, label = T, repel = T) + NoLegend() + NoAxes()
ggplot2::ggsave(r6, filename = paste0(name, "_UMAP_batch.pdf"), width = 6, height = 6)
r6=DimPlot(scmat, group.by = c("ident"), ncol = 2, label = T, repel = T) + NoLegend() + NoAxes()
ggplot2::ggsave(r6, filename = paste0(name, "_UMAP_clusters.pdf"), width = 6, height = 6)
}else{
r6=DimPlot(scmat, group.by = c("ident"), ncol = 1, label = T, repel = T) + NoLegend() + NoAxes()
ggplot2::ggsave(r6, filename = paste0(name, "_UMAP_clusters.pdf"), width = 6, height = 6)
}
}
# add immunoscores to FM format data matrix
fm.f <- GetAssayData(scmat)
# add gm based score:
add.scores=list(HLAIScore=c("B2M", "HLA-A", "HLA-B","HLA-C"), HLAIIScore=c("HLA-DMA","HLA-DMB","HLA-DPA1","HLA-DPB1","HLA-DRA","HLA-DRB1"), CytolyticScore=c("GZMA", "PRF1", "GNLY", "GZMH", "GZMM"))
if(!is.null(add.gm.score))add.scores=append(add.scores, add.gm.score)
gm.objects=do.call(rbind, lapply(seq(add.scores), function(i){
dat3=fm.f[rownames(fm.f)%in%add.scores[[i]],]
gm=t(apply(dat3, 2, gm_mean)) # done to normalized values
rownames(gm)=names(add.scores)[i]
return(gm)
}))
# also add to seurat object:
for(i in seq(add.scores)){
scmat[[names(add.scores)[i]]] <- gm.objects[i,]
}
# gene expression and scores to fm:
rownames(fm.f)=paste0("N:GEXP:", rownames(fm.f))
rownames(gm.objects)=paste0("N:SAMP:", rownames(gm.objects))
fm=rbind(gm.objects, fm.f)
# DE analysis:
markers.all=FindAllMarkers(scmat, only.pos = T)#, ...)
# go through each data matrix, add to seurat and fm if more matrices
if(!is.null(list.additional)){
for(i in seq(list.additional)){
mat.add=list.additional[[i]]
mat.add=subset(mat.add,cells = colnames(scmat))
scmat[[names(list.additional)[i]]] <- mat.add
scmat <- NormalizeData(scmat, assay = names(list.additional)[i], normalization.method = "CLR") # test scttransform too
scmat <- ScaleData(scmat, assay = names(list.additional)[i])
add.data.normalized <- data.matrix(GetAssayData(scmat, assay = names(list.additional)[i], slot = "data"))
rownames(add.data.normalized)=paste0("N:",names(list.additional)[i], ":", rownames(add.data.normalized))
fm=rbind(fm, add.data.normalized)
}
}
save(list=c("fm", "scmat", "markers.all"), file=paste0(name, "_scRNA.Rdata"))
return(paste0(name, "_scRNA.Rdata"))
}
# make sure data is on a linear non-log scale
# check 0 values, if a lot between 0-1 better to use add.one. If counts, remove works too quite well
# zero fix methods: https://arxiv.org/pdf/1806.06403.pdf
gm_mean=function(x, na.rm=TRUE, zero.handle=c("remove", "add.one", "zero.propagate")){
zero.handle=match.arg(zero.handle)
if(any(x < 0, na.rm = TRUE)){
warning("Negative values produced NaN - is the data on linear - non-log scale?")
return(NaN)
}
if(zero.handle=="remove"){
return(exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x)))
}
if(zero.handle=="add.one"){
return(exp(mean(log(x+1), na.rm = na.rm))-1)
}
if(zero.handle=="zero.propagate"){
if(any(x == 0, na.rm = TRUE)){
return(0)
}
return(exp(mean(log(x), na.rm = na.rm)))
}
}
get.cluster.label.singleR=function(singler){
# how to add this?
cluster=singler$singler[[1]]$SingleR.clusters$labels
cluster.main=singler$singler[[1]]$SingleR.clusters.main$labels
# order based on cell type:
# cat(paste0("'", colors.group[,1], "'"), sep=",")
lineage=c('HSC','MPP','CLP','CMP','GMP','MEP','Monocytes','DC','Macrophages','Macrophages M1','Macrophages M2','CD4+ T-cells','CD8+ T-cells','Tregs','T_cell:gamma-delta','T_cell:CD4+_Naive','T_cell:CD8+_naive','T_cell:Treg:Naive','CD4+ Tcm','CD4+ Tem','CD8+ Tcm','CD8+ Tem','NK cells','naive B-cells','Memory B-cells','B-cells','Class-switched memory B-cells','Plasma cells','Endothelial cells','Neutrophils','Eosinophils','Fibroblasts','Smooth muscle','Erythroblast','Erythrocytes','Megakaryocytes')
cluster=cluster[order(match(cluster[,1],lineage)),,drop=F]
# order and cluster identity;
cat(paste(rownames(cluster), collapse=","), paste(cluster, collapse=","), sep="\t")
cat("", sep="\n")
return(cluster)
}
# run cellphoneDB
run.cellphonedb=function(scmat, celltype, out=getwd(), cores=10, threshold=0.05){
if(!dir.exists(out))dir.create(out, recursive = T)
if(is.null(celltype)){
celltype=gsub(" ", "_", as.character(Idents(scmat)))
}
setwd(out)
# input1:
df=data.frame("Cell"=gsub("\\.", "_", colnames(scmat)), "cell_type"=celltype, stringsAsFactors = F)
df$cell_type[scmat[["SingleR.label"]]=="Tregs"]="Tregs"
meta=file.path(out, "meta.tsv")
write.table(df, meta, row.names = F, col.names = T, quote = F, sep="\t")
# input2:
counts <- data.matrix(GetAssayData(scmat, assay = "SCT", slot = "data"))
countsfile=tempfile()
# gene symbol was not working...
library("org.Hs.eg.db") # remember to install it if you don't have it already
ensambleid <- mapIds(org.Hs.eg.db, keys = rownames(counts), keytype = "SYMBOL", column="ENSEMBL")
rownames(counts)=ensambleid
counts=counts[!is.na(rownames(counts)),]
countsfile=file.path(out, "counts.tsv")
write.table(t(c("Gene", gsub("\\.", "_", colnames(counts)))), countsfile, row.names = F, col.names = F, quote = F, sep="\t")
write.table(counts, countsfile, row.names = T, col.names = F, quote = F, sep="\t", append=T)
system(paste0("cellphonedb method statistical_analysis ", meta, " ", countsfile, " --iterations=1000 --threads=", cores, " --threshold=", threshold))
# make heatmap of interactions:
system(paste0("cellphonedb plot heatmap_plot ", meta))
deconvoluted=data.table::fread(file.path(out, "out/deconvoluted.txt"), data.table = F)
pvalues=data.table::fread(file.path(out, "out/pvalues.txt"), data.table = F)
significant_means.txt=data.table::fread(file.path(out, "out/significant_means.txt"), data.table = F)
means=data.table::fread(file.path(out, "out/means.txt"), data.table = F)
unlink(countsfile)
unlink(meta)
return(list("deconvoluted"=deconvoluted, "pvalues"=pvalues, "significant_means"=significant_means.txt, "means"=means))
}
# found bug in reference list, was missing, also min genes was 200 not 0 causing errors...
CreateBigSingleRObject.custom=function (counts, annot = NULL, project.name, xy, clusters, N = 10000,
min.genes = 0, technology = "10X", species = "Human", citation = "",
ref.list = list(), normalize.gene.length = F, variable.genes = "de",
fine.tune = T, reduce.file.size = T, do.signatures = F, do.main.types = T,
temp.dir = getwd(), numCores = SingleR.numCores){
n = ncol(counts)
s = seq(1, n, by = N)
dir.create(paste0(temp.dir, "/singler.temp/"), showWarnings = FALSE)
for (i in s) {
print(i)
A = seq(i, min(i + N - 1, n))
singler = CreateSinglerObject(counts[, A], annot = annot[A], ref.list = ref.list,
project.name = project.name, min.genes = min.genes,
technology = technology, species = species, citation = citation,
do.signatures = do.signatures, clusters = NULL, numCores = numCores)
save(singler, file = paste0(temp.dir, "/singler.temp/",
project.name, ".", i, ".RData"))
}
singler.objects.file <- list.files(paste0(temp.dir, "/singler.temp/"),
pattern = "RData", full.names = T)
singler.objects = list()
for (i in 1:length(singler.objects.file)) {
load(singler.objects.file[[i]])
singler.objects[[i]] = singler
}
singler = SingleR.Combine.custom(singler.objects, order = colnames(counts),expr = counts, clusters = clusters, xy = xy)
singler
}
# bug also in this...
SingleR.Combine.custom=function (singler.list, order = NULL, clusters = NULL, expr = NULL, xy = NULL)
{
singler = c()
singler$singler = singler.list[[1]]$singler
for (j in 1:length(singler.list[[1]]$singler)) {
singler$singler[[j]]$SingleR.cluster = c()
singler$singler[[j]]$SingleR.cluster.main = c()
singler$singler[[j]]$SingleR.single$clusters = c()
}
singler$meta.data = singler.list[[1]]$meta.data
singler$meta.data$clusters = c()
singler$meta.data$xy = c()
singler$meta.data$data.sets = rep(singler$meta.data$project.name,
length(singler$meta.data$orig.ident))
for (i in 2:length(singler.list)) {
for (j in 1:length(singler$singler)) {
if (singler.list[[i]]$singler[[j]]$about$RefData !=
singler.list[[1]]$singler[[j]]$about$RefData) {
stop("The objects are not ordered by the same reference data.")
}
singler$singler[[j]]$about$Organism = c(singler$singler[[j]]$about$Organism,
singler.list[[i]]$singler[[j]]$about$Organism)
singler$singler[[j]]$about$Citation = c(singler$singler[[j]]$about$Citation,
singler.list[[i]]$singler[[j]]$about$Citation)
singler$singler[[j]]$about$Technology = c(singler$singler[[j]]$about$Technology,
singler.list[[i]]$singler[[j]]$about$Technology)
singler$singler[[j]]$SingleR.single$labels = rbind(singler$singler[[j]]$SingleR.single$labels,
singler.list[[i]]$singler[[j]]$SingleR.single$labels)
if (!is.null(singler$singler[[j]]$SingleR.single$labels1)) {
singler$singler[[j]]$SingleR.single$labels1 = rbind(singler$singler[[j]]$SingleR.single$labels1,
singler.list[[i]]$singler[[j]]$SingleR.single$labels1)
}
singler$singler[[j]]$SingleR.single$scores = rbind(singler$singler[[j]]$SingleR.single$scores,
singler.list[[i]]$singler[[j]]$SingleR.single$scores)
singler$singler[[j]]$SingleR.single.main$labels = rbind(singler$singler[[j]]$SingleR.single.main$labels,
singler.list[[i]]$singler[[j]]$SingleR.single.main$labels)
if (!is.null(singler$singler[[j]]$SingleR.single.main$labels1)) {
singler$singler[[j]]$SingleR.single.main$labels1 = rbind(singler$singler[[j]]$SingleR.single.main$labels1,
singler.list[[i]]$singler[[j]]$SingleR.single.main$labels1)
}
singler$singler[[j]]$SingleR.single.main$scores = rbind(singler$singler[[j]]$SingleR.single.main$scores,
singler.list[[i]]$singler[[j]]$SingleR.single.main$scores)
singler$singler[[j]]$SingleR.single$cell.names = c(singler$singler[[j]]$SingleR.single$cell.names,
singler.list[[i]]$singler[[j]]$SingleR.single$cell.names)
singler$singler[[j]]$SingleR.single.main$cell.names = c(singler$singler[[j]]$SingleR.single.main$cell.names,
singler.list[[i]]$singler[[j]]$SingleR.single.main$cell.names)
if (!is.null(singler$singler[[j]]$SingleR.single.main$pval)) {
singler$singler[[j]]$SingleR.single.main$pval = c(singler$singler[[j]]$SingleR.single.main$pval,
singler.list[[i]]$singler[[j]]$SingleR.single.main$pval)
}
if (!is.null(singler$singler[[j]]$SingleR.single$pval)) {
singler$singler[[j]]$SingleR.single$pval = c(singler$singler[[j]]$SingleR.single$pval,
singler.list[[i]]$singler[[j]]$SingleR.single$pval)
}
}
singler$meta.data$project.name = paste(singler$meta.data$project.name,
singler.list[[i]]$meta.data$project.name, sep = "+")
singler$meta.data$orig.ident = c(singler$meta.data$orig.ident,
singler.list[[i]]$meta.data$orig.ident)
singler$meta.data$data.sets = c(singler$meta.data$data.sets,
rep(singler.list[[i]]$meta.data$project.name, length(singler.list[[i]]$meta.data$orig.ident)))
}
for (j in 1:length(singler$singler)) {
if (!is.null(order)) {
singler$singler[[j]]$SingleR.single$labels = singler$singler[[j]]$SingleR.single$labels[order,
]
if (!is.null(singler$singler[[j]]$SingleR.single$labels1)) {
singler$singler[[j]]$SingleR.single$labels1 = singler$singler[[j]]$SingleR.single$labels1[order,
]
}
singler$singler[[j]]$SingleR.single$scores = singler$singler[[j]]$SingleR.single$scores[order,
]
singler$singler[[j]]$SingleR.single$cell.names = singler$singler[[j]]$SingleR.single$cell.names[order]
singler$singler[[j]]$SingleR.single.main$labels = singler$singler[[j]]$SingleR.single.main$labels[order,
]
if (!is.null(singler$singler[[j]]$SingleR.single.main$labels1)) {
singler$singler[[j]]$SingleR.single.main$labels1 = singler$singler[[j]]$SingleR.single.main$labels1[order,
]
}
singler$singler[[j]]$SingleR.single.main$scores = singler$singler[[j]]$SingleR.single.main$scores[order,
]
singler$singler[[j]]$SingleR.single.main$cell.names = singler$singler[[j]]$SingleR.single.main$cell.names[order]
if (!is.null(singler$singler[[j]]$SingleR.single$pval)) {
singler$singler[[j]]$SingleR.single$pval = singler$singler[[j]]$SingleR.single$pval[order]
singler$singler[[j]]$SingleR.single.main$pval = singler$singler[[j]]$SingleR.single.main$pval[order]
}
}
}
if (!is.null(clusters) && !is.null(expr)) {
for (j in 1:length(singler$singler)) {
if (is.character(singler$singler[[j]]$about$RefData)) {
ref = get(tolower(singler$singler[[j]]$about$RefData))
}
else {
ref = singler$singler[[j]]$about$RefData
}
singler$singler[[j]]$SingleR.clusters = SingleR("cluster",
expr, ref$data, types = ref$types, clusters = factor(clusters),
sd.thres = ref$sd.thres, genes = "de", fine.tune = T)
singler$singler[[j]]$SingleR.clusters.main = SingleR("cluster",
expr, ref$data, types = ref$main_types, clusters = factor(clusters),
sd.thres = ref$sd.thres, genes = "de", fine.tune = T)
}
singler$meta.data$clusters = clusters
if (!is.null(xy)) {
singler$meta.data$xy = xy
}
}
singler
}
Datasets
Models
