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

library(ArchR)
library(ComplexHeatmap)
library(DESeq2)
library(Seurat)
library(ggplot2)
library(dplyr)

###############################################
# Part 1: Read in TF motif prediction matrix 
###############################################
# Read in TF motif prediction matrix
enhancer.motifs<- read.delim("./motif_enhancers_scaled.txt",sep = "\t")
rownames(enhancer.motifs) <- enhancer.motifs[,1]
enhancer.motifs <- enhancer.motifs[,-1]
# Rename syntax of peaks
colnames(enhancer.motifs) <- sub("\\.",":",colnames(enhancer.motifs))
colnames(enhancer.motifs) <- sub("\\.","-",colnames(enhancer.motifs))


# Set labels of cell types of interest

# Only use clusters with 100% patient specificity 
labels <- c("31-Unciliated epithelia 1",
            "21-Unciliated epithelia 1",
            "19-Epithelial cell",
            "34-Epithelial cell",
            "3-Epithelial cell",
            "10-Epithelial cell",
            "16-Fibroblast",
            "17-Epithelial cell",
            "0-Fibroblast",
            "27-Fibroblast",
            "9-Epithelial cell")# Screen for "pure" patient clusters


###############################################
# Part 2: Pseudobulk TF expression 
###############################################

# Pseudobulk RNA TF expression
rna <- readRDS("./endo_ovar_All_scRNA_processed.rds")
rna.counts <- rna@assays$RNA@counts

rna.pseudobulk <- data.frame(rownames(rna))
for (i in labels){
  cells <- rownames(dplyr::filter(rna@meta.data,cell.type == i))
  
  rna.counts.sub <- rna.counts[,colnames(rna.counts) %in% cells]
  
  rna.counts.bulk <- rowSums(as.matrix(rna.counts.sub))
  
  rna.pseudobulk$i <- rna.counts.bulk
  colnames(rna.pseudobulk)[dim(rna.pseudobulk)[2]] <- i
  
}
rownames(rna.pseudobulk) <- rna.pseudobulk[,1]
rna.pseudobulk <- rna.pseudobulk[,-1]
dim(rna.pseudobulk)

# Remove any features/rows that have zero counts in ANY sample
dim(rna.pseudobulk)
rna.pseudobulk[rna.pseudobulk == 0] <- NA
rna.pseudobulk <- rna.pseudobulk[complete.cases(rna.pseudobulk),]
dim(rna.pseudobulk)


# rlog normalize data
library(DESeq2)
dds.rna <- DESeqDataSetFromMatrix(countData = rna.pseudobulk,
                              colData = data.frame(meta=colnames(rna.pseudobulk)),design = ~ 1)

dds.rna <- estimateSizeFactors(dds.rna)
sizeFactors(dds.rna)

rlog.rna <- rlog(dds.rna,blind = T)
saveRDS(dds.rna,"dds_rna.rds")
saveRDS(rlog.rna,"rlog_rna.rds")


###############################################
# Part 3: Pseudobulk enhancer activity
###############################################

# Pseudobulk ATAC enhancer matrix

atac <- readRDS("./final_archr_proj_archrGS.rds")

source("./Archr_Peak_Null_Permute.R")
source("./Archr_Peak_RawPval.R")
source("./getMatrixFromProject.R")
peaks.mat <- getMatrixFromProject.mod(atac,useMatrix = "PeakMatrix")
cell.names <- colnames(assay(peaks.mat))
peak.names <- paste0(seqnames(rowRanges(peaks.mat)),":", start(ranges(rowRanges(peaks.mat))),"-", end(ranges(rowRanges(peaks.mat))))


peaks.mat <- assay(peaks.mat)

dim(peaks.mat)
length(cell.names)
length(peak.names)

colnames(peaks.mat) <- cell.names
rownames(peaks.mat) <- peak.names
head(peaks.mat[,1:4])


# Construct pseudobulk peak/enhancer matrix
peaks.pseudobulk <- data.frame(rownames(peaks.mat))
for (i in labels){
  cells <- rownames(dplyr::filter(as.data.frame(atac@cellColData),predictedGroup_ArchR == i))
  
  peaks.sub <- peaks.mat[,colnames(peaks.mat) %in% cells]
  peaks.sub <- as(peaks.sub, "dgCMatrix")
  peaks.bulk <- Matrix::rowSums(peaks.sub)
  
  peaks.pseudobulk$i <- peaks.bulk
  colnames(peaks.pseudobulk)[dim(peaks.pseudobulk)[2]] <- i
  
  print("Iteration complete")
}
rownames(peaks.pseudobulk) <- peaks.pseudobulk$rownames.peaks.
peaks.pseudobulk <- peaks.pseudobulk[,-1]

# Remove features/peaks that have zero counts in ANY sample
dim(peaks.pseudobulk)
peaks.pseudobulk[peaks.pseudobulk == 0] <- NA
peaks.pseudobulk <- peaks.pseudobulk[complete.cases(peaks.pseudobulk),]
dim(peaks.pseudobulk)

#DESeq ATAC

library(DESeq2)
dds.atac <- DESeqDataSetFromMatrix(countData = peaks.pseudobulk,
                              colData = data.frame(meta=colnames(peaks.pseudobulk)),design = ~ 1)

dds.atac <- estimateSizeFactors(dds.atac)
sizeFactors(dds.atac)

rlog.atac <- rlog(dds.atac,blind = T)

saveRDS(rlog.atac,"rlog_atac.rds")
saveRDS(dds.atac,"dds_atac.rds")
rlog.rna <- readRDS("rlog_rna.rds")
rlog.atac <- readRDS("rlog_atac.rds")

###############################################
# Part 4: Subsetting and scaling from 0 to 1
###############################################

# Subset and reorder motif matrix, peak matrix and expression matrix
enhancer.mat <- assay(rlog.atac)
enhancer.mat <- as.data.frame(enhancer.mat)
enhancer.mat <- enhancer.mat[rownames(enhancer.mat) %in% colnames(enhancer.motifs),]

# Check order
motif.mat <- enhancer.motifs
motif.mat <- motif.mat[,colnames(motif.mat) %in% rownames(enhancer.mat)]
motif.mat <- motif.mat[ ,order(match(colnames(motif.mat), rownames(enhancer.mat))) ]
length(which(colnames(motif.mat) == rownames(enhancer.mat)))

expr.mat <- assay(rlog.rna)
expr.mat <- as.data.frame(expr.mat)

# Check order
tfs <- intersect(rownames(motif.mat),rownames(expr.mat))
expr.mat <- expr.mat[rownames(expr.mat) %in% tfs,]
motif.mat <- motif.mat[rownames(motif.mat) %in% tfs,]
expr.mat <- expr.mat[order(match(rownames(expr.mat), rownames(motif.mat))), ]

length(which(rownames(expr.mat)==rownames(motif.mat)))

# Rescale peak matrix and rna expr matrix
library(scales)
for (i in 1:ncol(expr.mat)){
  expr.mat[,i] <- rescale(expr.mat[,i],to = c(0,1))
}
library(scales)# Scale cell types
for (i in 1:ncol(enhancer.mat)){
  enhancer.mat[,i] <- rescale(enhancer.mat[,i],to = c(0,1))
}


dim(enhancer.mat)
dim(motif.mat)
dim(expr.mat)



length(which(rownames(motif.mat) == rownames(expr.mat)))

length(which(rownames(enhancer.mat) == colnames(motif.mat)))



length(which(rownames(motif.mat) == rownames(expr.mat)))

length(which(rownames(enhancer.mat) == colnames(motif.mat)))

dim(enhancer.mat)
dim(motif.mat)
dim(expr.mat)


###############################################
# Part 5: TFSEE matrix operations
###############################################

# Peform matrix multiplication of enhancer activity by motif prediction
int <- t(enhancer.mat) %*% t(motif.mat)

dim(int)
dim(expr.mat)

length(which(colnames(int)==rownames(expr.mat)))



# Multiply intermediate matrix by TF expression matrix 
tfsee <- int*t(expr.mat)
dim(tfsee)
head(tfsee)

saveRDS(tfsee,"tfsee_matrix.rds")
saveRDS(int,"int_matrix.rds")
saveRDS(expr.mat,"expr_matrix.rds")
saveRDS(motif.mat,"motif_matrix.rds")

tfsee <- readRDS("tfsee_matrix.rds")
# Write table of TFs for canSAR search:
write.csv(data.frame(colnames(tfsee)),"TFs_for_canSAR.csv")

##################################################################
# END OF TFSEE COMPUTATION
##################################################################

library(ArchR)
library(ComplexHeatmap)
library(DESeq2)
library(Seurat)
library(ggplot2)
library(dplyr)

# Read in tfsee matrix and write TFs out to csv for canSAR
tfsee <- readRDS("./tfsee_matrix.rds")
tfsee.t <- t(tfsee)
tfsee.t <- as.data.frame(tfsee.t)
tfsee.t$rowMeans <- rowMeans(as.matrix(tfsee.t))
tfsee.t <- dplyr::filter(tfsee.t,rowMeans > summary(tfsee.t$rowMeans)[4])
tfsee.t <- tfsee.t[,-12]# Remove rowMeans column
tfsee <- t(tfsee.t)
write.csv(data.frame(colnames(tfsee)),"TFs_for_canSAR.csv")
tfsee <- t(scale(t(tfsee)))

# Go to https://cansarblack.icr.ac.uk/cpat and input list of TFs

# Fold in druggability information for TFs and read in canSAR data
drug.scores <- read.csv("./cpat-results.csv")
drug.scores <- as.data.frame(cbind(drug.scores$gene_name,drug.scores$ligand_druggability_score))
colnames(drug.scores) <- c("TF","Score")
dim(drug.scores)

drug.scores <- drug.scores[duplicated(drug.scores) == FALSE,]
dim(drug.scores)
drug.scores <- dplyr::arrange(drug.scores,TF)
drug.scores <- drug.scores[-9,]# Remove duplicate DDIT3
drug.scores <- drug.scores[drug.scores$TF %in% colnames(tfsee),]
drug.scores.sub <- dplyr::filter(drug.scores,Score > 0)
length(which(drug.scores$TF== colnames(tfsee)))


# Rank frequency plot for subclones of endometrial metastasis
rownames(tfsee)
# Build TFSEE difference vector
tfsee.group1 <- tfsee[3,]
tfsee.group1 <- as.data.frame(tfsee.group1)
colnames(tfsee.group1) <- "group1"

tfsee.group2 <- tfsee[4,]
tfsee.group2 <- as.data.frame(tfsee.group2)
colnames(tfsee.group2) <- "group2"

print("Comparison: Subclone 1 v. Subclone 2")
print(rownames(tfsee)[3])
print(rownames(tfsee)[4])
print("Beginning comparison...")

length(which(rownames(tfsee.group1) == rownames(tfsee.group2)))
tfsee.group1 <- tfsee.group1
tfsee.group2 <- tfsee.group2
diff.tfsee <-tfsee.group1-tfsee.group2
hist(diff.tfsee)

colnames(diff.tfsee) <- "difference.tfsee"

total <- diff.tfsee

total$sum <- rowSums(total)

total.group1 <- dplyr::filter(total,sum > 0)
total.group2 <- dplyr::filter(total,sum < 0)

drug.group1 <- intersect(rownames(total.group1),drug.scores.sub$TF)
drug.group2 <- intersect(rownames(total.group2),drug.scores.sub$TF)

total.group1 <- total.group1[rownames(total.group1) %in% drug.group1,]
total.group2 <- total.group2[rownames(total.group2) %in% drug.group2,]

# Add druggability to DF:

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

drug.scores.new <- dplyr::filter(drug.scores,TF %in% rownames(total))

length(which(drug.scores.new$TF == rownames(total)))

total$drug.score <- drug.scores.new$Druggable

total <- dplyr::arrange(total,difference.tfsee)
total$rank <- 1:nrow(total)

p1<-ggplot(total,aes(x=rank,y = difference.tfsee))+
  geom_point(aes(color = drug.score),size = 2)+
  ylab("difference(Scaled TFSEE score)")+theme_bw()+
  xlab("Rank")+
  scale_color_manual(values = c("gray60","black"))+
  geom_hline(yintercept = 0.25,linetype="dashed")+
  geom_hline(yintercept = -0.25,linetype="dashed")+
  ggsave("TFSEE_Sublcone_Comparison.pdf",width =5,height = 4)



# Rank frequency plot for carcinosarcoma
rownames(tfsee)
# Build TFSEE fold change vector
tfsee.group1 <- tfsee[7,]
tfsee.group1 <- as.data.frame(tfsee.group1)
colnames(tfsee.group1) <- "group1"

tfsee.group2 <- tfsee[8,]
tfsee.group2 <- as.data.frame(tfsee.group2)
colnames(tfsee.group2) <- "group2"


print("Comparison: Sarcoma v. Carcinoma")
print(rownames(tfsee)[7])
print(rownames(tfsee)[8])
print("Beginning comparison...")

length(which(rownames(tfsee.group1) == rownames(tfsee.group2)))
tfsee.group1 <- tfsee.group1
tfsee.group2 <- tfsee.group2
diff.tfsee <-tfsee.group1-tfsee.group2
hist(diff.tfsee)

colnames(diff.tfsee) <- "difference.tfsee"

total <- diff.tfsee

total$sum <- rowSums(total)

total.group1 <- dplyr::filter(total,sum > 0)
total.group2 <- dplyr::filter(total,sum < 0)

drug.group1 <- intersect(rownames(total.group1),drug.scores.sub$TF)
drug.group2 <- intersect(rownames(total.group2),drug.scores.sub$TF)

total.group1 <- total.group1[rownames(total.group1) %in% drug.group1,]
total.group2 <- total.group2[rownames(total.group2) %in% drug.group2,]


# Add druggability to DF:

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

drug.scores.new <- dplyr::filter(drug.scores,TF %in% rownames(total))

length(which(drug.scores.new$TF == rownames(total)))

total$drug.score <- drug.scores.new$Druggable

total <- dplyr::arrange(total,difference.tfsee)
total$rank <- 1:nrow(total)

p2<-ggplot(total,aes(x=rank,y = difference.tfsee))+
  geom_point(aes(color = drug.score),size = 2)+
  ylab("difference(Scaled TFSEE score)")+theme_bw()+
  xlab("Rank")+
  scale_color_manual(values = c("gray60","black"))+
  geom_hline(yintercept = 0.25,linetype="dashed")+
  geom_hline(yintercept = -0.25,linetype="dashed")+
  ggsave("TFSEE_Carcinosarcoma.pdf",width =5,height = 4)




# Rank frequency plot for endometrial cancer of different histologies
rownames(tfsee)
# Build TFSEE fold change vector
tfsee.group1 <- colMeans(tfsee[3:4,])
tfsee.group1 <- as.data.frame(tfsee.group1)
colnames(tfsee.group1) <- "group1"

tfsee.group2 <- colMeans(tfsee[1:2,])
tfsee.group2 <- as.data.frame(tfsee.group2)
colnames(tfsee.group2) <- "group2"

print("Comparison: Serous v. Endometrioid EC")
print(paste0("Mean of:"))
print(rownames(tfsee[c(3,4),]))
print(paste0("Mean of:"))
print(rownames(tfsee[c(1,2),]))
print("Beginning comparison...")

length(which(rownames(tfsee.group1) == rownames(tfsee.group2)))
tfsee.group1 <- tfsee.group1
tfsee.group2 <- tfsee.group2
diff.tfsee <-tfsee.group1-tfsee.group2
hist(diff.tfsee)

colnames(diff.tfsee) <- "difference.tfsee"

total <- diff.tfsee

total$sum <- rowSums(total)

total.group1 <- dplyr::filter(total,sum > 0)
total.group2 <- dplyr::filter(total,sum < 0)

drug.group1 <- intersect(rownames(total.group1),drug.scores.sub$TF)
drug.group2 <- intersect(rownames(total.group2),drug.scores.sub$TF)

total.group1 <- total.group1[rownames(total.group1) %in% drug.group1,]
total.group2 <- total.group2[rownames(total.group2) %in% drug.group2,]


# Add druggability to DF:

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

drug.scores.new <- dplyr::filter(drug.scores,TF %in% rownames(total))

length(which(drug.scores.new$TF == rownames(total)))

total$drug.score <- drug.scores.new$Druggable

total <- dplyr::arrange(total,difference.tfsee)
total$rank <- 1:nrow(total)

p3<-ggplot(total,aes(x=rank,y = difference.tfsee))+
  geom_point(aes(color = drug.score),size = 2)+
  ylab("difference(Scaled TFSEE score)")+theme_bw()+
  xlab("Rank")+
  scale_color_manual(values = c("gray60","black"))+
  geom_hline(yintercept = 0.25,linetype="dashed")+
  geom_hline(yintercept = -0.25,linetype="dashed")+
  ggsave("TFSEE_SerousVEndometrioid_Endometrial.pdf",width =5,height = 4)





# Rank frequency plot for ovarian cancer of different histologies
rownames(tfsee)
# Build TFSEE fold change vector
tfsee.group1 <- colMeans(tfsee[c(6,11),])
tfsee.group1 <- as.data.frame(tfsee.group1)
colnames(tfsee.group1) <- "group1"

tfsee.group2 <- tfsee[5,]
tfsee.group2 <- as.data.frame(tfsee.group2)
colnames(tfsee.group2) <- "group2"

print("Comparison: Serous v. Endometrioid OC")
print(paste0("Mean of:"))
print(rownames(tfsee[c(6,11),]))
print(rownames(tfsee)[5])
print("Beginning comparison...")

length(which(rownames(tfsee.group1) == rownames(tfsee.group2)))
tfsee.group1 <- tfsee.group1
tfsee.group2 <- tfsee.group2
diff.tfsee <-tfsee.group1-tfsee.group2
hist(diff.tfsee)

colnames(diff.tfsee) <- "difference.tfsee"

total <- diff.tfsee

total$sum <- rowSums(total)

total.group1 <- dplyr::filter(total,sum > 0)
total.group2 <- dplyr::filter(total,sum < 0)

drug.group1 <- intersect(rownames(total.group1),drug.scores.sub$TF)
drug.group2 <- intersect(rownames(total.group2),drug.scores.sub$TF)

total.group1 <- total.group1[rownames(total.group1) %in% drug.group1,]
total.group2 <- total.group2[rownames(total.group2) %in% drug.group2,]


# Add druggability to DF:

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

drug.scores.new <- dplyr::filter(drug.scores,TF %in% rownames(total))

length(which(drug.scores.new$TF == rownames(total)))

total$drug.score <- drug.scores.new$Druggable

total <- dplyr::arrange(total,difference.tfsee)
total$rank <- 1:nrow(total)

p4<-ggplot(total,aes(x=rank,y = difference.tfsee))+
  geom_point(aes(color = drug.score),size = 2)+
  ylab("difference(Scaled TFSEE score)")+theme_bw()+
  xlab("Rank")+
  scale_color_manual(values = c("gray60","black"))+
  geom_hline(yintercept = 0.25,linetype="dashed")+
  geom_hline(yintercept = -0.25,linetype="dashed")+
  ggsave("TFSEE_SerousVEndometrioid_Ovarian.pdf",width =5,height = 4)



# Rank frequency plot for GIST versus HGSOC
rownames(tfsee)

# Build TFSEE fold change vector
tfsee.group1 <- colMeans(tfsee[c(9,10),])
tfsee.group1 <- as.data.frame(tfsee.group1)
colnames(tfsee.group1) <- "group1"

tfsee.group2 <- colMeans(tfsee[c(6,11),])
tfsee.group2 <- as.data.frame(tfsee.group2)
colnames(tfsee.group2) <- "group2"

print("Comparison: GIST v. HGSOC")
print(paste0("Mean of:"))
print(rownames(tfsee[c(9,10),]))
print(paste0("Mean of:"))
print(rownames(tfsee[c(6,11),]))
print("Beginning comparison...")

length(which(rownames(tfsee.group1) == rownames(tfsee.group2)))
tfsee.group1 <- tfsee.group1
tfsee.group2 <- tfsee.group2
diff.tfsee <-tfsee.group1-tfsee.group2
hist(diff.tfsee)

colnames(diff.tfsee) <- "difference.tfsee"

total <- diff.tfsee

total$sum <- rowSums(total)

total.group1 <- dplyr::filter(total,sum > 0)
total.group2 <- dplyr::filter(total,sum < 0)

drug.group1 <- intersect(rownames(total.group1),drug.scores.sub$TF)
drug.group2 <- intersect(rownames(total.group2),drug.scores.sub$TF)

total.group1 <- total.group1[rownames(total.group1) %in% drug.group1,]
total.group2 <- total.group2[rownames(total.group2) %in% drug.group2,]


# Add druggability to DF:

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

drug.scores.new <- dplyr::filter(drug.scores,TF %in% rownames(total))

length(which(drug.scores.new$TF == rownames(total)))

total$drug.score <- drug.scores.new$Druggable

total <- dplyr::arrange(total,difference.tfsee)
total$rank <- 1:nrow(total)

p5<-ggplot(total,aes(x=rank,y = difference.tfsee))+
  geom_point(aes(color = drug.score),size = 2)+
  ylab("difference(Scaled TFSEE score)")+theme_bw()+
  xlab("Rank")+
  scale_color_manual(values = c("gray60","black"))+
  geom_hline(yintercept = 0.25,linetype="dashed")+
  geom_hline(yintercept = -0.25,linetype="dashed")+
  ggsave("TFSEE_GISTvHGSOC.pdf",width =5,height = 4)


CombinePlots(list(p4,p3,p5),ncol=1)+ggsave("Freq_Distribtions_Suppl.pdf",
                                           width=5,height = 12)

CombinePlots(list(p1,p2),ncol=1)+ggsave("Freq_Distribtions_Main.pdf",
                                           width=5,height =8)


# Plot heatmap with druggability
# Make heatmap annotation
ha = HeatmapAnnotation(
  Site = factor(c(rep("1",5),rep("2",6))),
  Type = factor(c(rep("1",5),rep("2",5),rep("3",1))),
  Histology = factor(c(rep("1",5),rep("2",4),rep("3",1),rep("4",1))),
  Stage = factor(c(rep("1",1),rep("2",2),rep("3",1),rep("4",1),rep("5",4),rep("6",2))),
  Patient = factor(as.character(seq(1:11))))


library(ComplexHeatmap)

# Z-score the cell types
pdf(paste0("TFSEE_celltype_scaled.pdf"), width = 7, height =11)


drug.scores.sub <- dplyr::filter(drug.scores,Score > 0)

tf.names <- colnames(tfsee)[colnames(tfsee) %in% drug.scores.sub$TF]
idx <- match(tf.names,colnames(tfsee))

ha.2 = rowAnnotation(foo = anno_mark(at = idx, labels = tf.names))

drug.scores$Druggable <- ifelse(drug.scores$Score > 0, "Yes", "No")

activity=list(Druggability=c("Yes"="Black","No"="gray60")) 
ha.3 = HeatmapAnnotation(
  Druggability = drug.scores$Druggable,which = "row",col =activity)


Heatmap(scale(t(tfsee)),top_annotation = ha, 
        show_row_names =T,clustering_method_rows = "complete",clustering_method_columns = "complete",clustering_distance_rows = "euclidean",
        clustering_distance_columns = "euclidean",
        right_annotation = c(ha.2,ha.3))
dev.off()

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