-
Notifications
You must be signed in to change notification settings - Fork 0
/
00.cut_zonation_layer_and_pathway_module_score.R
19 lines (18 loc) · 6.35 KB
/
00.cut_zonation_layer_and_pathway_module_score.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
pathways<-list("acylglycerol metabolic process"=c("Apoa1","Apof","Cdk8","Cps1","Insig1","Ldlr","Pcsk9","Plb1","Plce1"),"ATP metabolic process"=c("Afg1l","Ak4","Aldob","Ampd3","Apoc3","Atp5a1","Atp5b","Atp5c1","Atp5e","Atp5g1","Atp5g3","Atp5h","Atp5j","Atp5j2","Atp5k","Atp5l","Cfh","Cox4i1","Cox5a","Cox5b","Cox7a2","Cox7a2l","Cox7c","Cyc1","Cycs","Ddit4","Dnm1l","Fbp1","Gm10358","Gm3839","Hk1","Hspa8","Igf1","Khk","Ndufa8","Ndufb6","Ndufb8","Ndufb9","Ndufc2","Ndufs2","Ndufs6","Ndufv2","Park7","Pklr","Prkag2","Sdhd","Slc25a25","Tkfc","Uqcr10","Uqcrfs1","Uqcrh"),"alpha-amino acid metabolic process"=c("Aass","Acmsd","Agxt","Amdhd1","Arg1","Asl","Aspg","Ass1","Cbs","Cth","Ftcd","Gcsh","Gldc","Gls2","Gmps","Gnmt","Got1","Hal","Hnf4a","Kyat1","Kyat3","Kynu","Mat1a","Mccc2","Nox4","Otc","Park7","Ppat","Pycr2","Qdpr","Sdsl","Sephs2","Slc7a7","Tat","Uroc1"),"electron transport chain"=c("Afg1l","Cox4i1","Cox5a","Cox5b","Cox7c","Cyb561","Cyc1","Cycs","Ndufa5","Ndufa8","Ndufb6","Ndufb8","Ndufb9","Ndufc2","Ndufs2","Ndufs6","Ndufv2","Park7","Sdhb","Sdhd","Slc25a12","Uqcr10","Uqcrfs1","Uqcrh"),"generation of precursor metabolites and energy"=c("Aco2","Chchd4","Cox7a1","Dlst","Eno1b","Etfrf1","Fh1","G6pc","G6pdx","Gpi1","Grb10","Idh3g","Lepr","Mdh1","Mybbp1a","Ndufs1","Nfatc4","Ogdh","Oxct1","Per2","Ppargc1a","Ppif","Sdha","Slc25a22","Slc37a2","Suclg1","Tpi1","Uqcrc1"),"organic anion transport"=c("Abat","Abcb1a","Ace","Apoa4","Apoc2","Arg2","Cyp4f18","G6pc","Gipc1","Kmo","Mfsd2a","Nfkbie","P2ry2","Per2","Plin2","Plscr1","Prelid2","Scp2","Slc16a5","Slc25a22","Slc26a1","Slc2a2","Slc37a2","Slc38a2","Slc51b","Slc6a8","Slc7a4","Slc9a3r1"),"cellular ketone metabolic process"=c("Apoa4","Apoc2","Cyp21a1","Cyp4f18","Fdxr","Fh1","Kmo","Mfsd2a","Mlycd","Oxct1","Pdss1","Ppargc1a","Scp2","Slc37a2","Srd5a1","Tdo2"),"arachidonic acid metabolic process"=c("Cyp2a22","Cyp2a4","Cyp2a5","Cyp2b10","Cyp2c37","Cyp2c38","Cyp2c40","Cyp2c69","Cyp2d12","Cyp2d22","Cyp2j6","Cyp4a12a","Cyp4a12b"),"fatty acid metabolic process"=c("Asah2","Cpt2","Cyp2a22","Cyp2a4","Cyp2a5","Cyp2b10","Cyp2c37","Cyp2c38","Cyp2c40","Cyp2c69","Cyp2d12","Cyp2d22","Cyp2j6","Cyp4a12a","Cyp4a12b","Etfa","Lipg","Pdk4","Pla2g10","Pparg","Slc25a17","Slc27a1","Them4","Trib3","Tysnd1"),"anion transmembrane transport"=c("Aacs","Abcd3","Abcd4","Abhd3","Acaa1a","Acaa1b","Acaa2","Acat1","Acat2","Acot1","Acot12","Acot2","Acot3","Acot4","Acot6","Acot8","Acox1","Acsf2","Acsl1","Acsm1","Acsm3","Acsm5","Acss2","Adtrp","Akr1c14","Akr1c20","Akr1c6","Aldh3a2","Avpr1a","Cyb5a","Cyp1a2","Cyp2c23","Cyp2c29","Cyp2c39","Cyp2c50","Cyp2c54","Cyp2c55","Cyp2c67","Cyp2c68","Cyp2d10","Cyp2d9","Cyp2e1","Cyp2g1","Cyp4a10","Cyp4a14","Cyp4a32","Cyp4f14","Cyp4f15","Decr1","Decr2","Dgat2","Ech1","Eci1","Eci2","Elovl1","Elovl3","Ephx2","Fabp1","Fabp2","Gcdh","Gstm4","Hacd1","Hacd2","Hacd3","Hacl1","Hadh","Hsd17b4","Lias","Lonp2","Mlxipl","Pex7","Phyh","Por","Rgn","Slc27a2","Slc27a5","Tmem189"),"organic acid catabolic process"=c("Abcd3","Abcd4","Abhd3","Acaa1a","Acaa1b","Acaa2","Acat1","Acat2","Acot2","Acot4","Acot8","Acox1","Adtrp","Agxt2","Akr1a1","Aldh4a1","Bckdhb","Blmh","Csad","Cyp26a1","Cyp26c1","Cyp4f14","Cyp4f15","Decr1","Eci1","Eci2","Fabp1","Fah","Gcdh","Glud1","Gstz1","Hacl1","Hadh","Hgd","Hibadh","Hmgcl","Hpd","Hsd17b4","Hyal1","Ldhd","Lonp2","Mtrr","Pex7","Phyh","Pon1","Pon3","Prodh","Shmt1","Slc27a2"),"sulfur compound metabolic process"=c("Abcc2","Acaa2","Acat1","Acot1","Acot12","Acot2","Acot3","Acot4","Acot6","Acot8","Acsl1","Acsm1","Acsm3","Acsm5","Acss2","Blmh","Comt","Cs","Csad","Dgat2","Dlat","Enpp1","Gcdh","Gclm","Ghr","Gsta3","Gstm1","Gstm3","Gstm4","Gstm6","Gstt1","Gstz1","Hmgcl","Hsd17b4","Idh1","Lias","Mat2a","Mgst1","Mtrr","Mvk","Nat8","Nfe2l2","Papss2","Pmvk","Sod1","Stat5a","Sult1b1","Sult1d1","Sult1e1","Ugdh"),"alcohol metabolic process"=c("Acer3","Adh1","Akr1a1","Akr1c14","Akr1c20","Akr1c6","Akr7a5","Aldh1a1","Aldh3a2","Apoc1","Apoe","Clcn2","Coq3","Cyp1a2","Cyp27a1","Dgat2","Dpm1","Ebp","Ephx2","Fdx1","Fech","Gba2","H6pd","Idh1","Lcat","Lmf1","Lpcat3","Lrp5","Mvk","Npc1","Nsdhl","P2ry1","Pctp","Plcb1","Pmvk","Pon1","Por","Prkg1","Rbp4","Rdh10","Soat2","Sod1","Sord","Srd5a3","Sult1b1","Sult1e1","Tm7sf2","Ttc39b"),"lipid homeostasis"=c("Abcd1","Abcg8","Acox2","Apoa2","Apoc4","Commd1","Cyp7a1","Gck","Lipc","Mia2","Nr1d1","Nr5a2"),"purine-containing compound metabolic process"=c("Abcd1","Acnat1","Acot7","Acsl5","Adk","Crot","Gamt","Gck","Mif","Mpc1","Mpc2","Nudt2","Pemt","Pipox","Ppara","Prps1l3","Slc25a13","Ttr"),"ER stress"=c("Dnajb11","Dnajb9","Hspa8","Hspa9","Xbp1","Manf","Atf4","Dnajc10"))
library(Seurat)
SeuObj<-readRDS("D0.rds") # Read the Seurat object which you want to deal with
SeuObj<-AddModuleScore(SeuObj,features=pathways) # Calculate the module score of each pathway
colnames(SeuObj@meta.data)=c(colnames(SeuObj@meta.data)[1:31],names(pathways))
library(ggplot2)
dir.create("pathways");
apply(as.matrix(colnames(SeuObj@meta.data)[32:46]),1,function(x){ # apply function runs faster than for loop
df=SeuObj@meta.data[SeuObj@meta.data[,x]<quantile(SeuObj@meta.data[,x],0.995)[[1]] & SeuObj@meta.data[,x]>quantile(SeuObj@meta.data[,x],0.005)[[1]],];
p<-ggplot(df,aes(x=coor_x,y=coor_y,fill=SeuObj@meta.data[,x][SeuObj@meta.data[,x]<quantile(SeuObj@meta.data[,x],0.995)[[1]] & SeuObj@meta.data[,x]>quantile(SeuObj@meta.data[,x],0.005)[[1]]]))+geom_tile()+coord_fixed()+theme_void()+scale_fill_viridis_c(option="B",direction=-1)+labs(fill="");
ggsave(filename=paste0("pathways/",x,".pdf"),plot=p)
})
CV=c("Alb","Aldh1b1","Aldob","Apoa4","Apoa5","Apoc2","Arg1","Asl","Ass1","C9","Cdh1","Cyp2f2","Etnppl","Fbp1","G6pc","Gls2","Gnmt","Hal","Hp","Hpx","Hsd17b6","Hsd17b13","Mfsd2a","Mup20","Pck1","Pigr","Rida","Sds","Selenbp2","Selenop","Serpina1c","Serpina1e","Tdo2","Trf","Ugt2b38","Uox")
PV=c("Akr1c6","Aldh1a1","Aldh3a2","Car3","Ces1c","Csad","Cyb5a","Cyp1a2","Cyp2a5","Cyp2c29","Cyp2c37","Cyp2c40","Cyp2c50","Cyp2c54","Cyp2e1","Cyp3a11","Cyp4a10","Cyp4a14","Glul","Gsta3","Gstm1","Gulo","Lect2","Mgst1","Mup11","Mup15","Mup16","Mup17","Mup18","Oat","Pon1","Rgn","Rnase4","Slc1a2","Slc22a1","Slco1b2","Sord","Ugt1a10","Ugt2b1")
SeuObj=AddModuleScore( SeuObj,features = list("CV"=sample(CV,size = x)),name = "CV")
SeuObj=AddModuleScore( SeuObj,features = list("PV"=sample(PV,size = x)),name = "PV")
SeuObj$rank=cut_number(SeuObj$PV1-SeuObj$CV1,n = 9,label=1:9) # cut_number function can divide numbers into even groups by their values.
saveRDS(SeuObj,"add.Pathway.layer.rds")