9步骤完成单细胞数据挖掘文章全部图表复现

library(Seurat)
library(stringr)
library(clustree)
library(dplyr)
library(GSVA)
library(msigdbr)
library(GSEABase)
library(pheatmap)
library(harmony)
rm(list = ls())
gc()

# 1. 循环读取每个样本的数据
filename <- paste('processed_data/',list.files('processed_data/'),sep = '')
sceList <- lapply(filename, function(x){
  obj <- CreateSeuratObject(counts = Read10X(x),
                            project = str_split(x,'/')[[1]][2])
})
names(sceList) <- list.files('processed_data/')
sce <- merge(sceList[[1]],sceList[-1],add.cell.ids = names(sceList),project='OC')
dim(sce) # 38224个gene， 67323个cell

# 2. 查看线粒体(MT开头)、核糖体(RPS/RPL开头)基因占比
grep('^MT',x=rownames(sce@assays$RNA@data),value = T)
grep('^RP[SL]',x=rownames(sce@assays$RNA@data),value = T)

sce <- PercentageFeatureSet(sce,pattern = '^MT',col.name = 'percent.MT')
sce <- PercentageFeatureSet(sce,pattern = '^RP[SL]',col.name = 'percent.RP')
VlnPlot(sce,features = 'percent.MT',pt.size = 0)
VlnPlot(sce,features = 'percent.RP',pt.size = 0)
VlnPlot(sce,features = 'nCount_RNA',pt.size = 0)
VlnPlot(sce,features = 'nFeature_RNA',pt.size = 0)
# 3. 过滤细胞
sce <- subset(sce,subset = nCount_RNA>3000 & nFeature_RNA>300 & percent.MT<40)
dim(sce) # 38224个gene， 29389个cell
# 4. 过滤gene
sce <- sce[rowSums(sce@assays$RNA@counts>0)>3,]
dim(sce) # 28074个gene， 29389个cell

# 5. 看下批次效应
sce <- NormalizeData(sce) # 默认用logNormalize
sce <- FindVariableFeatures(sce) # 默认选Top2000作为高变gene
sce <- ScaleData(sce) # 默认用variableFeature进行scale
sce <- RunPCA(sce) # 默认用variableFeature进行PCA降维，降维前一定要做scale
DimPlot(sce,reduction = 'pca',group.by = 'orig.ident')

# 6. SCTransform进行标准化
sce <- SCTransform(sce) # 这一步包含了NormalizeData 、ScaleData、FindVariableFeatures三步
DefaultAssay(sce)
sce <- RunPCA(sce)
DimPlot(sce,reduction = 'pca',group.by = 'orig.ident')
ElbowPlot(sce,ndims = 40)

# 7. 在PCA基础上进行umap降维,然后分群
sce <- RunUMAP(sce,reduction = 'pca',dims = 1:40)
sce <- FindNeighbors(sce,dims = 1:40)

sce_res <- sce
for (i in c(0.01, 0.05, 0.1, 0.15, 0.2, 0.3,0.4, 0.5,0.8,1)){
  sce_res <- FindClusters(sce_res,resolution = i)
}
clustree(sce_res,prefix = 'SCT_snn_res.')

sce <- FindClusters(sce,resolution = 0.4)
DimPlot(sce,reduction = 'umap',group.by = 'seurat_clusters',label = T)

# 8. 手动进行细胞注释

# 免疫细胞  5,6,8,9,15,20 
DotPlot(sce,features = c('PTPRC','CD45'))
# T细胞   0,[8],9,13
DotPlot(sce,features = c('CD3D','CD3E','CD8A')) 
# B细胞   14,16,18,[20]
DotPlot(sce,features = c('CD79A', 'CD37', 'CD19', 'CD79B', 'MS4A1','CD20'))
# 浆细胞 [14,16,18]
DotPlot(sce,features = c('IGHG1','MZB1','SDC1','CD79A'))
# NK细胞  5,[9,15]
DotPlot(sce,features = c('FGFBP2','FCG3RA','CX3CR1'))
# NK细胞  8,15
DotPlot(sce,features = c('CD160','NKG7','GNLY','CD247','CCL3','GZMB','CXCR1','TYOB','PRF1'))
# 内皮细胞 [11]
DotPlot(sce,features = c('PECAM1','VWF'))
# 成纤维细胞 [0,7],10,11
DotPlot(sce,features = c('FGF7','MME','DCN','LUM','GSN','PF4','PPBP'))
# 上皮细胞 1,2,[3],4,[12,13,17,19]
DotPlot(sce,features = c('EPCAM','KRT19','PROM1','ALDH1A1','CD24'))
# 单核细胞和巨噬细胞 5,6,9,15
DotPlot(sce,features = c('CD68','CD163','CD14','LYZ'))
# Ovarian somatic cell  [13]
DotPlot(sce,features = c('LGR5'))

marker <- data.frame(cluster = 0:20,cell = 0:20)
marker[marker$cluster %in% c(5,6,8,9,15,20),2] <- 'Immune cells'
marker[marker$cluster %in% c(0,8,9,13),2] <- 'T cells'
marker[marker$cluster %in% c(14,16,18,20),2] <- 'B cells'
marker[marker$cluster %in% c(14,16,18),2] <- 'plasma cells'
marker[marker$cluster %in% c(5,9,15),2] <- 'NK cells'
marker[marker$cluster %in% c(11),2] <- 'Endothelial cells'
marker[marker$cluster %in% c(5,6),2] <- 'myeloid cells'
marker[marker$cluster %in% c(1,2,3,4,12,13,17,19),2] <- 'epithelial cell'
marker[marker$cluster %in% c(0,7,10),2] <- 'fibroblast'
marker[marker$cluster %in% c(13),2] <- 'Ovarian somatic cell'

sce@meta.data$cell_type <- sapply(sce@meta.data$seurat_clusters,function(x){marker[x,2]})
DimPlot(sce,reduction = 'umap',group.by = 'cell_type',label = T)

# 9. marker gene 热图展示
heatmap_marker_gene <- c('CD3D','CD3E','CD8A',
                         'CD79A', 'CD19', 'CD79B', 'MS4A1','CD20',
                         'MZB1','SDC1','CD79A',
                         'CD68','CD163','CD14','LYZ')
heatmap_cells <- rownames(sce@meta.data[sce@meta.data$cell_type %in% c("myeloid cells","B cells","T cells"),])
sub_sce <- subset(sce,cells = heatmap_cells)
DoHeatmap(sub_sce,features = heatmap_marker_gene,group.by = 'cell_type',size = 3)

# 10. GSVA分析

Hallmark_gene_set <- msigdbr(species = 'Homo sapiens',category = 'H')
Hallmark_list <- split(Hallmark_gene_set$gene_symbol,Hallmark_gene_set$gs_name)

gsva_res <- gsva(expr = sce@assays$SCT@scale.data,
                 gset.idx.list = Hallmark_list,
                 method = 'ssgsea',
                 kcdf='Poisson',
                 parallel.sz = parallel::detectCores())
pheatmap(gsva_res,fontsize_row = 4,
         height = 10,
         width=18,
         annotation_col = data.frame(row.names = colnames(sce),group=sce@meta.data$orig.ident),
         show_colnames = F,
         show_rownames = T)

cell_type_gene_set <- msigdbr(species = 'Homo sapiens',category = 'C8')
cell_type_list <- split(cell_type_gene_set$gene_symbol,cell_type_gene_set$gs_name)
names(cell_type_list)

t_cell_marker <- read.table('CellMarker_T_cell.csv',sep = ',',header = 1)
b_cell_marker <- read.table('CellMarker_B_cell.csv',sep = ',',header = 1)
macrophage_marker <- read.table('CellMarker_Macrophage.csv',sep = ',',header = 1)

t_cell_marker <- t_cell_marker[t_cell_marker$Species=='Human',]
b_cell_marker <- b_cell_marker[b_cell_marker$Species=='Human',]
macrophage_marker <- macrophage_marker[macrophage_marker$Species=='Human',]

# T cell 挑了4种进行热图展示
cytotoxic <- unique(t_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(t_cell_marker$Cell.Type)),'cytotoxic')]
naive <- unique(t_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(t_cell_marker$Cell.Type)),'naive')]
regulatory <- unique(t_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(t_cell_marker$Cell.Type)),'regulatory')]
exhausted <- unique(t_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(t_cell_marker$Cell.Type)),'exhausted')]
t_cell_df <- t_cell_marker[t_cell_marker$Cell.Type %in% c(cytotoxic,naive,regulatory,exhausted),]
t_cell_df[t_cell_df$Cell.Type %in% naive,"Cell.Type"] <- 'naive'
t_cell_df[t_cell_df$Cell.Type %in% cytotoxic,"Cell.Type"] <- 'cytotoxic'
t_cell_df[t_cell_df$Cell.Type %in% regulatory,"Cell.Type"] <- 'regulatory'
t_cell_df[t_cell_df$Cell.Type %in% exhausted,"Cell.Type"] <- 'exhausted'
t_cell_list <- split(t_cell_df$Cell.Marker,t_cell_df$Cell.Type)

# B cell 挑了5种进行热图展示
unique(b_cell_marker$Cell.Type)
b_plasma <- unique(b_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(b_cell_marker$Cell.Type)),'plasma')]
b_memory <- unique(b_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(b_cell_marker$Cell.Type)),'memory')]
b_naive <- unique(b_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(b_cell_marker$Cell.Type)),'naive')]
b_activated <- unique(b_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(b_cell_marker$Cell.Type)),'activated')]
b_transitional <- unique(b_cell_marker$Cell.Type)[str_detect(str_to_lower(unique(b_cell_marker$Cell.Type)),'transitional')]

b_cell_df <- b_cell_marker[b_cell_marker$Cell.Type %in% c(b_plasma,b_memory,b_naive,b_activated,b_transitional),]
b_cell_df[b_cell_df$Cell.Type %in% b_plasma,"Cell.Type"] <- 'plasma'
b_cell_df[b_cell_df$Cell.Type %in% b_memory,"Cell.Type"] <- 'memory'
b_cell_df[b_cell_df$Cell.Type %in% b_naive,"Cell.Type"] <- 'naive'
b_cell_df[b_cell_df$Cell.Type %in% b_activated,"Cell.Type"] <- 'activated'
b_cell_df[b_cell_df$Cell.Type %in% b_transitional,"Cell.Type"] <- 'transitional'
b_cell_list <- split(b_cell_df$Cell.Marker,b_cell_df$Cell.Type)

# Macrophage cell
unique(macrophage_marker$Cell.Type)
M1 <- unique(macrophage_marker$Cell.Type)[str_detect(str_to_lower(unique(macrophage_marker$Cell.Type)),'m1')]
M2 <- unique(macrophage_marker$Cell.Type)[str_detect(str_to_lower(unique(macrophage_marker$Cell.Type)),'m2')]
macrophage_df <- macrophage_marker[macrophage_marker$Cell.Type %in% c(M1,M2),]
macrophage_list <- split(macrophage_df$Cell.Marker,macrophage_df$Cell.Type)

# 到这一步，gene list里的每一项，是一个字符串（gene symbol组成的）组成的vector，我写个函数，把每个字符串里的gene symbol拆出来，
# 最后组合成一个由gene symbol组成的vector，然后去个重
myfun <- function(x){
  res <- c()
  for (i in x){
    temp <- unlist(str_split(gsub(' ','',i),','))
    res <- c(res,temp)
  }
  unique(res)
}

t_cell_list <- lapply(t_cell_list, function(x){myfun(x)})
b_cell_list <- lapply(b_cell_list, function(x){myfun(x)})
macrophage_list <- lapply(macrophage_list, function(x){myfun(x)})

# 找到三种细胞分别对应的表达矩阵
t_cell <- colnames(sce[,sce@meta.data$cell_type=='T cells'])
b_cell <- colnames(sce[,sce@meta.data$cell_type=='B cells'])
macrophage_cell <- colnames(sce[,sce@meta.data$cell_type=='myeloid cells'])

t_cell_exp <- sce@assays$SCT@scale.data[,t_cell]
b_cell_exp <- sce@assays$SCT@scale.data[,b_cell]
macrophage_cell_exp <- sce@assays$SCT@scale.data[,macrophage_cell]

t_cell_gsva_res <- gsva(expr = t_cell_exp,
                           gset.idx.list = t_cell_list,
                           method = 'ssgsea',
                           kcdf='Poisson',
                           parallel.sz = parallel::detectCores())
pheatmap(t_cell_gsva_res,fontsize_row = 8,
         height = 10,
         width=18,
         annotation_col = data.frame(row.names = colnames(sce),group=sce@meta.data$orig.ident),
         show_colnames = F,show_rownames = T,cluster_cols = F)

b_cell_gsva_res <- gsva(expr = b_cell_exp,
                        gset.idx.list = b_cell_list,
                        method = 'ssgsea',
                        kcdf='Poisson',
                        parallel.sz = parallel::detectCores())

pheatmap(b_cell_gsva_res,fontsize_row = 8,
         height = 10,
         width=18,
         annotation_col = data.frame(row.names = colnames(sce),group=sce@meta.data$orig.ident),
         show_colnames = F,show_rownames = T,cluster_cols = F)

macrophage_cell_gsva_res <- gsva(expr = macrophage_cell_exp,
                        gset.idx.list = macrophage_list,
                        method = 'ssgsea',
                        kcdf='Poisson',
                        parallel.sz = parallel::detectCores())

pheatmap(macrophage_cell_gsva_res,fontsize_row = 8,
         height = 10,
         width=18,
         annotation_col = data.frame(row.names = colnames(sce),group=sce@meta.data$orig.ident),
         show_colnames = F,show_rownames = T,cluster_cols = F)

#加载需要的R包
library(TCGAbiolinks)
library(WGCNA)
library(stringr)
library(SummarizedExperiment)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tibble)
library(CIBERSORT)
library(tidyr)
library(ggplot2)
library(RColorBrewer)
library(readxl)
library(survival)
library(survminer)
library(WGCNA)
library(data.table)
library(dplyr)
library(IMvigor210CoreBiologies)
library(e1071)
library(pROC)
library(NMF)
library(pheatmap)
library(clusterProfiler)
library(ggsignif)
library(limma)
library(org.Hs.eg.db)
library(glmnet)
library(tibble)
library(ggrisk)

rm(list = ls())
gc()exp <- fread('TCGA-OV.htseq_fpkm.tsv',header = T)
exp <- column_to_rownames(exp,var = 'Ensembl_ID')

# 行名转换
exp <- as.data.frame(exp)
exp$ensembl_id <- str_split(rownames(exp),'\\.',simplify = T)[,1]
id_df <- bitr(exp$ensembl_id,fromType = 'ENSEMBL',toType = 'SYMBOL',OrgDb = org.Hs.eg.db)
id_df <- id_df[!duplicated(id_df$ENSEMBL),]
index <- match(exp$ensembl_id,id_df$ENSEMBL)
# 对symbol id这一列去重去空，设为行名
exp$symbol_id <- id_df[index,'SYMBOL']
exp <- exp[!is.na(exp$symbol_id),]
exp <- exp[!duplicated(exp$symbol_id),]
rownames(exp) <- exp$symbol_id
exp <- exp[,1:(ncol(exp)-2)]

save(exp,file = 'TCGA_OV.Rdata')

# 免疫浸润
load('TCGA_OV.Rdata')
lm22f <- system.file('extdata','LM22.txt',package = 'CIBERSORT')
TME.result <- cibersort(lm22f,as.matrix(exp),perm = 100,QN=T)
TME.result <- as.data.frame(TME.result)
save(result,file = 'TME_result.Rdata')

group <- ifelse(as.numeric(substr(colnames(exp),14,15))<10,'tumor','normal')
patient_exp <- exp[,which(group=='tumor')]
result <- TME.result[,-c(23:25)] %>% as.data.frame() %>% 
  rownames_to_column('Sample') %>% gather(key = cell_type,value = proportion,-Sample)

my_pallet <- colorRampPalette(brewer.pal(8,'Set1'))
ggplot(result,aes(Sample,proportion,fill=cell_type))+geom_bar(stat = 'identity')+
  scale_fill_manual(values = my_pallet(22))+
  labs(x='',y = "Estiamted Proportion",fill='cell_type')+theme_bw()+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "bottom")

# KM-plot
load('TCGA_OV.Rdata')
load('TME_result.Rdata')
# 准备生存数据，包含OS,OS.time,DFI,DFI.time这些信息
# https://gdc.cancer.gov/about-data/publications/pancanatlas
sur_data <- read_xlsx('TCGA-CDR-SupplementalTableS1.xlsx',sheet = 1)
sur_data <- column_to_rownames(sur_data,var = '...1')
sur_data <- sur_data[sur_data$type=='OV',c('bcr_patient_barcode','OS','OS.time','DFI','DFI.time')]
rownames(sur_data) <- sur_data$bcr_patient_barcode
save(sur_data,file = 'survival_data.Rdata')
# 建立sample和patient对应关系，把cibersort结果中Macrophages M1在每个样本中的信息转移过来
sample_patient_df <- data.frame(sample=colnames(exp),patient=substr(colnames(exp),1,12))
sample_patient_df$TAM_M1 <- TME.result$`Macrophages M1`
sample_patient_df$TAM_M1_type <- ifelse(sample_patient_df$TAM_M1>median(sample_patient_df$TAM_M1),'high','low')
sample_patient_df <- sample_patient_df[!duplicated(sample_patient_df$patient),]
# 下载的生存数据的临床信息中patient id 与 exp中patient id 取交集
sur_data <- sur_data[intersect(sample_patient_df$patient,sur_data$bcr_patient_barcode),]
# 把sample_patient_df中Macrophages M1的信息对应到sur_data中每个patient id上
index <- match(sur_data$bcr_patient_barcode,sample_patient_df$patient)
sur_data$TAM_M1 <- sample_patient_df[index,'TAM_M1_type']

M1_OS <- survfit(Surv(OS.time,OS)~TAM_M1,data = sur_data)
ggsurvplot(M1_OS,pval = T,risk.table = T,surv.median.line = 'hv',
           title='Overall survival',xlab='Days')
M1_DFI <- survfit(Surv(DFI.time,DFI)~TAM_M1,data = sur_data)
ggsurvplot(M1_DFI,pval = T,risk.table = T,surv.median.line = 'hv',
           title='Disease Free Interval',xlab='Days')
save(sur_data,file = 'sur_data.Rdata')

# WGCNA
load('TCGA_OV.Rdata')
load('TME_result.Rdata')
load('sur_data.Rdata')
exp <- as.data.frame(exp)

# 过滤异常值
exp <- exp[rowMads(as.matrix(exp))>0.01,]
dim(exp)

# 查找并删除异常值
temp <- exp
colnames(temp) <- 1:ncol(exp)
plot(hclust(dist(t(temp)))) # 无异常样本

# 构建基因共表达网络
exp <- as.data.frame(t(exp))

# 选择合适的power
sft <- pickSoftThreshold(exp,networkType = 'signed')

a <- ggplot(sft$fitIndices,aes(x=Power,y=-sign(slope)*SFT.R.sq))+
  geom_text(label=sft$fitIndices$Power,color='red')+
  geom_hline(yintercept = 0.9,color='red')+
  xlab("Soft Threshold (power)")+ylab("Scale Free Topology Model Fit,signed R^2")+
  ggtitle("Scale independence")+theme(plot.title = element_text(hjust = 0.5))

b <- ggplot(sft$fitIndices,aes(x=Power,y=mean.k.))+
  geom_text(label=sft$fitIndices$Power,color='red')+
  xlab("Soft Threshold (power)")+ylab("Mean Connectivity")+
  ggtitle('Mean connectivity')+theme(plot.title = element_text(hjust = 0.5))

a+b
power <- sft$powerEstimate
power <- 5

# 检验是否符合无尺度网络
k <- softConnectivity(exp,power = power)
scaleFreePlot(k)

# 构建邻接矩阵
adj <- adjacency(exp,power = power,type = 'signed')

# 计算TOM矩阵
TOM <- TOMsimilarity(adj,TOMType = 'signed')
save(TOM,file = 'TOM.Rdata')
load('TOM.Rdata')
# 计算gene间相异性
dissTOM <- 1-TOM

# 对gene进行聚类和可视化
geneTree <- hclust(as.dist(dissTOM),method = 'average')
sizeGrWindow(12,9)
plot(geneTree,labels = F,hang=0.04,xlab='',sub='')

# 将聚类树进行修剪，把gene归到不同模块里
dynamic_module <- cutreeDynamic(dendro = geneTree,distM = dissTOM,minClusterSize = 50,pamRespectsDendro = F)
table(dynamic_module)
module_color <- labels2colors(dynamic_module)
plotDendroAndColors(dendro = geneTree,colors = module_color,dendroLabels = F,hang=0.03,addGuide = T,
                    guideHang = 0.05,groupLabels='Dynamic color tree',
                    main='gene dendrogram and module colors')

# 随机选择400个基因画拓扑重叠热图
set.seed(10)
select = sample(5000, size = 400)
selectTOM = dissTOM[select, select]
selectTree = hclust(as.dist(selectTOM), method = "average") 
selectColors = module_color[select]

sizeGrWindow(9,9) 
plotDiss = selectTOM^power
diag(plotDiss) = NA
TOMplot(plotDiss, 
        selectTree, 
        selectColors, 
        main = "Network heatmap plot, selected genes") 

# 计算每个模块特征gene向量
MEList <- moduleEigengenes(exp,dynamic_module)
MEs <- MEList$eigengenes

# 计算特征基因的相异度，基于相异度对原有模块进行层次聚类
dissME <- 1-cor(MEs)
MEtree <- hclust(as.dist(dissME),method = 'average')

# 对模块进行合并
cut_height <- 0.6
sizeGrWindow(12,8)
plot(MEtree)
abline(h=cut_height,col='red')

merge_module <- mergeCloseModules(exprData = exp,colors=module_color,cutHeight = cut_height)
merged_color <- merge_module$colors
# merged_MEs每一行是一个病人ID，每一列是一个模块的特征gene向量
merged_MEs <- merge_module$newMEs

sizeGrWindow(12,8)
plotDendroAndColors(dendro = geneTree,colors = cbind(module_color,merged_color),dendroLabels = F,hang=0.03,addGuide = T,
                    guideHang = 0.05,groupLabels='Dynamic color tree',
                    main='gene dendrogram and module colors')
table(merged_color)

# 计算模块与TAM_M1，TAM_M2间的相关性
trait <- data.frame(row.names = rownames(exp),ID = substr(rownames(exp),1,12),
                    M1 = TME.result$`Macrophages M1`,M2 = TME.result$`Macrophages M2`)
index <- match(trait$ID,sur_data$bcr_patient_barcode)
trait$os <- sur_data[index,'OS']


module_trait_cor <- cor(merged_MEs,trait[,2:3],use = 'p')
module_trait_pvalue <- corPvalueStudent(module_trait_cor,ncol(exp))
trait_colors <- numbers2colors(trait[,2:3],signed = T,centered = T)
a <- paste(signif(module_trait_cor, 2), "\n(", signif(module_trait_pvalue, 1), ")", sep = "")
sizeGrWindow(8,6)
labeledHeatmap(module_trait_cor,xLabels = colnames(trait[,2:3]),yLabels = colnames(merged_MEs),colorLabels = FALSE, 
               colors = blueWhiteRed(50),textMatrix = a,
               setStdMargins = FALSE,cex.text = 0.65,zlim = c(-1,1), 
               main = paste("Module-trait relationships"))

# 找到M1相关的关键模块
module <- 'MEblue'
module_membership <- signedKME(exp,merged_MEs)
# 计算表达矩阵exp和OS之间相关性
gene_os_sig <- apply(exp,2,function(x){cor(x,trait$os)})

verboseScatterplot(abs(module_membership$kMEblue),
                   abs(gene_os_sig),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for OS",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2)

# 找出MEblue模块内所有基因
module_membership <- module_membership[!is.na(module_membership$kMEblue),]
# 筛选出MM>0.5 和 exp-OS相关性>0.1 的作为hub gene
MM_gene <- rownames(module_membership[abs(module_membership$kMEblue)>0.5,])
OS_sig_gene <- names(gene_os_sig[abs(gene_os_sig)>0.1])
hub_gene <- intersect(MM_gene,OS_sig_gene)

save(hub_gene,file = 'WGCNA_hub_gene.Rdata')

# SVM预测免疫类型
# DEseq安装 http://bioconductor.org/packages/3.8/bioc/html/DESeq.html
# IMvigor210CoreBiologies安装：http://research-pub.gene.com/IMvigor210CoreBiologies/packageVersions/
data(cds)
expreSet <- as.data.frame(counts(cds))
annoData <- fData(cds)
phenoData <- pData(cds)

# 清洗数据
expreSet <- expreSet[rowSums(expreSet>0)>3,]

phenoData$`Immune phenotype` <- as.character(phenoData$`Immune phenotype`)
phenoData <- phenoData[!is.na(phenoData$`Immune phenotype`),]
phenoData$`Immune phenotype` <- as.factor(phenoData$`Immune phenotype`)

expreSet <- expreSet[,rownames(phenoData)]
# count转TPM
index <- match(rownames(expreSet),annoData$entrez_id)
expreSet$gene_length <- annoData[index,'length']
kb <- expreSet$gene_length/1000
rpk <- expreSet[,-ncol(expreSet)] / kb
tpm <- t(t(rpk) / colSums(rpk) * 1E6)
exp <- tpm

# 取log
exp <- log2(exp+1)
exp <- as.data.frame(exp)

# 划分训练集和测试集
test_size <- 0.2
test_sample <- sample(colnames(exp),size = as.integer(ncol(exp)*test_size))
train_sample <- setdiff(colnames(exp),test_sample)
# 训练集和测试集的表达矩阵
test_exp <- as.data.frame(t(exp[,test_sample]))
train_exp <- as.data.frame(t(exp[,train_sample]))
# 训练集和测试集的标签  
test_label <- phenoData[test_sample,"Immune phenotype"]
train_label <- phenoData[train_sample,"Immune phenotype"]
# 训练集标签转换为二分类
train_label_desert <- as.factor(ifelse(train_label=='desert','desert','not desert'))
train_label_inflamed <- as.factor(ifelse(train_label=='inflamed','inflamed','not inflamed'))
train_label_excluded <- as.factor(ifelse(train_label=='excluded','excluded','not excluded'))
# 测试集标签转换为二分类
test_label_desert <- as.factor(ifelse(test_label=='desert','desert','not desert'))
test_label_inflamed <- as.factor(ifelse(test_label=='inflamed','inflamed','not inflamed'))
test_label_excluded <- as.factor(ifelse(test_label=='excluded','excluded','not excluded'))
# 测试集标签转换为数字
test_label_desert <- ifelse(test_label=='desert',1,0)
test_label_inflamed <- ifelse(test_label=='inflamed',1,0)
test_label_excluded <- ifelse(test_label=='excluded',1,0)

# 训练模型
desert_model <- svm(x = train_exp, 
                    y = train_label_desert,
                    type = 'C',kernel = 'radial',probability = T)
inflamed_model <- svm(x = train_exp, 
                    y = train_label_inflamed,
                    type = 'C',kernel = 'radial',probability = T)
excluded_model <- svm(x = train_exp, 
                    y = train_label_excluded,
                    type = 'C',kernel = 'radial',probability = T)
# 进行预测
test_predict_desert <- predict(desert_model,test_exp,prob = T)
# 得出预测结果的概率
desert_prob <- attr(test_predict_desert,'probabilities')[,2]
test_predict_inflamed <- predict(inflamed_model,test_exp)
inflamed_prob <- attr(test_predict_desert,'probabilities')[,2]
test_predict_excluded <- predict(excluded_model,test_exp)
excluded_prob <- attr(test_predict_desert,'probabilities')[,2]

desert_roc_obj <- roc(response = as.numeric(test_label_desert),
                      predictor = as.numeric(desert_prob),na.rm=T)
desert_auc <- round(auc(test_label_desert,desert_prob),2)
inflamed_roc_obj <- roc(response = as.numeric(test_label_inflamed),
                        predictor = as.numeric(inflamed_prob),na.rm=T)
inflamed_auc <- round(auc(test_label_inflamed,inflamed_prob),2)
excluded_roc_obj <- roc(response = as.numeric(test_label_excluded),
                        predictor = as.numeric(excluded_prob),na.rm=T)
excluded_auc <- round(auc(test_label_excluded,excluded_prob),2)

desert <- ggroc(desert_roc_obj)
inflamed <- ggroc(inflamed_roc_obj)
excluded <- ggroc(excluded_roc_obj)

roc_data <- list(desert_roc_obj,inflamed_roc_obj,excluded_roc_obj)
names(roc_data) <- c('desert','inflamed','excluded')

## sensitivity(敏感性): 也称recall或TPR(真阳性率，实际为正例并且模型预测为正例 占 所有所有正例 的比例)，越接近1越好
## specificity(特异性)：也称TNR(真阴性率，实际为负例并且模型预测为负例 占 所有所有负例 的比例)，越接近1越好
ggroc(roc_data)+geom_segment(aes(x = 1, xend = 0, y = 0, yend = 1), color="red",linetype=6)+
  annotate('text',x=c(0.12,0.12,0.12),y=c(0.1,0.06,0.02),label=c(paste0('desert_auc: ',desert_auc),
                                                              paste0('inflamed_auc: ',inflamed_auc),
                                                              paste0('excluded_auc: ',excluded_auc))

load('TCGA_OV.Rdata')
load('WGCNA_hub_gene.Rdata')

hubgene_exp <- exp[hub_gene,]

estimate <- nmf(hubgene_exp,rank = 2:10,method = 'brunet')
save(estimate,file = 'estimate.Rdata')
plot(estimate)

rank <- 6
nmf.rank <- nmf(hubgene_exp,rank = rank,method = 'brunet')
group <- predict(nmf.rank)
index <- extractFeatures(nmf.rank,'max')
index <- unlist(index)

nmf.exp <- hubgene_exp[index,]
nmf.exp <- na.omit(nmf.exp)

jco <- c("#2874C5","#EABF00","#C6524A","#868686")
consensusmap(nmf.rank,labRow=NA,labCol=NA,
             annCol = data.frame("cluster"=group[colnames(nmf.exp)]))

anno_col <- data.frame(sample = colnames(nmf.exp),group = group)
anno_col <- anno_col[order(anno_col$group),]
pheatmap(nmf.exp[,rownames(anno_col)],show_colnames = F,annotation_col = select(anno_col,group),
         cluster_cols = F)

# 挑出两个组分析
group12 <- group[group==2 | group==1]
group12 <- sort(group12)
exp_temp <- nmf.exp[,names(group12)]
pheatmap(exp_temp,show_colnames = F,annotation_col = data.frame(row.names = colnames(exp_temp),group = group12),
         cluster_cols = F)

# KM-plot
load('survival_data.Rdata')
load('TME_result.Rdata')
temp <- data.frame(row.names = colnames(exp_temp),bcr_patient_barcode = substr(colnames(exp_temp),1,12),
                   group=group12)
temp <- left_join(temp,sur_data,by='bcr_patient_barcode')
group_OS <- survfit(Surv(OS.time,OS)~group,data = temp)
ggsurvplot(group_OS,pval = T,risk.table = T,surv.median.line = 'hv',
           title='Overall survival',xlab='Days')

# M1细胞的小提琴图
vln_df <- data.frame(row.names = names(group12),M1=TME.result[names(group12),'Macrophages M1'],group = group12)
ggplot(vln_df,aes(x=group,y=M1))+geom_violin(aes(fill=group))+geom_boxplot(width=0.2)+
    geom_signif(comparisons = list(c('1','2')),map_signif_level = F)

# KEGG富集
kegg_exp <- exp[,names(group12)]
design <- model.matrix(~factor(as.character(group12),levels = c('1','2')))
fit <- lmFit(kegg_exp,design=design)
fit <- eBayes(fit)
deg <- topTable(fit,coef = 2,number = Inf)
deg$change <- ifelse((deg$logFC>log2(2))&(deg$adj.P.Val<0.05),'up',
                     ifelse((deg$logFC<log2(0.5))&(deg$adj.P.Val<0.05),'down','nochange'))
save(deg,file = 'TCGA-deg.Rdata')
up_1_id <- bitr(rownames(deg[deg$change=='up',]),fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = org.Hs.eg.db)
kegg <- enrichKEGG(up_1_id$ENTREZID,organism = 'hsa',keyType = 'ncbi-geneid')
dotplot(kegg)

up_2_id <- bitr(rownames(deg[deg$change=='down',]),fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = org.Hs.eg.db)
kegg <- enrichKEGG(up_2_id$ENTREZID,organism = 'hsa',keyType = 'ncbi-geneid')
dotplot(kegg)

load('TCGA_OV.Rdata')
load('survival_data.Rdata')
load('TCGA-deg.Rdata')
# 数据处理
deg_gene <- rownames(deg[deg$change!='nochange',])
exp <- as.data.frame(t(exp))
exp <- exp[,deg_gene]
exp$ID <- substr(rownames(exp),1,12)
exp <- exp[!duplicated(exp$ID),]
index <- match(rownames(sur_data),exp$ID)
index <- index[!is.na(index)]

sur_data <- merge(sur_data,exp,by.x='bcr_patient_barcode',by.y='ID')

# 批量单因素COX回归
gene <- c()
p_value <- c()
HR <- c()
lower95 <- c()
upper95 <- c()
for (i in 6:ncol(sur_data)){
  res <- coxph(Surv(OS.time,OS)~sur_data[,i],data=sur_data)
  sum_res <- summary(res)
  p <- sum_res$coefficients[,'Pr(>|z|)']
  if (p<0.05){
    gene <- c(gene,colnames(sur_data)[i])
    p_value <- c(p_value,p)
    HR <- c(HR,sum_res$conf.int[,'exp(coef)'])
    lower95 <- c(lower95,sum_res$conf.int[,'lower .95'])
    upper95 <- c(upper95,sum_res$conf.int[,'upper .95'])
  }
}
cox_df <- data.frame(row.names = gene,pvalue=p_value,HR=HR,lower.95=lower95,upper.95=upper95)

# LASSO进一步筛选gene
sur_data <- column_to_rownames(sur_data,'bcr_patient_barcode')
sur_data <- select(sur_data,c(OS,OS.time,rownames(cox_df)))
sur_data <- sur_data[!is.na(sur_data$OS.time),]
x <- as.matrix(sur_data[,5:ncol(sur_data)])
y <- Surv(sur_data$OS.time,sur_data$OS)
glm.fit <- glmnet(x,y,family = 'cox',alpha = 1)
plot(glm.fit)
cv.fit <- cv.glmnet(x,y,alpha = 1,nfolds = 10,family="cox")
plot(cv.fit)
lambda <- cv.fit$lambda.min

glm.fit <- glmnet(x,y,family = 'cox',alpha = 1,lambda = lambda)
lasso_filter_gene <- names(glm.fit$beta[glm.fit$beta[,1]!=0,1])
coef <- glm.fit$beta[glm.fit$beta[,1]!=0,1]

# 筛选后的gene用多因素cox回归建模
sur_data_temp <- select(sur_data,c(OS,OS.time,lasso_filter_gene))
# gene名IGKV4-1，在ggrisk时会报错，改成.
colnames(sur_data_temp) <- c("OS","OS.time","LAMP3","STAT1","BATF2","CXCL10","CXCL11","IGKV4.1")
cox.res <- coxph(Surv(OS.time,OS)~.,data = sur_data_temp)
riskScore <- predict(cox.res,type = 'risk',newdata=sur_data_temp)
riskScore <- as.data.frame(riskScore)
riskScore$ID <- rownames(riskScore)
riskScore$risk <- ifelse(riskScore$riskScore>median(riskScore$riskScore),'high','low')
riskScore$OS <- sur_data$OS
riskScore$OS.time <- sur_data$OS.time

ggrisk(cox.res)

# riskscore做生存分析
surv_fit <- survfit(Surv(OS.time,OS)~risk,data = riskScore)
ggsurvplot(surv_fit,data = riskScore,pval = T,risk.table = T,surv.median.line = 'hv',
           legend.title = 'RiskScore',title = 'Overall survival',
           ylab='Cummulative survival',xlab='Time(Days)')