preprocess_integrate_present_study_with_Bassler_control.Rmd

---
title: "ILD project seurat analysis"
author: "Jiangyan Yu (jiangyan.yu@uni-bonn.de)"
date: "`r Sys.Date()`"
output:
  html_document: 
    code_folding: hide
    number_sections: yes
    theme: united
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
---
notes: 
1) integrate all cells from IPF project with control cells from Kevin's COPD paper
2) script is executed in docker: docker run –rm -p 8888:8787 -e USER=jyu -e PASSWORD=jyu -e ROOT=TRUE -v /home/agschlitzer/sciebo/Projects/IPF_ownerJY:/home/jyu/rstudio jiangyanyu/jyu_r4.1.2:1.0
3) report generated by: Rscript -e "rmarkdown::render('integrate_persent_study_with_Bassler_control.Rmd')"
# general steps

```{r global_options}
knitr::opts_chunk$set(warning=FALSE, messgae=FALSE, fig.path='Figs/', results = "hide")
## fig.width=4, fig.height=4
```

## load library

```{r include=FALSE}
rm(list=ls())
gc()
#CRAN packages
list.of.packages <- c("readr",
                      "cowplot",
                      "useful",
                      "stringr",
                      "umap",
                      "ggplot2",
                      "reshape2",
                      "dplyr",
                      "tidyr",
                      "Matrix.utils",
                      "VGAM",
                      "plotly",
                      "future",
                      "data.table"
)

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)>0) install.packages(new.packages)

#BioconductoR packages
list.of.bioc.packages <- c("tximport",
                           "DESeq2",
                           "Seurat",
                           "slingshot",
                           "flowCore",
                           "biomaRt",
                           "clusterProfiler",
                           "org.Hs.eg.db",
                           "GSEABase",
                           "DOSE",
                           "BiocGenerics",
                           "DelayedArray",
                           "DelayedMatrixStats",
                           "limma",
                           "S4Vectors",
                           "SingleCellExperiment",
                           "SummarizedExperiment",
                           "batchelor")
# "annotables")
new.packages.bioc <- list.of.bioc.packages[!(list.of.bioc.packages %in% installed.packages()[,"Package"])]
# 
if(length(new.packages.bioc)>0)if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install(new.packages.bioc, update = FALSE)

lapply(c(list.of.packages,list.of.bioc.packages), require, character.only = TRUE)
```

## parallel computing

multicore does not work in rstudio. better to use plan::multiprocess

```{r include=FALSE}

# library(parallel)
# plan("multiprocess", workers = detectCores()-1)
# options(future.globals.maxSize = 20000 * 1024^2)
```
## working directory

```{r}
working.dir = "/home/jyu/rstudio/"
```

# upload IPF objects

```{r}
print("prepare seurat object list...")
ILD_pheno = readxl::read_xlsx(path = paste0(working.dir,"/IPF_manuscript/meta_data/ILD_FACs_SeqWell_phenotype_overview.xlsx"),sheet = "Sheet1")
ILD_pheno = ILD_pheno[,c("patient","Sampling_date","Seqwell","Age","Gender","Disease","DiseaseGroup","hospital")]
seurat_list = list()

for (i in c(1:4)){

  ## read into saved seurat objects
  tmp_seurat = readRDS(file = paste0(working.dir,"/scRNA_processed/4batches_rds/ILD_batch",i,"_seurat.rds"))
  tmp_seurat$batch = paste0("batch",i)

  ## add clinical information to seurat objects
  ILD_pheno_tmp = ILD_pheno[ILD_pheno$patient %in% tmp_seurat$patient,]
  pheno_meta = merge(tmp_seurat@meta.data,ILD_pheno_tmp,by.x="patient",by.y="patient")
  rownames(pheno_meta) = pheno_meta$cell
  tmp_seurat= AddMetaData(tmp_seurat,metadata = pheno_meta[,c("Sampling_date","Age","Gender","Disease","DiseaseGroup","hospital")])
  selected_cells = tmp_seurat@meta.data[tmp_seurat@meta.data$DiseaseGroup %in% c("Control","IPF"),]
  tmp_seurat = tmp_seurat[, rownames(selected_cells)]
  tmp_seurat = CreateSeuratObject(counts = tmp_seurat@assays$RNA,meta.data = tmp_seurat@meta.data[,c("patient","Age","Gender","Disease","hospital")])
  ## save tmp seurat to
  seurat_list[[i]] = tmp_seurat
}

rm(tmp_seurat,pheno_meta,ILD_pheno_tmp,i, selected_cells)
```

# upload copd dataset

```{r}
copd_bal = readRDS(file = paste0(working.dir,"/scRNA_processed/Kevin_copd_bal_analyzedJY20220825.rds"))
copd_bal$batch = paste0("copd_ck")
copd_meta = read.csv(paste0(working.dir,"/IPF_manuscript/meta_data/COPD_meta_data.csv"))
pheno_meta = merge(copd_bal@meta.data,copd_meta,by.x="orig.ident",by.y="patient1")
rownames(pheno_meta) = rownames(copd_bal@meta.data)
tmp_seurat= AddMetaData(copd_bal,metadata = pheno_meta[,c("patient","Sampling_date","Age","Gender","Disease","DiseaseGroup","hospital")])
selected_cells = tmp_seurat@meta.data[tmp_seurat@meta.data$Disease %in% c("Control"),]
tmp_seurat = tmp_seurat[, rownames(selected_cells)]
tmp_seurat = CreateSeuratObject(counts = tmp_seurat@assays$RNA,meta.data = tmp_seurat@meta.data[,c("patient","Age","Gender","Disease","hospital")])
seurat_list[[5]] = tmp_seurat

rm(copd_bal,copd_meta,tmp_seurat,pheno_meta, selected_cells)
```

# integrate

```{r}
ILD_combined_sct_list = merge(x=seurat_list[[1]],y=list(seurat_list[[2]],seurat_list[[3]],seurat_list[[4]],seurat_list[[5]]))

rm(seurat_list)

batch_list = SplitObject(ILD_combined_sct_list,split.by = "patient")

print("find integration features...")
## use SCT transformation
## Note that this single command replaces NormalizeData(), ScaleData(), and FindVariableFeatures().
batch_list <- lapply(X = batch_list, FUN = function(x) {
  x <- NormalizeData(x, verbose = FALSE)
  x <- FindVariableFeatures(x, verbose = FALSE)
})

features <- SelectIntegrationFeatures(object.list = batch_list)
# saveRDS(features,"features.rds")
batch_list <- lapply(X = batch_list, FUN = function(x) {
  x <- ScaleData(x, features = features, verbose = FALSE)
  x <- RunPCA(x, features = features, verbose = FALSE)
})

# saveRDS(batch_list,"batch_list.rds")

### find anchors using rpca
print("find anchors and integrate data...")
batch_anchors <- FindIntegrationAnchors(object.list = batch_list, reduction = "rpca",
                                        dims = 1:50)
# saveRDS(batch_anchors,"batch_anchors.rds")

# batch_anchors = readRDS("/home/jyu/rstudio/batch_anchors.rds")

ILD_combined_sct <- IntegrateData(anchorset = batch_anchors, dims = 1:50)


rm(batch_list, features, batch_anchors)
```

# UMAP

```{r}
DefaultAssay(ILD_combined_sct) = "integrated"

print("UMAP")
ILD_combined_sct <- ScaleData(ILD_combined_sct, verbose = FALSE)
ILD_combined_sct <- RunPCA(ILD_combined_sct, verbose = FALSE)
ILD_combined_sct <- RunUMAP(ILD_combined_sct,
                            reduction = "pca",
                            dims = 1:30,
                            n.neighbors = 30,
                            min.dist = 0.000001,
                            n.components = 3,
                            spread = 4)


DefaultAssay(ILD_combined_sct) <- "integrated"
ILD_combined_sct <- FindNeighbors(ILD_combined_sct, reduction = "pca", dims = 1:30)
ILD_combined_sct <- FindClusters(ILD_combined_sct, resolution = 1)

# pdf(paste0(working.dir,"umap.pdf"))
DimPlot(ILD_combined_sct, reduction = "umap", label = TRUE,dims = c(1,2),raster=TRUE)
DimPlot(ILD_combined_sct, reduction = "umap", label = TRUE,dims = c(1,3),raster=TRUE)
DimPlot(ILD_combined_sct, reduction = "umap", label = TRUE,dims = c(2,3),raster = TRUE)

# dev.off()

# ILD_combined_sct = readRDS("/home/jyu/rstudio/IPF_manuscript/data/ILD_combined_sct_rpca20230720.rds")
```


# find subcluster

```{r}
# ILD_combined_sct = readRDS("/home/jyu/rstudio/IPF_manuscript/data/ILD_combined_sct_rpca20230720.rds")
Idents(ILD_combined_sct) = "seurat_clusters"
ILD_combined_sct = FindSubCluster(ILD_combined_sct,
                                   cluster = 4,
                                   graph.name = "integrated_snn",
                                   subcluster.name = "subcluster",
                                   resolution = 0.16,
                                   algorithm = 1)

DimPlot(ILD_combined_sct,group.by = "subcluster",dims = c(1,3))
```

# QC

## mt% 

```{r}
DefaultAssay(ILD_combined_sct) <- "RNA"
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
ILD_combined_sct[["percent.mt"]] <- PercentageFeatureSet(ILD_combined_sct, pattern = "^MT-")
VlnPlot(ILD_combined_sct, features = c("percent.mt"), ncol = 1,pt.size = 0)

```

## find doublets

using two approaches: 1) in silico estimation of doublet score using computeDoubletDensity; 2) check # of degs of suspicious clusters. if low # of deg, then doublets

```{r}
# BiocManager::install("scDblFinder")
library(scDblFinder)
## convert RNA assay which contains counts
sce = Seurat::as.SingleCellExperiment(ILD_combined_sct)
## specify cluster information for each cell
dbl.out <- findDoubletClusters(sce,clusters = ILD_combined_sct$integrated_snn_res.1)

dbl.out
sce_inte = Seurat::as.SingleCellExperiment(ILD_combined_sct,assay = "integrated")
dbl.dens <- computeDoubletDensity(sce, d=ncol(reducedDim(sce_inte)))
summary(dbl.dens)

sce_inte$DoubletScore <- dbl.dens
ILD_combined_sct$DoubletScore = dbl.dens
VlnPlot(ILD_combined_sct,"DoubletScore",pt.size = 0)+
  ylim(c(0,10))

# plotUMAP(sce_inte, colour_by="DoubletScore")
# dbl.calls <- doubletThresholding(data.frame(score=dbl.dens),
#     method="griffiths", returnType="call")
# summary(dbl.calls)
# 
# tmp_dbl = ILD_combined_sct@meta.data[,c("DoubletScore","seurat_clusters")]
# rownames(tmp_dbl) = rownames(ILD_combined_sct@meta.data)

# write.csv(tmp_dbl, file = paste0(working.dir,"/ILD_combined_sct_doubletscore.csv") )

# plotColData(sce_inte, x="integrated_snn_res.1", y="DoubletScore", colour_by=I(dbl.calls))
```

# select single live cells

```{r}
ILD_combined_sct = subset(ILD_combined_sct,idents=c(0,4:17,19:22,24))

ILD_combined_sct = subset(ILD_combined_sct,subset= patient != "pnt2718174")
```

## change patient name

```{r}
tmp_meta = ILD_combined_sct@meta.data
tmp_meta$cell = rownames(tmp_meta)
tmp_meta$order = c(1:nrow(tmp_meta))
# ,"subcluster"
## remove some columns
tmp_meta = tmp_meta[,c("orig.ident","nCount_RNA","nFeature_RNA","patient","Age","Gender","Disease","hospital","percent.mt","integrated_snn_res.1","cell","order")]

pnt_meta = readxl::read_xlsx(path = paste0(working.dir,"/manuscript/figures_20230524/data/pnt_number_in_paper.xlsx"),sheet = "Sheet1")

tmp_meta = merge(tmp_meta,pnt_meta,by.x = "patient",by.y="Patient")
tmp_meta = tmp_meta[order(tmp_meta$order),]

# tmp_meta = tmp_meta[,c(3:ncol(tmp_meta))]

tmp1 = tmp_meta$cell %>% stringr::str_split_fixed(pattern = "_",2) %>% as.data.frame()
tmp2 = tmp1$V1 %>% stringr::str_split_fixed(pattern = "-",3) %>% as.data.frame()
tmp2$V2 = ifelse(tmp2$V2 == "",rownames(tmp2),tmp2$V2)
tmp2$patient = tmp_meta$patient
tmp2$order = c(1:nrow(tmp2))

tmp2 = merge(tmp2,pnt_meta,by.x = "patient",by.y = "Patient")
tmp2 = tmp2[order(tmp2$order),]

tmp_meta$cellID = paste0(tmp2$pnt_in_paper,"-",tmp2$V2,"-",tmp2$V3,"_",tmp1$V2)

rm(tmp1,tmp2)


rownames(tmp_meta) = tmp_meta$cell

ILD_combined_sct = AddMetaData(ILD_combined_sct,
                                       metadata = tmp_meta[,c("pnt_in_paper","Stage","cellID","cell"),],
                               col.name = c("pnt_in_paper","IPFgroup","cellID","cell"))

ILD_combined_sct = RenameCells(ILD_combined_sct,new.names = ILD_combined_sct$cellID)

rm(tmp_meta,pnt_meta)


```

# save seurat

```{r}
ILD_combined_sct$orig.ident = NULL
ILD_combined_sct$patient = NULL
# ILD_combined_sct$integrated_snn_res.0.3 = NULL
# ILD_combined_sct$defined = NULL
# ILD_combined_sct$cell = NULL

saveRDS(ILD_combined_sct,paste0(working.dir,"/IPF_manuscript/data/ILD_combined_sct_rpca20230720.rds"))
```

# DEG

```{r,fig.width=8,fig.height=20}
print("find DEGs...")
DefaultAssay(ILD_combined_sct) = "RNA"
ILD_combined_sct.DEG = FindAllMarkers(ILD_combined_sct)

selected_cells = ILD_combined_sct@meta.data
selected_cells$cell = rownames(selected_cells)
selected_cells = selected_cells %>% group_by(seurat_clusters) %>% sample_n(100)


# pdf(file = paste0(working.dir,"DEG_heatmap.pdf"),width = 10,height = 15)
DoHeatmap(ILD_combined_sct[,selected_cells$cell],
          features = ILD_combined_sct.DEG %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) %>% .$gene,
          slot = "scale.data",
          assay = "integrated")
# dev.off()


write.csv(ILD_combined_sct.DEG,"/home/jyu/rstudio/ILD_combined_sct_rpca.DEG.csv")
```

# session info

```{r}
sessionInfo()
```