getwd()path="G:/silicosis/geo/GSE104154_scRNA-seq_fibrotic MC_bleomycin"#空间转录组dir.create(path)setwd(path)getwd()list.files()raw_counts=read.csv("G:/silicosis/geo/GSE104154_scRNA-seq_fibrotic MC_bleomycin/GSE104154_d0_d21_sma_tm_Expr_raw/GSE104154_d0_d21_sma_tm_Expr_raw.csv")head(raw_counts)[1:4,1:4]counts=raw_counts[,-1]head(counts)[1:4,1:4]rownames(counts)=counts$symbolhead(raw_counts)[1:4,1:4]counts=raw_counts[,-2]head(counts)[1:4,1:4]rownames(counts)=counts$idcounts=counts[,-1]library(Seurat)#https://zhuanlan.zhihu.com/p/385206713rawdata=CreateSeuratObject(counts=counts,project="blem",assay="RNA")ids=raw_counts[,1:2]head(ids)colnames(ids)=c('ENSEMBL','SYMBOL')head(ids)dim(ids)#[1]16428ids=na.omit(ids)dim(ids)#[1]15504length(unique(ids$SYMBOL))#[1]15494#这里的关系超级乱，互相之间都不是一对一#凡是混乱的ID一律删除即可ids=ids[!duplicated(ids$SYMBOL),]ids=ids[!duplicated(ids$ENSEMBL),]dim(ids)pos=match(ids$ENSEMBL,rownames(rawdata))hp_sce=rawdata[pos,]hp_sce#rownames(hp_sce)=ids$SYMBOL#RenameGenesSeurat-----------------------------------------------#创建函数改名字RenameGenesSeurat<-function(obj,newnames){#ReplacegenenamesindifferentslotsofaSeuratobject.Runthisbeforeintegration.Runthisbeforeintegration.#Itonlychangesobj@assays$RNA@counts,@dataand@scale.data.print("Run this before integration. It only changes obj@assays$RNA@counts, @data and @scale.data.")RNA<-obj@assays$RNAif(nrow(RNA)==length(newnames)){if(length(RNA@counts))RNA@counts@Dimnames[[1]]<-newnamesif(length(RNA@data))RNA@data@Dimnames[[1]]<-newnamesif(length(RNA@scale.data))RNA@scale.data@Dimnames[[1]]<-newnames}else{"Unequal gene sets: nrow(RNA) != nrow(newnames)"}obj@assays$RNA<-RNAreturn(obj)}hp_sce=RenameGenesSeurat(obj=hp_sce,newnames=ids$SYMBOL)getwd()#save(hp_sce,file='first_sce.Rdata')hp_scerownames(hp_sce)[grepl('^mt-',rownames(hp_sce))]rownames(hp_sce)[grepl('^Rp[sl]',rownames(hp_sce))]hp_sce[["percent.mt"]]<-PercentageFeatureSet(hp_sce,pattern="^mt-")fivenum(hp_sce[["percent.mt"]][,1])rb.genes<-rownames(hp_sce)[grep("^Rp[sl]",rownames(hp_sce))]C<-GetAssayData(object=hp_sce,slot="counts")percent.ribo<-Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100hp_sce<-AddMetaData(hp_sce,percent.ribo,col.name="percent.ribo")getwd()plot1<-FeatureScatter(hp_sce,feature1="nCount_RNA",feature2="percent.mt")plot2<-FeatureScatter(hp_sce,feature1="nCount_RNA",feature2="nFeature_RNA")CombinePlots(plots=list(plot1,plot2))VlnPlot(hp_sce,features=c("percent.ribo","percent.mt"),ncol=2)VlnPlot(hp_sce,features=c("nFeature_RNA","nCount_RNA"),ncol=2)VlnPlot(hp_sce,features=c("percent.ribo","nCount_RNA"),ncol=2)hp_scehp_sce1<-subset(hp_sce,subset=nFeature_RNA>200&nCount_RNA>1000&percent.mt<20)hp_sce1sce=hp_sce1scecolnames(sce)grep(colnames(sce),pattern=".1")grep(colnames(sce),pattern=".2")sce@meta.data$stim<-c(rep("PBS",length(grep("1$",sce@assays$RNA@counts@Dimnames[[2]]))),rep("PBS",length(grep("2$",sce@assays$RNA@counts@Dimnames[[2]]))),rep("PBS",length(grep("3$",sce@assays$RNA@counts@Dimnames[[2]]))),rep("Bleomycin",length(grep("4$",sce@assays$RNA@counts@Dimnames[[2]]))),rep("Bleomycin",length(grep("5$",sce@assays$RNA@counts@Dimnames[[2]]))),rep("Bleomycin",length(grep("6$",sce@assays$RNA@counts@Dimnames[[2]]))))##8186,7947;table(sce$stim)library(dplyr)sce[["RNA"]]@meta.features<-data.frame(row.names=rownames(sce[["RNA"]]))All=sce%>%Seurat::NormalizeData(verbose=FALSE)%>%FindVariableFeatures(selection.method="vst",nfeatures=2000)%>%ScaleData(verbose=FALSE)All=RunPCA(All,npcs=50,verbose=FALSE)pdf("2_ElbowPlot.pdf")ElbowPlot(All,ndims=50)dev.off()library(cowplot)#All@meta.data$stim<-c(rep("case",length(grep("1$",All@assays$RNA@counts@Dimnames[[2]]))),rep("ctrl",length(grep("2$",All@assays$RNA@counts@Dimnames[[2]]))))##8186,7947;pdf("2_pre_harmony_harmony_plot.pdf")options(repr.plot.height=5,repr.plot.width=12)p1<-DimPlot(object=All,reduction="pca",pt.size=.1,group.by="stim")p2<-VlnPlot(object=All,features="PC_1",group.by="stim",pt.size=.1)plot_grid(p1,p2)dev.off()##########################runharmonyAll<-All%>%RunHarmony("stim",plot_convergence=TRUE)harmony_embeddings<-Embeddings(All,'harmony')pdf("2_after_harmony_harmony_plot.pdf")options(repr.plot.height=5,repr.plot.width=12)p3<-DimPlot(object=All,reduction="harmony",pt.size=.1,group.by="stim")p4<-VlnPlot(object=All,features="harmony_1",group.by="stim",pt.size=.1)plot_grid(p3,p4)dev.off()#############cluster#library(harmony)All<-All%>%RunUMAP(reduction="harmony",dims=1:30)%>%RunTSNE(reduction="harmony",dims=1:30)%>%FindNeighbors(reduction="harmony",dims=1:30)All<-All%>%FindClusters(resolution=3)%>%identity()options(repr.plot.height=4,repr.plot.width=10)pdf("3_after_harmony_umap_two_group.pdf")DimPlot(All,reduction="umap",group.by="stim",pt.size=.1)dev.off()pdf("3_after_harmony_cluster_UMAP.pdf")DimPlot(All,reduction="umap",label=TRUE,pt.size=.1)dev.off()pdf("3_umap_samples_split.pdf")DimPlot(All,reduction="umap",pt.size=.1,split.by="stim",label=T)dev.off()pdf("3_after_harmony_tsne_two_group.pdf")DimPlot(All,reduction="tsne",group.by="stim",pt.size=.1)dev.off()pdf("3_after_harmony_cluster_tSNE.pdf")DimPlot(All,reduction="tsne",label=TRUE,pt.size=.1)dev.off()pdf("3_tSNE_samples_split.pdf")DimPlot(All,reduction="tsne",pt.size=.1,split.by="stim",label=T)dev.off()getwd()#save(All,file="G:/silicosis/geo/GSE104154_scRNA-seq_fibrotic MC_bleomycin/All_for_clustering.rds")load("G:/silicosis/geo/GSE104154_scRNA-seq_fibrotic MC_bleomycin/All_for_clustering.rds")