```{r}
library(Seurat)
library(ggplot2)
library(dplyr)
#library(Matrix)
library(ggpubr)
library(gplots)

library(data.table)


library(future)
plan("multisession", workers = 32)
```

```{r}
feature=readRDS("1488_feature.RDS")
raw.lipid=readRDS("1488_all.RDS")
fake.lipid.pos=read.csv("diaoyu-IMSML-coords.csv")
#lipid.xpos.old=lipid.xpos$x
#lipid.ypos.old=lipid.ypos$x

lipid.xpos=fake.lipid.pos$x_original
lipid.ypos=fake.lipid.pos$y_original
lipid.registered=read.csv("A1-IMSML-coords.csv") 
lipid.xpos.registered=lipid.registered$x_micro_physical
lipid.ypos.registered=lipid.registered$y_micro_physical

raw.lipid=raw.lipid %>% t()

rownames(raw.lipid)=feature$mzLow %>% factor()
colnames(raw.lipid)=seq(1,ncol(raw.lipid)) %>% as.character()
```



```{r}
lipid.fake.spatial=data.frame(x=lipid.xpos,y=lipid.ypos)
ggplot(lipid.fake.spatial,aes(x,y))+
  scale_y_reverse()+
  geom_point()
coords=data.frame(x_new=-(lipid.xpos-min(lipid.xpos))+max(lipid.xpos),y_new=-(lipid.ypos-min(lipid.ypos))+max(lipid.ypos)

)
ggplot(coords,aes(x_new,y_new))+
  scale_y_reverse()+
  geom_point()+
  geom_vline(xintercept = 1390)

coords$drug=drug_h
ggplot(coords, aes(x = x_new, y = y_new, color = drug)) +
  geom_point(size=0.0000001,alpha=0.5) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Drug Concentration Scatter Plot", 
       x = "X Coordinate", 
       y = "Y Coordinate", 
       color = "Drug Concentration") +
  theme_minimal()+
  scale_y_reverse()
```
```{r}


filtered.lung.lipid.noct=raw.lipid[,coords$y_new[colnames(raw.lipid) %>% as.integer()]>200]

filtered.lung.lipid.noct.xpos=coords$x_new[coords$y_new[colnames(raw.lipid) %>% as.integer()]>200]

filtered.lung.lipid.noct.ypos=coords$y_new[coords$y_new[colnames(raw.lipid) %>% as.integer()]>200]


filtered.lung.lipid.bg=raw.lipid[,coords$y_new[colnames(raw.lipid) %>% as.integer()]<200]
filtered.lung.lipid.bg.xpos=coords$x_new[coords$y_new[colnames(raw.lipid) %>% as.integer()]<200]
filtered.lung.lipid.bg.ypos=coords$y_new[coords$y_new[colnames(raw.lipid) %>% as.integer()]<200]




lipid.fake.spatial=data.frame(x=filtered.lung.lipid.noct.xpos,y=filtered.lung.lipid.noct.ypos)
ggplot(lipid.fake.spatial,aes(x,y))+
  geom_tile()
lipid.fake.spatial=data.frame(x=filtered.lung.lipid.bg.xpos,y=filtered.lung.lipid.bg.ypos)
ggplot(lipid.fake.spatial,aes(x,y))+
  geom_tile()

```
```{r}


spec.ct=data.frame(mz=feature$mzLow%>% as.numeric(),lung.ct=filtered.lung.lipid.bg %>% rowMeans()) 
spec.mean.plot <- melt(spec.ct, id = "mz")
ggplot(spec.mean.plot, aes(x = mz, y = value, color = variable)) + 
  geom_line(size = 1)  + 
  scale_color_brewer(palette = "Set1")

spec.mean=data.frame(mz=feature$mzLow %>% as.numeric(),tissue=filtered.lung.lipid.noct %>% rowMeans(),lung.ct=filtered.lung.lipid.bg %>% rowMeans())

spec.mean.plot <- melt(spec.mean, id = "mz")

ggplot(spec.mean, aes(x = mz, y =tissue)) + 
  geom_line(size = 1)  + 
  scale_color_brewer(palette = "Set1")

ggplot(spec.mean, aes(x = mz, y =lung.ct)) + 
  geom_line(size = 1)  + 
  scale_color_brewer(palette = "Set1")

ggplot(spec.mean.plot, aes(x = mz, y = value, color = variable)) + 
  geom_line(size = 0.5) +
  scale_color_brewer(palette = "Set1") +
  theme_minimal()+
  ylab("Intensity")+
  xlab("m/z")


```
```{r}
biggerthanback=spec.mean$tissue>10*spec.mean$lung.ct
nrow(spec.mean[biggerthanback,])
```

```{r}
filtered.lung.lipid.noct.threholded=filtered.lung.lipid.noct[biggerthanback,]
rm(raw.lipid)
rm(filtered.lung.lipid.noct)
filtered.lung.lipid.noct.threholded.rowSum.log=rowMeans(filtered.lung.lipid.noct.threholded) %>% log10() %>% data.frame()
length(filtered.lung.lipid.noct.threholded.rowSum.log[filtered.lung.lipid.noct.threholded.rowSum.log>0.6])
hist(filtered.lung.lipid.noct.threholded.rowSum.log$.,breaks=100)

```


```{r}
drug=readRDS("/gpfs/data/yanailab/1488_drug.RDS")
drug_name=readRDS("1488_drug_feature.RDS")
drug %>% colMeans()
#drug_h=drug[,3]+drug[,5]+drug[,1]+drug[,6]+drug[,2]
drug_h=drug[,2]+drug[,6]
hist(drug_h %>% log10(),breaks=100)
#drug_ct=drug_h[coords$y_new[colnames(filtered.lung.lipid.bg) %>% as.integer()]]
drug_ct=drug_h[colnames(filtered.lung.lipid.bg) %>% as.integer()]
drug_1=drug_h[coords$y_new>200&coords$x_new<500]
drug_2=drug_h[coords$y_new>200&coords$x_new>500&coords$x_new<1000]
drug_3=drug_h[coords$y_new>200&coords$x_new>1000&coords$x_new<1390]
drug_4=drug_h[coords$y_new>200&coords$x_new>1390]

data <- data.frame(
  value = c(drug_ct, drug_1, drug_2, drug_3, drug_4),
  group = factor(rep(c("drug_ct", "drug_1", "drug_2", "drug_3", "drug_4"), 
                     times = c(length(drug_ct), length(drug_1), length(drug_2), length(drug_3), length(drug_4)))
))

# Create the histogram
ggplot(data, aes(x = value, fill = group)) +
  geom_histogram(bins = 30, position = "identity", alpha = 0.5) +
  theme_minimal() +
  labs(x = "Value", y = "Count", fill = "Drug Group") +
  scale_fill_brewer(palette = "Set1")  # You can choose different palettes

ggplot(data, aes(x = value)) +
  geom_histogram(bins = 100, fill = "blue", color = "black") +
  facet_wrap(~ group, scales = "fixed") +  # Fixed scales for all facets
  theme_minimal() +
  labs(x = "Value", y = "Count") +
  theme(strip.text.x = element_text(size = 12))  # Adjust facet label size
```

```{r}
lipid_1=filtered.lung.lipid.noct.threholded[,filtered.lung.lipid.noct.xpos<500]
lipid_2=filtered.lung.lipid.noct.threholded[,filtered.lung.lipid.noct.xpos>500&filtered.lung.lipid.noct.xpos<1000]
lipid_3=filtered.lung.lipid.noct.threholded[,filtered.lung.lipid.noct.xpos>1000&filtered.lung.lipid.noct.xpos<1390]
lipid_4=filtered.lung.lipid.noct.threholded[,filtered.lung.lipid.noct.xpos>1390]
feature_afterfilter=feature[biggerthanback,]
```
#correlation between lipids and drug
```{r}
library(doParallel)
library(foreach)

# Register the parallel backend to use several cores
numCores <- 32 # You mentioned having 32 cores
registerDoParallel(cores=numCores)

# Define the function to compute correlations, adapted for parallel execution
compute_correlations_parallel <- function(lipid, drug) {
  correlations <- foreach(i=1:nrow(lipid), .combine=c) %dopar% {
    cor(lipid[i, ], drug, use="complete.obs")
  }
  return(correlations)
}

# Assuming lipid_1, lipid_2, lipid_3, lipid_4, drug_1, drug_2, drug_3, drug_4 are already defined
# Compute correlations in parallel
correlations_lipid_1 <- compute_correlations_parallel(lipid_1, drug_1)
correlations_lipid_2 <- compute_correlations_parallel(lipid_2, drug_2)
correlations_lipid_3 <- compute_correlations_parallel(lipid_3, drug_3)
correlations_lipid_4 <- compute_correlations_parallel(lipid_4, drug_4)

# Stop the parallel backend when done
stopImplicitCluster()
correlations_lipid_2 %>% max()
```
```{r}
# Assuming correlations_lipid_1, correlations_lipid_2, correlations_lipid_3, correlations_lipid_4 
# are the vectors containing the correlation coefficients you've computed

# Plot histogram for correlations of lipid_1 and drug_1
hist(correlations_lipid_1, main="Histogram of Correlations for Lipid 1 and Drug 1", xlab="Correlation Coefficient", col="blue", border="black",breaks=100)

# Plot histogram for correlations of lipid_2 and drug_2
hist(correlations_lipid_2, main="Histogram of Correlations for Lipid 2 and Drug 2", xlab="Correlation Coefficient", col="green", border="black",breaks=100)

# Plot histogram for correlations of lipid_3 and drug_3
hist(correlations_lipid_3, main="Histogram of Correlations for Lipid 3 and Drug 3", xlab="Correlation Coefficient", col="red", border="black",breaks=100)

# Plot histogram for correlations of lipid_4 and drug_4
hist(correlations_lipid_4, main="Histogram of Correlations for Lipid 4 and Drug 4", xlab="Correlation Coefficient", col="yellow", border="black",breaks=100)

```
```{r}

```


```{r}
all.lipid=CreateSeuratObject(counts=filtered.lung.lipid.noct.threholded,min.cells = 1000,min.features = 100)
plan(sequential)
all.lipid <- NormalizeData(all.lipid, verbose = FALSE, assay = "RNA")

all.lipid <- ScaleData(all.lipid, features = rownames(filtered.lung.lipid.noct.threholded))
plan("multisession", workers = 24)
#all.lipid =FindVariableFeatures(all.lipid)
#all.lipid <- RunPCA(all.lipid)
all.lipid <- RunPCA(all.lipid,features = rownames(filtered.lung.lipid.noct.threholded))
VizDimLoadings(all.lipid, dims = 1:2, reduction = "pca")
DimPlot(all.lipid, reduction = "pca")
ElbowPlot(all.lipid)
all.lipid <- RunUMAP(all.lipid, dims = 1:20)
all.lipid<- FindNeighbors(all.lipid, dims = 1:20)
all.lipid <- FindClusters(all.lipid, resolution = 0.4)
DimPlot(all.lipid, reduction = "umap")
DimPlot(all.lipid, reduction = "pca")
```

```{r}

filtered.lipid.fake.spatial=data.frame(x=lipid.xpos.registered,y=lipid.ypos.registered)
filtered.lipid.fake.spatial$iden=Idents(B1.lipid)
ggplot(filtered.lipid.fake.spatial,aes(x,y))+
  geom_point(aes(color=iden))+
#  geom_point(data=D1.ST.tissuepos,size=1.4,shape=16,aes(fill="black"))+ 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),plot.background=element_blank())+
  coord_fixed()+
  scale_y_reverse()
```

```{r}
rm(all.lipid)
load("1488.RData")
#a=filtered.lung.lipid.noct.xpos>500&filtered.lung.lipid.noct.xpos<1000
a=coords$y_new[colnames(filtered.lung.lipid.noct.threholded) %>% as.integer()]>200&coords$x_new[colnames(filtered.lung.lipid.noct.threholded) %>% as.integer()]>500&coords$x_new[colnames(filtered.lung.lipid.noct.threholded) %>% as.integer()]<1000
B1.lipid=filtered.lung.lipid.noct.threholded[,a]
B1.lipid.count=B1.lipid
B1.xpos=filtered.lung.lipid.noct.xpos[a]
B1.ypos=filtered.lung.lipid.noct.ypos[a]
```


```{r}


#filtered.lipid.fake.spatial=data.frame(x=B1.xpos,y=B1.ypos)
filtered.lipid.fake.spatial=data.frame(x = coords$x_new[coords$y_new>200&coords$x_new>500&coords$x_new<1000], y =  coords$y_new[coords$y_new>200&coords$x_new>500&coords$x_new<1000])
filtered.lipid.fake.spatial$drug=drug_2


ggplot(filtered.lipid.fake.spatial, aes(x=x,y=y, color = drug)) +
  geom_point(size=0.01,alpha=0.5) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "Drug Concentration Scatter Plot", 
       x = "X Coordinate", 
       y = "Y Coordinate", 
       color = "Drug Concentration") +
  theme_minimal()+
  scale_y_reverse()
```


```{r}
B1.lipid=CreateSeuratObject(counts=B1.lipid.count,min.cells =0,min.features = 0)
plan(sequential)
B1.lipid <- NormalizeData(B1.lipid, verbose = FALSE, assay = "RNA")

B1.lipid <- ScaleData(B1.lipid, features = rownames(filtered.lung.lipid.noct.threholded))
plan("multisession", workers = 14)
#B1.lipid =FindVariableFeatures(B1.lipid)
#B1.lipid <- RunPCA(B1.lipid)
B1.lipid <- RunPCA(B1.lipid,features = rownames(filtered.lung.lipid.noct.threholded))
VizDimLoadings(B1.lipid, dims = 1:2, reduction = "pca")
DimPlot(B1.lipid, reduction = "pca")
ElbowPlot(B1.lipid)
B1.lipid <- RunUMAP(B1.lipid, dims = 1:20)
B1.lipid<- FindNeighbors(B1.lipid, dims = 1:20)
B1.lipid <- FindClusters(B1.lipid, resolution = 0.3)
DimPlot(B1.lipid, reduction = "umap")
DimPlot(B1.lipid, reduction = "pca")
```

```{r}
filtered.lipid.fake.spatial=data.frame(x=B1.xpos,y=B1.ypos)
filtered.lipid.fake.spatial$iden=Idents(B1.lipid)
ggplot(filtered.lipid.fake.spatial,aes(x,y))+
  geom_point(size=0.5,aes(color=iden))+
#  geom_point(data=D1.ST.tissuepos,size=1.4,shape=16,aes(fill="black"))+ 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),plot.background=element_blank())+
  coord_fixed()+
  scale_y_reverse()
```

```{r}
a=filtered.lung.lipid.noct.xpos>1000&filtered.lung.lipid.noct.xpos<1390
C1.lipid=filtered.lung.lipid.noct.threholded[,a]
C1.lipid.count=C1.lipid
C1.xpos=filtered.lung.lipid.noct.xpos[a]
C1.ypos=filtered.lung.lipid.noct.ypos[a]
```

```{r}
C1.lipid=CreateSeuratObject(counts=C1.lipid.count,min.cells =0,min.features = 0)
plan(sequential)
C1.lipid <- NormalizeData(C1.lipid, verbose = FALSE, assay = "RNA")

C1.lipid <- ScaleData(C1.lipid, features = rownames(filtered.lung.lipid.noct.threholded))
plan("multisession", workers = 32)
#C1.lipid =FindVariableFeatures(C1.lipid)
#C1.lipid <- RunPCA(C1.lipid)
C1.lipid <- RunPCA(C1.lipid,features = rownames(filtered.lung.lipid.noct.threholded))
VizDimLoadings(C1.lipid, dims = 1:2, reduction = "pca")
DimPlot(C1.lipid, reduction = "pca")
ElbowPlot(C1.lipid)
C1.lipid <- RunUMAP(C1.lipid, dims = 1:20)
C1.lipid<- FindNeighbors(C1.lipid, dims = 1:20)
C1.lipid <- FindClusters(C1.lipid, resolution = 0.3)
DimPlot(C1.lipid, reduction = "umap")
DimPlot(C1.lipid, reduction = "pca")
```

```{r}
filtered.lipid.fake.spatial=data.frame(x=C1.xpos,y=C1.ypos)
filtered.lipid.fake.spatial$iden=Idents(C1.lipid)
ggplot(filtered.lipid.fake.spatial,aes(x,y))+
  geom_point(size=0.5,aes(color=iden))+
#  geom_point(data=D1.ST.tissuepos,size=1.4,shape=16,aes(fill="black"))+ 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),plot.background=element_blank())+
  coord_fixed()+
  scale_y_reverse()
```

```{r}
a=filtered.lung.lipid.noct.xpos>1390
D1.lipid=filtered.lung.lipid.noct.threholded[,a]
D1.lipid.count=D1.lipid
D1.xpos=filtered.lung.lipid.noct.xpos[a]
D1.ypos=filtered.lung.lipid.noct.ypos[a]
```

```{r}
D1.lipid=CreateSeuratObject(counts=D1.lipid.count,min.cells =0,min.features = 0)
plan(sequential)
D1.lipid <- NormalizeData(D1.lipid, verbose = FALSE, assay = "RNA")

D1.lipid <- ScaleData(D1.lipid, features = rownames(filtered.lung.lipid.noct.threholded))
plan("multisession", workers = 32)
#D1.lipid =FindVariableFeatures(D1.lipid)
#D1.lipid <- RunPCA(D1.lipid)
D1.lipid <- RunPCA(D1.lipid,features = rownames(filtered.lung.lipid.noct.threholded))
VizDimLoadings(D1.lipid, dims = 1:2, reduction = "pca")
DimPlot(D1.lipid, reduction = "pca")
ElbowPlot(D1.lipid)
D1.lipid <- RunUMAP(D1.lipid, dims = 1:20)
D1.lipid<- FindNeighbors(D1.lipid, dims = 1:20)
D1.lipid <- FindClusters(D1.lipid, resolution = 0.3)
DimPlot(D1.lipid, reduction = "umap")
DimPlot(D1.lipid, reduction = "pca")
```

```{r}
filtered.lipid.fake.spatial=data.frame(x=D1.xpos,y=D1.ypos)
filtered.lipid.fake.spatial$iden=Idents(D1.lipid)
ggplot(filtered.lipid.fake.spatial,aes(x,y))+
  geom_point(size=0.5,aes(color=iden))+
#  geom_point(data=D1.ST.tissuepos,size=1.4,shape=16,aes(fill="black"))+ 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),plot.background=element_blank())+
  coord_fixed()+
  scale_y_reverse()
```

```{r}

torgb=col2rgb(hue_pal()(11))
lookup_table <- data.frame(
  iden = c("0", "1", "2", "3", "4","5","6","7","8","9","10"),
 # iden = c("0", "1", "2", "3", "4","5"),
  R=torgb[1,],
  G=torgb[2,],
  B=torgb[3,]
)
filtered.lipid.fake.spatial<- filtered.lipid.fake.spatial %>%
  left_join(lookup_table, by = "iden")
write.csv(filtered.lipid.fake.spatial,"/scratch/jw7078/D1_lipid_regis.csv")
```

