单细胞 芯片 转录组测序的数据挖掘文章一比一复现
自己造轮子
(生信技能树优秀学员“细心的网友”)
马拉松授课线上直播课程，最近的一期是5月8号（今晚七点半）开课哦：
生信入门&数据挖掘线上直播课5月班
文章复现：Integration of Single-Cell RNA Sequencing and Bulk RNA Sequencing Data to Establish and Validate a Prognostic Model for Patients With Lung Adenocarcinoma
全文数据分析，没有湿实验。在GEO数据库下载了一份scRNA的数据，两份基因芯片的数据，还下载了TCGA的数据。研究肺腺癌（LUAD），整合scRNA-seq和传统的RNA-seq数据来构建LUAD患者的预后模型，并采用两个外部验证队列来验证其风险分层能力。
文章整体思路如下：
单细胞降维分群注释，找出tumor和normal间关键细胞类型。单细胞功能富集，拟时序分析，细胞通讯分析TCGA数据进行差异分析，GO，KEGG功能富集TCGA数据进行WGCNA分析，找到hub gene后，取hub gene和TCGA DEGs的交集，进行后续分析单因素Cox回归鉴定潜在的预后DEGs，非负矩阵分解进行sample clustering（亚型），不同cluster的生存分析，免疫浸润分析单因素Cox回归筛选出的gene用LASSO模型进一步筛选，筛选完后进行多因素Cox回归，计算每个样本的risk score，根据risk score分成两类，进行生存分析，将构建的模型应用于两个GEO芯片数据，进行同样分析高危低危人群与临床信息相关性分析，分析risk score是否可作为预后指标，基因突变分析ssGSEA分析，免疫检查点与肿瘤突变负荷分析
本次复现只是为了复现图表，数据不一定具有生物学意义（大家可以理解为凑图或者灌水）
Step 1.单细胞数据处理
先下载好数据，一共4个样本就新建4个文件夹，名字命名成样本名，文件夹内带'GSM'开头的数据是下载的原始数据，重新复制一份，把文件名改成标准的read10x函数读入格式，一定注意名字不要写错。
一定要保证这3个文件同时存在，而且在同一个文件夹下面。
示例代码是：
rm(list=ls())options(stringsAsFactors = F)library(Seurat)sce1 <- CreateSeuratObject(Read10X('../10x-results/WT/'), "wt")
重点就是 Read10X 函数读取 文件夹路径，比如：../10x-results/WT/ ，保证文件夹下面有3个文件。
1.1 QC
library(Seurat)library(harmony)
library(dplyr)
library(stringr)
rm(list = ls())
# 多个数据读取与合并
rawdata_path <- './rawdata'
filename <- list.files(rawdata_path)
rawdata_path <- paste(rawdata_path,filename,sep = '/')
sceList <- lapply(rawdata_path, function(x){
obj <- CreateSeuratObject(counts = Read10X(x),
project = str_split(x,'/')[[1]][3])
})
names(sceList) <- filename
这里sceList是一个列表，列表里共4个元素，每个元素都是一个seurat对象，每个元素的名字就是上面的样本名。rawdata_path就是上面下载的.gz文件 所在文件夹的路径，Read10X函数只需要传入文件夹路径就可以自动读入数。
project参数传入的是一个字符串，一般用样本名即可，这里用str_split函数分割路径字符串得到样本名
sce <- merge(sceList[[1]],sceList[-1],add.cell.ids = names(sceList),project = 'luad')
sce@meta.data$group <- str_split(sce@meta.data$orig.ident,'_',simplify = T)[,2]
把4个seurat对象用merge函数合并后，看一下每个细胞的基本信息，这里我根据作者给出的信息，把疾病和正常分成两组，又增加了一列，ade1是疾病，nor是正常
# 查看线粒体(MT开头)、核糖体(RPS/RPL开头)、血红细胞所占比例
grep('^MT',x=rownames(sce@assays$RNA@data),value = T)
grep('^RP[SL]',x=rownames(sce@assays$RNA@data),value = T)
grep('^HB[^(P)]',x=rownames(sce@assays$RNA@data),value = T)
sce <- PercentageFeatureSet(sce,'^MT',col.name = 'percent_MT')
sce <- PercentageFeatureSet(sce,'^RP[SL]',col.name = 'percent_RP')
sce <- PercentageFeatureSet(sce,'^HB[^(P)]',col.name = 'percent_HB')
VlnPlot(sce,features = "nCount_RNA",pt.size = 0,y.max = 10000)
VlnPlot(sce,features = "nFeature_RNA",pt.size = 0,y.max = 2500)
VlnPlot(sce,features = "percent_MT",pt.size = 0)
VlnPlot(sce,features = "percent_RP",pt.size = 0)
VlnPlot(sce,features = "percent_HB",pt.size = 0,y.max = 0.1)
VlnPlot(sce,features = c("nCount_RNA","nFeature_RNA","percent_MT"),pt.size = 0,group.by = 'orig.ident')
对于每个细胞，我又添加了三列信息：线粒体基因占比，核糖体基因占比，红细胞基因占比。这几列用来过滤数据，这里重新给sce赋值是因为，PercentageFeatureSet函数的返回值还是一个seurat对象，和原来的sce相比， 就只是在meta data里多了一列而已，所以直接重新赋值，把原来的seurat对象sce给更新掉
过滤前先看下数据，32538个gene，13060个细胞
# 过滤细胞
sce <- subset(sce,subset = nCount_RNA>1000 & nFeature_RNA>300 & percent_MT<25)
# 过滤基因
sce <- sce[rowSums(sce@assays$RNA@counts>0)>3,]
# 过滤线粒体、核糖体、血红细胞进占比高的细胞
sce <- subset(sce,subset = percent_MT<25 & percent_RP<30 & percent_HB<0.1)
过滤完后22249个gene，7892个细胞
# 细胞周期评分
# S.Score较高的为S期，G2M.Score较高的为G2M期，都比较低的为G1期
s_feature <- cc.genes.updated.2019$s.genes
g2m_feature <- cc.genes.updated.2019$g2m.genes
sce <- CellCycleScoring(sce,
s.features = s_feature,
g2m.features = g2m_feature,
set.ident = T)
VlnPlot(sce,features = c('S.Score','G2M.Score'),group.by = 'orig.ident',pt.size = 0)
saveRDS(sce,'sce_qc.rds')
再看一下细胞周期，把细胞周期两个评分加到meta data里，更新sce
1.2 用harmony去除批次效应
library(harmony)library(dplyr)
library(Seurat)
library(clustree)
rm(list = ls())
sce <- readRDS('sce_qc.rds')
sce <- NormalizeData(sce,
normalization.method = 'LogNormalize',
scale.factor = 10000)
sce <- FindVariableFeatures(sce,
selection.method = "vst",
nfeatures = 2000)
# 默认用variableFeature做scale
sce <- ScaleData(sce)
sce <- RunPCA(sce,features = VariableFeatures(sce))
DimPlot(sce,reduction = 'pca',group.by = 'group')
从PCA降维的结果看，两个组的批次效应不大，均匀的混合在一起，但是这里还是用harmony试试
sce <- RunHarmony(sce,group.by.vars = 'orig.ident')
ElbowPlot(sce,reduction = 'harmony')
sce <- RunUMAP(sce,dims = 1:10,reduction = 'harmony')
sce <- RunTSNE(sce,dims = 1:10,reduction = 'harmony')
DimPlot(sce,reduction = 'umap',label = T,group.by = 'group')
DimPlot(sce,reduction = 'tsne',label = T,group.by = 'group')
选择orig.ident这一列，去除4个样本之间的批次效应。根据肘部图确定降维后的维度。
sce <- FindNeighbors(sce,reduction = 'harmony',dims = 1:10)
sce_res <- FindClusters(sce_res,resolution = 0.8)
# 设置不同的分辨率查看分组情况
sce_res <- sce
for (i in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)){
sce_res <- FindClusters(sce_res,resolution = i)
}
clustree(sce_res,prefix = 'RNA_snn_res.')
FindNeighbors函数计算细胞间的距离，FindClusters函数确定分群结果。上面的UMAP看着大概分了7坨，但实际上并不是7个分群，实际分群数量还是要看FindClusters的结果。FindClusters函数分辨率不同，分群数量会不同，具体分多少需要摸索。这里分辨率选了0.5，对应第六行蓝色，13个分群，如果选0.8的话，感觉分群数量会突然增加好多，不太合适。
sce <- FindClusters(sce,resolution = 0.5)
saveRDS(sce,file = 'step2_harmony.rds')
DimPlot(sce,reduction = 'umap',group.by = 'orig.ident')
DimPlot(sce,reduction = 'umap',group.by = 'seurat_clusters')
确定好分辨率后，分群赋值给sce，在meta data里会多出来两列，这两列是一样的，表示每个细胞分在哪个cluster里，最后保存一下更新后的sce对象，再来看下分群结果
1.3 细胞注释
library(SingleR)library(celldex)
library(Seurat)
library(dplyr)
library(stringr)
library(pheatmap)
library(ReactomeGSA)
library(monocle)
library(limma)
library(ggplot2)
library(msigdbr)
library(singleseqgset)
rm(list = ls())
gc()
sce <- readRDS('step2_harmony.rds')
# singleR注释
hpca.se <- HumanPrimaryCellAtlasData()
bpe.se <- BlueprintEncodeData()
anno <- SingleR(sce@assays$RNA@data,
ref = list(BP=bpe.se,HPCA=hpca.se),
labels = list(bpe.se$label.main,hpca.se$label.main),
clusters = sce@meta.data$seurat_clusters
)
用SingleR自动注释，HumanPrimaryCellAtlasData和BlueprintEncodeData是SingleR自带的两个数据库，用两个数据库注释效果好像比一个会好点。得到的anno也是一个对象，每个细胞群的注释信息就藏在里面，13个分群被注释成7种细胞，其中有6个分群都是上皮细胞，也许上皮细胞还能继续细分。想继续细分的话，对sce取子集，把标签为上皮细胞的细胞取出来，降维分群，丢到SingleR()函数里再注释
plotScoreHeatmap(anno,clusters = anno@rownames,show_colnames = T)
再看下注释结果的打分，不满意的话也可以手动注释
sce@meta.data$singleR_label <- unlist(lapply(sce@meta.data$seurat_clusters, function(x){anno$labels[x]}))
DimPlot(sce,reduction = 'umap',group.by = 'singleR_label',label = T)
DimPlot(sce,reduction = 'tsne',group.by = 'singleR_label',label = T)
anno里面有每个分群对应的注释信息，把这个信息添加到sce的meta data里，然后就可以在分群的图上添加上注释信息了
1.4 鉴定关键细胞类型
cell_list <- split(colnames(sce), sce$singleR_label)deg <- c()
for (i in names(cell_list)){
sce_temp <- sce[,colnames(sce) %in% cell_list[[i]]]
sce_temp <- CreateSeuratObject(counts = sce_temp@assays$RNA@counts,
meta.data = sce_temp@meta.data,
min.cells = 3,
min.features = 200)
sce_temp <- NormalizeData(sce_temp)
sce_temp <- FindVariableFeatures(sce_temp)
sce_temp <- ScaleData(sce_temp,vars.to.regress = c('nFeature_RNA','"percent_MT"'))
Idents(sce_temp) <- sce_temp$group
sce_temp_markers <- FindAllMarkers(object=sce_temp,
only.pos = T,
logfc.threshold = 2,
min.pct = 0.1)
sce_temp_markers <- sce_temp_markers[order(sce_temp_markers$cluster,sce_temp_markers$avg_log2FC),]
deg_gene <- sce_temp_markers[sce_temp_markers$p_val_adj<0.05,'gene']
deg <- c(deg,length(deg_gene))
}
df <- data.frame(cell_type = names(cell_list),deg_gene = deg)
把每种类型的细胞分别拿出来，用FindAllMarkers找出ade1和nor组之间的差异基因
1.5 单细胞通路富集
temp <- sceIdents(temp) <- temp@meta.data$singleR_label
reactome_obj <- analyse_sc_clusters(temp,create_reactome_visualization=T,use_interactors=F)
path_ways <- pathways(reactome_obj)
path_ways$diff <- rowMax(as.matrix(path_ways[,2:ncol(path_ways)])) - rowMin(as.matrix(path_ways[,2:ncol(path_ways)]))
path_ways <- path_ways[order(path_ways$diff,decreasing = T),]
plot_gsva_heatmap(reactome_obj,rownames(path_ways)[1:10],margins = c(10,20))
saveRDS(sce,file = 'step3_sce.rds')
用pathway()提取富集出来的结果，每行是一个通路，每列是一种细胞，计算每个通路在7种细胞中表达差异的diff，根据diff排序，画出前10个通路
1.6 拟时序分析
library(Seurat)library(monocle)
sce <- readRDS('step3_sce.rds')
# 1. 将seurat对象转化为monocle的CDS对象
sparse_data <- as(as.matrix(sce@assays$RNA@counts),'sparseMatrix')
mdata <- new('AnnotatedDataFrame',data=sce@meta.data) # 行为cell
fData <- data.frame(gene_short_name=row.names(sparse_data),row.names = row.names(sparse_data))
fd <- new('AnnotatedDataFrame',data=fData) # 行为gene 列为gene description
monocle_cds <- newCellDataSet(cellData = sparse_data,
phenoData = mdata,
featureData = fd,
lowerDetectionLimit = 0.5,
expressionFamily = negbinomial.size())
# 2. 计算size factor和离散度
monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)
# 3.过滤低质量细胞
monocle_cds <- detectGenes(monocle_cds,min_expr = 0.1)
# 4. 找高变基因、降维分群
# 用seurat找的高变基因
sce_var_gene <- VariableFeatures(sce)
monocle_cds <- setOrderingFilter(monocle_cds,sce_var_gene)
monocle_cds <- reduceDimension(monocle_cds,num_dim=10,norm_method = 'log',reduction_method = 'tSNE')
monocle_cds <- clusterCells(monocle_cds,num_clusters = 10)
plot_cell_clusters(monocle_cds,color_by = 'singleR_label')
diff_test_gene <- differentialGeneTest(monocle_cds[sce_var_gene,],fullModelFormulaStr = '~singleR_label')
diff_gene <- row.names(subset(diff_test_gene,qval<0.01))
monocle_cds <- setOrderingFilter(monocle_cds,diff_gene)
# 5. 用reversed graph embedding降维
monocle_cds <- reduceDimension(monocle_cds,reduction_method = 'DDRTree')
# 6. 计算细胞拟时间
monocle_cds <- orderCells(monocle_cds)
这里运行到第六步，orderCells()函数大概率会报错，详细的解决办法可以参考知乎上的贴子https://zhuanlan.zhihu.com/p/589134519 ，这里说一下大概解决步骤
下载一个包，'igraph’下载一个R文件，order_cells.R，链接在知乎贴子里，然后source调用把orderCells()中用DDRTree计算的那部分代码单独拿出来，写成一个函数，后面直接调用，记得给出返回值。或者不封装成函数的形式也行，那就不需要返回值了# 6. 计算细胞拟时间
source('order_cells.R')
library('igraph')
my_ordercell <- function(cds){
root_state = NULL
num_paths = NULL
reverse = NULL
root_cell <- select_root_cell(cds, root_state, reverse)
cds@auxOrderingData <- new.env(hash = TRUE)
if (cds@dim_reduce_type == "DDRTree") {
if (is.null(num_paths) == FALSE) {
message("Warning: num_paths only valid for method 'ICA' in reduceDimension()")
}
cc_ordering <- extract_ddrtree_ordering(cds, root_cell)
pData(cds)$Pseudotime <- cc_ordering[row.names(pData(cds)), ]$pseudo_time
K_old <- reducedDimK(cds)
old_dp <- cellPairwiseDistances(cds)
old_mst <- minSpanningTree(cds)
old_A <- reducedDimA(cds)
old_W <- reducedDimW(cds)
cds <- project2MST(cds, project_point_to_line_segment)
minSpanningTree(cds) <- cds@auxOrderingData[[cds@dim_reduce_type]]$pr_graph_cell_proj_tree
root_cell_idx <- which(V(old_mst)$name == root_cell, arr.ind = T)
cells_mapped_to_graph_root <- which(cds@auxOrderingData[["DDRTree"]]$pr_graph_cell_proj_closest_vertex == root_cell_idx)
if (length(cells_mapped_to_graph_root) == 0) {
cells_mapped_to_graph_root <- root_cell_idx
}
cells_mapped_to_graph_root <- V(minSpanningTree(cds))[cells_mapped_to_graph_root]$name
tip_leaves <- names(which(degree(minSpanningTree(cds)) == 1))
root_cell <- cells_mapped_to_graph_root[cells_mapped_to_graph_root %in% tip_leaves][1]
if (is.na(root_cell)) {
root_cell <- select_root_cell(cds, root_state, reverse)
}
cds@auxOrderingData[[cds@dim_reduce_type]]$root_cell <- root_cell
cc_ordering_new_pseudotime <- extract_ddrtree_ordering(cds, root_cell)
pData(cds)$Pseudotime <- cc_ordering_new_pseudotime[row.names(pData(cds)), ]$pseudo_time
if (is.null(root_state) == TRUE) {
closest_vertex <- cds@auxOrderingData[["DDRTree"]]$pr_graph_cell_proj_closest_vertex
pData(cds)$State <- cc_ordering[closest_vertex[, 1], ]$cell_state
}
}
cds
}
monocle_cds <- my_ordercell(monocle_cds)
# 可以画图了
plot_cell_trajectory(monocle_cds,color_by = 'Pseudotime')
plot_cell_trajectory(monocle_cds,color_by = 'State')
plot_cell_trajectory(monocle_cds,color_by = 'singleR_label')
plot_cell_trajectory(monocle_cds,color_by = 'singleR_label') facet_wrap(~singleR_label,nrow = 3)
saveRDS(monocle_cds,file = 'monocle.rds')
1.7 细胞通讯
library(CellChat)library(patchwork)
library(Seurat)
library(dplyr)
library(ggalluvial)
library(NMF)
options(stringsAsFactors = FALSE)
rm(list = ls())
gc()
sce <- readRDS('step3_sce.rds')
# Step1.创建cellchat对象
data.input <- sce@assays$RNA@data
meta.data <- sce@meta.data
cellchat <- createCellChat(object=data.input,
meta = meta.data,
group.by='singleR_label')
cellchat <- addMeta(cellchat,meta = meta.data)
# Step2.加载CellChatDB数据库
cellchatDB <- CellChatDB.human
cellchat@DB <- cellchatDB
# Step3.处理表达数据
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# Step4.计算通讯概率，推断细胞通讯网络 这一步很耗时
cellchat <- computeCommunProb(cellchat,population.size = F)
cellchat <- filterCommunication(cellchat,min.cells = 10)
# Step5. 提取预测的细胞通讯网络为data frame
df.net <- subsetCommunication(cellchat)
df.pathway <- subsetCommunication(cellchat,slot.name = 'netP')
# Step6. 在信号通路水平推断细胞通讯
cellchat <- computeCommunProbPathway(cellchat)
# Step7. 计算加和的cell-cell通讯网络
cellchat <- aggregateNet(cellchat)
par(mfrow = c(1,2))
netVisual_circle(cellchat@net$count,vertex.weight = groupSize,
weight.scale = T,label.edge = F,
title.name = 'number of Interaction')
netVisual_circle(cellchat@net$weight,vertex.weight = groupSize,
weight.scale = T,label.edge = F,
title.name = 'Interaction Weight')
mat <- cellchat@net$weight
par(mfrow = c(3,3),xpd=T)
for (i in 1:nrow(mat)){
mat2 <- matrix(0,nrow = nrow(mat),ncol = ncol(mat),
dimnames = dimnames(mat))
mat2[i,] <- mat[i,]
netVisual_circle(mat2,vertex.weight = groupSize,
weight.scale = T,edge.weight.max = max(mat),
title.name = rownames(mat)[i])
}
# Step8. 使用层次图（Hierarchical plot），圆圈图（Circle plot）或和弦图（Chord diagram）可视化每个信号通路
cellchat@netP$pathways
pathway.show <- c('SEMA4')
vertex.receiver <- c(1,2,3,4)
netVisual_aggregate(cellchat,signaling = pathway.show,
vertex.receiver = vertex.receiver,
layout = 'hierarchy')
par(mfrow=c(1,2),xpd=T)
netVisual_aggregate(cellchat,signaling = pathway.show,
layout = 'circle')
netVisual_aggregate(cellchat,signaling = pathway.show,
layout = 'circle',label.edge=T)
pairLR.CXCL <- extractEnrichedLR(cellchat,signaling = pathway.show,
geneLR.return = F)
LR.show <- pairLR.CXCL[1,]
netVisual_individual(cellchat,signaling = pathway.show,
pairLR.use = 'SEMA4A_PLXNB2',
layout = 'circle')
netVisual_aggregate(cellchat,signaling = pathway.show,layout = 'chord',
title.name = "Chord diagram 1")
Step 2. TCGA差异分析
library(TCGAbiolinks)library(SummarizedExperiment)
library(AnnoProbe)
library(stringr)
library(limma)
library(tinyarray)
library(data.table)
library(ggplot2)
library(dplyr)
library(tibble)
library(pheatmap)
library(clusterProfiler)
library(org.Hs.eg.db)
rm(list = ls())
# 下载LUAD表达矩阵
query <- GDCquery(project = 'TCGA-LUAD',
data.category = 'Transcriptome Profiling',
data.type = 'Gene Expression Quantification',
workflow.type = 'STAR - Counts')
GDCdownload(query)
dat <- GDCprepare(query)
exp <- assay(dat)
save(exp,file='TCGA_LUAD_exp.Rdata')
# 下载LUAD临床信息
query <- GDCquery(project = 'TCGA-LUAD',
data.category = 'Clinical',
data.type = 'Clinical Supplement',
file.type = 'xml')
GDCdownload(query)
dat <- GDCprepare_clinic(query,clinical.info = 'patient')
k = apply(dat, 2, function(x){!all(is.na(x))});table(k)
clinical <- dat[,k]
save(clinical,file = 'TCGA_LUAD_Clinical.Rdata')
注意一下TCGA的样本命名含义：
以这个为例：TCGA-A6-6650-01A-11R-1774-07
A6：Tissue source site，组织来源编码
6650：Participant, 参与者编号
01：Sample, 编号01~09表示肿瘤，10~19表示正常对照
A：Vial, 在一系列患者组织中的顺序，绝大多数样本该位置编码都是A;
很少数的是B，表示福尔马林固定石蜡包埋组织，已被证明用于测序分析的效果不佳，所以不建议使用-01B的样本数据
11：Portion, 同属于一个患者组织的不同部分的顺序编号，同一组织会分割为100-120mg的部分，分别使用
R：RNA, 分析的分子类型，对应关系如下所示：
1774：Plate, 在一系列96孔板中的顺序，值大表示制板越晚
07：Center, 测序或鉴定中心编码
# count 转换成TPM
gene_length <- fread('gene_len.txt')
colnames(gene_length) <- c('id','length')
gene_length <- as.data.frame(gene_length)
gene_length <- gene_length[!duplicated(gene_length$id),]
exp <- as.data.frame(exp)
exp$id <- unlist(lapply(rownames(exp),function(x){str_split(x,'\\.',simplify = T)[,1]}))
merge <- left_join(x=exp,y=gene_length,by='id')
merge <- na.omit(merge)
merge <- merge[!duplicated(merge$id),]
rownames(merge) <- merge$id
merge <- merge[,-601]
# 计算TPM
kb <- merge$length / 1000
rpk <- merge[,-601] / kb
tpm <- t(t(rpk) / colSums(rpk) * 1E6)
exp <- tpm
采用limma做差异分析，先把counts转换成TPM。不转也可以用DEseq2做差异分析。gene_length是每个基因的长度，这个文件是把gtf文件中，基因的终止位置-起始位置得到的。
# 更改行名：ENSEMBL -> SYMBOL
a <- annoGene(str_split(rownames(exp),'\\.',simplify = T)[,1],'ENSEMBL')
a <- a[!duplicated(a$ENSEMBL),]
rownames(a) <- a$ENSEMBL
rownames(exp) <- a[str_split(rownames(exp),'\\.',simplify = T)[,1],'SYMBOL']
# 给样本分组 tumor/normal
group <- ifelse(sapply(str_split(colnames(exp),'-',simplify = T)[,4],function(x){as.integer(substr(x,1,2))})>=10,
'normal','tumor')
sample_group <- data.frame(sample=colnames(exp),group=group)
a是annoGene函数返回的一个ID转换的结果，是一个dataframe，多个ENSEMBL ID可能对应同一个SYMBOL ID，对symbol去重后，再映射成exp的行名。
按照TCGA名字中第四部分，把样本分成tumor和normal
# 过滤
# 1. 去除FFPE样本
sample_group$sample_type <- ifelse(sapply(str_split(colnames(exp),'-',simplify = T)[,4],function(x){substr(x,3,3)})=='A',
'tissue','ffpe')
exp <- exp[,sample_group$sample_type!='ffpe']
sample_group <- sample_group[sample_group$sample_type=='tissue',]
# 2. 保留在一半以上样本中有值的gene
exp <- exp[rowSums(exp>0)>as.integer(ncol(exp)/2),]
dim(exp)
save(exp,file = 'TCGA_LUAD_TPM_exp.Rdata')
FFPE样本无论做测序还是蛋白组，代谢组，效果都不怎么好
# 对表达矩阵取log2
for (i in 1:ncol(exp)){
exp[,i] <- log2(exp[,i] 0.0001)
}
# 差异筛选
design <- model.matrix(~factor(sample_group$group,levels = c('tumor','normal')))
fit <- lmFit(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'))
deg <- na.omit(deg)
table(deg$change)
做完差异筛选，有10个gene含空值，删去有空值的行即可。更新下exp并保存好
# deg比exp少了几个gene
exp <- exp[intersect(rownames(exp),deg$ID),]
save(exp,file = 'TCGA_exp.Rdata')
save(deg,file = 'TCGA_DEGS.Rdata')
# volcano plot
ggplot(deg,aes(x=deg$logFC,y=-log10(deg$adj.P.Val),color=change))
geom_point()
geom_hline(yintercept = -log10(0.05),linetype=2,linewidth=1)
geom_vline(xintercept = log2(0.5),linetype=2,linewidth=1)
geom_vline(xintercept = log2(2),linetype=2,linewidth=1)
scale_color_manual(values=c('blue','gray','red'))
# heatmap
deg_heatmap <- deg[order(deg$logFC,decreasing = T),]
deg_gene <- deg_heatmap[deg_heatmap$change!='nochange','ID']
heatmap_deg_gene <- as.data.frame(exp)[deg_gene,]
heatmap_deg_gene <- na.omit(heatmap_deg_gene)
anno_col <- data.frame(group=sample_group$group)
rownames(anno_col) <- sample_group$sample
pheatmap(heatmap_deg_gene,show_rownames = F,cluster_rows = F,scale = 'row',
show_colnames = F,annotation_col=anno_col ,
breaks = seq(-2.5,2.5,length.out=100))
dim(heatmap_deg_gene)
# pca plot
exp <- as.data.frame(exp)
pca <- prcomp(t(exp))
ggplot(data.frame(pca$x),aes(x=PC1,y=PC2,color=sample_group$group))
geom_point(size = 3) stat_ellipse(level = 0.95, show.legend = F)
labs(title = 'PCA Analysis')
theme(plot.title = element_text(hjust = 0.5))
guides(color = guide_legend(title = 'Group'))
# GO enrich
go_bp <- enrichGO(deg_gene,org.Hs.eg.db,'SYMBOL','BP')
go_mf <- enrichGO(deg_gene,org.Hs.eg.db,'SYMBOL','MF')
go_cc <- enrichGO(deg_gene,org.Hs.eg.db,'SYMBOL','CC')
dotplot(go_bp)
barplot(go_bp,showCategory = 4)
# KEGG enrich
entrez_id <- bitr(deg_gene,fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = org.Hs.eg.db)
kegg <- enrichKEGG(entrez_id$ENTREZID,organism = 'hsa',keyType = 'ncbi-geneid')
dotplot(kegg)
barplot(kegg,showCategory = 4)
Step 3. WGCNA分析
library(WGCNA)library(reshape2)
library(stringr)
library(dplyr)
library(patchwork)
library(ggplot2)
rm(list = ls())
gc()
load('TCGA_DEGS.Rdata')
load('TCGA_exp.Rdata')
load('TCGA_LUAD_Clinical.Rdata')
options(stringsAsFactors = F) # 开启多线程
enableWGCNAThreads()
# 1. 数据处理
# >>> 选取表达量前10000的gene，用全部gene的话，计算后面的邻接矩阵和TOM矩阵太慢了
temp <- as.data.frame(exp)
temp$avg_exp <- rowMeans(temp)
temp <- temp[order(temp$avg_exp,decreasing = T),]
top10000_gene <- rownames(temp[1:10000,])
exp <- as.data.frame(exp[top10000_gene,])
# >>> 查找并删除异常值
temp <- exp
colnames(temp) <- 1:ncol(exp)
plot(hclust(dist(t(temp)))) # 没有异常样本
# 2. 构建基因共表达网络
exp <- as.data.frame(t(exp))
# >>> 2.1 选择合适的power
power_vec <- c(seq(1, 10, by = 1), seq(12, 20, by = 2))
sft <- pickSoftThreshold(exp,powerVector = power_vec)
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=12，达到0.8即可
power <- sft$powerEstimate
# >>> 2.2 检验所选power是否符合无尺度网络
k <- softConnectivity(exp,power = power)
scaleFreePlot(k)
power=12确实也符合无尺度网络，R2较高
# >>> 2.3 计算gene在sample间的corr或者distance，构建临近矩阵。默认用dist函数计算（欧氏距离）
adj <- adjacency(exp,power = power)
# >>> 2.4 根据临近矩阵计算TOM（拓扑重叠矩阵）
TOM <- TOMsimilarity(adj)
# >>> 2.5 计算gene间相异性
dissTOM <- 1-TOM
# 3. 对gene进行聚类并可视化
geneTree <- hclust(as.dist(dissTOM),method = 'average')
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)
# >>> 3.2 将gene聚类树进行修剪后，把gene归到不同的模块里
dynamic_modules <- cutreeDynamic(dendro = geneTree,distM = dissTOM,minClusterSize = 50,
deepSplit = 2,pamRespectsDendro=F,minSplitHeight = NULL)
table(dynamic_modules)
module_color <- labels2colors(dynamic_modules)
plotDendroAndColors(geneTree,module_color,dendroLabels = F,hang=0.03,addGuide = T,
guideHang = 0.05,groupLabels='Dynamic color tree',
main='gene dendrogram and module colors')
# >>> 3.3 计算每个模块的特征基因向量(对每个模块的gene表达量降维后，只保留PC1，即MEs)
MEList <- moduleEigengenes(exp,colors = dynamic_modules)
MEs <- MEList$eigengenes
# >>> 3.4 计算特征基因的相异度，基于相异度对原有模块进行层次聚类
dissME <- 1-cor(MEs)
MEtree <- hclust(as.dist(dissME),method = 'average')
plotEigengeneNetworks(MEs,
"Eigengene adjacency heatmap",
marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE,
xLabelsAngle = 90)
# >>> 3.5 如果模块太多，可以对模块进行合并.设置剪切高度,剪切高度以下的模块进行合并
cut_height <- 0.2
sizeGrWindow(8,6)
plot(MEtree,main = "Clustering of module eigengenes",
xlab = "", sub = "")
abline(h=cut_height,col='red')
把高度0.2以下的模块进行合并，再看看合并前后模块的变化
module_merge <- mergeCloseModules(exp,module_color,cutHeight = cut_height)
merged_color <- module_merge$colors
merged_MEs <- module_merge$newMEs
sizeGrWindow(8,6)
plotDendroAndColors(geneTree,cbind(module_color,merged_color),dendroLabels = F,hang=0.03,
addGuide = T,groupLabels = c("Dynamic Tree Cut", "Merged dynamic"))
# 4. 计算每个模块与临床信息的相关性、pvalue。找出与临床信息相关性高的模块
nGenes <- ncol(exp)
nSamples <- nrow(exp)
################################准备临床信息####################################
clinical <- clinical[,c('bcr_patient_barcode','gender','residual_tumor')]
traits <- data.frame(row.names = rownames(exp),ID=rownames(exp))
traits$ID <- substr(traits$ID,1,12)
index <- match(traits$ID,clinical$bcr_patient_barcode)
traits$gender <- as.numeric(as.factor(clinical[index,'gender']))
traits$residual_tumor <- as.numeric(as.factor(clinical[index,'residual_tumor']))
traits$type <- ifelse(sapply(str_split(rownames(traits),'-',simplify = T)[,4],
function(x){as.integer(substr(x,1,2))})>=10,'normal','tumor')
traits$type <- as.numeric(as.factor(traits$type))
traits <- traits[,-1]
################################准备临床信息####################################
clinical是前面下载的TCGA临床数据，这里随便取了两个比较关心的临床信息。clinical的每一行是一个病人的ID，而表达矩阵每一行是一个样本，一个病人身上可能有多个样本。在把临床信息和表达矩阵一一对应的时候，如果直接把exp的行名取前12位，就会产生重复值。在做的时候，先准备好样本和病人ID的对应关系，再把clinical里病人的信息给匹配过来
最终的临床信息是这样，要把字符串变成数值型变量
module_trait_cor <- cor(merged_MEs,traits,use = 'p')
module_trait_pvalue <- corPvalueStudent(module_trait_cor,nSamples)
trait_colors <- numbers2colors(traits,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(traits),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"))
模块与tumor/normal之间的相关性并不是很高。。。。。就拿粉色模块当做感兴趣的模块吧
# 5. 找出感兴趣的模块内的枢纽gene
# >>> 5.1 计算gene表达量与临床性状的相关性
gs <- apply(exp,2,function(x){cor(x,traits$type)})
# >>> 5.2 基于网络计算gene与临床性状的相关性
ns <- networkScreening(y=traits$type,datME = merged_MEs,datExpr = exp)
# >>> 5.3 计算每个gene与特征gene的相关性
kme <- signedKME(datExpr = exp,datME = merged_MEs)
# >>> 5.4 根据gs, ns, kme 找出hub gene
hubgene_filter <- (gs>0.2 & ns$q.Weighted<0.01 & abs(kme$kMEmagenta)>0.4)
hub_gene <- colnames(exp)[hubgene_filter]
plotNetworkHeatmap(exp,
plotGenes = hub_gene,
networkType = "unsigned",
useTOM = TRUE,
power=8,
main="unsigned correlations")
看下kme，其中粉色模块里，gene的相关性并不是很高，这里筛选了>0.4的，最后找到25个hub gene。如果阈值再设高一点，hub gene就全卡掉了。。。这跟没有用所有的3w多个gene去计算有关，也跟一些参数有关，我没有一个个去调整优化，只是展示一下如何做图和分析思路
hub_gene <- intersect(hub_gene,deg$ID)
save(hub_gene,file = 'TCGA_WGCNA_hub_gene.Rdata')
最后找出hub gene和前面deg的交集部分，作为hub gene，存档
Step 4. 非负矩阵分解算法
4.1 NMF算法对样本分类
library(NMF)library(survival)
library(stringr)
library(data.table)
library(survminer)
library(dplyr)
library(AnnoProbe)
library(CIBERSORT)
library(tidyr)
library(tibble)
library(RColorBrewer)
library(ggalluvial)
library(ggplot2)
options(stringsAsFactors = F)
rm(list = ls())
gc()
load('TCGA_LUAD_TPM_exp.Rdata')
load('TCGA_DEGS.Rdata')
load('TCGA_WGCNA_hub_gene.Rdata')
dim(exp)
# 对表达矩阵取log2
exp <- log2(exp 1)
ranks <- seq(2,10)
hubgene_exp <- exp[hub_gene,]
estimate <- nmf(hubgene_exp,2:10,nrun=50,method = 'brunet')
plot(estimate)
注意hubgene_exp里不能有负数，不然会报错，rank为2时，波动最大
rank <- 2
nmf.rank4 <- nmf(hubgene_exp,rank,nrun=50,method = 'brunet')
index <- extractFeatures(nmf.rank4,'max')
index <- unlist(index)
nmf.exp <- hubgene_exp[index,]
nmf.exp <- na.omit(nmf.exp)
group <- predict(nmf.rank4)
jco <- c("#2874C5","#EABF00","#C6524A","#868686")
consensusmap(nmf.rank4,labRow=NA,labCol=NA,
annCol = data.frame("cluster"=group[colnames(nmf.exp)]),
annColors = list(cluster=c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])))
cluster_df <- data.frame(row.names = colnames(exp),ID=substr(colnames(exp),1,12),cluster=group)
NMF算法把所有样本分成两类，顺便看下最后的clueter_df
4.2 生存分析
load('TCGA_LUAD_TPM_exp.Rdata')load('TCGA_LUAD_Clinical.Rdata')
# 1. 取log（TPM）
for (i in 1:ncol(exp)){
exp[,i] <- log2(exp[,i] 0.0001)
}
# 2. 以病人为中心，去重
exp <- as.data.frame(exp)
exp <- exp[,sort(colnames(exp))]
index <- !duplicated(substr(colnames(exp),1,12))
exp <- exp[,index]
# 3.读取从xena下载的生存数据
sur <- fread('TCGA-LUAD.survival.tsv')
sur <- as.data.frame(sur)
patient_id <- substr(colnames(exp),1,12)
clin <- data.frame(row.names = patient_id,ID=patient_id)
index <- match(clin$ID,sur[,'_PATIENT'])
clin <- cbind(clin,sur[index,c('OS','OS.time')])
clin <- na.omit(clin)
index2 <- match(clin$ID,clinical$bcr_patient_barcode)
clin <- cbind(clin,clinical[index2,'gender'])
colnames(clin) <- c(colnames(clin)[1:ncol(clin)-1],'gender')
clin$gender <- as.numeric(as.factor(clin$gender))
下载下来的生存数据长这样
之前从TCGA下载的临床数据长这样
合并后的clin是这样，所有变量类型都要转成数值型
# 合并NMF结果
cluster_type_group <- ifelse(sapply(str_split(rownames(cluster_df),'-',simplify = T)[,4],
function(x){as.integer(substr(x,1,2))}
)>=10,
'normal','tumor')
cluster_df <- cluster_df[which(cluster_type_group=='tumor'),]
clin <- left_join(clin,cluster_df,'ID',multiple='first')
再把NMF分类的结果添加进来，就是cluster那一列
# 7. 过滤生存数据，不要<30天的，再把天数转换成月
clin <- clin[clin$OS.time>30,]
clin$OS.time <- clin$OS.time / 30
# 8. exp的列名和临床信息的行名一一对应。NMF将样本分成若干个亚型，画出亚型的生存曲线
colnames(exp) <- substr(colnames(exp),1,12)
index <- match(clin$ID,colnames(exp))
exp <- exp[,index]
sfit <- survfit(Surv(OS.time,OS)~cluster,data=clin)
ggsurvplot(sfit,pval = T)
save(clin,file = 'clin.Rdata')
4.3 免疫浸润
my_pallet <- colorRampPalette(brewer.pal(8,'Set1'))lm22f = system.file("extdata", "LM22.txt", package = "CIBERSORT")
cluster1 <- exp[,rownames(cluster_df[cluster_df$cluster==1,])]
cluster2 <- exp[,rownames(cluster_df[cluster_df$cluster==2,])]
TME.result1 <- cibersort(lm22f,as.matrix(cluster1),perm = 100,QN=T) # 耗时很久
TME.result2 <- cibersort(lm22f,as.matrix(cluster2),perm = 100,QN=T)
原文中免疫浸润用MCPcounter做的，这里我用CIBERSORT做，运行完的结果就是各种免疫细胞在每个样本中的含量
result1 <- TME.result1[,-c(23:25)] %>% as.data.frame() %>%
rownames_to_column('Sample') %>% gather(key = cell_type,value = proportion,-Sample)
result2 <- TME.result2[,-c(23:25)] %>% as.data.frame() %>%
rownames_to_column('Sample') %>% gather(key = cell_type,value = proportion,-Sample)
result1$cluster <- rep('cluster1', times=nrow(result1))
result2$cluster <- rep('cluster2', times=nrow(result2))
dat <- rbind(result1,result2)
ggplot(dat,aes(axis1 = cluster, axis2 = cell_type))
geom_alluvium(aes(fill = as.factor(cell_type)),width = 2/5, discern = FALSE)
geom_stratum(width = 2/5, discern = FALSE)
geom_text(stat = "stratum", discern = FALSE,aes(label = after_stat(stratum)))
theme(legend.position = "none",#去除刻度线和背景颜色
panel.background = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(size =15,colour = "black"),#坐标轴名
axis.title = element_blank())
scale_x_discrete(position = "top")
把上面的透视表转换成长表，添加上cluster的信息，再把cluster1和cluster2对应的长表合并起来，就可以画图了
cell <- 'Plasma cells'
ggviolin(dat[dat$cell_type==cell,],x='cluster',y='proportion',fill = 'cluster',
add = 'boxplot',add.params = list(fill = 'white',width=0.1,linetype=1),title = cell)
stat_compare_means(method = 't.test',label = 'p.format',aes(label = 'p.format'),
label.x.npc = 'center',label.y = 0.5,size=5)
4.4 批量单因素Cox回归
load('TCGA_WGCNA_hub_gene.Rdata')load('TCGA_LUAD_TPM_exp.Rdata')
load('clin.Rdata')
hubgene_exp <- exp[hub_gene,]
hubgene_exp <- as.data.frame(t(hubgene_exp))
hubgene_exp$ID <- substr(rownames(hubgene_exp),1,12)
clin <- left_join(clin,hubgene_exp,by='ID')
把WGCNA中找到的hub gene挑出来，拿到他们在所有样本中的表达值。把表达值和生存信息合并到一起，最后的clin长这样。在合并之前要注意，删掉normal样本，我这里没有做。
res.cox <- coxph(Surv(OS.time,OS)~SAPCD2,data=clin)
sumres.cox <- summary(res.cox)
挑一个基因SAPCD2做单因素cox回归，最后需要的信息都藏在sumres.cox里
# 批量单因素cox回归
gene <- c()
p_value <- c()
HR <- c()
lower95 <- c()
upper95 <- c()
for (i in 6:ncol(clin)){
res <- coxph(Surv(OS.time,OS)~clin[,i],data=clin)
sum_res <- summary(res)
p <- sum_res$coefficients[,'Pr(>|z|)']
if (p<0.05){
gene <- c(gene,colnames(clin)[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_input_gene <- rownames(cox_df[cox_df$pvalue<0.01,])
save(lasso_input_gene,file = 'Lasso_input_gene.Rdata')
批量单因素cox回归得到每个gene的pvalue，挑出pvalue<0.01的再进行后续筛选
Step 5. 预后模型构建和验证
load('Lasso_input_gene.Rdata')load('TCGA_LUAD_TPM_exp.Rdata')
load('clin.Rdata')
# 对表达矩阵取log2
exp <- log2(exp 0.00001)
exp <- as.data.frame(exp)
# Lasso回归进一步筛选gene
exp2 <- exp[lasso_input_gene,]
exp2 <- as.data.frame(t(exp2))
exp2$ID <- substr(rownames(exp2),1,12)
exp2 <- exp2[!duplicated(exp2$ID),]
clin <- clin[,1:3]
exp3 <- left_join(clin,exp2,by='ID')
rownames(exp3) <- exp3$ID
exp3 <- exp3[,-1]
exp3 <- na.omit(exp3)
exp3是生存数据和上一步批量单因素Cox回归筛选出来的基因的表达量合并起来的结果
x <- as.matrix(exp3[,3:(ncol(exp3)-1)])
y <- Surv(exp3$OS.time,exp3$OS)
fit1 <- glmnet(x,y,alpha = 1,family = 'cox')
plot(fit1)
cv_fit <- cv.glmnet(x,y,alpha = 1,nfolds = 10,family="cox")
plot(cv_fit)
best_lambda <- cv_fit$lambda.min
fit2 <- glmnet(x,y,lambda = best_lambda,alpha = 1,family = 'cox')
lasso_filter_gene <- names(fit2$beta[fit2$beta[,1]!=0,1])
由于之前筛选出来的gene效果不怎么好， 经过lasso过滤后，只剩2个gene了，就用这两个来做多因素cox回归
# 多因素cox回归构建预后模型
exp4 <- exp3 %>% dplyr::select(OS,OS.time,lasso_filter_gene)
multicox <- coxph(Surv(OS.time,OS)~.,data = exp4)
sum_multicox <- summary(multicox)
riskScore <- predict(multicox,type = 'risk',newdata = exp4)
riskScore <- as.data.frame(riskScore)
riskScore$ID <- rownames(riskScore)
riskScore$risk <- ifelse(riskScore$riskScore>median(riskScore$riskScore),'high','low')
做完多因素回归后，计算每个样本的riskscore，用riskScore把样本分成两类
# KM plot
km_data <- left_join(riskScore,clin,by='ID')
fit <- survfit(Surv(OS.time,OS)~risk,data=km_data)
ggsurvplot(fit,data = km_data,pval = T,risk.table = T,surv.median.line = 'hv',
legend.title = 'RiskScore',title = 'Overall survival',
ylab='Cummulative survival',xlab='Time(Days)')
# 绘制ROC曲线
roc_data <- data.frame(x = roc$FP[,1],
y = roc$TP[,1],
time = roc$times[1])
for (i in 2:length(roc$times)) {
temp <- data.frame(x = roc$FP[,i],
y = roc$TP[,i],
time = roc$times[i])
roc_data <- rbind(roc_data, temp)
}
ggplot(roc_data, aes(x = x, y = y, color = as.factor(time)))
geom_line()
scale_color_manual(values = c("gray", "blue", "black"))
geom_abline(slope = 1, intercept = 0, linetype = "dashed")
xlab("1 - Specificity")
ylab("Sensitivity")
ggtitle("ROC Curve")
theme_minimal()
# 绘制风险图
ggrisk(multicox)
# risk score与临床病理特征的相关性
load('TCGA_LUAD_Clinical.Rdata')
clinical_df <- clinical[,c('bcr_patient_barcode','gender','stage_event_pathologic_stage')]
clinical_df <- merge(km_data,clinical_df,by.x = 'ID',by.y = 'bcr_patient_barcode')
clinical_df <- clinical_df[!duplicated(clinical_df$ID),]
clinical_df$stage_event_pathologic_stage <- case_when(clinical_df$stage_event_pathologic_stage=='Stage IA'~'Stage I',
clinical_df$stage_event_pathologic_stage=='Stage IB'~'Stage I',
clinical_df$stage_event_pathologic_stage=='Stage IIA'~'Stage II',
clinical_df$stage_event_pathologic_stage=='Stage IIB'~'Stage II',
clinical_df$stage_event_pathologic_stage=='Stage IIIA'~'Stage III',
clinical_df$stage_event_pathologic_stage=='Stage IIIB'~'Stage III',
)
ggviolin(clinical_df,x='gender',y='riskScore',fill = 'gender',
add = 'boxplot',add.params = list(fill = 'white',width=0.1,linetype=1),)
stat_compare_means(method = 'wilcox.test',label = 'p.format',aes(label = 'p.format'),
label.x.npc = 'center',label.y = 20,size=5)
ggviolin(clinical_df,x='stage_event_pathologic_stage',y='riskScore',fill = 'stage_event_pathologic_stage',
add = 'boxplot',add.params = list(fill = 'white',width=0.1,linetype=1),)
stat_compare_means(method = 'wilcox.test',label = 'p.signif',aes(label = 'p.format'),
label.x = 1,label.y = 20,size=5)
Step 6. riskScore和临床信息的相关性
# 森林图# 批量单因素回归
df_forest <- clinical_df
df_forest$gender <- as.numeric(as.factor(df_forest$gender))
df_forest$stage_event_pathologic_stage <- as.numeric(as.factor(df_forest$stage_event_pathologic_stage))
gene <- c()
p_value <- c()
HR <- c()
lower95 <- c()
upper95 <- c()
means <- c()
for (i in c('riskScore','gender','stage_event_pathologic_stage')){
res <- coxph(Surv(OS.time,OS)~df_forest[,i],data=df_forest)
sum_res <- summary(res)
p <- sum_res$coefficients[,'Pr(>|z|)']
gene <- c(gene,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'])
means <- c(means,res$means)
}
cox_df <- data.frame(row.names = gene,pvalue=p_value,HR=HR,lower.95=lower95,upper.95=upper95,means=means)
cox_df$pvalue <- round(cox_df$pvalue,3)
cox_df$HR <- round(cox_df$HR,3)
cox_df <- rownames_to_column(cox_df,var = 'Variable')
cox_df <- rbind(colnames(cox_df),cox_df)
forestplot(labeltext=as.matrix(cox_df)[,1:3],lower=as.numeric(cox_df[,4]),
upper=as.numeric(cox_df[,5]),mean=as.numeric(cox_df[,6]),
is.summary=c(T,F,F,F),
zero=0,lineheight=unit(0.5,'cm'),graphwidth=unit(35,'mm'),
colgap=unit(2,'mm'),lwd.zero=2,lwd.ci=3,
col=fpColors(box = '#7AC5CD',lines = 'black',zero = 'purple'),
lwd.xaxis=2,lty.ci='solid',ci.vertices.height=0.05,graph.pos=2,
xlim=c(0,2.5),box.size=0.5,
hrzl_lines=list('2'=gpar(lwd=3,columns=1:4,col='dark green'),
'5'=gpar(lwd=3,columns=1:4,col='dark green'))
)
# 多因素回归
res <- coxph(Surv(OS.time,OS)~riskScore gender stage_event_pathologic_stage,data=df_forest)
sum_res <- summary(res)
df <- data.frame(row.names = rownames(sum_res$conf.int),pvalue=sum_res$coefficients[,'Pr(>|z|)'],
HR=sum_res$conf.int[,'exp(coef)'],lower=sum_res$conf.int[,'lower .95'],
upper=sum_res$conf.int[,'upper .95'],means=res$means)
df$pvalue <- round(df$pvalue,3)
df$HR <- round(df$HR,3)
df <- rownames_to_column(df,var='Variable')
df <- rbind(colnames(df),df)
forestplot(labeltext=as.matrix(df)[,1:3],lower=as.numeric(df[,4]),
upper=as.numeric(df[,5]),mean=as.numeric(df[,6]),
is.summary=c(T,F,F,F),
zero=0,lineheight=unit(0.5,'cm'),graphwidth=unit(35,'mm'),
colgap=unit(2,'mm'),lwd.zero=2,lwd.ci=3,
col=fpColors(box = '#7AC5CD',lines = 'black',zero = 'purple'),
lwd.xaxis=2,lty.ci='solid',ci.vertices.height=0.05,graph.pos=2,
xlim=c(0,2.5),box.size=0.5,
hrzl_lines=list('2'=gpar(lwd=3,columns=1:4,col='dark green'),
'5'=gpar(lwd=3,columns=1:4,col='dark green'))
)
# GSEA
df_exp <- as.data.frame(t(exp))
df_exp <- rownames_to_column(df_exp,var = 'sample_ID')
df_exp$sample_ID <- substr(df_exp$sample_ID,1,12)
df_exp <- df_exp[!duplicated(df_exp$sample_ID),]
rownames(df_exp) <- df_exp$sample_ID
df_exp <- as.data.frame(t(df_exp[,-1]))
df_exp <- df_exp[,km_data$ID]
group <- as.factor(km_data$risk)
design <- model.matrix(~km_data$risk)
fit <- lmFit(df_exp,design)
fit <- eBayes(fit)
df_deg <- topTable(fit,coef = 2,number = Inf)
df_deg$change <- ifelse(df_deg$logFC>log2(2) & df_deg$adj.P.Val<0.05,'up',
ifelse(df_deg$logFC<log2(0.5) & df_deg$adj.P.Val<0.05,'down','nochange'))
df_deg <- df_deg[order(df_deg$logFC,decreasing = T),]
symbol_2_entrez <- bitr(rownames(df_deg),fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = org.Hs.eg.db)
symbol_2_entrez <- symbol_2_entrez[!duplicated(symbol_2_entrez$SYMBOL),]
rownames(symbol_2_entrez) <- symbol_2_entrez$SYMBOL
df_deg$entrezid <- symbol_2_entrez[rownames(df_deg),'ENTREZID']
df_deg <- df_deg[!duplicated(df_deg$entrezid),]
df_deg <- na.omit(df_deg)
a <- df_deg[(df_deg$change=='up'),'logFC']
names(a) <- df_deg[(df_deg$change=='up'),'entrezid']
res_a <- gseKEGG(a,'hsa')
b <- df_deg[(df_deg$change=='down'),'logFC']
names(b) <- df_deg[(df_deg$change=='down'),'entrezid']
res_b <- gseKEGG(b,'hsa')
gseaplot2(res_a,geneSetID = 1:5)
gseaplot2(res_b,geneSetID = 1:5)
先做差异分析，找到high risk vs low risk的差异基因，上调和下调的gene分别拿去做GSEA再画图
# waterfall plot
query <- GDCquery(project = 'TCGA-LUAD',
data.category = "Simple Nucleotide Variation",
data.type = "Masked Somatic Mutation",
access = 'open')
GDCdownload(query)
GDCprepare(query,save = T,save.filename = 'TCGA-LUAD-SNP.Rdata')
load('TCGA-LUAD-SNP.Rdata')
luad_snp <- read.maf(maf=data)
plotmafSummary(luad_snp,rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
oncoplot(luad_snp,top = 15)
Step 7. 高低风险组间的免疫浸润
library(GSVA)library(GSEABase)
library(ggplot2)
library(tidyr)
rm(list = ls())
#1. 获取geneSets
# 下载28种免疫细胞的参考基因集 </TISIDB/data/download/CellReports.txt>
geneSet <- read.csv("cellReport.txt",header = F,sep = "\t",)
geneSet[1:5,1:5]
geneSet <- t(tibble::column_to_rownames(geneSet,var = 'V1'))
gene_list <- list()
for (i in colnames(geneSet)){
gene_list[[i]] <- geneSet[,i]
gene_list[[i]] <-gene_list[[i]][gene_list[[i]]!='']
}
直接读取的geneSet长这样，有空值，需要转变成gene_list，list长度为28，其中每个元素是某个细胞所有基因的向量，名字是对应的细胞
# 载入表达矩阵,把表达矩阵分成high risk 和low risk组
load('TCGA_LUAD_TPM_exp.Rdata')
load('risk.Rdata')
# 对表达矩阵取log2
exp <- log2(exp 0.00001)
exp <- as.data.frame(exp)
colnames(exp) <- substr(colnames(exp),1,12)
exp <- exp[,!duplicated(colnames(exp))]
high_risk_sample <- km_data[km_data$risk=='high','ID']
low_risk_sample <- km_data[km_data$risk=='low','ID']
high_risk_exp <- exp[,high_risk_sample]
low_risk_exp <- exp[,low_risk_sample]
# ssGSEA
high_risk_gsea <- gsva(high_risk_exp,gene_list,method='ssgsea',kcdf='Gaussian',abs.ranking=T)
low_risk_gsea <- gsva(low_risk_exp,gene_list,method='ssgsea',kcdf='Gaussian',abs.ranking=T)
要注意gsva函数的输入得是matrix，不能是data.frame，不然会报错
hr_gsea_norm <- high_risk_gsea
for (i in colnames(hr_gsea_norm)) {
diff <- max(hr_gsea_norm[,i] )-min(hr_gsea_norm[,i] )
hr_gsea_norm[,i] <- (hr_gsea_norm[,i] -min(hr_gsea_norm[,i]))/diff
}
hr_gsea_norm <- as.data.frame(t(hr_gsea_norm))
hr_gsea_norm$type <- rep('high',times=nrow(hr_gsea_norm))
lr_gsea_norm <- low_risk_gsea
for (i in colnames(lr_gsea_norm)) {
diff <- max(lr_gsea_norm[,i] )-min(lr_gsea_norm[,i])
lr_gsea_norm[,i] <- (lr_gsea_norm[,i] -min(lr_gsea_norm[,i]))/diff
}
lr_gsea_norm <- as.data.frame(t(lr_gsea_norm))
lr_gsea_norm$type <- rep('low',times=nrow(lr_gsea_norm))
gsea_norm <- rbind(hr_gsea_norm,lr_gsea_norm)
longer <- pivot_longer(gsea_norm,cols = !type,names_to = 'cell_type',values_to = 'Score')
ggplot(longer,aes(x=cell_type,y=Score,color=type)) geom_boxplot()
theme(axis.text.x = element_text(angle = 45,hjust = 1))
stat_compare_means(label = "p.signif",size=3,method = 'wilcox.test')
先把high risk gsea和low risk gsea的数据分别进行min max scale，新增一列标记上high和low，然后合并起来，行是样本，列是细胞。再把透视表转成长表，就可以画图了
# 免疫检查点&TBM
risk_exp <- as.data.frame(t(exp))
risk_exp <- risk_exp[km_data$ID,]
risk_exp$type <- km_data$risk
risk_exp <- risk_exp[,c('CD274','PDCD1','CTLA4','LAG3','TIGIT','type')]
risk_long <- pivot_longer(risk_exp,cols = !type,names_to = 'ICIs',values_to = 'Expression')
ggplot(risk_long,aes(x=ICIs,y=Expression,fill=type)) geom_boxplot()
stat_compare_means(label = "p.format",size=3,method = 'wilcox.test')
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
生物信息学马拉松授课（买一得五） ，你的生物信息学入门课144线程640Gb内存服务器共享一年仍然是仅需800千呼万唤始出来的独享生物信息学云服务器
0000- 0000
0000
0000- 0000
