单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析11

library(GenomicRanges)
##基因组版本转换：Liftover
library(liftOver)
## rtracklayer读取gtf文件
library(rtracklayer)
library(genomation)
library(plyranges)

# Make Ovarian epithelial peakset
files <- list.files(pattern = "iOSE")
##读取narrowpeak文件
grObject.1 <- readNarrowPeak(files[1])
grObject.2 <- readNarrowPeak(files[2])
grObject.3 <- readNarrowPeak(files[3])
grObject.4 <- readNarrowPeak(files[4])

##读取基因组的位置
grObject <- unlist(GenomicRanges::reduce(GRangesList(grObject.1,grObject.2,
                                  grObject.3,grObject.4)))

# Liftover:
chainObject = import.chain("hg19ToHg38.over.chain")
results <- as.data.frame(liftOver(grObject, chainObject))
saveRDS(results,"Ovarian_Epithelial_Cell_line_Peaks.rds")

###################################################
# Matt Regner
# Franco Lab 
# 2020-2021
# Description: plot figures for Fig. 2
###################################################

source("P2G_Heatmap_Distal.R")
source("Archr_Peak_RawPval.R")
source("Archr_Peak_Null_Permute.R")
source("ArchRBrowser.R")
library(ggplot2)
library(Seurat)
library(scales)
library(forcats)
library(RColorBrewer)
library(ArchR)
library(stringr)
library(stringi)
library(dplyr)
library(tidyr)
addArchRThreads(32)
h5disableFileLocking()
####################################################################
# PART 1: summary distribution histograms (correlation and p-values)
# 1) Observed data
# 2) Null data example
####################################################################
# Read in P2G data
p2g <- readRDS("./All_P2G_Observed.rds")
p2g$RawPVal <- as.numeric(p2g$RawPVal)
p2g$Correlation <- as.numeric(p2g$Correlation)

# Plot correlation histogram
p1 <- ggplot(p2g,aes(x=Correlation))+geom_histogram(color="black",fill="gray90",size=0.3,boundary = 0,bins=15)+
  theme_bw()+
  scale_y_continuous(expand = c(0,0),limits=c(0,1500000)) +
  scale_x_continuous(expand= c(0,0),limits = c(-1,1))

# Plot p-value histogram
p2 <- ggplot(p2g,aes(x=RawPVal))+geom_histogram(color="black",fill="gray90",size=0.3,boundary = 0,bins=15)+
  theme_bw()+
  scale_y_continuous(expand = c(0,0),
                     limits = c(0,1.75e+06)) +
  scale_x_continuous(expand = c(0,0),
                     limits = c(0,1))


# Read in P2G data (null example)
p2g <- readRDS("./All_P2G_Null_example.rds")
p2g$RawPVal <- as.numeric(p2g$RawPVal)
p2g$Correlation <- as.numeric(p2g$Correlation)

# Plot correlation histogram
p3 <- ggplot(p2g,aes(x=Correlation))+geom_histogram(color="black",fill="gray90",size=0.3,boundary = 0,bins=15)+
  theme_bw()+
  scale_y_continuous(limits=c(0,1500000),expand=c(0,0)) +
  scale_x_continuous(expand = c(0,0),limits = c(-1,1))

# Plot p-value histogram
p4 <- ggplot(p2g,aes(x=RawPVal))+geom_histogram(color="black",fill="gray90",size=0.3,boundary = 0,bins=15)+
  theme_bw()+
  scale_y_continuous(expand = c(0,0),
                     limits = c(0,1.75e+06)) +
  scale_x_continuous(expand = c(0,0),limits = c(0,1))
# Plot side-by-side
CombinePlots(list(p1,p2,p3,p4),ncol=2)+ggsave("Histograms_obs_null.pdf",height =3,width = 4)
##上面这段代码不错，一张图实现多图叠加



############################################################################
# PART 1.5.0.5: How many cancer specific peaks are deferentially accessible?
############################################################################

# Marker Peaks for malignant clusters with 100% patient specificity 
atac <- readRDS("final_archr_proj_archrGS-P2Gs.rds")
markersPeaks <- getMarkerFeatures(
  ArchRProj = atac, 
  useMatrix = "PeakMatrix", 
  groupBy = "predictedGroup_ArchR",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)
levels(factor(atac$predictedGroup_ArchR))
labels <- c("0-Fibroblast",
            "10-Epithelial cell" ,
            "11-Unciliated epithelia 1",
            "16-Fibroblast",
            "17-Epithelial cell" ,
            "19-Epithelial cell" ,
            "20-Ciliated",
            "21-Unciliated epithelia 1",
            "22-Unciliated epithelia 2",
            "3-Epithelial cell",
            "31-Unciliated epithelia 1",
            "34-Epithelial cell",
            "9-Epithelial cell",
            "27-Fibroblast" 
            
)

markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.1 & Log2FC >= 0.5")
markerList.sub <- markerList[names(markerList) %in% labels]
markers.peaks <- unlist(markerList.sub)
#paste0：默认以空格作为连接字符
markers.peaks <- paste0(markers.peaks$seqnames,":",markers.peaks$start,"-",markers.peaks$end)
markers.peaks <- unique(markers.peaks)

p2g.cancer <- readRDS("./Cancer_specific_P2G_table.rds")
p2g.cancer.diff <- p2g.cancer[p2g.cancer$peakName %in% markers.peaks,]
##Round函数返回一个数值，该数值是按照指定的小数位数进行四舍五入运算的结果
print(paste0(round(nrow(p2g.cancer.diff)/nrow(p2g.cancer)*100,3),"% of cancer-specifc P2Gs have differentially accessible peaks in the malignant cell populations! (FDR <= 0.1 & Log2FC >= 0.5)"))
print(paste0(round(length(unique(p2g.cancer.diff$peakName))/length(unique(p2g.cancer$peakName))*100,3),"% of cancer-specific peaks are differentially accessible in the malignant cell populations! (FDR <= 0.1 & Log2FC >= 0.5)"))

####################################################################
# PART 1.5.1: Investigate the near neighboring gene rule in P2Gs:
####################################################################
p2g <- readRDS("./All_P2G_Observed.rds")
p2g <- p2g[p2g$Correlation >= 0.45,]
p2g <- p2g[p2g$RawPVal <= 1e-12,]
p2g <- p2g[p2g$peakType == "Distal",]# Subset to distal P2Gs
p2g$idx <- paste0(p2g$idxATAC,"-",p2g$idxRNA)
p2g.cancer <- readRDS("./Cancer_specific_P2G_table.rds")
p2g$Cancer.Specific <- ifelse(p2g$idx %in% p2g.cancer$idx,"Cancer","Normal")

p2g$p2g.pair <- paste0(p2g$peakName,"::",p2g$geneName)

atac <- readRDS("final_archr_proj_archrGS-P2Gs.rds")
peak.info <- getPeakSet(atac)
p2g.nearest.pair <- paste0(peak.info@seqnames,":",peak.info@ranges,"::",peak.info$nearestGene)

p2g$Nearest <- ifelse(p2g$p2g.pair %in% p2g.nearest.pair,"Yes","No")

p2g.norm <- p2g[p2g$Cancer.Specific == "Normal",]
p2g.cancer <- p2g[p2g$Cancer.Specific == "Cancer",]

p2g.yes <- p2g[p2g$Nearest == "Yes",]
print(length(unique(p2g.yes$peakName))/length(unique(p2g$peakName)))
print(nrow(p2g.yes)/nrow(p2g))

ggplot(p2g) +
  aes(x = Cancer.Specific, fill = factor(Nearest),boundary=0) +
  geom_bar(position = "fill")+
  scale_y_continuous(expand = c(0,0))+
  scale_x_discrete(expand=c(0,0))+
  theme_classic()+ggsave("ProportionofP2Gs_predicted_by_NearestRule.pdf",width = 4,height = 3)

# Proportions differ significantly by fisher exact test:
cancer.nearest <- p2g.cancer[p2g.cancer$Nearest == "Yes",]
norm.nearest <- p2g.norm[p2g.norm$Nearest == "Yes",]
##Fisher's exact test：是用于分析列联表（ contingency tables ）统计显著性检验方法，它用于检验两个分类的关联（association)。虽然实际中常常使用于小数据情况，但同样适用于大样本的情况。
res <- fisher.test(matrix(c(nrow(cancer.nearest), nrow(norm.nearest), 
                            nrow(p2g.cancer)-nrow(cancer.nearest), nrow(p2g.norm)-nrow(norm.nearest)),
                          ncol = 2))
print(res)
print(paste0("The proportion of cancer-specific P2Gs predictable by nearest neighboring rule is significantly less relative to the normal P2Gs (p=",res$p.value,")"))
####################################################################

##################################################################
# PART 2: Plot proportion of peaks per number of peaks and 
# average number of target genes per cancer and normal groups
####################################################################
# Number of peaks per number of genes
p2g.cancer <- readRDS("./Cancer_specific_P2G_table.rds")
p2g.normal <- readRDS("./All_P2G_Observed.rds")
##我比较好奇这里的阈值标准，相关性的不是很严格，但是后面的pvalue设的很严格
p2g.normal <- p2g.normal[p2g.normal$Correlation >= 0.45,]
p2g.normal <- p2g.normal[p2g.normal$RawPVal <= 1e-12,]
p2g.normal <- p2g.normal[p2g.normal$peakType == "Distal",]
p2g.normal <- p2g.normal[p2g.normal$peakName %ni% p2g.cancer$peakName,]

df.cancer <- data.frame(num.genes = table(p2g.cancer$peakName))
df.cancer$cat <- ifelse(df.cancer$num.genes.Freq < 3,"1-2",df.cancer$num.genes.Freq)
df.cancer$cat <- ifelse(df.cancer$num.genes.Freq >2,"3+",df.cancer$cat)

df.normal <- data.frame(num.genes = table(p2g.normal$peakName))
df.normal$cat <- ifelse(df.normal$num.genes.Freq < 3,"1-2",df.normal$num.genes.Freq)
df.normal$cat <- ifelse(df.normal$num.genes.Freq >2,"3+",df.normal$cat)

head(df.normal)
head(df.cancer)

df.normal$type <- "normal"
df.cancer$type <- "cancer"

##rbind这个函数单细胞中经常用，因为经常要做数据整合
df<- rbind(df.normal,df.cancer)

p1 <- ggplot(df, aes(x=type, y=num.genes.Freq, fill = type)) +
  stat_summary(geom = "bar", fun.y = mean, position = "dodge",width=0.3) +
  stat_summary(geom = "errorbar", fun.data = mean_se, position = "dodge",width=0.15)+
  theme_classic()+  coord_cartesian(ylim=c(0.00,2.0))+NoLegend()+
  scale_y_continuous(expand = c(0,0)) 
##mutate的衍生函数主要是按列对数据赋予function
df.cancer %>% 
  count(cat) %>% 
  mutate(perc = (n / nrow(df.cancer)*100)) -> cancer
df.normal %>% 
  count(cat) %>% 
  mutate(perc = (n / nrow(df.normal)*100)) -> normal

cancer$type <- "cancer"
normal$type <- "normal"

comb <- rbind(cancer,normal)
print(comb)
p2 <- ggplot(comb,aes(x=cat,y=perc,fill=type))+
  geom_bar(stat="identity",position = position_dodge(width=0.75),width=0.6)+
  theme_classic()+
  geom_text(aes(label=round(perc,3)), vjust=-0.4, size=3.5,position = position_dodge(width=0.75))+
  scale_y_continuous(expand = c(0,0),limits=c(0,100))+NoLegend()


CombinePlots(list(p2,p1),ncol=2)+ggsave("Barcharts_P2G_plots.pdf",width=4,height=2)

# Cancer specific peaks link to more genes on average with statistical significance:
test <- wilcox.test(df.cancer$num.genes.Freq,df.normal$num.genes.Freq,correct = F)
print(test)
print(test$p.value)
print(paste0("Cancer-specific peaks link to more genes on average (1.57 v. 1.44 genes, p=",test$p.value,")"))

# Proportion of 3+ peaks is greater in cancer-specific peaks relative to normal peaks:
res <- fisher.test(matrix(c(414, 1979, 
                            3274, 20187),
                          ncol = 2))
print(res)
print(res$p.value)
print(paste0("The proportion of cancer-specific peaks linking to 3 or more genes is significantly greater relative to the normal peaks (p=",res$p.value,")"))

####################################################################
# PART 3: browser track for RHEB enhancers
# 1) Plot browser track
# 2) Verify cancer-specific enhancers have differential accessibility
####################################################################

# Read in other annotation features:
encode.all <- read.delim("./GRCh38-ccREs.bed",header =F)
colnames(encode.all)[1:3] <- c("seqnames","start","end")
##makeGRangesFromDataFrame 将类似数据框的对象作为输入，并尝试自动查找描述基因组范围的列。
encode.all <- makeGRangesFromDataFrame(encode.all)

ft.peaks <- readRDS("./Fallopian_Tube_Cell_line_Peaks.rds")
ft.peaks <- ft.peaks[,3:5]
colnames(ft.peaks)[1:3] <- c("seqnames","start","end")
ft.peaks <- makeGRangesFromDataFrame(ft.peaks)

ov.peaks <- readRDS("./Ovarian_Epithelial_Cell_line_Peaks.rds")
ov.peaks <- ov.peaks[,3:5]
colnames(ov.peaks)[1:3] <- c("seqnames","start","end")
ov.peaks <- makeGRangesFromDataFrame(ov.peaks)

atac <- readRDS("final_archr_proj_archrGS-P2Gs.rds")
# ATAC
levels(factor(atac$predictedGroup_ArchR))
my_levels <- as.character(c(3,9,10,16,17,
                            11,20,21,22,31,19,34,
                            0,27,
                            6,8,12,14,15,18,24,25,26,29,7,23,
                            1,33,
                            2,4,30,
                            5,13,
                            32,
                            28,35))

for ( i in levels(factor(atac$predictedGroup_ArchR))){
  num <-  gsub("-.*","",i)
  idx <- match(num,my_levels)
  atac$predictedGroup_ArchR <- str_replace(atac$predictedGroup_ArchR,pattern = i,replacement = paste0(idx,"_",atac$predictedGroup_ArchR))
  print("iter complete")
}
test <- as.data.frame(atac@cellColData)
test <- test[test$Sample == "38FE7L",]
levels(factor(test$predictedGroup_ArchR))

new.idents <- setdiff(levels(factor(atac$predictedGroup_ArchR)),
                      c("16_8-Stromal fibroblasts",
                        "17_12-Stromal fibroblasts",
                        "18_14-Stromal fibroblasts",
                        "19_15-Stromal fibroblasts",
                        "20_18-Fibroblast",
                        "21_24-Fibroblast",
                        "23_215_23_26-Fibroblast",
                        "24_29-Fibroblast",
                        "26_23-Stromal fibroblasts",
                        "31_30-T cell",
                        "30_4-Lymphocytes",
                        "34_32-Mast cell",
                        "33_13-Macrophages",
                        "36_35-B cell",
                        "32_5-Macrophage",
                        "35_28-B cell" ))

##BiocGenerics:: 用于确定重复元素。 showMethods 用于显示为给定通用函数定义的方法的摘要。
idxSample <- BiocGenerics::which(atac$predictedGroup_ArchR %in% new.idents)
cellsSample <- atac$cellNames[idxSample]
atac.sub <- atac[cellsSample, ]

#########################################################################
# Make modified getP2G function:
##从 ArchRProject 获取峰到基因的链接
getPeak2GeneLinks.mod <- function(
  ArchRProj = NULL, 
  corCutOff = 0.45, 
  PValCutOff = 0.0001,
  varCutOffATAC = 0.25,
  varCutOffRNA = 0.25,
  resolution = 1, 
  returnLoops = TRUE
){
  
  .validInput(input = ArchRProj, name = "ArchRProj", valid = "ArchRProject")
  .validInput(input = corCutOff, name = "corCutOff", valid = "numeric")
  .validInput(input = PValCutOff, name = "PValCutOff", valid = "numeric")
  .validInput(input = varCutOffATAC, name = "varCutOffATAC", valid = "numeric")
  .validInput(input = varCutOffRNA, name = "varCutOffRNA", valid = "numeric")
  .validInput(input = resolution, name = "resolution", valid = c("integer", "null"))
  .validInput(input = returnLoops, name = "returnLoops", valid = "boolean")
  
  if(is.null(ArchRProj@peakSet)){
    return(NULL)
  }
  
  if(is.null(metadata(ArchRProj@peakSet)$Peak2GeneLinks)){
    
    return(NULL)
    
  }else{
    
    p2g <- metadata(ArchRProj@peakSet)$Peak2GeneLinks
    p2g <- p2g[which(p2g$Correlation >= corCutOff & p2g$RawPVal <= PValCutOff), ,drop=FALSE]
    
    if(!is.null(varCutOffATAC)){
      p2g <- p2g[which(p2g$VarQATAC > varCutOffATAC),]
    }
    
    if(!is.null(varCutOffRNA)){
      p2g <- p2g[which(p2g$VarQRNA > varCutOffRNA),]
    }
    
    if(returnLoops){
      
      peakSummits <- resize(metadata(p2g)$peakSet, 1, "center")
      geneStarts <- resize(metadata(p2g)$geneSet, 1, "start")
      
      if(!is.null(resolution)){
        summitTiles <- floor(start(peakSummits) / resolution) * resolution + floor(resolution / 2)
        geneTiles <- floor(start(geneStarts) / resolution) * resolution + floor(resolution / 2)
      }else{
        summitTiles <- start(peakSummits)
        geneTiles <- start(geneTiles)
      }
      
      loops <- .constructGR(
        seqnames = seqnames(peakSummits[p2g$idxATAC]),
        start = summitTiles[p2g$idxATAC],
        end = geneTiles[p2g$idxRNA]
      )
      mcols(loops)$value <- p2g$Correlation
      mcols(loops)$FDR <- p2g$FDR
      
      loops <- loops[order(mcols(loops)$value, decreasing=TRUE)]
      loops <- unique(loops)
      loops <- loops[width(loops) > 0]
      loops <- sort(sortSeqlevels(loops))
      
      loops <- SimpleList(Peak2GeneLinks = loops)
      
      return(loops)
      
    }else{
      
      return(p2g)
      
    }
    
  }
  
}
##这里只能说写的一手好循环嵌套


#########################################################################
# Color rows:

epithelial.cols <- colorRampPalette(c("#A0E989", "#245719"))
epithelial.cols <- epithelial.cols(14)

fibro.cols <-colorRampPalette(c("#FABFD2", "#B1339E"))
fibro.cols <- fibro.cols(10)

smooth.cols <- c("#b47fe5","#d8b7e8")

endo.cols <- c("#93CEFF","#4A99FF")

t.cols <- c("gray60","gray20","gray40")

macro.cols <- c("#ff6600","#ff9d5c")

mast.cols <- "gold3"

b.cols <- c("#B22222","#CD5C5C")


cols <- c(epithelial.cols,fibro.cols,smooth.cols,endo.cols,t.cols,macro.cols,mast.cols,b.cols)
cols <- cols[-c(16:21,23:24,26,30:36)]
cols <- cols[c(8:12,1:7,13:20)]
##plotBrowserTrack：此函数将以浏览器轨迹的样式绘制输入区域的覆盖率。它允许对信号进行归一化，从而实现跨样本的直接比较。
plot <- plotBrowserTrack(atac.sub,geneSymbol ="RHEB", groupBy = "predictedGroup_ArchR",
                         features = GRangesList(TrackA = encode.all,TrackB = ft.peaks,TrackC = ov.peaks), 
                         loops = getPeak2GeneLinks.mod(atac,corCutOff = 0.45,
                                                       PValCutOff = 1e-12,varCutOffATAC = 0,
                                                       varCutOffRNA = 0),upstream = 6000,downstream = 35000,
                         pal=cols)

pdf("RHEB_final.pdf",width = 6,height = 8)
#grid.draw(x, recording= TRUE) Arguments x 类 "grob" 或 NULL 的对象。 recording 一个逻辑值,表示是否应该在网格显示列表中记录绘图操作
grid::grid.draw(plot[[1]])
dev.off()


names <- gsub(".*_","",atac.sub$predictedGroup_ArchR)
saveRDS(names,"names.rds")
rm(atac.sub)
rm(atac)


# Plot MUC16 and CD117 examples:
atac <- readRDS("final_archr_proj_archrGS-P2Gs.rds")
# ATAC
levels(factor(atac$predictedGroup_ArchR))
my_levels <- as.character(c(11,20,21,22,31,19,34,
                            3,9,10,16,17,
                            0,27,
                            6,8,12,14,15,18,24,25,26,29,7,23,
                            1,33,
                            2,4,30,
                            5,13,
                            32,
                            28,35))

for ( i in levels(factor(atac$predictedGroup_ArchR))){
  num <-  gsub("-.*","",i)
  idx <- match(num,my_levels)
  atac$predictedGroup_ArchR <- str_replace(atac$predictedGroup_ArchR,pattern = i,replacement = paste0(idx,"_",atac$predictedGroup_ArchR))
  print("iter complete")
}

epithelial.cols <- colorRampPalette(c("#A0E989", "#245719"))
epithelial.cols <- epithelial.cols(14)

fibro.cols <-colorRampPalette(c("#FABFD2", "#B1339E"))
fibro.cols <- fibro.cols(10)

smooth.cols <- c("#b47fe5","#d8b7e8")

endo.cols <- c("#93CEFF","#4A99FF")

t.cols <- c("gray60","gray20","gray40")

macro.cols <- c("#ff6600","#ff9d5c")

mast.cols <- "gold3"

b.cols <- c("#B22222","#CD5C5C")


cols <- c(epithelial.cols,fibro.cols,smooth.cols,endo.cols,t.cols,macro.cols,mast.cols,b.cols)
plot <- plotBrowserTrack(atac,geneSymbol ="MUC16", groupBy = "predictedGroup_ArchR",
                         features = GRangesList(TrackA = encode.all,TrackB = ft.peaks,TrackC = ov.peaks), 
                         loops = getPeak2GeneLinks.mod(atac,corCutOff = 0.45,
                                                       PValCutOff = 1e-12,varCutOffATAC = 0,
                                                       varCutOffRNA = 0),upstream = 75000,downstream = 110000,
                         pal=cols,
                         ylim=.995)

pdf("MUC16_final.pdf",width = 6,height = 8)
grid::grid.draw(plot[[1]])
dev.off()

plot <- plotBrowserTrack(atac,geneSymbol ="KIT", groupBy = "predictedGroup_ArchR",
                         features = GRangesList(TrackA = encode.all,TrackB = ft.peaks,TrackC = ov.peaks), 
                         loops = getPeak2GeneLinks.mod(atac,corCutOff = 0.45,
                                                       PValCutOff = 1e-12,varCutOffATAC = 0,
                                                       varCutOffRNA = 0),upstream = 85000,downstream = 30000,
                         pal=cols)

pdf("KIT_final.pdf",width = 6,height = 8)
grid::grid.draw(plot[[1]])
dev.off()


plot <- plotBrowserTrack(atac,geneSymbol = "LAPTM4B", groupBy = "predictedGroup_ArchR",
                         features = GRangesList(TrackA = encode.all,TrackB = ft.peaks,TrackC = ov.peaks), 
                         loops = getPeak2GeneLinks.mod(atac,corCutOff = 0.45,
                                                       PValCutOff = 1e-12,varCutOffATAC = 0,
                                                       varCutOffRNA = 0),
                         pal=cols,
                         upstream = 41500,
                         downstream = 10000
)


pdf("LAPTM4B_final.pdf",width = 6,height = 8)
grid::grid.draw(plot[[1]])
dev.off()


####################################################################
# PART 4: plot matching violin plots for RHEB expression and mTOR
####################################################################
names <- readRDS("./names.rds")
rna <- readRDS("./endo_ovar_All_scRNA_processed.rds")
rna.sub <- rna[,rna$cell.type %in% intersect(levels(factor(rna$cell.type)),levels(factor(names)))]

my_levels <- as.character(c(3,9,10,16,17,
                            11,20,21,22,31,19,34,
                            0,27,
                            6,25,7,
                            1,33,
                            2))

# Make violin plots for RHEB expression and mTOR pathway expression
# Relevel object@ident
rna.sub@active.ident <- factor(x =rna.sub$RNA_snn_res.0.7, levels = rev(my_levels))
p1 <- VlnPlot(rna.sub,features = "RHEB",pt.size = 0)+coord_flip()+NoLegend()
p1 <- ggplot(p1$data,aes(y=ident,x=RHEB))+geom_boxplot(aes(fill=ident),lwd=0.45,outlier.size = 0.95,fatten = 0.95)+NoLegend()+
  theme_classic()+NoLegend()+ylab("Cluster #")+xlab("RHEB expression")


# Compute mTOR pathway member expression
gset <- read.delim("./BIOCARTA_MTOR_PATHWAY.txt",header = T)
gset <- as.character(gset$BIOCARTA_MTOR_PATHWAY[2:length(gset$BIOCARTA_MTOR_PATHWAY)])
rna.sub <- AddModuleScore(rna.sub,features = list(gset),name = "mTOR_members",search = T)


p2 <- VlnPlot(rna.sub,features = "mTOR_members1",pt.size = 0)+coord_flip()+NoLegend()
p2 <- ggplot(p2$data,aes(y=ident,x=mTOR_members1))+geom_boxplot(aes(fill=ident),lwd=0.45,outlier.size = 0.95,fatten = 0.95)+NoLegend()+
  theme_classic()+NoLegend()+ylab("Cluster #")+xlab("mTOR member expression")
CombinePlots(list(p1,p2),ncol=2)+ggsave("VlnPlots.pdf",width = 6,height = 8)

# Differential enrichment of mTOR pathway members
res <- kruskal.test(data=p2$data,mTOR_members1 ~ ident)
print(res)
print(res$p.value)
# Differential expression of RHEB in cluster 3
grouped.markers <- FindMarkers(rna.sub,ident.1 = "3",
                               ident.2 = as.character(c(9,10,16,17,
                                                        11,20,21,22,31,19,34,
                                                        0,27,
                                                        6,25,7,
                                                        1,33,
                                                        2)),only.pos = T)
grouped.markers$gene <- rownames(grouped.markers)
grouped.markers.RHEB <- grouped.markers[grouped.markers$gene == "RHEB",]
print(grouped.markers.RHEB)# Significant up-regulation in cluster 3


# make boxplots for MUC16 and KIT expression
rna <- readRDS("./endo_ovar_All_scRNA_processed.rds")

my_levels <- as.character(c(11,20,21,22,31,19,34,
                            3,9,10,16,17,
                            0,27,
                            6,8,12,14,15,18,24,25,26,29,7,23,
                            1,33,
                            2,4,30,
                            5,13,
                            32,
                            28,35))

# Make violin plots for MUC16 expression and mTOR pathway expression
# Relevel object@ident
rna@active.ident <- factor(x =rna$RNA_snn_res.0.7, levels = rev(my_levels))
p1 <- VlnPlot(rna,features = "MUC16",pt.size = 0)+coord_flip()+NoLegend()
p1 <- ggplot(p1$data,aes(y=ident,x=MUC16))+geom_boxplot(aes(fill=ident),lwd=0.45,outlier.size = 0.95,fatten = 0.95)+NoLegend()+
  theme_classic()+NoLegend()+ylab("Cluster #")+xlab("MUC16 expression")+ggsave("Vln-MUC16.pdf",width=3,height = 8)

# Make violin plots for KIT expression and mTOR pathway expression
# Relevel object@ident
rna@active.ident <- factor(x =rna$RNA_snn_res.0.7, levels = rev(my_levels))
p1 <- VlnPlot(rna,features = "KIT",pt.size = 0)+coord_flip()+NoLegend()
p1 <- ggplot(p1$data,aes(y=ident,x=KIT))+geom_boxplot(aes(fill=ident),lwd=0.45,outlier.size = 0.95,fatten = 0.95)+NoLegend()+
  theme_classic()+NoLegend()+ylab("Cluster #")+xlab("KIT expression")+ggsave("Vln-KIT.pdf",width=3,height = 8)

# Make violin plots for LAPTM4B expression and mTOR pathway expression
# Relevel object@ident
rna@active.ident <- factor(x =rna$RNA_snn_res.0.7, levels = rev(my_levels))
p1 <- VlnPlot(rna,features = "LAPTM4B",pt.size = 0)+coord_flip()+NoLegend()
p1 <- ggplot(p1$data,aes(y=ident,x=LAPTM4B))+geom_boxplot(aes(fill=ident),lwd=0.45,outlier.size = 0.95,fatten = 0.95)+NoLegend()+
  theme_classic()+NoLegend()+ylab("Cluster #")+xlab("LAPTM4B expression")+ggsave("Vln-LAPTM4B.pdf",width=3,height = 8)


writeLines(capture.output(sessionInfo()), "sessionInfo.txt")