单细胞转录组（scRNA）

getwd()
setwd("~/data_analysis/Scanpy分析/")

#method1----------------------------------------------------------------------------
library(SeuratDisk)
Convert('./pbmc3k.h5ad', dest = "h5seurat", overwrite = TRUE)
seuratObj <- LoadH5Seurat("./pbmc3k.h5seurat", meta.data = FALSE, misc = FALSE)

meta <- read.csv('./meta.csv',header = T)
seuratObj@meta.data <- meta
rownames(seuratObj@meta.data) <- meta$X

Idents(seuratObj) <- 'leiden'
DimPlot(seuratObj)

df_markers <- FindAllMarkers(seuratObj, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)


#method2----------------------------------------------------------------------------
devtools::install_github("scverse/anndataR")
library(anndataR)
# adata.raw.to_adata().write("./pbmc3k_withoutX.h5ad")
adata <- anndataR::read_h5ad("./pbmc3k_withoutX.h5ad", to = "HDF5AnnData")
seuratObj2 <- adata$to_Seurat()


a = adata$obsm$X_umap
colnames(a) <- c("umap_1","umap_2")
rownames(a) <- adata$obs_names
A <- Seurat::CreateDimReducObject(as.matrix(a), key="umap",assay = 'RNA')

seuratObj2@reductions[['umap']] <- A


Idents(seuratObj2) <- 'leiden'
DimPlot(seuratObj2)

# df_markers <- FindAllMarkers(seuratObj2, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)

#method3----------------------------------------------------------------------------
#https://github.com/cellgeni/schard
devtools::install_github("cellgeni/schard")


seuratObj3 = schard::h5ad2seurat('pbmc3k.h5ad',use.raw=T)

Idents(seuratObj3) <- 'leiden'
DimPlot(seuratObj3)
df_markers <- FindAllMarkers(seuratObj3, logfc.threshold = 0.2, min.pct = 0.2, only.pos = T)


#method4----------------------------------------------------------------------------
#如果可以就用，不行就放弃使用别的方法
# devtools::install_github("cellgeni/sceasy")
# library(sceasy)
# library(reticulate)
# use_condaenv('sceasy')
# loompy <- reticulate::import('loompy')
# 
# ad_path <- "./pbmc3k.h5ad"
# sceasy::convertFormat(ad_path, from="anndata", to="seurat", outFile="file.rds")

#详细参考：https://bioconductor.org/packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.html

###GSVA codes###
library(Seurat)
library(RColorBrewer)
library(GSVA)
library(GSEABase)
setwd("D:\\KS项目\\公众号文章\\GSVA_reanalyiss")
human_data <- readRDS("D:/KS项目/公众号文章/human_data.rds")

##1.GSVA分析gene set可在GSEA官网下载，也可以使用msigdbr包获取
#msigdbr包的可以参考我们之前的帖子https://mp.weixin.qq.com/s/Ic7v6gpENJ1GDTFMAqSVoA
#比如这里我想要分析human KEGG在GSEA官网下载人的gene set gmt文件：https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp
genesets <- getGmt('./c2.cp.v2023.2.Hs.symbols.gmt')

##2.获取单细胞matrix,这里用counts或者normalized都行
gene_matrix <- GetAssayData(human_data, layer = 'data')

##3. GSVA analysis,如果矩阵室counts kcdf="Poisson"，normalized用默认
GSVA <- gsva(as.matrix(gene_matrix), genesets, min.sz=10, max.sz=1000,verbose=TRUE)
head(rownames(GSVA))
# [1] "SA_B_CELL_RECEPTOR_COMPLEXES" "SA_CASPASE_CASCADE"          
# [3] "SA_G1_AND_S_PHASES"           "SA_MMP_CYTOKINE_CONNECTION"  
# [5] "SA_PROGRAMMED_CELL_DEATH"     "SA_PTEN_PATHWAY"  
rownames(GSVA) <- tolower(rownames(GSVA))#将大写换为小写
rownames(GSVA) <- gsub('_',' ',rownames(GSVA))#将_转化为空格
library(Hmisc)
rownames(GSVA) <- capitalize(rownames(GSVA))#首字母大写

##4. 差异分析
meta=human_data@meta.data[,c("orig.ident","celltype","group")]
row.names(meta)=colnames(GSVA)
GSVA_Seurat <- CreateSeuratObject(counts = GSVA, meta.data = meta, project = "GSVA_singleCell")
#celltype
Idents(GSVA_Seurat) <- "celltype"
GSVA_celltype=FindAllMarkers(GSVA_Seurat, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.1)

#group
Macrophage <- subset(GSVA_Seurat, celltype=='Macrophage')
GSVA_DE = FindMarkers(Macrophage, ident.1 = "GM",ident.2 = "BM",
                      group.by = "group", min.pct = 0, logfc.threshold = 0)

#火山图
library(ggplot2)
library(ggrepel)
ggplot(GSVA_DE, aes(x=avg_log2FC, y=-log10(p_val_adj))) +
  geom_point(color = "#B3DFB2",size=2)+
  geom_vline(aes(xintercept=0), color = "#999999", linetype="dashed", size=1) + #添加纵线
  geom_point(data = GSVA_DE[GSVA_DE$p_val_adj <0.01&(GSVA_DE$avg_log2FC >0.25),],  size=3, shape=16, color="#AE1A41") +
  geom_point(data = GSVA_DE[GSVA_DE$p_val_adj <0.01&(GSVA_DE$avg_log2FC < -0.25),], size=3, shape=16,color='#3C91BB') + 
  labs(x = "Log2 fold change",y = "-Log10(P-value)", title = "") + #标题横纵轴名称
  theme_classic() 

#Dotplot
DotPlot(GSVA_Seurat, features = c("Wp circadian rhythm genes",
                                  "Reactome pi3k akt signaling in cancer",
                                  "Kegg jak stat signaling pathway",
                                  "Wp aerobic glycolysis",
                                  "Wp glycolysis in senescence",
                                  "Kegg medicus reference translation initiation"),group.by = "orig.ident")+coord_flip()+NoLegend()+
  theme(axis.title = element_blank())

#heatmap
GSVA_exp <- AverageExpression(GSVA_Seurat,features = rownames(GSVA),group.by = 'celltype', assays = 'RNA', slot = "counts") 
GSVA_exp <- as.data.frame(GSVA_exp$RNA)

pathways <- c("Biocarta thelper pathway",
              "Biocarta tcytotoxic pathway",
              "Reactome digestion and absorption",
              "Reactome aspirin adme",
              "Biocarta il12 pathway",
              "Biocarta no2il12 pathway",
              "Reactome response to metal ions",
              "Kegg medicus reference translation initiation",
              "Reactome pd 1 signaling",
              "Kegg asthma",
              "Kegg arachidonic acid metabolism")

gsva_plot <- GSVA_exp[pathways,]

library(ComplexHeatmap)
pheatmap(gsva_plot,cluster_rows = F,cluster_cols = F,scale = 'none',
         colorRampPalette(c("#2166AC",'#478ABF','#90C0DC', "white",'#EF8C65','#CF4F45',"#B2182B"))(100),
         border=FALSE,cellwidth = 25, cellheight = 25,heatmap_legend_param = list(title="GSVA score"))

#================================================================================

#                        Seurat V5 

#================================================================================
#install packages
install.packages('Seurat')
library(Seurat)
#安装一些额外的包
setRepositories(ind = 1:3, addURLs = c('https://satijalab.r-universe.dev', 'https://bnprks.r-universe.dev/'))
install.packages(c("BPCells", "presto", "glmGamPoi"))
remotes::install_github("satijalab/seurat-data", quiet = TRUE)
remotes::install_github("satijalab/azimuth", quiet = TRUE)
remotes::install_github("satijalab/seurat-wrappers", quiet = TRUE)
# If users encounter any errors related to the Matrix package, please resolve by re-installing the TFBSTools package using the command below and opening a fresh R session:
BiocManager::install("TFBSTools", type = "source", force = TRUE)

#================================================================================

#                        Seurat V5 read data and QC

#================================================================================

setwd("~/data_analysis/seuratV5_test/scRNA_analysis/")

#read data and creat seurat obj
WT <- Read10X("./scRNA_data/WT_E18/")
WT <- WT$`Gene Expression`
WT <- CreateSeuratObject(counts = WT, project = "WT", min.cells = 3, min.features = 200)

GO <- Read10X("./scRNA_data/GO_E18/")
GO <- GO$`Gene Expression`
GO <- CreateSeuratObject(counts = GO, project = "GO", min.cells = 3, min.features = 200)


#Standard pre-processing workflow
#线粒体比例
WT[["percent.mt"]] <- PercentageFeatureSet(WT, pattern = "^mt-")
GO[["percent.mt"]] <- PercentageFeatureSet(GO, pattern = "^mt-")

#血红蛋白基因
WT[["percent.hb"]] <- PercentageFeatureSet(WT, pattern = "^Hb[^(p)]")
GO[["percent.hb"]] <- PercentageFeatureSet(GO, pattern = "^Hb[^(p)]")

#核糖体基因
WT[["percent.rb"]] <- PercentageFeatureSet(WT, pattern = "^Rbs|Rpl")
GO[["percent.rb"]] <- PercentageFeatureSet(GO, pattern = "^Rbs|Rpl")

#细胞周期
s.genes=Seurat::cc.genes.updated.2019
g2m.genes=Seurat::cc.genes.updated.2019$g2m.genes
test.seu <- CellCycleScoring(test.seu, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)




options(warn = 1)
# VlnPlot(WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), layer = 'counts')
p1= VlnPlot(WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
# Error: (converted from warning) Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
p2 = VlnPlot(GO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p1/p2

#QC质控
WT <- subset(WT, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 30000 &  percent.mt < 5)
GO <- subset(GO, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 30000 &  percent.mt < 5)

p3= VlnPlot(WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p4 = VlnPlot(GO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.hb","percent.rb"),pt.size = 0.1, ncol = 5)
p3/p4

#================================================================================

#                        Seurat V5 Integrated Data

#================================================================================
library(Seurat)
library(clustree)
library(ggplot2)
library(SeuratWrappers)
library(rlang)

setwd("~/data_analysis/seuratV5_test/scRNA_analysis/")
# save(WT, GO, file='./scRNA_data/sample_seurat_obj.RData')
#merge data
sce <- merge(WT, y=GO,add.cell.ids = c("WT", "GO"))


#NormalizeData & ScaleData
sce <- NormalizeData(sce)
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce, vars.to.regress = c("percent.mt"))
sce <- RunPCA(sce, verbose=F)

# methods of Integration
# CCA integration (method=CCAIntegration)
# RPCA integration (method=RPCAIntegration)
# Harmony (method=HarmonyIntegration)
# JointPCA (method= JointPCAIntegration)

# FastMNN (method= FastMNNIntegration)
# scVI (method=scVIIntegration)

###Integrated with CCA
sce_cca <- IntegrateLayers(object = sce, 
                           method = CCAIntegration, 
                           orig.reduction = "pca", 
                           new.reduction = "integrated.cca",verbose = FALSE)

# re-join layers after integration
sce_cca[["RNA"]] <- JoinLayers(sce_cca[["RNA"]])

###Integrated with harmony

sce_harm <- IntegrateLayers(object = sce, 
                            method = HarmonyIntegration, 
                           orig.reduction = "pca", 
                           new.reduction = "integrated.Harmony",
                           verbose = FALSE)
# re-join layers after integration
sce_harm[["RNA"]] <- JoinLayers(sce_harm[["RNA"]])

###Integrated with scVI
#first, install scVI-tool before IntegrateLayers if you did not install previously
# conda create -n scvi-env python=3.8
# conda activate scvi-env
# pip install scvi-tools -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com
# pip install scanpy -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com

sce_scvi <- IntegrateLayers(object = sce, 
                            method = scVIIntegration, 
                            orig.reduction = "pca", 
                            new.reduction = "integrated.scvi",
                            conda_env="~/miniconda3/envs/scvi-env",
                            verbose = FALSE)

# re-join layers after integration
sce_scvi[["RNA"]] <- JoinLayers(sce_scvi[["RNA"]])

#================================================================================
#Perform cca reduction
Seurat::ElbowPlot(sce_cca, ndims = 50)
sce_cca <- FindNeighbors(sce_cca, reduction = "integrated.cca", dims = 1:20)
sce_cca <- FindClusters(sce_cca, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_cca <- RunUMAP(sce_cca, dims = 1:20, reduction = "integrated.cca")
# clustree(sce_cca)


#Perform charmony reduction
Seurat::ElbowPlot(sce_harm, ndims = 50)
sce_harm <- FindNeighbors(sce_harm, reduction = "integrated.Harmony", dims = 1:20)
sce_harm <- FindClusters(sce_harm, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_harm <- RunUMAP(sce_harm, dims = 1:20, reduction = "integrated.Harmony")
# clustree(sce_harm)


#Perform scVI reduction
Seurat::ElbowPlot(sce_scvi, ndims = 50)
sce_scvi <- FindNeighbors(sce_scvi, reduction = "integrated.scvi", dims = 1:20)
sce_scvi <- FindClusters(sce_scvi, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_scvi <- RunUMAP(sce_scvi, dims = 1:20, reduction = "integrated.scvi")
# clustree(sce_scvi)


#cluster plot
DimPlot(sce_cca, reduction = "umap", group.by = "orig.ident")+
  ggtitle("CCA")
DimPlot(sce_harm, reduction = "umap", group.by = "orig.ident")+
  ggtitle("Harmony")
DimPlot(sce_scvi, reduction = "umap", group.by = "orig.ident")+
  ggtitle("scvi")


DimPlot(sce_cca, reduction = "umap", label = T)+
  ggtitle("CCA")
DimPlot(sce_harm, reduction = "umap", label = T)+
  ggtitle("Harmony")
DimPlot(sce_scvi, reduction = "umap", label = T)+
  ggtitle("scvi")

#================================================================================
#SCTransform-normalized
#merge data
sce <- merge(WT, y=GO,add.cell.ids = c("WT", "GO"))
sce <- SCTransform(sce)
sce <- RunPCA(sce, verbose=F)
sce_sct_cca <- IntegrateLayers(object = sce, method = CCAIntegration, normalization.method = "SCT", verbose = F)
sce_sct_cca[["RNA"]] <- JoinLayers(sce_sct_cca[["RNA"]])

sce_sct_cca <- FindNeighbors(sce_sct_cca, reduction = "integrated.dr", dims = 1:25)
sce_sct_cca <- FindClusters(sce_sct_cca, resolution = seq(from = 0.1, to = 1.0, by = 0.1))
sce_sct_cca <- RunUMAP(sce_sct_cca, dims = 1:25, reduction = "integrated.dr")
save(sce_sct_cca, file = "sce_sct_cca.RData")


#plot
DimPlot(sce_sct_cca, reduction = "umap", group.by = "orig.ident")
DimPlot(sce_sct_cca, reduction = "umap", label = T)

#================================================================================

#                        Cellannotation & labels

#================================================================================
library(Seurat)
library(ggplot2)

Allmarkers <- FindAllMarkers(sce_cca, logfc.threshold = 0.3, min.pct = 0.3)
write.csv(Allmarkers, file = 'Allmarkers.csv')
#================================================================================
#Manual annotation, reference to published articles
markers <- c("Pparg", "Myh11", "Mrc1", "Flt1", "Col11a1", "Mymk", "Pax7", "Pdgfra","Ttn","Sox2")
DotPlot(sce_cca, features = markers, col.min = 0)+coord_flip()

#20 Adipocytes
#19 SMC
#13 Macrophages
#14 Endothelial
#7,10,21 Tenocytes
#9 Myoblasts
#11,12 MuSCs
#0,1,3,22 Mesenchymal
#2,4,5,6,8,15,16,18 Myonuclei
#17 NPCs

#================================================================================
#singleR,这里仅仅是演示一下操作，这个效果一般
library(SingleR)
# load("~/data_analysis/Zhou_2024/re-analysis/ref_Mouse_all.RData")
# SingleR_data <- GetAssayData(sce_cca, layer="data")
# SingleR_cluster <- sce_cca$seurat_clusters
# #annotation
# All_singleR <- SingleR::SingleR(test = SingleR_data, 
#                           ref = ref_Mouse, 
#                           labels = ref_Mouse$label.fine, 
#                           method = "cluster", 
#                           clusters = SingleR_cluster, 
#                           assay.type.test = "logcounts", 
#                           assay.type.ref = "logcounts")
# 
# All_celltype = data.frame(ClusterID=rownames(All_singleR),
#                            celltype=All_singleR$labels, stringsAsFactors = FALSE)
# write.csv(All_celltype,"All_celltype.csv",row.names = FALSE)


#================================================================================
#cell annotation

sce_cca$celltype <- "no"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('20'))] <- "Adipocytes"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('19'))] <- "SMC"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('13'))] <- "Macrophages"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('14'))] <- "Endothelial"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('7','10','21'))] <- "Tenocytes"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('9'))] <- "Myoblasts"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('17'))] <- "NPCs"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('11','12'))] <- "MuSCs"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('0','1','3','22'))] <- "Mesenchymal"
sce_cca$celltype[which(sce_cca$seurat_clusters %in% c('2','4','5','6','8','15','16','18'))] <- "Myonuclei"

Idents(sce_cca) <- 'celltype'
DimPlot(sce_cca, reduction = "umap", label = T)

#save files 
save(sce_cca, file = 'sce_cca_cellltype.RData')


#================================================================================
#marker plot
markers <- c("Pparg", "Myh11", "Mrc1", "Flt1", "Col11a1", "Mymk", "Pax7", "Pdgfra","Ttn","Sox2")
DotPlot(sce_cca, features = markers, col.min = 0)+
  theme_bw()+theme(axis.title = element_blank(),
                   axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0),
                   panel.grid.major = element_blank(),
                   panel.grid.minor = element_blank())


VlnPlot(sce_cca, features = markers, stack=T)+
  theme(axis.title = element_blank())+NoLegend()

#================================================================================

#                        Cell-Cycle

#================================================================================
library(Seurat)

# For human data, Cell-Cycle genes 
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

#For mouse data, Cell-Cycle genes
s.genes <- c("Mcm4", "Exo1", "Slbp", "Gmnn", "Cdc45","Msh2", "Mcm6",
             "Rrm2", "Pold3","Blm",  "Ubr7", "Mcm5", "Clspn","Hells",   
             "Nasp", "Rpa2", "Rad51ap1", "Tyms", "Rrm1", "Rfc2", "Prim1",   
             "Brip1","Usp1", "Ung",  "Pola1","Mcm2", "Fen1", "Tipin" ,  
             "Pcna", "Cdca7","Uhrf1","Casp8ap2","Cdc6", "Dscc1","Wdr76",   
             "E2f8", "Dtl",  "Ccne2","Atad2","Gins2","Chaf1b","Pcna-ps2")

g2m.genes <- c("Nuf2", "Psrc1","Ncapd2","Ccnb2","Smc4", "Lbr",  "Tacc3",  
               "Cenpa","Kif23","Cdca2","Anp32e","G2e3", "Cdca3","Anln",
               "Cenpe","Gas2l3","Tubb4b","Cenpf","Dlgap5","Hjurp","Cks1brt",
               "Gtse1","Bub1", "Birc5","Ube2c","Rangap1","Hmmr", "Ect2",
               "Tpx2", "Ckap5","Cbx5", "Nek2", "Ttk",  "Cdca8","Nusap1", 
               "Ctcf", "Cdc20","Cks2", "Mki67","Tmpo", "Ckap2l","Aurkb",  
               "Kif2c","Cdk1", "Kif20b","Top2a","Aurka","Ckap2","Hmgb2",  
               "Cdc25c","Ndc80","Kif11")

#load data 
# load("~/data_analysis/seuratV5_test/scRNA_analysis/sce_cca_cellltype.RData")
#Assign Cell-Cycle Scores
sce_cca <- CellCycleScoring(sce_cca, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
# plot 
DimPlot(sce_cca, group.by = 'Phase')
DimPlot(sce_cca, group.by = 'Phase', reduction = 'pca')

#If cell cycle affect the cluster in your data 
#Please Regress out cell cycle scores during data scaling

#================================================================================

#                        Differentially expressed genes analysis

#================================================================================
# load("~/data_analysis/seuratV5_test/scRNA_analysis/sce_cca_cellltype.RData")

#for normal scaled data
#所有细胞类型中两组的差异
celltypes <- unique(sce_cca$celltype)
DEGs_celltype <- list()

for (i in 1:length(celltypes)) {
  
  data = subset(sce_cca, celltype==celltypes[i])
  deg = FindMarkers(data,
                    group.by="orig.ident",
                    ident.1 = 'GO',
                    ident.2 = "WT",
                    logfc.threshold=0.25,
                    min.pct = 0.25)
  
  DEGs_celltype[[i]] <- deg
  names(DEGs_celltype)[i] <- celltypes[i]
  
}


#for SCT data
# load("~/data_analysis/seuratV5_test/scRNA_analysis/sce_sct_cca.RData")

# SCT_deg <- FindMarkers(sce_sct_cca,
#                        group.by="orig.ident",
#                        ident.1 = 'GO',
#                        ident.2 = "WT",
#                        logfc.threshold=0.25,
#                        min.pct = 0.25)
# Error in FindMarkers.SCTAssay(object = data.use, slot = slot, cells.1 = cells$cells.1,  : 
#                                 Object contains multiple models with unequal library sizes. Run `PrepSCTFindMarkers()` before running `FindMarkers()`.

# perform differential expression
sce_sct_cca <- PrepSCTFindMarkers(sce_sct_cca)
sce_sct_cca$newgroup <- paste(sce_sct_cca$orig.ident, sce_sct_cca$seurat_clusters, sep = "_")
Idents(sce_sct_cca) <- "newgroup"
unique(sce_sct_cca$newgroup)

DEGs_0 <- FindMarkers(sce_sct_cca,
                      ident.1 = "GO_0", 
                      ident.2 = "WT_0",
                      logfc.threshold=0.25,
                      min.pct = 0.25)