Skip to content
Snippets Groups Projects
Commit 99423e5c authored by Gervaise Henry's avatar Gervaise Henry :cowboy:
Browse files

Update analysis, cutoff for nGenes increase due to deeper seq, don't recluster... 

Update analysis, cutoff for nGenes increase due to deeper seq, don't recluster subcluster, change genesets to not use readxl
parent 0f46f0cd
Branches
Tags
1 merge request!2Merge develop into master
Expand all Hide whitespace changes
Inline Side-by-side
genesets/DEG_C2.CGP.M10970.txt
+ 1
− 2
View file @ 99423e5c
CHUANG_OXIDATIVE_STRESS_RESPONSE_UP CHUANG_OXIDATIVE_STRESS_RESPONSE_UP
> Genes up-regulated in MCF7 cells (breast cancer) after treatment with the oxydants: hydrogen peroxyde, menadione, and t-butyl hydroperoxyde [PubChem=784;4055;6410].
CTAGE5 CTAGE5
CYP1B1 CYP1B1
DDIT3 DDIT3
...@@ -27,4 +26,4 @@ PSMD3 ...@@ -27,4 +26,4 @@ PSMD3
RCBTB1 RCBTB1
RPL38 RPL38
SLC35B3 SLC35B3
ZBTB4 ZBTB4
\ No newline at end of file
genesets/SupTab3_Consensus_Sigs.csv 0 → 100644
+ 288
− 0
View file @ 99423e5c
Supplementary Table 3,,,,,
"Consensus marker genes for cell-type clusters from both initial 3' droplet-based and full-length plate-based scRNA-seq data (7,193 + 301 cells, in Sup. Tabs 1 and 2, respectively)",,,,,
"Thresholds: Fisher's combined FDR: 0.001 for plate-based and maximum FDR: 0.05 for droplet-based, Minimum log2 fold-change:0.25",,,,,
,,,,,
,,,,,
,,,,,
Basal,Club,Ciliated,Tuft,Neuroendocrine,Ionocyte
Aqp3,Scgb1a1,AU040972,Lrmp,Chga,Stap1
Icam1,Bpifa1,Tmem212,Gnat3,Calca,Gm933
Krt17,Scgb3a2,Dynlrb2,Gnb3,Cxcl13,Ascl3
Krt5,Sftpd,Ccdc153,Plac8,Scg2,Foxi1
Krt15,Alox15,Tppp3,Trpm5,Nov,Atp6v0d2
Sfn,Pon1,1110017D15Rik,Gng13,Scg5,Moxd1
Perp,Cyp4a12b,Fam183b,Ltc4s,Pcsk1,Gnas
Fxyd3,Wfdc1,Sec14l3,Rgs13,Dnajc12,Ldhb
Sdc1,Mgst1,Pltp,Hck,Uchl1,Rasd1
Gstm2,Akr1c18,Rsph1,Alox5ap,Ddc,P2ry14
F3,Sftpa1,Mlf1,Avil,Snca,Tfcp2l1
Epas1,Lypd2,2410004P03Rik,Alox5,Snap25,Cd81
Abi3bp,Ces1f,Tekt1,Ptpn6,Cib3,Cftr
Adh7,Gsto1,Cdh26,Atp2a3,Bex2,Slc12a2
Id1,Tst,Cdhr3,Plk2,Guk1,Asgr1
Upk1b,Ffar4,Nme5,Zfp428,Chgb,Parm1
Hspa1a,Hp,Ccdc17,Ethe1,Ascl1,Atp6v1e1
Dapl1,Chad,Adam8,Mctp1,Ngf,Iqgap1
Hspb1,Selenbp1,Cyp2s1,Fxyd6,Cplx2,Rbpms
Zfp296,Mgat3,1700007K13Rik,1810046K07Rik,Igfbp5,
Dcn,Gabrp,Tm4sf1,Pik3r5,Pkib,
Rplp0,,Elof1,Cd300lf,Lrp11,
Fam107b,,Ly6a,Rassf6,Tcerg1l,
Tacstd2,,Sntn,Dclk1,Rbp4,
Rpl13,,Ppil6,Spib,Mthfd2,
Sat1,,Ifitm1,Pik3cg,Cd9,
Sult1d1,,Tctex1d4,Pde2a,Meis2,
Hpgd,,Aldh3b1,Matk,Ly6h,
S100a14,,Lrrc48,Vav1,Tmem158,
Rps18,,Foxj1,Sh2d7,Cacna2d1,
Emp1,,1700026L06Rik,Bmx,Syt7,
Sgk1,,Csrp2,Fyb,Ptp4a3,
Gsta4,,Pifo,Ccdc28b,,
Rps4x,,Acot1,Ptpn18,,
Rps5,,Dnali1,Bpgm,,
Oat,,Ccdc113,Itgb7,,
Rps8,,Pcp4l1,Ccdc109b,,
S100a16,,Spag6,Anxa4,,
Phlda3,,Capsl,Rgs2,,
Dusp1,,Lrrc51,Tspan6,,
Crlf1,,2610028H24Rik,Prss53,,
Sh3bgrl3,,Ldlrad1,Ly6g6f,,
Rps2,,Stmnd1,Inpp5d,,
Junb,,1700024G13Rik,Ivns1abp,,
Pcolce,,Dmkn,Nkd1,,
Cotl1,,Ccdc78,Sh3bgrl,,
Sord,,Chchd6,Pou2f3,,
Igfbp7,,1700016K19Rik,Sox9,,
Rpl8,,Calml4,Chn2,,
Rpl4,,Vpreb3,Krt18,,
,,Odf3b,Calm2,,
,,Ak7,Kctd12,,
,,1700003M02Rik,Etv1,,
,,1700026D08Rik,Coprs,,
,,Lrrc23,Cd47,,
,,BC051019,Espn,,
,,Tubb4b,Tmem245,,
,,Enkur,Pitpnc1,,
,,Bphl,Hap1,,
,,Spata18,Reep5,,
,,Fam216b,H2afj,,
,,Riiad1,Selm,,
,,1700088E04Rik,NA,,
,,1700007G11Rik,NA,,
,,Rsph4a,NA,,
,,Fam47e,NA,,
,,Agr3,NA,,
,,Fbxo36,NA,,
,,1700001C02Rik,NA,,
,,Cmbl,NA,,
,,Spag16,NA,,
,,Bok,NA,,
,,Ift43,NA,,
,,Ubxn10,NA,,
,,Tekt4,NA,,
,,Ccdc108,NA,,
,,Erich2,NA,,
,,6820408C15Rik,NA,,
,,Prdx5,NA,,
,,Fam92b,NA,,
,,Smim5,NA,,
,,Hdc,NA,,
,,Serpina9,NA,,
,,Ly6c1,NA,,
,,Aldh1a1,NA,,
,,Fhl1,NA,,
,,4833427G06Rik,NA,,
,,Rfk,NA,,
,,Dnaic2,NA,,
,,Ppp1r36,NA,,
,,Ak8,NA,,
,,B9d1,NA,,
,,Chchd10,NA,,
,,Zmynd10,NA,,
,,Rsph9,NA,,
,,Nek5,NA,,
,,Gm973,NA,,
,,Ccdc39,NA,,
,,Efhc1,NA,,
,,Hsph1,NA,,
,,Tuba1b,NA,,
,,Ubxn11,NA,,
,,Slc23a2,NA,,
,,Spa17,NA,,
,,Nat9,NA,,
,,Efcab10,NA,,
,,Cd177,NA,,
,,Hsp90aa1,NA,,
,,Iqca,NA,,
,,Dync2li1,NA,,
,,Ruvbl1,NA,,
,,Ppp1r32,NA,,
,,Spef1,NA,,
,,Ccdc81,NA,,
,,Dnajb13,NA,,
,,Acox2,NA,,
,,Fhad1,NA,,
,,Spag8,NA,,
,,Nudc,NA,,
,,Ccdc146,NA,,
,,Pih1d2,NA,,
,,Vwa3a,NA,,
,,Fam179a,NA,,
,,Fth1,NA,,
,,Ccdc151,NA,,
,,Lrrc36,NA,,
,,Cxcl17,NA,,
,,Kif9,NA,,
,,Ccdc37,NA,,
,,Mapk15,NA,,
,,Lhb,NA,,
,,Cd82,NA,,
,,Fam161a,NA,,
,,Rabl2,NA,,
,,Fank1,NA,,
,,Wdr63,NA,,
,,Glipr1,NA,,
,,Zfp474,NA,,
,,Lgals3bp,NA,,
,,Dpcd,NA,,
,,Ttc29,NA,,
,,Ttc16,NA,,
,,Mdh1b,NA,,
,,Ruvbl2,NA,,
,,Gipc2,NA,,
,,Arc,NA,,
,,Mak,NA,,
,,Dcdc2b,NA,,
,,Ccdc30,NA,,
,,1700028P14Rik,NA,,
,,Mycbp,NA,,
,,Oscp1,NA,,
,,Morn5,NA,,
,,Iqcg,NA,,
,,Eno4,NA,,
,,Ift74,NA,,
,,Lrrc34,NA,,
,,Cetn4,NA,,
,,Ankrd66,NA,,
,,Dnajc15,NA,,
,,Ulk4,NA,,
,,Tbata,NA,,
,,Traf3ip1,NA,,
,,Akr1b8,NA,,
,,6430531B16Rik,NA,,
,,Armc4,NA,,
,,Maats1,NA,,
,,Plekhb1,NA,,
,,Meig1,NA,,
,,Morn3,NA,,
,,Lrrc71,NA,,
,,1110004E09Rik,NA,,
,,Ift88,NA,,
,,Timp4,NA,,
,,Fam166b,NA,,
,,Ccdc121,NA,,
,,Usp18,NA,,
,,Ccdc33,NA,,
,,Slc23a1,NA,,
,,Gm166,NA,,
,,1810037I17Rik,NA,,
,,Syt5,NA,,
,,Dnaaf1,NA,,
,,Mns1,NA,,
,,Dusp14,NA,,
,,Wdr78,NA,,
,,Dynll1,NA,,
,,Gas8,NA,,
,,Ccdc114,NA,,
,,Crip2,NA,,
,,Ctxn1,NA,,
,,1700040L02Rik,NA,,
,,Tcp11,NA,,
,,Kif19a,NA,,
,,Ccdc65,NA,,
,,Tmem107,NA,,
,,Adgb,NA,,
,,Cib1,NA,,
,,Dnaic1,NA,,
,,Odc1,NA,,
,,Dnaaf3,NA,,
,,Cdkl4,NA,,
,,Mrps17,NA,,
,,Wdr34,NA,,
,,Slc9a3r1,NA,,
,,Kcnmb2,NA,,
,,Swi5,NA,,
,,Plin5,NA,,
,,Zc2hc1a,NA,,
,,Jam3,NA,,
,,Lrrc6,NA,,
,,Sptlc3,NA,,
,,Hydin,NA,,
,,Fam216a,NA,,
,,4932443I19Rik,NA,,
,,Ankrd42,NA,,
,,Wdr60,NA,,
,,Lrrc45,NA,,
,,Spef2,NA,,
,,Ezr,NA,,
,,Wnt11,NA,,
,,Fam213a,NA,,
,,Gtsf1l,NA,,
,,Lbp,NA,,
,,Wdr95,NA,,
,,Syne1,NA,,
,,Abhd3,NA,,
,,Tctex1d2,NA,,
,,Ttc39a,NA,,
,,Lrrc46,NA,,
,,Stk33,NA,,
,,Slc2a12,NA,,
,,Mipep,NA,,
,,Mst1,NA,,
,,Myl4,NA,,
,,Oaz2,NA,,
,,Ttc12,NA,,
,,Ccdc176,NA,,
,,Ift57,NA,,
,,Wdr54,NA,,
,,4930444P10Rik,NA,,
,,Tuba1a,NA,,
,,1700013F07Rik,NA,,
,,Mxra8,NA,,
,,D430036J16Rik,NA,,
,,Spag17,NA,,
,,Ddo,NA,,
,,Sri,NA,,
,,Ctsk,NA,,
,,Calm1,NA,,
,,4430402I18Rik,NA,,
,,2310007B03Rik,NA,,
,,Psmc3ip,NA,,
,,Nudt4,NA,,
,,Ttll9,NA,,
,,Dusp18,NA,,
,,Rfc2,NA,,
,,1700029J07Rik,NA,,
,,Naga,NA,,
,,Akap14,NA,,
,,Atp2a2,NA,,
,,Pttg1,NA,,
,,Casc1,NA,,
,,Ccdc60,NA,,
,,Dydc2,NA,,
,,Lrrc10b,NA,,
,,Mycbpap,NA,,
,,1700001L19Rik,NA,,
,,Dlec1,NA,,
,,Kif6,NA,,
,,Ropn1l,NA,,
,,Eml2,NA,,
,,Ech1,NA,,
,,Slc1a4,NA,,
,,Sftpc,NA,,
,,1700101E01Rik,NA,,
,,Efcab1,NA,,
,,Sirt3,NA,,
,,Pacrg,NA,,
,,Cdc14a,NA,,
,,Itln1,NA,,
genesets/SupTab6_Krt13_Hillock.csv 0 → 100644
+ 5327
− 0
View file @ 99423e5c
This diff is collapsed.
r.scripts/sc-TissueMapper.R
+ 74
− 56
View file @ 99423e5c
...@@ -491,14 +491,27 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){ ...@@ -491,14 +491,27 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){
number.clusters <- length(unique(sc10x@ident)) number.clusters <- length(unique(sc10x@ident))
labels <- paste0("Cluster_",as.vector(factor(sc10x@ident))) labels <- paste0("Cluster_",as.vector(factor(sc10x@ident)))
if ((ncol(sc10x@data)>ds & ds!=0)){
rnd <- sample(1:ncol(sc10x@data),ds) #if ((ncol(sc10x@data)>ds & ds!=0)){
data <- sc10x@data[,rnd] # rnd <- sample(1:ncol(sc10x@data),ds)
labels <- labels[rnd] # data <- sc10x@data[,rnd]
} else { # labels <- labels[rnd]
data <- sc10x@data #} else {
# data <- sc10x@data
#}
cell.sample <- NULL
for (i in 1:number.clusters){
cell <- names(sc10x@ident[sc10x@ident==i])
if (length(cell)>ds & ds!=0){
rnd <- sample(1:length(cell),ds)
cell <- cell[rnd]
}
cell.sample <- c(cell.sample,cell)
} }
data <- sc10x@data[,colnames(sc10x@data) %in% cell.sample]
data <- as.data.frame(as.matrix(data)) data <- as.data.frame(as.matrix(data))
labels <- labels[colnames(sc10x@data) %in% cell.sample]
#Make labels for QuSAGE #Make labels for QuSAGE
clust <- list() clust <- list()
...@@ -541,7 +554,7 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){ ...@@ -541,7 +554,7 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){
if (max(results.cor[results.cor[,4]==i,][,2])>=0){ if (max(results.cor[results.cor[,4]==i,][,2])>=0){
results.clust.id <- rbind(results.clust.id,results.cor[results.cor[,4]==i,][which.max(results.cor[results.cor[,4]==i,][,2]),]) results.clust.id <- rbind(results.clust.id,results.cor[results.cor[,4]==i,][which.max(results.cor[results.cor[,4]==i,][,2]),])
} else { } else {
results.clust.id <- rbind(results.clust.id,c("Unknown",0,0,i)) results.clust.id <- rbind(results.clust.id,data.frame(pathway.name="Unknown",log.fold.change=0,FDR=0,Cluster=i))
} }
} }
rownames(results.clust.id) <- NULL rownames(results.clust.id) <- NULL
...@@ -601,7 +614,7 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){ ...@@ -601,7 +614,7 @@ scQuSAGE <- function(sc10x,gs,res.use=0.1,ds=25000,nm="Lin",folder="lin"){
merge.cluster <- c(merge.cluster,"Unknown") merge.cluster <- c(merge.cluster,"Unknown")
}} }}
sc10x@ident <- plyr::mapvalues(x=sc10x@ident,from=1:number.clusters,to=merge.cluster) sc10x@ident <- plyr::mapvalues(x=sc10x@ident,from=1:number.clusters,to=merge.cluster)
sc10x@ident <- factor(sc10x@ident,levels=names(gs)) sc10x@ident <- factor(sc10x@ident,levels=c(names(gs),"Unknown"))
sc10x <- StashIdent(object=sc10x,save.name=nm) sc10x <- StashIdent(object=sc10x,save.name=nm)
postscript(paste0("./analysis/tSNE/",folder,"/tSNE_",nm,".eps")) postscript(paste0("./analysis/tSNE/",folder,"/tSNE_",nm,".eps"))
...@@ -638,14 +651,29 @@ scQuSAGEsm <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){ ...@@ -638,14 +651,29 @@ scQuSAGEsm <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){
clusters <- unique(sc10x@ident) clusters <- unique(sc10x@ident)
labels <- as.vector(factor(sc10x@ident)) labels <- as.vector(factor(sc10x@ident))
if ((ncol(sc10x@data)>ds & ds!=0)){
rnd <- sample(1:ncol(sc10x@data),ds) #if ((ncol(sc10x@data)>ds & ds!=0)){
data <- sc10x@data[,rnd] # rnd <- sample(1:ncol(sc10x@data),ds)
labels <- labels[rnd] # data <- sc10x@data[,rnd]
} else { # labels <- labels[rnd]
data <- sc10x@data #} else {
# data <- sc10x@data
#}
#data <- as.data.frame(as.matrix(data))
cell.sample <- NULL
for (i in clusters){
cell <- names(sc10x@ident[sc10x@ident==i])
if (length(cell)>ds & ds!=0){
rnd <- sample(1:length(cell),ds)
cell <- cell[rnd]
}
cell.sample <- c(cell.sample,cell)
} }
data <- sc10x@data[,colnames(sc10x@data) %in% cell.sample]
data <- as.data.frame(as.matrix(data)) data <- as.data.frame(as.matrix(data))
labels <- labels[colnames(sc10x@data) %in% cell.sample]
#Make labels for QuSAGE #Make labels for QuSAGE
clust <- list() clust <- list()
...@@ -744,14 +772,28 @@ scQuSAGElg <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){ ...@@ -744,14 +772,28 @@ scQuSAGElg <- function(sc10x,gs,ds=25000,nm="Lin",folder="lin"){
clusters <- unique(sc10x@ident) clusters <- unique(sc10x@ident)
labels <- as.vector(factor(sc10x@ident)) labels <- as.vector(factor(sc10x@ident))
if ((ncol(sc10x@data)>ds & ds!=0)){
rnd <- sample(1:ncol(sc10x@data),ds) #if ((ncol(sc10x@data)>ds & ds!=0)){
data <- sc10x@data[,rnd] # rnd <- sample(1:ncol(sc10x@data),ds)
labels <- labels[rnd] # data <- sc10x@data[,rnd]
} else { # labels <- labels[rnd]
data <- sc10x@data #} else {
# data <- sc10x@data
#}
#data <- as.data.frame(as.matrix(data))
cell.sample <- NULL
for (i in clusters){
cell <- names(sc10x@ident[sc10x@ident==i])
if (length(cell)>ds & ds!=0){
rnd <- sample(1:length(cell),ds)
cell <- cell[rnd]
}
cell.sample <- c(cell.sample,cell)
} }
data <- sc10x@data[,colnames(sc10x@data) %in% cell.sample]
data <- as.data.frame(as.matrix(data)) data <- as.data.frame(as.matrix(data))
labels <- labels[colnames(sc10x@data) %in% cell.sample]
#Make labels for QuSAGE #Make labels for QuSAGE
clust <- list() clust <- list()
...@@ -1071,37 +1113,13 @@ scCCA <- function(sc10x.1,sc10x.2,nm.1="D17",nm.2="D27",cc=FALSE){ ...@@ -1071,37 +1113,13 @@ scCCA <- function(sc10x.1,sc10x.2,nm.1="D17",nm.2="D27",cc=FALSE){
} }
sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=FALSE,cca=25,acca=15){ sc3CCA <- function(sc10x.1,sc10x.2,sc10x.3,nm.1="D17",nm.2="D27",nm.3="D35",cc=FALSE,cca=25,acca=15,lx=0.15,hx=3.5,ly=0.75){
gc() sc10x.1 <- FindVariableGenes(object=sc10x.1,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=lx,x.high.cutoff=hx,y.cutoff=ly,do.plot=FALSE)
if (cc==TRUE){ sc10x.2 <- FindVariableGenes(object=sc10x.2,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=lx,x.high.cutoff=hx,y.cutoff=ly,do.plot=FALSE)
sc10x.1 <- ScaleData(object=sc10x.1,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45) sc10x.3 <- FindVariableGenes(object=sc10x.3,mean.function=ExpMean,dispersion.function=LogVMR,x.low.cutoff=lx,x.high.cutoff=hx,y.cutoff=ly,do.plot=FALSE)
} else { genes.hvg.1 <- sc10x.1@var.genes
sc10x.1 <- ScaleData(object=sc10x.1,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45) genes.hvg.2 <- sc10x.2@var.genes
} genes.hvg.3 <- sc10x.3@var.genes
gc()
gc()
if (cc==TRUE){
sc10x.2 <- ScaleData(object=sc10x.2,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
} else {
sc10x.2 <- ScaleData(object=sc10x.2,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45)
}
gc()
gc()
if (cc==TRUE){
sc10x.3 <- ScaleData(object=sc10x.3,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
} else {
sc10x.3 <- ScaleData(object=sc10x.3,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45)
}
gc()
sc10x.1 <- FindVariableGenes(sc10x.1,do.plot=FALSE)
sc10x.2 <- FindVariableGenes(sc10x.2,do.plot=FALSE)
sc10x.3 <- FindVariableGenes(sc10x.3,do.plot=FALSE)
genes.hvg.1 <- head(rownames(sc10x.1@hvg.info),2000)
genes.hvg.2 <- head(rownames(sc10x.2@hvg.info),2000)
genes.hvg.3 <- head(rownames(sc10x.3@hvg.info),2000)
genes.hvg.Comb <- unique(c(genes.hvg.1,genes.hvg.2,genes.hvg.3)) genes.hvg.Comb <- unique(c(genes.hvg.1,genes.hvg.2,genes.hvg.3))
genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.1@scale.data)) genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.1@scale.data))
genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.2@scale.data)) genes.hvg.Comb <- intersect(genes.hvg.Comb,rownames(sc10x.2@scale.data))
...@@ -1238,10 +1256,10 @@ scPseudotime <- function(sc10x,i,ds){ ...@@ -1238,10 +1256,10 @@ scPseudotime <- function(sc10x,i,ds){
scMonocle <- reduceDimension(scMonocle,method='DDRTree') scMonocle <- reduceDimension(scMonocle,method='DDRTree')
scMonocle <- orderCells(scMonocle) scMonocle <- orderCells(scMonocle)
postscript(paste0("./analysis/pseudotime/Trajectory.State.eps"),paper="special",width=10,height=10,horizontal=FALSE) postscript(paste0("./analysis/pseudotime/Trajectory.State.eps"),paper="special",width=10,height=10,horizontal=FALSE)
plot_cell_trajectory(scMonocle,color_by="State") plot(plot_cell_trajectory(scMonocle,color_by="State"))
dev.off() dev.off()
postscript(paste0("./analysis/pseudotime/Trajectory.",i,".eps"),paper="special",width=10,height=10,horizontal=FALSE) postscript(paste0("./analysis/pseudotime/Trajectory.",i,".eps"),paper="special",width=10,height=10,horizontal=FALSE)
plot_cell_trajectory(scMonocle,color_by=i) plot(plot_cell_trajectory(scMonocle,color_by=i))
dev.off() dev.off()
my_pseudotime_de <- differentialGeneTest(scMonocle,fullModelFormulaStr = "~sm.ns(Pseudotime)",cores=50) my_pseudotime_de <- differentialGeneTest(scMonocle,fullModelFormulaStr = "~sm.ns(Pseudotime)",cores=50)
...@@ -1249,17 +1267,17 @@ scPseudotime <- function(sc10x,i,ds){ ...@@ -1249,17 +1267,17 @@ scPseudotime <- function(sc10x,i,ds){
my_pseudotime_gene <- my_pseudotime_gene$gene_short_name my_pseudotime_gene <- my_pseudotime_gene$gene_short_name
my_pseudotime_gene <- as.character(my_pseudotime_gene) my_pseudotime_gene <- as.character(my_pseudotime_gene)
postscript(paste0("./analysis/pseudotime/Genes.State.eps"),paper="special",width=10,height=10,horizontal=FALSE) postscript(paste0("./analysis/pseudotime/Genes.State.eps"),paper="special",width=10,height=10,horizontal=FALSE)
plot_genes_branched_pseudotime(scMonocle[my_pseudotime_gene,],color_by="State") plot(plot_genes_branched_pseudotime(scMonocle[my_pseudotime_gene,],color_by="State"))
dev.off() dev.off()
postscript(paste0("./analysis/pseudotime/Genes.",i,".eps"),paper="special",width=10,height=10,horizontal=FALSE) postscript(paste0("./analysis/pseudotime/Genes.",i,".eps"),paper="special",width=10,height=10,horizontal=FALSE)
plot_genes_branched_pseudotime(scMonocle[my_pseudotime_gene,],color_by=i) plot(plot_genes_branched_pseudotime(scMonocle[my_pseudotime_gene,],color_by=i))
dev.off() dev.off()
my_pseudotime_de %>% arrange(qval) %>% head(50) %>% select(gene_short_name) -> gene_to_cluster my_pseudotime_de %>% arrange(qval) %>% head(50) %>% select(gene_short_name) -> gene_to_cluster
gene_to_cluster <- gene_to_cluster$gene_short_name gene_to_cluster <- gene_to_cluster$gene_short_name
postscript(paste0("./analysis/pseudotime/Heatmap.eps"),paper="special",width=10,height=10,horizontal=FALSE) postscript(paste0("./analysis/pseudotime/Heatmap.eps"),paper="special",width=10,height=10,horizontal=FALSE)
my_pseudotime_cluster <- plot_pseudotime_heatmap(scMonocle[gene_to_cluster,],show_rownames=TRUE,return_heatmap=TRUE) my_pseudotime_cluster <- plot_pseudotime_heatmap(scMonocle[gene_to_cluster,],show_rownames=TRUE,return_heatmap=TRUE)
dev.off() dev.off()
#BEAM_res <- BEAM(scMonocle,branch_point=1,cores=50) #BEAM_res <- BEAM(scMonocle,branch_point=1,cores=50)
......
r.scripts/sc-TissueMapper_RUN.Pd.R
+ 130
− 66
View file @ 99423e5c
...@@ -3,7 +3,7 @@ library(methods) ...@@ -3,7 +3,7 @@ library(methods)
library(optparse) library(optparse)
library(Seurat) library(Seurat)
library(readr) library(readr)
library(readxl) #library(readxl)
library(fBasics) library(fBasics)
library(pastecs) library(pastecs)
library(qusage) library(qusage)
...@@ -94,28 +94,28 @@ option_list=list( ...@@ -94,28 +94,28 @@ option_list=list(
make_option("--p",action="store",default="DPrF",type='character',help="Project Name"), make_option("--p",action="store",default="DPrF",type='character',help="Project Name"),
make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"), make_option("--g",action="store",default="ALL",type='character',help="Group To analyze"),
make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"), make_option("--lg",action="store",default=500,type='integer',help="Threshold for cells with minimum genes"),
make_option("--hg",action="store",default=2500,type='integer',help="Threshold for cells with maximum genes"), make_option("--hg",action="store",default=3000,type='integer',help="Threshold for cells with maximum genes"),
make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"), make_option("--lm",action="store",default=0,type='numeric',help="Threshold for cells with minimum %mito genes"),
make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"), make_option("--hm",action="store",default=0.1,type='numeric',help="Threshold for cells with maximum %mito genes"),
make_option("--lx",action="store",default=0.0125,type='numeric',help="x low threshold for hvg selection"), make_option("--lx",action="store",default=0.2,type='numeric',help="x low threshold for hvg selection"),
make_option("--hx",action="store",default=3.5,type='numeric',help="x high threshold for hvg selection"), make_option("--hx",action="store",default=5,type='numeric',help="x high threshold for hvg selection"),
make_option("--ly",action="store",default=0.5,type='numeric',help="y low threshold for hvg selection"), make_option("--ly",action="store",default=1,type='numeric',help="y low threshold for hvg selection"),
make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"), make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"),
make_option("--cca",action="store",default=50,type='integer',help="Number of CCAs to cacluate"), make_option("--cca",action="store",default=50,type='integer',help="Number of CCAs to cacluate"),
make_option("--acca",action="store",default=30,type='integer',help="Number of CCAs to align"), make_option("--acca",action="store",default=30,type='integer',help="Number of CCAs to align"),
make_option("--pc",action="store",default=25,type='integer',help="Number of PCs to cacluate"), make_option("--pc",action="store",default=50,type='integer',help="Number of PCs to cacluate"),
make_option("--hpc",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, pre-stress"), #make_option("--hpc",action="store",default=0.9,type='numeric',help="Max variance cutoff for PCs to use, pre-stress"),
make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"), make_option("--res.prestress",action="store",default=1,type='numeric',help="Resolution to cluster, pre-stress"),
make_option("--st",action="store",default=TRUE,type='logical',help="Remove stressed cells?"), make_option("--st",action="store",default=TRUE,type='logical',help="Remove stressed cells?"),
make_option("--stg",action="store",default="dws",type='character',help="Geneset to use for stress ID"), make_option("--stg",action="store",default="go",type='character',help="Geneset to use for stress ID"),
make_option("--cut.stress",action="store",default=0.9,type='numeric',help="Cutoff for stress score"), make_option("--cut.stress",action="store",default=0.9,type='numeric',help="Cutoff for stress score"),
#make_option("--hpc.poststress",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, post-stress"), #make_option("--hpc.poststress",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, post-stress"),
make_option("--res.poststress",action="store",default=0.15,type='numeric',help="Resolution to cluster, post-stress"), make_option("--res.poststress",action="store",default=0.2,type='numeric',help="Resolution to cluster, post-stress"),
make_option("--ds",action="store",default=10000,type='integer',help="Number of cells to downsample"), #make_option("--ds",action="store",default=250,type='integer',help="Number of cells to downsample"),
#make_option("--hpc.epi",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, Epi"), #make_option("--hpc.epi",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, Epi"),
make_option("--res.epi",action="store",default=0.15,type='numeric',help="Resolution to cluster, Epi"), #make_option("--res.epi",action="store",default=0.2,type='numeric',help="Resolution to cluster, Epi"),
#make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"), #make_option("--hpc.st",action="store",default=0.85,type='numeric',help="Max variance cutoff for PCs to use, St"),
make_option("--res.st",action="store",default=0.15,type='numeric',help="Resolution to cluster, St"), #make_option("--res.st",action="store",default=0.15,type='numeric',help="Resolution to cluster, St"),
make_option("--cut.ne",action="store",default=0.999,type='numeric',help="Cutoff for NE score") make_option("--cut.ne",action="store",default=0.999,type='numeric',help="Cutoff for NE score")
) )
opt=parse_args(OptionParser(option_list=option_list)) opt=parse_args(OptionParser(option_list=option_list))
...@@ -143,50 +143,68 @@ counts.cell.filtered <- results[[4]] ...@@ -143,50 +143,68 @@ counts.cell.filtered <- results[[4]]
counts.gene.filtered <- results[[5]] counts.gene.filtered <- results[[5]]
rm(results) rm(results)
gc()
if (opt$cc==TRUE){
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito","S.Score","G2M.Score"),display.progress=FALSE,do.par=TRUE,num.cores=45)
} else {
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nUMI","percent.mito"),display.progress=FALSE,do.par=TRUE,num.cores=45)
}
gc()
sc10x.D17 <- scSubset(sc10x,"D17","D17") sc10x.D17 <- scSubset(sc10x,"D17","D17")
sc10x.D27 <- scSubset(sc10x,"D27","D27") sc10x.D27 <- scSubset(sc10x,"D27","D27")
sc10x.D35 <- scSubset(sc10x,"D35","D35") sc10x.D35 <- scSubset(sc10x,"D35","D35")
results <- sc3CCA(sc10x.D17,sc10x.D27,sc10x.D35,"D17","D27","D35",cc=opt$cc,cca=opt$cca,acca=opt$acca) results <- sc3CCA(sc10x.D17,sc10x.D27,sc10x.D35,"D17","D27","D35",cc=opt$cc,cca=opt$cca,acca=opt$acca,lx=opt$lx,hx=opt$hx,ly=opt$ly)
sc10x <- results[[1]] sc10x <- results[[1]]
genes.hvg.cca <- results[[2]] genes.hvg.cca <- results[[2]]
rm(results) rm(results)
#results <- scPC(sc10xlx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="pre.stress",cca=TRUE) rm(sc10x.D17)
rm(sc10x.D27)
rm(sc10x.D35)
#results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="pre.stress",cca=TRUE)
#sc10x <- results[[1]] #sc10x <- results[[1]]
#genes.hvg.prestress <- results[[2]] #genes.hvg.prestress <- results[[2]]
#pc.use.prestress <- results[[3]] #pc.use.prestress <- results[[3]]
#rm(results) #rm(results)
sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned") sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.prestress,folder="pre.stress",red="cca.aligned")
#sc10x <- scCluster(sc10x,pc.use=pc.use.prestress,res.use=opt$res.prestress,folder="pre.stress",red="pca")
if (opt$st==TRUE){ if (opt$st==TRUE){
results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,pc.use=pc.use.prestress,cut=opt$cut.stress) results <- scStress(sc10x,stg=opt$stg,res.use=opt$res.prestress,cut=opt$cut.stress)
sc10x <- results[[1]] sc10x <- results[[1]]
counts.cell.filtered.stress <- results[[2]] counts.cell.filtered.stress <- results[[2]]
sc10x.Stress <- results[[3]] sc10x.Stress <- results[[3]]
rm(results) rm(results)
#results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="post.stress",cca=TRUE) #results <- scPC(sc10x,lx=opt$lx,hx=opt$hx,ly=opt$ly,cc=opt$cc,pc=opt$pc,hpc=opt$hpc,file="post.stress",cca=FALSE)
#sc10x <- results[[1]] #sc10x <- results[[1]]
#genes.hvg.poststress <- results[[2]] #genes.hvg.poststress <- results[[2]]
#pc.use.poststress <- results[[3]] #pc.use.poststress <- results[[3]]
#rm(results) #rm(results)
#sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=opt$res.poststress,folder="post.stress",red="pca")
sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned") sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=opt$res.poststress,folder="post.stress",red="cca.aligned")
#sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=0.15,folder="post.stress",red="cca.aligned")
#sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=0.25,folder="post.stress",red="cca.aligned")
#sc10x <- scCluster(sc10x,pc.use=opt$acca,res.use=0.5,folder="post.stress",red="cca.aligned")
} }
gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1) gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Epi" names(gene.set1) <- "Epi"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set,gene.set1) gene.set <- c(gene.set,gene.set1)
rm(gene.set1) rm(gene.set1)
gc() gc()
results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=opt$ds,nm="Lin",folder="lin") min.all <- min(table(sc10x@meta.data[,paste0("res",opt$res.poststress)]))
results <- scQuSAGE(sc10x,gs=gene.set,res.use=opt$res.poststress,ds=min.all,nm="Lin",folder="lin")
sc10x <- results[[1]] sc10x <- results[[1]]
results.cor.Lin <- results[[2]] results.cor.Lin <- results[[2]]
results.clust.Lin.id <- results[[3]] results.clust.Lin.id <- results[[3]]
...@@ -199,8 +217,27 @@ if (any(levels(sc10x@ident)=="Unknown")){ ...@@ -199,8 +217,27 @@ if (any(levels(sc10x@ident)=="Unknown")){
} else { } else {
sc10x.St <- scSubset(sc10x,i="Lin",g="St") sc10x.St <- scSubset(sc10x,i="Lin",g="St")
} }
sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",opt$res.poststress))
sc10x.Epi <- BuildClusterTree(sc10x.Epi,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
sc10x.Epi <- StashIdent(object=sc10x.Epi,save.name=paste0("res",opt$res.poststress))
sc10x.St <- SetAllIdent(object=sc10x.St,id=paste0("res",opt$res.poststress))
sc10x.St <- BuildClusterTree(sc10x.St,do.reorder=TRUE,reorder.numeric=TRUE,do.plot=FALSE)
sc10x.St <- StashIdent(object=sc10x.St,save.name=paste0("res",opt$res.poststress))
sc10x.Epi <- scCluster(sc10x.Epi,pc.use=opt$acca,res.use=opt$res.epi,folder="epi",red="cca.aligned") sc10x.Epi <- RunTSNE(object=sc10x.Epi,reduction.use="cca.aligned",dims.use=1:opt$acca,do.fast=TRUE)
postscript(paste0("./analysis/tSNE/epi/tSNE_Sample.eps"))
plot <- TSNEPlot(object=sc10x.Epi,group.by="samples",pt.size=2.5,do.return=TRUE,vector.friendly=FALSE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
plot(plot)
dev.off()
postscript(paste0("./analysis/tSNE/epi/tSNE_res",opt$res.poststress,".eps"))
plot <- TSNEPlot(object=sc10x.Epi,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE,vector.friendly=FALSE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
plot(plot)
dev.off()
rm(plot)
gene.set1 <- read_delim("./genesets/DEG_BE_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE) gene.set1 <- read_delim("./genesets/DEG_BE_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1) gene.set1 <- as.list(gene.set1)
...@@ -216,14 +253,15 @@ names(gene.set1) <- "OE" ...@@ -216,14 +253,15 @@ names(gene.set1) <- "OE"
gene.set <- c(gene.set,gene.set1) gene.set <- c(gene.set,gene.set1)
rm(gene.set1) rm(gene.set1)
gc() gc()
results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=0.15,ds=opt$ds,nm="Epi.dws",folder="epi") min.epi <- min(table(sc10x.Epi@meta.data[,paste0("res",opt$res.poststress)]))
results <- scQuSAGE(sc10x.Epi,gs=gene.set,res.use=opt$res.poststress,ds=min.epi,nm="Epi.dws",folder="epi")
sc10x.Epi <- results[[1]] sc10x.Epi <- results[[1]]
results.cor.Epi.dws <- results[[2]] results.cor.Epi.dws <- results[[2]]
results.clust.Epi.dws.id <- results[[3]] results.clust.Epi.dws.id <- results[[3]]
rm(results) rm(results)
rm(gene.set) rm(gene.set)
sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",opt$res.epi)) sc10x.Epi <- SetAllIdent(object=sc10x.Epi,id=paste0("res",opt$res.poststress))
OE <- levels(factor(sc10x.Epi@ident[sc10x.Epi@meta.data$Epi.dws=="OE"])) OE <- levels(factor(sc10x.Epi@ident[sc10x.Epi@meta.data$Epi.dws=="OE"]))
OE.cells <- NULL OE.cells <- NULL
for (i in 1:length(OE)){ for (i in 1:length(OE)){
...@@ -243,6 +281,7 @@ dev.off() ...@@ -243,6 +281,7 @@ dev.off()
rm(plot) rm(plot)
rm(OE) rm(OE)
rm(OE.cells) rm(OE.cells)
rm(i)
gene.set1 <- read_csv("./genesets/Basal cells-signature-genes.csv") gene.set1 <- read_csv("./genesets/Basal cells-signature-genes.csv")
gene.set1 <- gene.set1[,2] gene.set1 <- gene.set1[,2]
...@@ -261,11 +300,71 @@ names(gene.set1) <- "CGC" ...@@ -261,11 +300,71 @@ names(gene.set1) <- "CGC"
gene.set<- c(gene.set,gene.set1) gene.set<- c(gene.set,gene.set1)
rm(gene.set1) rm(gene.set1)
gc() gc()
results.cor.Epi.lgea <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=opt$ds,nm="Epi.dws.sub",folder="lgea") results.cor.Epi.lgea <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=min.epi,nm="Epi.dws.sub",folder="lgea")
rm(gene.set) rm(gene.set)
#sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=opt$res.st,folder="st",red="cca.aligned")
#sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=0.2,folder="st",red="cca.aligned")
#sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=0.25,folder="st",red="cca.aligned")
#sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=0.5,folder="st",red="cca.aligned")
#sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=0.1,folder="st",red="cca.aligned")
sc10x.St <- RunTSNE(object=sc10x.St,reduction.use="cca.aligned",dims.use=1:opt$acca,do.fast=TRUE)
postscript(paste0("./analysis/tSNE/st/tSNE_Sample.eps"))
plot <- TSNEPlot(object=sc10x.St,group.by="samples",pt.size=2.5,do.return=TRUE,vector.friendly=FALSE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
plot(plot)
dev.off()
postscript(paste0("./analysis/tSNE/st/tSNE_res",opt$res.poststress,".eps"))
plot <- TSNEPlot(object=sc10x.St,pt.size=5,do.label=TRUE,label.size=10,do.return=TRUE,vector.friendly=FALSE)
plot <- plot+theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),legend.text=element_text(size=20))
plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
plot(plot)
dev.off()
rm(plot)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Endo"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "SM"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Fib"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M10124.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Leu"
gene.set <- c(gene.set,gene.set1)
rm(gene.set1)
gc()
min.st <- min(table(sc10x.Epi@meta.data[,paste0("res",opt$res.poststress)]))
results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=0.2,ds=min.st,nm="St.go",folder="st")
sc10x.St <- results[[1]]
results.cor.St.go <- results[[2]]
results.clust.St.go.id <- results[[3]]
rm(results)
rm(gene.set)
sc10x.Epi.NE <- scNE(sc10x.Epi,neg="EurUro",cut=opt$cut.ne)
sc10x.Epi <- scMergeSubClust(sc10x.Epi,i="Epi.dws.sub",g=c("BE","LE"),nm="Merge")
sc10x.St <- scMergeSubClust(sc10x.St,i="St.go",g=c("Endo","SM","Fib","Leu"),nm="Merge")
sc10x <- scMerge(sc10x,sc10x.Epi,sc10x.St,i.1="Merge",i.2="Merge",nm="Merge_Epi.dws_St.go")
sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws_St.go",i.2="NE",nm="Merge_Epi.dws_St.go_NE")
sc10x.Epi <- scMerge(sc10x.Epi,sc10x.Epi,sc10x.Epi.NE,i.1="Epi.dws.sub",i.2="NE",nm="Epi.dws.sub_NE")
gene.orthog <- read.delim("./genesets/Ensemble.mus-hum.txt") gene.orthog <- read.delim("./genesets/Ensemble.mus-hum.txt")
gene.set1 <- read_excel("./genesets/SupTab3_Consensus_Sigs.xlsx",skip=5) #gene.set1 <- read_excel("./genesets/SupTab3_Consensus_Sigs.xlsx",skip=5)
gene.set1 <- read_csv("./genesets/SupTab3_Consensus_Sigs.csv",skip=6)
gene.set2 <- as.data.frame(gene.set1$Basal[!is.na(gene.set1$Basal)]) gene.set2 <- as.data.frame(gene.set1$Basal[!is.na(gene.set1$Basal)])
colnames(gene.set2) <- "genes" colnames(gene.set2) <- "genes"
gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4]) gene.set2 <- as.data.frame(merge(gene.set2,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
...@@ -303,7 +402,8 @@ gene.set2 <- as.list(gene.set2) ...@@ -303,7 +402,8 @@ gene.set2 <- as.list(gene.set2)
names(gene.set2) <- "Ionocyte" names(gene.set2) <- "Ionocyte"
gene.set <- c(gene.set,gene.set2) gene.set <- c(gene.set,gene.set2)
rm(gene.set2) rm(gene.set2)
gene.set1 <- read_excel("./genesets/SupTab6_Krt13_Hillock.xlsx",skip=5) #gene.set1 <- read_excel("./genesets/SupTab6_Krt13_Hillock.xlsx",skip=5)
gene.set1 <- read_csv("./genesets/SupTab6_Krt13_Hillock.csv",skip=6)
gene.set1 <- gene.set1[gene.set1$FDR<=0.05 & gene.set1$'log2 fold-change (MAST)'>=1.5,1] gene.set1 <- gene.set1[gene.set1$FDR<=0.05 & gene.set1$'log2 fold-change (MAST)'>=1.5,1]
colnames(gene.set1) <- "genes" colnames(gene.set1) <- "genes"
gene.set1 <- as.data.frame(merge(gene.set1,gene.orthog,by.x="genes",by.y="Gene.name")[,4]) gene.set1 <- as.data.frame(merge(gene.set1,gene.orthog,by.x="genes",by.y="Gene.name")[,4])
...@@ -313,48 +413,11 @@ gene.set <- c(gene.set,gene.set1) ...@@ -313,48 +413,11 @@ gene.set <- c(gene.set,gene.set1)
rm(gene.set1) rm(gene.set1)
rm(gene.orthog) rm(gene.orthog)
gene.set <- lapply(gene.set,droplevels) gene.set <- lapply(gene.set,droplevels)
results.cor.Epi.MusLungHierarchy <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=opt$ds,nm="Epi.dws.sub",folder="MusLungHierarchy") results.cor.Epi.MusLungHierarchy <- scQuSAGEsm(sc10x.Epi,gs=gene.set,ds=min.epi,nm="Epi.dws.sub_NE",folder="MusLungHierarchy")
rm(gene.set)
sc10x.St <- scCluster(sc10x.St,pc.use=opt$acca,res.use=opt$res.st,folder="st",red="cca.aligned")
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M11704.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Endo"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M10794.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "SM"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M13024.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Fib"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/DEG_C5.BP.M10124.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- as.list(gene.set1[-1,])
names(gene.set1) <- "Leu"
gene.set <- c(gene.set,gene.set1)
rm(gene.set1)
gc()
results <- scQuSAGE(sc10x.St,gs=gene.set,res.use=0.15,ds=opt$ds,nm="St.go",folder="st")
sc10x.St <- results[[1]]
results.cor.St.go <- results[[2]]
results.clust.St.go.id <- results[[3]]
rm(results)
rm(gene.set) rm(gene.set)
sc10x.Epi.NE <- scNE(sc10x.Epi,neg="EurUro",cut=opt$cut.ne)
sc10x.Epi <- scMergeSubClust(sc10x.Epi,i="Epi.dws.sub",g=c("BE","LE"),nm="Merge")
sc10x.St <- scMergeSubClust(sc10x.St,i="St.go",g=c("Endo","SM","Fib","Leu"),nm="Merge")
sc10x <- scMerge(sc10x,sc10x.Epi,sc10x.St,i.1="Merge",i.2="Merge",nm="Merge_Epi.dws_St.go")
sc10x <- scMerge(sc10x,sc10x,sc10x.Epi.NE,i.1="Merge_Epi.dws_St.go",i.2="NE",nm="Merge_Epi.dws_St.go_NE")
gene.set.c2.all <- read.gmt("./genesets/c2.all.v6.1.symbols.gmt") gene.set.c2.all <- read.gmt("./genesets/c2.all.v6.1.symbols.gmt")
results.cor.c2.ALL <- scQuSAGElg(sc10x,gs=gene.set.c2.all,ds=opt$ds,nm="Merge_Epi.dws_St.go",folder="c2.all") results.cor.c2.ALL <- scQuSAGElg(sc10x,gs=gene.set.c2.all,ds=min.all,nm="Merge_Epi.dws_St.go",folder="c2.all")
rm(gene.set.c2.all) rm(gene.set.c2.all)
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go") sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws_St.go")
...@@ -476,6 +539,7 @@ for (g in c("D17","D27","D35")){ ...@@ -476,6 +539,7 @@ for (g in c("D17","D27","D35")){
} }
sctSNEHighlight(sc10x,i="Pz",g="Pz",file="D") sctSNEHighlight(sc10x,i="Pz",g="Pz",file="D")
sctSNEHighlight(sc10x,i="Tz",g="Tz",file="D") sctSNEHighlight(sc10x,i="Tz",g="Tz",file="D")
rm(i)
rm(g) rm(g)
scCustHeatmap(sc10x.Epi,i="Epi.dws.sub",gs=c(genes.deg.BE.unique,genes.deg.LE.unique,genes.deg.OE1.unique,genes.deg.OE2.unique),g=c("BE","LE","OE1","OE2")) scCustHeatmap(sc10x.Epi,i="Epi.dws.sub",gs=c(genes.deg.BE.unique,genes.deg.LE.unique,genes.deg.OE1.unique,genes.deg.OE2.unique),g=c("BE","LE","OE1","OE2"))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment