a1<-intersect(rownames(GSE20680_anno),rownames(GSE42148_anno)) a2<-intersect(rownames(GSE25724_anno),a1) GSE20966_anno<-GSE20966_anno[a2,] GSE38642_anno<-GSE38642_anno[a2,] GSE25724_anno<-GSE25724_anno[a2,] target<-c(rep("GSE20966",20),rep("GSE38642",63),rep("GSE25724",13)) heji<-cbind(GSE20966_anno,GSE38642_anno,GSE25724_anno) expr <- heji boxplot(expr, las = 2, outline = F, col = as.factor(target)) library(factoextra) library(RColorBrewer) boxplot(expr, las = 2, outline = F, col = as.factor(GSE2_targets$group)) #PCA plot eset1<-eset[int,] eset1<-as.data.frame(t(GSE2_combat)) esetzong<-cbind(eset1,pdata) combot<-as.data.frame(t(fenxing)) iris<-cbind(combot,aim$group) library(factoextra) library(ggplot2) library("FactoMineR") library("factoextra") identical(colnames(eset),rownames(pdata)) fpkmToTpm <- function(fpkm) { exp(log(fpkm) - log(sum(fpkm)) + log(1e6)) } tpms <- apply(combot,2,fpkmToTpm) eset1<-as.data.frame(t(eset1)) esetzong<-cbind(eset1,pdata) combot<-normalize(combot) res.pca <- prcomp(eset1, scale = TRUE) fviz_screeplot(res.pca, addlabels = TRUE) pdf("C:/Users/admin/Desktop/guanxinbing/PCApca1.pdf",width = 8,height = 8, onefile=FALSE) fviz_pca_ind(res.pca, label="none", habillage= target, addEllipses=T, ellipse.level=0.95, palette = c("#91612D", "#E69F00", "#56B4E9")) dev.off() ## PCA PCA_new(expr, ntop = nrow(expr), group = target) pdf("C:/Users/admin/Desktop/guanxinbing/pca.pdf", height = 10, width = 10) dev.off() ## hclust dim(expr) colnames(expr) <- paste(target, 1:96, sep = "_") dist <- dist(t(expr)) hc = hclust(dist) library(factoextra) library(RColorBrewer) pdf("C:/Users/admin/Desktop/guanxinbing/clust4.pdf", height = 10, width = 30) fviz_dend(hc, k = 4, cex = 0.5, k_colors = brewer.pal(4, "Dark2"), color_labels_by_k = TRUE, ggtheme = theme_classic()) dev.off() library(dendextend) #### 6 merge batch <- GSE2_targets$GSE_num BiocManager::install("sva") library(sva) load("D:/Program Files//RawdatSEzong.rdata") mod <- model.matrix(~ as.factor(group), data= GSE2_targets) GSE2_combat <- ComBat(dat = expr, batch = target, # mod = mod, par.prior = TRUE, prior.plots = FALSE) GSE2_combat1 <- ComBat(dat = zonggene, batch = target, # mod = mod, par.prior = TRUE, prior.plots = FALSE) #### 7 ---- ## boxplot boxplot(GSE2_combat, las = 2, outline = F, col = as.factor(target)) ## PCA pdf("C:/Users/admin/Desktop/guanxinbing/box hou3.pdf", height = 10, width = 10) PCA_new(GSE2_combat, ntop = nrow(GSE2_combat), group = target) dev.off() ## hclust expr <- GSE2_combat colnames(expr) <- paste(target, 1:476, sep = "_") dist <- dist(t(expr)) hc = hclust(dist) library(factoextra) library(RColorBrewer) pdf("D:/Program Files/cafs/clust4.pdf", height = 10, width = 30) fviz_dend(hc, k = 6, cex = 0.5, k_colors = brewer.pal(4, "Set2"), color_labels_by_k = TRUE, ggtheme = theme_classic()) dev.off() # library(limma) library(tidyverse) aw<-GSE2_combat aw<-as.data.frame(aw) group <- factor(pdatazongt$group, levels = c("yes","non")) design <- model.matrix( ~ 0 + group) colnames(design) <- levels(group) design contrast.matrix <- makeContrasts("yes-non", levels = design) contrast.matrix expr expr <- aw fit <- lmFit(expr, design) fit <- contrasts.fit(fit, contrast.matrix) fit <- eBayes(fit) GSE89632_DEG <- topTable(fit, coef = 1, n = Inf) %>% rownames_to_column(var = "symbol") GSE89632_DEG_p <- GSE89632_DEG %>% dplyr::filter( abs(logFC) > 1,P.Value < 0.05) # #### volancae ---- library(tidyverse) library(ggrepel) GSE89632_DEG<-GSE89632_DEG[!duplicated(GSE89632_DEG$symbol), ] DEG_all <- GSE89632_DEG %>% dplyr::select(symbol, logFC, Pvalue = adj.P.Val) %>% mutate(direction = factor(ifelse(Pvalue < 0.05 & abs(logFC) > 0.5, ifelse(logFC > 0.5, "Up", "Down"),"NS"), levels=c('Up','Down','NS'))) #### ---- pdf('C:/Users/admin/Desktop/guanxinbing/volcanoplot.pdf',onefile = F,width = 8,height = 6) print(p1) dev.off() p1<-ggplot(data = DEG_all, aes(x = logFC, y = -log10(Pvalue), colour = direction)) + #数据映射 geom_point(alpha = 0.6) +# scale_color_manual(values = c("#DC143C", "#00008B", "#808080")) + # geom_text_repel(data = DEG_all %>% filter(Pvalue < 0.05, abs(logFC) > 1),# aes(label = symbol),# size = 3, segment.color = "black", # show.legend = FALSE) + theme_bw() +#设定主题 theme(legend.title = element_blank()) +# ylab(expression(-log[10]("Adjusted P Value"))) + xlab(expression(log[2]("Fold Change"))) + xlim(-2, 2) + ylim(0, 8) + geom_vline(xintercept = c(-1, 1), lty = 2, col = "black", lwd = 0.6) + geom_hline(yintercept = -log10(0.05), lty = 2, col = "black", lwd = 0.6) #WGCNA tangxiangguan<-intersect(GSE89632_DEG_p$symbol,azonggene) # GSE20680_pdata<-GSE20680_pdata[,c(34,8)] GSE20681_pdata<-GSE20681_pdata[,c(39,8)] GSE42148_pdata<-GSE42148_pdata[,c(36,38)] GSE12288_pdata<-GSE12288_pdata[,c(43,45)] colnames(GSE12288_pdata)<-c("group","tissue") pdataguan<-rbind(GSE20680_pdata,GSE20681_pdata,GSE42148_pdata,GSE12288_pdata) target<-c(rep("GSE20680",195),rep("GSE20681",198),rep("GSE42148",24),rep("GSE12288",222)) hejixin<-cbind(GSE20680_anno,GSE20681_anno,GSE42148_anno,GSE12288_anno) expr <- hejixin boxplot(expr, las = 2, outline = F, col = as.factor(target)) library(factoextra) library(RColorBrewer) boxplot(expr, las = 2, outline = F, col = as.factor(GSE2_targets$group)) #PCA2 eset1<-eset[int,] eset1<-as.data.frame(t(GSE2_combatxin )) esetzong<-cbind(eset1,pdata) combot<-as.data.frame(t(fenxing)) iris<-cbind(combot,aim$group) library(factoextra) library(ggplot2) library("FactoMineR") library("factoextra") identical(colnames(eset),rownames(pdata)) fpkmToTpm <- function(fpkm) { exp(log(fpkm) - log(sum(fpkm)) + log(1e6)) } tpms <- apply(combot,2,fpkmToTpm) eset1<-as.data.frame(t(eset1)) esetzong<-cbind(eset1,pdata) combot<-normalize(combot) res.pca <- prcomp(eset1, scale = TRUE) fviz_screeplot(res.pca, addlabels = TRUE) pdf("C:/Users/admin/Desktop/guanxinbing/pcaxinqian.pdf",width = 8,height = 8, onefile=FALSE) fviz_pca_ind(res.pca, label="none", habillage= target, addEllipses=T, ellipse.level=0.95, palette = c("#91612D", "#E69F00", "#56B4E9",'#029149')) dev.off() ## PCA PCA_new(expr, ntop = nrow(expr), group = target) pdf("C:/Users/admin/Desktop/guanxinbing/pca.pdf", height = 10, width = 10) dev.off() ## hclust dim(expr) colnames(expr) <- paste(target, 1:96, sep = "_") dist <- dist(t(expr)) hc = hclust(dist) library(factoextra) library(RColorBrewer) pdf("C:/Users/admin/Desktop/guanxinbing/clust4.pdf", height = 10, width = 30) fviz_dend(hc, k = 4, cex = 0.5, k_colors = brewer.pal(4, "Dark2"), color_labels_by_k = TRUE, ggtheme = theme_classic()) dev.off() library(dendextend) #### 6 ---- batch <- GSE2_targets$GSE_num BiocManager::install("sva") library(sva) load("D:/Program Files/RawdatSEzong.rdata") mod <- model.matrix(~ as.factor(group), data= GSE2_targets) GSE2_combatxin <- ComBat(dat = expr, batch = target, # mod = mod, par.prior = TRUE, prior.plots = FALSE) GSE2_combat1 <- ComBat(dat = zonggene, batch = target, # mod = mod, par.prior = TRUE, prior.plots = FALSE) #### 7 ---- ## boxplot boxplot(GSE2_combat, las = 2, outline = F, col = as.factor(target)) #kmeans library(ConsensusClusterPlus) library(limma) library(tidyverse) aw<-GSE2_combatxin group <- factor(pdatazongt$group, levels = c("yes","non")) design <- model.matrix( ~ 0 + group) colnames(design) <- levels(group) design contrast.matrix <- makeContrasts("yes-non", levels = design) contrast.matrix expr expr <- aw fit <- lmFit(expr, design) fit <- contrasts.fit(fit, contrast.matrix) fit <- eBayes(fit) GSE89632_DEG <- topTable(fit, coef = 1, n = Inf) %>% rownames_to_column(var = "symbol") GSE89632_DEG_p <- GSE89632_DEG %>% dplyr::filter( abs(logFC) > 1,P.Value < 0.05) library(tidyverse) library(ggrepel) GSE89632_DEG<-GSE89632_DEG[!duplicated(GSE89632_DEG$symbol), ] DEG_all <- GSE89632_DEG %>% dplyr::select(symbol, logFC, Pvalue = adj.P.Val) %>% mutate(direction = factor(ifelse(Pvalue < 0.05 & abs(logFC) > 0.5, ifelse(logFC > 0.5, "Up", "Down"),"NS"), levels=c('Up','Down','NS'))) fzong<-c(cxcl,ccl,CCR,CXCR,il,chkine,ifterion) GSE2_combatxin1<-as.data.frame((GSE2_combatxin)) GSE2_combatxin1<-GSE2_combatxin1[fzong,] GSE2_combatxin1<-na.omit(GSE2_combatxin1) #single cell analysis folders=c("G1") sceList = lapply(folders,function(folder){ CreateSeuratObject(counts = Read10X(folder), project = folder ) }) aw<-sceList[1] sced <- merge(x=sceList[[1]],y=Seurat_object) #merge rm(list = ls()[match(c("ScRNAdata","pd","monocle.matrix"),ls())]) sce[["percent.mt"]] <- PercentageFeatureSet(sce, pattern = "^MT-")# sce[["percent.Ribo"]] <- PercentageFeatureSet(sce, pattern = "^RP[SL]")# sce<- CreateSeuratObject(counts = Read10X("G1"), project = "G1" ) sce <- merge(x=Seurat_object,y=sceList[1:length(sceList)],add.cell.ids = folders) #merge raw_meta=sce@meta.data raw_count <- table(raw_meta$orig.ident) raw_count sum(raw_count)#10232 cols <- c("#a696c8","#F38181","#95E1D3", '#80B1D3', '#FDB462', '#B3DE69', '#FCCDE5', '#BC80BD', '#CCEBC5', '#FFED6F', '#E41A1C', '#377EB8', '#4DAF4A', '#984EA3', '#FF7F00', '#A65628','#F0027F') pearplot_befor<-VlnPlot(sce,group.by ='orig.ident', features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.Ribo"), pt.size = 0, ncol = 4) pearplot_befor #featurescatter Feature_ber1<-FeatureScatter(sce,feature1 = 'nFeature_RNA',feature2 = 'nCount_RNA',group.by = 'orig.ident') Feature_ber2<-FeatureScatter(sce,feature1 = 'percent.mt',feature2 = 'nCount_RNA',group.by = 'orig.ident' ) Feature_ber3<-FeatureScatter(sce,feature1 = 'percent.mt',feature2 = 'nFeature_RNA',group.by = 'orig.ident') Feature_ber1=Feature_ber1+theme(legend.position = 'none') Feature_ber2=Feature_ber2+theme(legend.position = 'none') Feature_ber<-ggarrange(Feature_ber1,Feature_ber2,Feature_ber3,ncol = 3,nrow = 1,widths = c(1,1,1.2)) Feature_ber sce=subset(x =sce, subset = nFeature_RNA > 100 & nFeature_RNA < 5000 & percent.mt < 25 & nCount_RNA > 100) #merge clean_meta=sce@meta.data clean_count <- table(clean_meta$orig.ident) clean_count sum(clean_count)#9292 pearplot_after <- VlnPlot(sce,group.by ='orig.ident', features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.Ribo"), pt.size = 0, ncol = 4) pearplot_after gc() #Normalizing the data sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000) #Identification of highly variable features (feature selection) sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000, mean.cutoff=c(0.0125,3), dispersion.cutoff =c(1.5,Inf)) #PCA sce <- ScaleData(sce, features = rownames(sce)) sce <- RunPCA(sce, features = VariableFeatures(sce)) sce <- RunTSNE(sce, features = VariableFeatures(sce)) sce <- RunUMAP(sce, features = VariableFeatures(sce)) dimplot1 <- DimPlot(sce, reduction = "pca",group.by = 'orig.ident') elbowplot1 <- ElbowPlot(sce, ndims=50, reduction="pca") sc_pca <- dimplot1+elbowplot1 figs1=ggarrange(Feature_ber,pearplot_befor,pearplot_after,sc_pca, nrow = 4,labels = c('A','B','C','D'),ncol = 1,widths = c(1,1,1,2)) ggsave(filename = 'FigS11.pdf',plot = figs1,he=12,width = 10) ggsave(filename = 'FigS1.jpg',plot = figs1,he=12,width = 10,dpi = 300) sce <- JackStraw(object = sce, num.replicate = 100) sce <- ScoreJackStraw(object = sce, dims = 1:15) pdf(file="02.pcaJackStraw.pdf",width=8,height=6) JackStrawPlot(object = sce, dims = 1:20) dev.off() Dims <- 15 sce <- RunTSNE(sce, dims=1:Dims, reduction="pca", perplexity=30, max_iter=1000) sce <- RunUMAP(sce, dims=1:Dims, reduction="pca") figs1=ggarrange(Feature_ber,pearplot_befor,pearplot_after,sc_pca, nrow = 4,labels = c('A','B','','C'),ncol = 1,widths = c(1,1,1,2)) ggsave(filename = 'FigS1.pdf',plot = figs1,he=12,width = 10) ggsave(filename = 'rFigS1.jpg',plot = figs1,he=12,width = 10,dpi = 300) library(clustree) sce <- FindNeighbors(sce, dims = 1:Dims) sce <- FindClusters( object = sce, resolution = c(seq(.1,1,.1)) ) colnames(sce@meta.data) clustree(sce@meta.data, prefix = "integrated_snn_res.") pdf('clust.snn_res.pdf',he=15,wi=15) clustree(sce@meta.data, prefix = "RNA_snn_res.") dev.off() colnames(sce@meta.data) Resolution <- 0.6 sce <- FindNeighbors(object = sce, dims = 1:15) sce <- FindClusters(object = sce, resolution = Resolution) #clustree(sce) DimPlot(sce,reduction = "tsne",label = T) DimPlot(sce,reduction = 'tsne',group.by = 'orig.ident') ###################################04.SingleR ################################### counts<-sce@assays$RNA@counts clusters<-sce@meta.data$seurat_clusters ann=sce@meta.data$orig.ident #ref=get(load("ref_Human_all.RData")) ref=celldex::HumanPrimaryCellAtlasData() singler=SingleR(test=counts, ref =ref, labels=ref$label.main, clusters = clusters) clusterAnn=as.data.frame(singler) clusterAnn=cbind(id=row.names(clusterAnn), clusterAnn) clusterAnn=clusterAnn[,c("id", "labels")] write.table(clusterAnn,file="04.clusterAnn.txt",quote=F,sep="\t", row.names=F) library(SingleR) #cluster newLabels=singler$labels names(newLabels)=levels(pbmc) pbmc=RenameIdents(pbmc, newLabels) pdf(file="04.TSNE.pdf",width=6.5,height=6) TSNEPlot(object = pbmc, pt.size = 2, label = TRUE) #TSNE可视化 dev.off() ##cluster pbmc.markers=FindAllMarkers(object = pbmc, only.pos = FALSE, min.pct = 0.25, logfc.threshold = logFCfilter) sig.cellMarkers=pbmc.markers[(abs(as.numeric(as.vector(pbmc.markers$avg_log2FC)))>logFCfilter & as.numeric(as.vector(pbmc.markers$p_val_adj))