require(Seurat)
require(magrittr)
require(decoupleR)
require(reshape2)
require(tidyverse)
require(tidyr)
require(broom)
require(limma)
require(ggsci)
require(fossil)
require(clusterProfiler)
require(msigdbr)
require(plotly)
require(destiny)
require(tidyverse)
require(ggplot2)
require(pryr)
source("script/20240910.R")
#source("../bryony/manuscript/script.R")


find.npc <- function(sd){
  k <- length(sd)
  s <- abs(diff(sd))
  mu <- mean(s[(0.8*k):k-1])
  sigma <- sd(s[(0.8*k):k-1])
  sk <- mu + 6*sigma
  npcs <- max(which(s > sk))+1
  npcs
}
################################# Load data and Seurat pre-processing #################################
tmp.obj <- readRDS("data/sclc/1ba7d495-c1a8-4809-b56d-548fbea77c8a.rds")
meta.features <- tmp.obj@assays$RNA@meta.features
rownames(meta.features)<-  as.character(meta.features$feature_name)
tmp.obj@assays$RNA@counts@Dimnames[[1]] <- rownames(meta.features)
tmp.obj@assays$RNA@data@Dimnames[[1]] <-  rownames(meta.features)
obj <- CreateSeuratObject(counts = tmp.obj@assays$RNA@counts, meta.data = tmp.obj@meta.data)
LayerData(obj, "data") <- tmp.obj@assays$RNA$data
obj@reductions <- tmp.obj@reductions
obj@misc <- tmp.obj@misc
rm(tmp.obj)

# Seurat pipeline
obj.mac <- subset(obj, cell_type == "macrophage")
obj.mac <- obj.mac[-grep("ENSG", rownames(obj.mac)),]
obj.mac <- FindVariableFeatures(obj.mac)
obj.mac <- ScaleData(obj.mac)
obj.mac <- RunPCA(obj.mac, npcs = 100)
obj.mac <- RunUMAP(obj.mac, dims = 1:10)

# Save the result
saveRDS(obj.mac, "result/obj.sclc.mac.rds")

################################# Supplementary Figure A Dimplot full data #################################
cols <- rep("lightgrey", 11)
names(cols) <- levels(obj$cell_type)
cols["macrophage"] <- "darkred"

# Supplementary Figure A
p <- DimPlot(obj, group.by = "cell_type") + theme(text = element_text(size = 13, face = "bold"),
                                                  legend.position = "none",
                                                  axis.text = element_blank(), axis.ticks = element_blank()
) + xlab("UMAP1")+ylab("UMAP2")+
  scale_color_manual(values = cols)+
  ggtitle("SCLC atlas")

p <- LabelClusters(p, id = "cell_type", fontface = "bold", color = "black",size = 3.5, box = 1.3, alpha=0.8)
p

################################# Figure B Dimplot macrophage subset #################################
obj.mac <- readRDS("result/obj.sclc.mac.rds")
p1 <- DimPlot(obj.mac, group.by = "disease", pt.size = 0.6) + theme(text = element_text(face = "bold", size = 13),
                                                                    legend.position = "none",
                                                                    axis.text = element_blank(), axis.ticks = element_blank()
) + xlab("UMAP1")+ylab("UMAP2")+
  scale_color_manual(values = c(ggsci::pal_bmj()(9)[1], "red2", "lightgrey"))+
  ggtitle("Macrophge subset")

p1 <- LabelClusters(p1, id = "disease", fontface = "bold", color = "black",size = 4, box = 1.3, alpha=0.8)
p1

################################# Figure C Differential expression analysis #################################
obj.mac <- readRDS("result/obj.sclc.mac.rds")
Idents(obj.mac) <- "disease"
fc <- FoldChange(obj.mac, ident.1 = "small cell lung carcinoma", ident.2 = "normal" )
fc <- subset(fc, (pct.1 > 0.1) & (pct.1 - pct.2 > 0))
test.result <- FindMarkers(obj.mac, ident.1 = "small cell lung carcinoma", ident.2 = "normal", features = rownames(fc), test.use = "MAST")
test.result$pi <- -log10(test.result$p_val_adj) * (test.result$avg_log2FC)
test.result <- arrange(test.result, desc(pi))
write.table(test.result, "result/sclc.de.txt")

test.result$gene <- rownames(test.result)
ggplot(data = test.result) + geom_point(aes(x = avg_log2FC, y = pi, col = p_val_adj < 0.05)) +
  geom_point(data = head(test.result, 50), aes(x = avg_log2FC, y = pi),col = "red3", size = 1)+
  geom_text_repel(data = head(test.result, 50), aes(x = avg_log2FC, y = pi, label = gene),
                  size = 3.5, fontface = "italic")  +
  scale_color_manual("Adj p-val < 0.05", values = c("lightgrey", "red3"))+
  theme_classic()+
  theme(text = element_text(face = "bold", size =13), legend.position = "top") +
  xlab("Average Log2 FC") + ylab("Pi")

################################# NNet analysis #################################
cells <- readLines("result/sclc.cells")
test.result <- read.table("result/sclc.de.txt")
obj.mac <- readRDS("result/obj.sclc.mac.rds")
rt.ppr <- get.ppr()
genes <- select.gene(obj.mac, min.cells = 20)
obj.mac <- prepare.seurat(obj.mac, genes = genes$genes)

mem_change({
  obj.mac <- prepare.graph(obj.mac)
  # obj.mac <- select.cell(obj.mac, n = 1500)
  obj.mac@misc$setting$cells <- as.numeric(cells)
  names(obj.mac@misc$setting$cells) <- colnames(obj.mac)[as.numeric(cells)]
  t1 <- system.time({
    obj.mac <- prepare.reg(obj.mac, predictors = genes$tfs,
                           responses = genes$targets, check.expressed = T)
  })
  responses <- head(rownames(test.result), 50)
  t2 <- system.time({
    obj.mac <- run.nn.reg(obj.mac, responses = responses, return.p.val = T)
  })
  #obj.mac <- set.defaults(obj.mac, defaults = list(cutoff = 0.99))
  t3 <- system.time({
    obj.mac <- build.meta.network(obj.mac)
  })
})
t1
t2
t3

# Select central genes to visualise co-expression
central.genes  <- select.central.genes(obj.mac, keep.responses = T, n.per.component = 20, n.net = 3)
obj.mac <- prepare.visualise(obj.mac, central.genes = central.genes, n.clu = 6)
cells <- obj.mac@misc$setting$cells

# Save cells used in the analysis
# writeLines(obj.mac@misc$setting$cells %>% as.character, "result/sclc.cells")

################################# Figure D Featureplot with cells being colored by their nPCA weighting #################################
ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = "MN 3"), size = 0.8, alpha = 10*obj.mac@misc$mod$meta.network$npca.loadings[,1])+
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = "MN 1"), size = 0.8, alpha = 10*obj.mac@misc$mod$meta.network$npca.loadings[,2])+
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = "MN 2"), size = 0.8, alpha = 10*obj.mac@misc$mod$meta.network$npca.loadings[,3])+
  theme_classic()+
  scale_color_manual("", values = c("#3182BD", "#31A354", "darkred"))+
  theme(text = element_text(face = "bold", size = 13),
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5),
        legend.position = "bottom"
  )+ xlab("")+ylab("") +ggtitle("NMF weighting")+guides(colour = guide_legend(override.aes = list(size=4)))

RColorBrewer::brewer.pal(n = 5, "Purples")[4]
RColorBrewer::brewer.pal(n = 5, "Greens")[4]
RColorBrewer::brewer.pal(n = 5, "Blues")[4]

# Cell selection (discarded)
plot_ly() %>% add_markers(data = data.frame(obj.mac@reductions$umap@cell.embeddings),
                          x = ~umap_1, y = ~umap_2) %>%
  add_markers(data = data.frame(obj.mac@reductions$umap@cell.embeddings[cells,]),
              x = ~umap_1, y = ~umap_2, color = obj.mac@misc$mod$meta.network$npca.loadings[,3],  text = cells)

################################# Figure E Network plot #################################
# Meta network (discarded)
visualise.network(obj.mac,2, meta.network = T, cutoff = 0.25, swap.layers = F, show.pathways = T, n.clu = 6,
                  radius = c(0.4,0.7,0.9,1),pie.radius = 0.03, text.size = 4)

# Cell-specific network of the most representative cell of meta-network 1
cells[which.max(obj.mac@misc$mod$meta.network$npca.loadings[,3])]
visualise.network(obj.mac, 8700, meta.network = F, cutoff = 0.5, swap.layers = F, show.pathways = F, n.clu = 6,pie.radius = 0.04)
visualise.network(obj.mac, 8700  , meta.network = F, cutoff = 0.25, swap.layers = F, show.pathways = T, n.clu = 6,
                  radius = c(0.4,0.7,0.9,1),pie.radius = 0.035, text.size = 4.5, n.per.component = 20)

################################# Supplementary Figure C1,2 Network plot #################################
# Cell-specific network of the most representative cell of meta-network 2
cells[which.max(obj.mac@misc$mod$meta.network$npca.loadings[,1])]
visualise.network(obj.mac, 8566, meta.network = F, cutoff = 0.25, swap.layers = F, show.pathways = F, n.clu = 6,pie.radius = 0.04)
visualise.network(obj.mac, 8566, meta.network = F, cutoff = 0.5, swap.layers = F, show.pathways = T, n.clu = 6,
                  radius = c(0.4,0.7,0.9,1),pie.radius = 0.035, text.size = 4.5, n.per.component = 20)

# Cell-specific network of the most representative cell of meta-network 3
cells[which.max(obj.mac@misc$mod$meta.network$npca.loadings[,2])]
visualise.network(obj.mac, 5554, meta.network = F, cutoff = 0.25, swap.layers = F, show.pathways = F, n.clu = 6,pie.radius = 0.04)
visualise.network(obj.mac, 5554   , meta.network = F, cutoff = 0.5, swap.layers = F, show.pathways = T, n.clu = 6,
                  radius = c(0.4,0.7,0.9,1),pie.radius = 0.035, text.size = 4, n.per.component = 20)

################################# Figure F Ternary plot of receptor activity #################################
act <- receptor.activity(obj.mac, meta.network = F, cutoff = 0.84)
# Extract receptor activity in individual cells
receptor.act <- scale(as.matrix(act$receptor.act), center =F, scale = T)
receptor.act[is.na(receptor.act)] <- 0
# Aggregate receptor activity
receptor.act <-  receptor.act %*% obj.mac@misc$mod$meta.network$npca.loadings[,c(1,2,3)]
receptor.act <- data.frame(receptor.act)
names(receptor.act) <- c("MN3","MN1","MN2")
receptor.act$receptors <- rownames(act$receptor.act)
ggtern::geom_point_swap()
ggtern::ggtern(receptor.act, aes(x=`MN1`, y=`MN3`, z=`MN2`, label = `receptors`)) +
  geom_point(col = "lightgrey", size = 1)+
  geom_point(data=receptor.act[c("NOTCH1", "NOTCH2", "NOTCH3","NOTCH4"),,drop=F], col = "red", size =2) +
  ggtern::theme_rgbw() + theme(text = element_text(size = 13, face = "bold"))



# ggtern::geom_(aes(x=`MN1`, y=`MN3`, z=`MN2`), col = "lightgrey", size = 1) +
# geom_text(aes(x=`MN1`, y=`MN3`, z=`MN2`, label = `receptors`),
#           check_overlap = T, hjust = -0.1, vjust = -0.2, col = "darkblue", size = 3.2) +
# geom_point(data=receptor.act[c("NOTCH1", "NOTCH2", "NOTCH3","NOTCH4"),,drop=F],
#            aes(x=`MN1`, y=`MN3`, z=`MN2`), col = "red", size =2) +
# ggtern::theme_rgbw() + theme(text = element_text(size = 13, face = "bold"))

################################# Figure G, H NOTCH2 and RBPJ activity #################################
# Extract receptor activity
receptor.act <- scale(as.matrix(act$receptor.act), center =F, scale=T)
receptor.act[is.na(receptor.act)] <- 0

# Extract TF activity
tf.act <- scale(as.matrix(act$tf.act), center =F, scale=T)
tf.act[is.na(tf.act)] <- 0
target.act <- scale(as.matrix(act$target.act), center =F, scale=T)
target.act[is.na(target.act)] <- 0

# Calculate correlation between receptor activity and TF activity
top.tfs <- cor(t(act$receptor.act["NOTCH2",]), t(act$tf.act))
top.tfs[is.na(top.tfs)] <- 0
top.tfs <- top.tfs[,order(top.tfs, decreasing = T)][1:20]

# Calculate correlation between receptor activity and target activity (discarded)
top.targets <- cor(t(act$receptor.act["NOTCH2",]), t(act$target.act))
top.targets[is.na(top.targets)] <- 0
top.targets <- top.targets[,order(top.targets, decreasing = T)][1:20]

################################# Figure G NOTCH2 activity #################################
# Featureplot of NOTCH2 activity
ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = as.numeric(receptor.act["NOTCH2",] )), size = 0.8)+
  theme_classic()+
  scale_color_viridis()+
  theme(text = element_text(face = "bold", size = 13),
        legend.position = "none",
        plot.title = ggtext::element_markdown(hjust = 0.5),
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank() ,
  )+ xlab("")+ylab("") + ggtitle("*NOTCH2* (receptor) activity")

################################# Figure H RBPJ activity #################################
# Featureplot of RBPJ activity
ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = as.numeric(tf.act["RBPJ",] )), size = 1)+
  theme_classic()+
  scale_color_viridis()+
  theme(text = element_text(face = "bold", size = 13),
        legend.position = "none",
        plot.title = ggtext::element_markdown(hjust = 0.5),
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank() ,
  )+ xlab("")+ylab("") + ggtitle("*RBPJ* (TF) activity")

################################# NNet with RBPJ as response, TFs as predictors ################################
obj.mac <- prepare.reg(obj.mac, predictors = genes$targets,
                       responses = genes$tfs)
obj.mac <- run.nn.reg(obj.mac, responses = "RBPJ", return.p.val = T)
obj.mac <- build.meta.network(obj.mac)

top.targets <- cor(receptor.act["NOTCH2",] %>% matrix, get.network(obj.mac) %>% t)[1,]
top.targets <- sort(top.targets, decreasing = T) %>% head(., n=30)
top.targets.act <- get.gr.adj(t = 2)["RBPJ",names(top.targets)]
plot.df <- data.frame(gene = names(top.targets), cor = top.targets, activity = top.targets.act, act.size = abs(top.targets.act))
plot.df$activity <-  sapply(plot.df$activity, function(x) if(x > 0) "stimualte" else if(x < 0) "suppress" else "none")
plot.df$gene <- factor(plot.df$gene, levels=plot.df$gene %>% rev)

ggplot(plot.df) + geom_point(aes(cor, gene, col = activity, size = act.size)) +
  theme_bw()+
  theme(text = element_text(face = "bold", size = 13),
        plot.title = element_text(hjust = 0.5),
        axis.text.y = element_text(face = "bold.italic"),
  )+ xlab("Correlation")+ylab("Targets") +
  ggtitle("NOTCH2 activity vs RBPJ-target coexpression")+
  scale_color_manual("Act direction", values = c("orange","red2","blue2")) +
  scale_size_continuous("Act strength")

################################# Deprecated analysis showing meta-network  ################################
central.genes  <- select.central.genes(obj.mac, keep.responses = T, n.per.component = 30, n.net = 1)
central.genes$central.predictors <- tail(top.targets, n = 30) %>% names
central.genes$predictor.module <- rep("M1", 30) %>% as.factor
obj.mac <- prepare.visualise(obj.mac, central.genes = central.genes, n.clu = 1, as.g2 = "response")
visualise.network(obj.mac, 1, meta.network = T, cutoff = 0.25, swap.layers = F, show.pathways = F, n.clu = 1)
visualise.network(obj.mac, 5302, meta.network = F, cutoff = 0.25, swap.layers = F, show.pathways = F, n.clu = 6,pie.radius = 0.04)

################################# Supplementary Figure D Dotplot ################################
features <- c("DLL1", "DLL4", "JAG1", "JAG2", "CCR2", "CCR3", "CCR5", "CCR4")
names(features) <- c(rep("NOTCH signaling",4), rep("CCL13 signaling",4))
Seurat::DotPlot(subset(obj, disease != "normal"), features, group.by = "cell_type", ) +
  theme(text = element_text(size = 13, face = "bold"),
        axis.text.x = element_text(angle = 90),
  ) +
  scale_color_gradientn(colours = gg.color.spec(11))+ylab("")

################################# Deprecated test code ################################
FeaturePlot(obj, "IL1A")+
  scale_color_gradientn(colors =gg.color.spec(11))

FeaturePlot(obj.mac, "HES1")+
  scale_color_gradientn(colors =gg.color.spec(11))

obj.mac <- run.nn.reg(obj.mac, responses = "SPP1", return.p.val = T)

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[as.numeric(cells),],
             aes(umap_1, umap_2, col = obj.mac@misc$mod$p.val["SPP1","NFKB1",]), size = 1)+
  theme_classic()+
  scale_color_gradientn(colors =gg.color.spec(11))+
  theme(text = element_text(face = "bold", size = 15),
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj.mac@misc$mod$p.val["SPP1","KLF4",] > 0.95), size = 1)+
  theme_classic()+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj.mac@misc$setting$nn.scale.gene["SPP1",]), size = 1)+
  theme_classic()+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj.mac@misc$setting$lra[cells,"RBPJ"]), size = 1)+
  theme_classic()+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(aes(obj.mac@misc$setting$lra[cells,"SPP1"], obj.mac@misc$setting$lra[cells,"RPBJ"], col = obj.mac@misc$mod$p.val["SPP1","RPBJ",]))

FeaturePlot(obj.mac, "HIF1A")+
  scale_color_gradientn(colors =gg.color.spec(11))

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = get.network(obj.mac, responses = "RBPJ",assay = "p.val", remove.self.loops = F)["SELENOP",]), size = 1)+
  theme_classic()+
  scale_color_gradientn(colors =gg.color.spec(11), limits = c(0,1))+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj.mac@misc$mod$meta.network$npca.loadings[,1]), size = 1)+
  theme_classic()+
  scale_color_gradientn(colors =gg.color.spec(11))+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

FeaturePlot(obj.mac, "KLF4",pt.size = 2)+
  scale_color_gradientn(colors =gg.color.spec(11))



p <- "GAPDH"
new.obj <- run.nn.reg(obj.mac, p, return.p.val = T)
mu <- new.obj@misc$mod$mus[p]
sigma <- new.obj@misc$mod$sigmas[p]

FeaturePlot(new.obj, p, pt.size = 1, order = T)+xlab("UMAP1")+ylab("UMAP2")+
  theme_classic()+
  scale_color_gradientn(colors =gg.color.spec(11))+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() + geom_point(data = obj.mac@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj.mac@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = new.obj@misc$mod$p.val[p,p,] > 0.95), size = 1)+
  theme_classic()+
  theme(text = element_text(face = "bold", size = 15),
        legend.position = "none",
        axis.line = element_blank(),
        axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(hjust = 0.5)
  )+ xlab("")+ylab("")

ggplot() +
  geom_jitter(aes(x= "Response", y= log(abs(new.obj@misc$mod$effect[p,p,]))), size= 1) +
  geom_errorbar(aes(x = "Null distribution", ymin = mu-1.64*sigma, ymax =mu + 1.64*sigma), color = "darkred",size = 2)+
  theme_classic()+
  theme(text = element_text(size = 15, face = "bold"),
        axis.ticks = element_blank(), legend.position = "none") + xlab("Co-expression type") + ylab("Log absolute co-expression")

data <- cbind(t(obj.mac@assays$RNA$scale.data), matrix(rnorm(9251), nrow = 9251, ncol = 1))
data <- cbind(t(obj.mac@assays$RNA$scale.data))
#simulated.loading <- RSpectra::svds(data, k = 67)


data <- cbind(t(obj.mac@assays$RNA$scale.data), matrix(rnorm(9251), nrow = 9251, ncol = 1))
simulated.loading <- RSpectra::eigs(crossprod(data)/9251, k = 67)
plot(simulated.loading$vectors[5245,])
plot(simulated.loading$vectors[5245,] * sqrt((simulated.loading$values-1)^2/simulated.loading$values))

simulated.loading <- RSpectra::svds(data, k = 67)
(simulated.loading$d)^2/9251

#simulated.loading$d

plot(simulated.loading[,2], obj.mac@misc$setting$pcs[,2])

sds <- apply(obj.mac@misc$setting$pcs, 2, var)
plot(simulated.loading[5245,]*(sds-1))
RSpectra::svds(simulated.loading, k = 67)$u

plot(simulated.loading[1,])
plot(rowSums(simulated.loading^2))
rank(rowSums(simulated.loading^2)) %>% last
plot(simulated.loading[5245,]*sds)
hist(simulated.loading[,1])
hist(simulated.loading[,65])
dim(simulated.loading)

ks.test(simulated.loading[,65], simulated.loading[,64])
dim(simulated.loading)
sds <- apply(obj.mac@misc$setting$pcs, 2, sd)

