library(dplyr)
library(Seurat)

# load data
pbmc.data <- Read10X(data.dir="PBMC-10x/")
pbmc <- CreateSeuratObject(counts=pbmc.data, project="pbmc3k", min.cells=3, min.features=200)
pbmc

# basic filtering
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern="^MT-")
pbmc[["reliable"]] <- pbmc[["nCount_RNA"]]<=7500 & pbmc[["percent.mt"]]<=7 # number of genes was imposed at loading
pbmc <- subset(pbmc, subset=reliable)
pbmc <- NormalizeData(pbmc)

# basic clustering
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features=all.genes)
pbmc <- FindVariableFeatures(pbmc, nfeatures=2000)
pbmc <- RunPCA(pbmc, features=VariableFeatures(pbmc), ndims.print=1:3, nfeatures.print=10)
pbmc <- FindNeighbors(pbmc, dims=1:10)
pbmc <- FindClusters(pbmc, resolution=0.5) # resolution controls the size of the communities
pbmc <- RunUMAP(pbmc, dims=1:10)
DimPlot(pbmc, reduction="umap")
pbmc <- RunTSNE(pbmc, dims=1:10)
DimPlot(pbmc, reductio="tsne")

# population markers
pbmc.markers <- FindAllMarkers(pbmc, only.pos=TRUE, min.pct=0.25, logfc.threshold=0.25)
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n=2, order_by=avg_log2FC)

# cell type calling
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
                     "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction="tsne", label=TRUE, pt.size=0.5) + NoLegend()


# SingleCellSignalR is in Bioconductor and it can be installed with BiocManager::install
# Since we are making some changes to adapt to an underlying package (circos) evolution,
# it is better to install from GitHub.com directly. This is also something that you should
# be able to do for installing "just published" packages.

# install devtools before loading it
library(devtools)
# for windows users, install Rtools, a separate application that is not a R package
devtools::install_github(repo="https://github.com/SCA-IRCM/SingleCellSignalR")

library(SingleCellSignalR)
counts <- as.matrix(GetAssayData(object=pbmc, slot="counts"))
counts[1:10,1:5]
q <- apply(counts, 2, function(x) quantile(x[x>0],0.99))
ncounts <- log(1+sweep(counts, 2, q/median(q), "/"))
ncounts[1:10,1:5]
cluster.num <- as.numeric(Idents(pbmc))
LR.inter <- cell_signaling(ncounts, rownames(ncounts), cluster.num, c.names=new.cluster.ids)
visualize_interactions(LR.inter)
visualize_interactions(LR.inter, show.in=c("B-NK"))