Skip to content

Commit 9c943a1

Browse files
author
Baohong Zhang
committed
2 parents e049486 + 7f41b57 commit 9c943a1

9 files changed

Lines changed: 244 additions & 69 deletions

File tree

README.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
# What's new
22
* Upgrade to cellxgene v 1.1.1. After pull, please run "config.sh"
3+
* Install new modules if needed, e.g., pip install anndata==0.8.0
34

45
# Cellxgene VIP unleashes full power of interactive visualization and integrative analysis of scRNA-seq, spatial transcriptomics, and multiome data
56

@@ -64,7 +65,7 @@ R -q -e 'if(!require(devtools)) install.packages("devtools",repos = "http://cran
6465
R -q -e 'if(!require(Cairo)) devtools::install_version("Cairo",version="1.5-12",repos = "http://cran.us.r-project.org")'
6566
R -q -e 'if(!require(foreign)) devtools::install_version("foreign",version="0.8-76",repos = "http://cran.us.r-project.org")'
6667
R -q -e 'if(!require(ggpubr)) devtools::install_version("ggpubr",version="0.3.0",repos = "http://cran.us.r-project.org")'
67-
R -q -e 'if(!require(ggrastr)) devtools::install_version("ggrastr",version="0.1.9",repos = "http://cran.us.r-project.org")'
68+
R -q -e 'if(!require(ggrastr)) devtools::install_version("ggrastr",version="0.2.1",repos = "http://cran.us.r-project.org")'
6869
R -q -e 'if(!require(arrow)) devtools::install_version("arrow",version="2.0.0",repos = "http://cran.us.r-project.org")'
6970
R -q -e 'if(!require(Seurat)) devtools::install_version("Seurat",version="3.2.3",repos = "http://cran.us.r-project.org")'
7071
R -q -e 'if(!require(rmarkdown)) devtools::install_version("rmarkdown",version="2.5",repos = "http://cran.us.r-project.org")'

VIP.macOS.yml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@ dependencies:
8181
- zlib=1.2.11
8282
- pip:
8383
- aniso8601==8.1.1
84-
- anndata==0.7.4
84+
- anndata==0.8.0
8585
- attrs==20.3.0
8686
- batchglm==0.7.4
8787
- bleach==3.2.3
@@ -109,7 +109,7 @@ dependencies:
109109
- fsspec==0.7.4
110110
- get-version==2.1
111111
- gunicorn==20.0.4
112-
- h5py==2.10.0
112+
- h5py==3.7.0
113113
- idna==2.10
114114
- iniconfig==1.1.1
115115
- itsdangerous==1.1.0

VIP.yml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -85,7 +85,7 @@ dependencies:
8585
- zlib=1.2.11
8686
- pip:
8787
- aniso8601==8.1.1
88-
- anndata==0.7.4
88+
- anndata==0.8.0
8989
- attrs==20.3.0
9090
- batchglm==0.7.4
9191
- bleach==3.2.3
@@ -113,7 +113,7 @@ dependencies:
113113
- fsspec==0.7.4
114114
- get-version==2.1
115115
- gunicorn==20.0.4
116-
- h5py==2.10.0
116+
- h5py==3.7.0
117117
- idna==2.10
118118
- iniconfig==1.1.1
119119
- itsdangerous==1.1.0

VIPInterface.py

Lines changed: 21 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ def getObs(data):
127127
anno = []
128128
sel = list(set(selAnno)&set(dAnno))
129129
if len(sel)>0:
130-
tmp = scD.data.obs.loc[selC,sel].astype('str')
130+
tmp = scD.data.obs.iloc[selC,:].loc[:,sel].astype('str')
131131
tmp.index = cNames
132132
anno += [tmp]
133133
sel = list(set(selAnno)-set(dAnno))
@@ -770,7 +770,7 @@ def getGSEA(data):
770770
return json.dumps(sorted([os.path.basename(i).replace(".symbols.gmt","") for i in glob.glob(strGSEA+"*.symbols.gmt")]))
771771

772772
def DEG(data):
773-
adata = None;
773+
adata = None
774774
genes = data['genes']
775775
data['genes'] = []
776776
comGrp = 'cellGrp'
@@ -839,8 +839,8 @@ def DEG(data):
839839
raise ValueError('Unknown DE methods:'+data['DEmethod'])
840840
#res = de.test.two_sample(adata,comGrp,test=data['DEmethod'],noise_model=nm)
841841
deg = res.summary()
842-
deg = deg.sort_values(by=['qval']).loc[:,['gene','log2fc','pval','qval']]
843842
deg['log2fc'] = -1 * deg['log2fc']
843+
deg = deg.sort_values(by=['qval']).loc[:,['gene','log2fc','pval','qval']]
844844
#del adata
845845
## plot in R
846846
#strF = ('/tmp/DEG%f.csv' % time.time())
@@ -1638,7 +1638,7 @@ def getBWinfo(data):
16381638
strD = re.sub(".h5ad$","/",data["h5ad"])
16391639
if os.path.isdir(strD):
16401640
for one in os.listdir(strD):
1641-
if not re.search("bw$",one)==None:
1641+
if one.endswith("bw"):#not re.search("bw$",one)==None:
16421642
BWinfo["BWfile"].append(one)
16431643
elif one=="annotation.rds":
16441644
BWinfo["BWannotation"]="annotation.rds"
@@ -1647,42 +1647,39 @@ def getBWinfo(data):
16471647
elif one=="links.rds":
16481648
BWinfo["BWlink"]="links.rds"
16491649
elif one=="bw.cluster":
1650-
BWinfo["BWcluster"]="bw.cluster"
1650+
BWinfo["BWcluster"]=pd.read_csv(strD+one,sep="\t",header=0) #"bw.cluster" .to_csv(index=False)
1651+
if len(BWinfo["BWcluster"])>0:
1652+
BWinfo["BWcluster"]=BWinfo["BWcluster"][BWinfo["BWcluster"]['Wig'].isin(BWinfo["BWfile"])].to_csv(index=False)
16511653
return json.dumps(BWinfo)
16521654

16531655
def plotBW(data):
16541656
strD = re.sub(".h5ad$","/",data["h5ad"])
16551657
strCSV = ('%s/BW%f.csv' % (data["CLItmp"],time.time()))
16561658
## select all cells
16571659
strType = strD + 'bw.cluster'
1658-
data['bw']=['%s.bw'%one for one in data['bw']]
1659-
grpFlag = False
16601660
if os.path.isfile(strType) and len(data['genes'])>0:
1661-
with open(strType,"r") as f:
1662-
grp = f.readline().strip()
1663-
#with app.get_data_adaptor(url_dataroot=data['url_dataroot'],dataset=data['dataset']) as scD:
1661+
clusterInfo = pd.read_csv(strType,sep="\t",header=0,index_col=0)
1662+
clusterInfo = clusterInfo.loc[data['bw'],:]
1663+
grp = []
16641664
scD=data['data_adapter']
16651665
if scD is not None:
16661666
dAnno = list(scD.get_obs_keys())
1667-
if grp in dAnno:
1668-
grpFlag = True
1669-
if grpFlag:
1670-
data['grp'] = [grp]
1667+
grp = [i for i in clusterInfo.columns if i in dAnno]
1668+
if len(grp)>0:
1669+
data['grp'] = grp
16711670
adata = createData(data)
1672-
if len(adata)==0:
1673-
grpFlag = False
1674-
else:
1675-
cluster = pd.read_csv(strType,sep="\t",header=None,index_col=1,skiprows=1)#delimiter="\n",
1676-
cluster = cluster[cluster[0].isin(data['bw'])]
1677-
adata = adata[adata.obs[grp].isin(list(cluster.index)),:]
1678-
obsCluster = pd.DataFrame(list(cluster.loc[adata.obs[grp],:][0]),index=adata.obs.index,columns=[grp])
1679-
pd.concat([obsCluster,adata.to_df()], axis=1, sort=False).to_csv(strCSV)
1671+
if len(adata)>0:
1672+
selC=adata.obs[grp[0]].isin(list(clusterInfo[grp[0]]))
1673+
for oneG in grp:
1674+
selC = selC | adata.obs[oneG].isin(list(clusterInfo[oneG]))
1675+
adata = adata[selC,:]
1676+
pd.concat([adata.obs,adata.to_df()], axis=1, sort=False).to_csv(strCSV)
16801677
## plot in R
1681-
#strCMD = ' '.join([strExePath+'/browserPlot.R',strD,data['region'],str(data['exUP']),str(data['exDN']),strCSV,str(data['cutoff']),data['figOpt']['img'],str(data['figOpt']['fontsize']),str(data['figOpt']['dpi']),data['Rlib']])
1678+
#strCMD = ' '.join([strExePath+'/browserPlot.R',strD,data['region'],','.join(data['bw']),str(data['exUP']),str(data['exDN']),strCSV,str(data['cutoff']),data['figOpt']['img'],str(data['figOpt']['fontsize']),str(data['figOpt']['dpi']),data['Rlib']])
16821679
#ppr.pprint(strCMD)
16831680
res = subprocess.run([strExePath+'/browserPlot.R',strD,data['region'],','.join(data['bw']),str(data['exUP']),str(data['exDN']),strCSV,str(data['cutoff']),data['figOpt']['img'],str(data['figOpt']['fontsize']),str(data['figOpt']['dpi']),data['Rlib']],capture_output=True)#
16841681
img = res.stdout.decode('utf-8')
1685-
if grpFlag:
1682+
if len(adata)>0:
16861683
os.remove(strCSV)
16871684
if 'Error' in res.stderr.decode('utf-8'):
16881685
raise SyntaxError("in R: "+res.stderr.decode('utf-8'))

VIP_conda_R.yml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,7 @@ dependencies:
131131
- zstd=1.4.5
132132
- pip:
133133
- aniso8601==8.1.1
134-
- anndata==0.7.4
134+
- anndata==0.8.0
135135
- async-generator==1.10
136136
- attrs==20.3.0
137137
- batchglm==0.7.4
@@ -159,7 +159,7 @@ dependencies:
159159
- fsspec==0.7.4
160160
- get-version==2.1
161161
- gunicorn==20.0.4
162-
- h5py==2.10.0
162+
- h5py==3.7.0
163163
- idna==2.10
164164
- iniconfig==1.1.1
165165
- itsdangerous==1.1.0

VIPlight.yml

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
channels:
2+
- conda-forge
3+
- defaults
4+
- bioconda
5+
dependencies:
6+
- scanpy=1.9.1
7+
- matplotlib=3.6.2
8+
- seaborn=0.12.1
9+
- plotly=5.11.0
10+
- pip=23.0.1
11+
- pyarrow=10.0.1
12+
- r-base=4.2.2
13+
- r-ggplot2=3.4.0
14+
- r-reshape2=1.4.4
15+
- r-ggforce=0.4.1
16+
- r-patchwork=1.1.2
17+
- r-tidyverse=1.3.2
18+
- r-ggpubr=0.5.0
19+
- r-cairo=1.6_0
20+
- r-dplyr=1.0.10
21+
- r-png=0.1_8
22+
- r-mass=7.3_58.1
23+
- r-data.table=1.14.6
24+
- r-rcolorbrewer=1.1_3
25+
- r-glue=1.6.2
26+
- r-gridextra=2.3
27+
- r-ggrepel=0.9.2
28+
- r-ggrastr=1.0.1
29+
- r-viridis=0.6.2
30+
- bioconductor-rtracklayer=1.58.0
31+
- bioconductor-complexheatmap=2.14.0
32+
- bioconductor-fgsea=1.24.0
33+
- pip:
34+
- diffxpy==0.7.4

browserPlot.R

Lines changed: 40 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ main <- function(){
1616
## process input ----
1717
args <- commandArgs(trailingOnly=TRUE)
1818
suppressMessages(suppressWarnings(load(tail(args,1))))
19-
strPath <- args[1]
19+
strPath <- paste0(normalizePath(args[1]),"/")
2020
region <- args[2]
2121
bwList <- unlist(strsplit(args[3],","))
2222
extend.upstream <- as.numeric(args[4])
@@ -29,24 +29,25 @@ main <- function(){
2929

3030
## obtain the region -----
3131
annotations <- NULL
32-
if(file.exists(paste0(strPath,"/annotation.rds"))) annotations <- readRDS(paste0(strPath,"/annotation.rds"))
32+
if(file.exists(paste0(strPath,"annotation.rds"))) annotations <- readRDS(paste0(strPath,"annotation.rds"))
3333
region <- customFindRegion(
3434
region = region,
3535
annotations = annotations,
3636
extend.upstream = extend.upstream,
3737
extend.downstream = extend.downstream
3838
)
3939

40-
## plot the bigwig -----
41-
strBW <- list.files(strPath,"bw$",full.names=T)
42-
strBW <- strBW[basename(strBW)%in%bwList]
40+
## plot the bigwig -----
41+
bwList <- bwList[bwList%in%list.files(strPath,"bw$")]
42+
strBW <- paste0(strPath,bwList)
4343
AllPlots <- suppressWarnings(suppressMessages(customBigwigTrack(strBW,region,fontsize)))
4444
#h <- rep(1,length(AllPlots))
4545
h <- length(strBW)
4646
w <- 6
4747

4848
## plot gene expression ----
49-
expPlot <- customExpressionPlot(strExp,strBW,expCutoff,fontsize)
49+
strCluster <- paste0(strPath,"/bw.cluster")
50+
expPlot <- customExpressionPlot(strExp,strCluster,bwList,expCutoff,fontsize)
5051
if(!is.null(expPlot)){
5152
w <- c(w,nlevels(expPlot$data$gene))
5253
}
@@ -113,7 +114,7 @@ customFindRegion <- function(region,annotations=NULL,sep = c("-", "-"),extend.up
113114
}
114115
)
115116
if (is.null(x = region)) {
116-
stop("Gene not found")
117+
stop("Gene/Region is not found")
117118
}
118119
}
119120
region <- suppressWarnings(expr = Extend(
@@ -124,6 +125,17 @@ customFindRegion <- function(region,annotations=NULL,sep = c("-", "-"),extend.up
124125
)
125126
return(region)
126127
}
128+
StringToGRanges <- function(regions, sep = c("-", "-"), ...) {
129+
ranges.df <- data.frame(ranges = regions)
130+
ranges.df <- tidyr::separate(
131+
data = ranges.df,
132+
col = "ranges",
133+
sep = paste0(sep[[1]], "|", sep[[2]]),
134+
into = c("chr", "start", "end")
135+
)
136+
granges <- makeGRangesFromDataFrame(df = ranges.df, ...)
137+
return(granges)
138+
}
127139
customLookupGeneCoords <- function(annotations, gene) {
128140
isgene <- annotations$gene_name == gene
129141
isgene <- !is.na(x = isgene) & isgene
@@ -137,20 +149,30 @@ customLookupGeneCoords <- function(annotations, gene) {
137149
return(gr)
138150
}
139151
}
140-
customExpressionPlot <- function(strExp,strBW,geneCutoff=-0.1,fontsize=9) {
152+
customExpressionPlot <- function(strExp,strCluster,bwList,geneCutoff=-0.1,fontsize=9) {
141153
if(!file.exists(strExp)) return(NULL)
142154
Exp <- read.csv(strExp,row.names=1,as.is=T)
143-
if(sum(!Exp[,1]%in%basename(strBW))>0) return(NULL)
155+
clusterInfo <- read.table(strCluster,row.names=1,header=T,as.is=T,sep="\t")[bwList,,drop=F]
156+
grp <- intersect(colnames(Exp),colnames(clusterInfo))
157+
if(length(grp)==0) return(NULL)
158+
features <- colnames(Exp)[!colnames(Exp)%in%grp][-1]
144159
# get data
145-
data.plot <- t(Exp[,-1,drop=F])
146-
obj.groups <- Exp[,1,drop=F]
147-
features <- row.names(data.plot)
160+
df <- NULL
161+
for(oneBW in bwList){
162+
selC <- rep(T,nrow(Exp))
163+
for(oneG in grp){
164+
if(nchar(clusterInfo[oneBW,oneG])>0){
165+
selC <- selC & Exp[,oneG]==clusterInfo[oneBW,oneG]
166+
}
167+
}
168+
if(sum(selC)==0) stop(paste("No cell for",oneBW))
169+
oneD = setNames(melt(Exp[selC,features,drop=F],id.vars=NULL),c("gene","expression"))
170+
oneD$group <- oneBW
171+
df <- rbind(df,oneD)
172+
}
173+
df$group <- factor(df$group,levels=bwList)
148174
# set the color
149-
levels.use <- basename(strBW)
150-
colors_all <- setNames(gg_color_hue(length(levels.use)),levels.use)
151-
# construct data frame
152-
df <- setNames(melt(data.plot),c("gene","Cell","expression"))
153-
df <- cbind(df,group=obj.groups[df$Cell,])
175+
colors_all <- setNames(gg_color_hue(length(bwList)),bwList)
154176

155177
p.list <- list()
156178
for (i in seq_along(along.with = features)) {
@@ -167,6 +189,7 @@ customExpressionPlot <- function(strExp,strBW,geneCutoff=-0.1,fontsize=9) {
167189
scale_y_discrete(position = "top") +
168190
scale_x_continuous(position = "bottom", limits = c(0, NA)) +
169191
geom_text(data=df.rate,aes(xpos,ypos,label=text,color=group),hjust="right",vjust=1.1, size=fontsize-7)+
192+
scale_fill_manual(values = colors_all)+
170193
theme(
171194
axis.text.x.top = element_text(size=fontsize, face="bold", angle=0), # Gene names @ top
172195
axis.text.x.bottom = element_text(size=fontsize-1, angle=60, vjust=0.6),
@@ -175,9 +198,6 @@ customExpressionPlot <- function(strExp,strBW,geneCutoff=-0.1,fontsize=9) {
175198
strip.text.y = element_blank(),
176199
legend.position = "none"
177200
)
178-
if (!is.null(x = levels.use)) {
179-
p <- p + scale_fill_manual(values = colors_all)
180-
}
181201
p.list[[i]] <- p
182202
}
183203
p <- patchwork::wrap_plots(p.list, ncol = length(x = p.list))

0 commit comments

Comments
 (0)