单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析1：https://cloud.tencent.com/developer/article/2055573
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析2：https://cloud.tencent.com/developer/article/2072069
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析3：https://cloud.tencent.com/developer/article/2078159
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析4：https://cloud.tencent.com/developer/article/2078348
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析5：https://cloud.tencent.com/developer/article/2084580
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析6：https://cloud.tencent.com/developer/article/2085385
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析7：https://cloud.tencent.com/developer/article/2085705
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析8：https://cloud.tencent.com/developer/article/2086805
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析9：https://cloud.tencent.com/developer/article/2087563
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析10：https://cloud.tencent.com/developer/article/2090290
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析11：https://cloud.tencent.com/developer/article/2093123
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析12：https://cloud.tencent.com/developer/article/2093208
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析13：https://cloud.tencent.com/developer/article/2093567
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析14: https://cloud.tencent.com/developer/article/2094034
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析15：https://cloud.tencent.com/developer/article/2102671
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析16:https://cloud.tencent.com/developer/article/2103030
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析17：https://cloud.tencent.com/developer/user/8402865
单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析18：https://cloud.tencent.com/developer/article/2120405
image.pngPart III: Assemble inputs and run TFSEE:
image.pnglibrary(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
##################################################################Part IV: Plotting Figure 6
image.pnglibrary(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")
