批次效應強時脑奠，用CCA進行校正

先用PCA圖查看不整合前的批次效應
原理在每個批次中找到相同的一組表達差異低的基因，用這組基因在每個亞群中的差異代表整體基因的差異箱锐，按照這個差異調(diào)整整體基因表達水平

# 該代碼用于理解并分析批次差異比較大的樣本躏尉，通過CCA算法將批次效應進行校正
# 并基于校正之后的數(shù)據(jù)進行下游分析

#加載分析使用的包
library(Seurat)
library(monocle)
library(devtools)
devtools::install_github("cole-trapnell-lab/monocle-release@develop")

library(monocle)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)

#讀入原始表達數(shù)據(jù)
exp1<-read.table("./CCA_monocle2_data/01.CRC01-CRC02_500sc_Read_Counts.txt",
                 sep="\t",header=T,row.names=1)
colnames(exp1)<-paste(colnames(exp1),"read",sep="_")#第一個樣本叫read

exp2<-read.table("./CCA_monocle2_data/02.CRC03-04-06-09-10-11_688sc_UMI_Counts.txt",
                 sep="\t",header=T,row.names=1)
colnames(exp2)<-paste(colnames(exp2),"umi",sep="_") #第二個樣本叫umi

#不使用CCA去批次效應蓖救，直接合并分析
combined<-cbind(exp1,exp2)

#創(chuàng)建一個文件夾用于寫分析結(jié)果
result.name <- "cca_monocle_result"
if(!dir.exists(result.name)){
  dir.create(result.name)
}

#合并的樣本分析看批次效應
combined.seurat<- CreateSeuratObject(
  combined,
  project = "combined", 
  min.cells = 5,
  min.features = 200,
  names.field = 4,
  names.delim = "_")
combined.seurat[["percent.mt"]] <- PercentageFeatureSet(combined.seurat,pattern = "^mt-")
combined.seurat<-NormalizeData(combined.seurat,verbose = FALSE)
combined.seurat<- FindVariableFeatures(combined.seurat, selection.method = "vst",nfeatures = 2000)
combined.seurat <- ScaleData(combined.seurat, verbose = FALSE)
combined.seurat <- RunPCA(combined.seurat, npcs = 50, verbose = FALSE)
pdf(paste0("./",result.name,"/PCA-DimPlot-noCCA.pdf"),width = 5,height = 4)
DimPlot(combined.seurat, reduction = "pca", )
dev.off()

接下來用CCA校正洪规，需要校正幾個就創(chuàng)建幾個seurat對象

#創(chuàng)建Seurat對象，標準化及尋找高表達變異基因-exp1
exp1.seurat<- CreateSeuratObject(
  exp1,
  project = "exp1", 
  min.cells = 5,
  min.features = 200,
  names.field = 4,
  names.delim = "_")
exp1.seurat[["percent.mt"]] <- PercentageFeatureSet(exp1.seurat,pattern = "^mt-")
exp1.seurat<-NormalizeData(exp1.seurat,verbose = FALSE)
exp1.seurat<- FindVariableFeatures(exp1.seurat, selection.method = "vst",nfeatures = 2000)

#創(chuàng)建Seurat對象循捺，標準化尋找高表達變異基因-exp2
exp2.seurat<- CreateSeuratObject(
  exp2,
  project = "exp2", 
  min.cells = 5,
  min.features = 200,
  names.field = 4,
  names.delim = "_")
exp2.seurat[["percent.mt"]] <- PercentageFeatureSet(exp2.seurat,pattern = "^mt-")
exp2.seurat<-NormalizeData(exp2.seurat,verbose = FALSE)
exp2.seurat<- FindVariableFeatures(exp2.seurat, selection.method = "vst",nfeatures = 2000)

#整合2批數(shù)據(jù)(k.filter默認為200斩例，一定要小于最小細胞數(shù)，不然會報錯)
exp.anchors<- FindIntegrationAnchors(object.list = c(exp1.seurat,exp2.seurat), dims = 1:30)
exp.seurat <- IntegrateData(anchorset = exp.anchors, dims = 1:30)

#設置當前使用的Assay(不運行也行从橘，整合完以后默認使用整合后的Assay)
DefaultAssay(exp.seurat)<-"integrated"

#均一化
exp.seurat <- ScaleData(exp.seurat, verbose = FALSE)

#運行PCA
exp.seurat <- RunPCA(exp.seurat, npcs = 50, verbose = FALSE)

#PCA結(jié)果展示
pdf(paste0("./",result.name,"/PCA-DimPlot-CCA.pdf"),width = 5,height = 4)
DimPlot(exp.seurat, reduction = "pca")
dev.off()

#PCA結(jié)果展示 - 碎石圖
pdf(paste0("./",result.name,"/PCA-ElbowPlot.pdf"),width = 5,height = 4)
ElbowPlot(exp.seurat, ndim=50,reduction="pca")
dev.off()

#增加樣本以及批次的信息到metadata中
cell.name<-row.names(exp.seurat@meta.data)
cell.name<-strsplit(cell.name,"_")
sample<-c()
batch=c()
for (i in 1:length(cell.name)){sample<-c(sample,substring(cell.name[[i]][2],1,2));batch<-c(batch,cell.name[[i]][4])}
exp.seurat@meta.data[["sample"]]<-sample
exp.seurat@meta.data[["batch"]]<-batch
rm(cell.name,sample,batch,i)
View(exp.seurat@meta.data)

#確定用于細胞分群的PC
dim.use <- 1:20

#細胞分群TSNE算法
exp.seurat <- FindNeighbors(exp.seurat, dims = dim.use)
exp.seurat <- FindClusters(exp.seurat, resolution = 0.5)
exp.seurat <- RunTSNE(exp.seurat, dims = dim.use, do.fast = TRUE)

#展示TSNE分類結(jié)果
pdf(paste0("./",result.name,"/CellCluster-TSNEPlot_",max(dim.use),"PC.pdf"),width = 5,height = 4)
DimPlot(object = exp.seurat, pt.size=1,label = T)
dev.off()

#按照數(shù)據(jù)來源分組展示細胞異同-batch
pdf(paste0("./",result.name,"/CellCluster-TSNEPlot_batch_",max(dim.use),"PC.pdf"),width = 5,height = 4)
DimPlot(object = exp.seurat, group.by="batch", pt.size=1)
dev.off()

#按照數(shù)據(jù)來源分組展示細胞異同-sample
pdf(paste0("./",result.name,"/CellCluster-TSNEPlot_sample_",max(dim.use),"PC.pdf"),width = 5,height = 4)
DimPlot(object = exp.seurat, group.by="sample", pt.size=1)
dev.off()

#plot在同一張圖中
plot1 <- DimPlot(object = exp.seurat, pt.size=1,label = T)
plot2 <- DimPlot(object = exp.seurat, group.by="batch", pt.size=1)
plot3 <- DimPlot(object = exp.seurat, group.by="sample", pt.size=1)
pdf(paste0("./",result.name,"/Combined_plot_exp.seurat_",max(dim.use),"PC.pdf"),width = 14,height = 4)
CombinePlots(plots = list(plot1,plot2,plot3),ncol=3)
dev.off()
rm(plot1,plot2,plot3)

#計算marker基因
all.markers<-FindAllMarkers(exp.seurat, only.pos = TRUE, 
                            min.pct = 0.25, logfc.threshold = 0.25)

# 遍歷每一個cluster然后展示其中前4個基因
marker.sig <- all.markers %>% filter(p_val_adj <= 0.05)
for(cluster in unique(marker.sig$cluster)){
  cluster.markers <- FindMarkers(exp.seurat, ident.1 = cluster, min.pct = 0.25)
  cluster.markers <- as.data.frame(cluster.markers) %>% 
    mutate(Gene = rownames(cluster.markers))
  cl4.genes <- cluster.markers %>% arrange(desc(avg_logFC))
  cl4.genes <- cl4.genes[1:min(nrow(cl4.genes),4),"Gene"]
  
  #RidgePlot
  pvn<-RidgePlot(exp.seurat, features = cl4.genes, ncol = 2)
  pdf(paste0("./",result.name,"/MarkerGene-RidgePlot_cluster",cluster,"_tsne_",max(dim.use),"PC.pdf"),width = 7,height = 6)
  print(pvn)
  dev.off()
  
  #VlnPlot
  pvn <- VlnPlot(exp.seurat, features = cl4.genes)
  pdf(paste0("./",result.name,"/MarkerGene-VlnPlot_cluster",cluster,"_tsne_",max(dim.use),"PC.pdf"),width = 7,height = 6)
  print(pvn)
  dev.off()
  
  #Feather plot 
  pvn <- FeaturePlot(exp.seurat,features=cl4.genes)
  pdf(paste0("./",result.name,"/MarkerGene-FeaturePlot_cluster",cluster,"_tsne_",max(dim.use),"PC.pdf"),width = 7,height = 6)
  print(pvn)
  dev.off()
  
}
rm(cl4.genes,cluster.markers,pvn)

#熱圖展示Top marker基因
top10 <- marker.sig %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)

pdf(paste0("./",result.name,"/MarkerGene-Heatmap_all_cluster_tsne_",max(dim.use),"PC.pdf"),width = 8,height = 5)
DoHeatmap(exp.seurat, features = top10$gene,size = 2)
dev.off()
write.table(top10,file=paste0("./",result.name,"_top10_marker_genes_tsne_",max(dim.use),"PC.txt"),sep="\t",quote = F)

pdf(paste0("./",result.name,"/MarkerGene-DotPlot_all_cluster_tsne_",max(dim.use),"PC.pdf"),width = 100,height = 10)
DotPlot(exp.seurat, features = unique(top10$gene))+RotatedAxis()
dev.off()