文章复现：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数据进行WGCNA分析，找到hub gene后，取hub gene和TCGA DEGs的交集，进行后续分析 单因素Cox回归鉴定潜在的预后DEGs，非负矩阵分解进行sample clustering（亚型），不同cluster的生存分析，免疫浸润分析 单因素Cox回归筛选出的gene用LASSO模型进一步筛选，筛选完后进行多因素Cox回归，计算每个样本的risk score，根据risk score分成两类，进行生存分析，将构建的模型应用于两个GEO芯片数据，进行同样分析 高危低危人群与临床信息相关性分析，分析risk score是否可作为预后指标，基因突变分析 本次复现只是为了复现图表，数据不一定具有生物学意义（大家可以理解为凑图或者灌水）
Step 1.单细胞数据处理 先下载好数据，一共4个样本就新建4个文件夹，名字命名成样本名，文件夹内带'GSM'开头的数据是下载的原始数据，重新复制一份，把文件名改成标准的read10x函数读入格式，一定注意名字不要写错。
一定要保证这3个文件同时存在，而且在同一个文件夹下面。
示例代码是：
rm(list=ls())F )library (Seurat)'../10x-results/WT/' ),"wt" )重点就是 Read10X 函数读取 文件夹路径，比如：../10x-results/WT/ ，保证文件夹下面有3个文件。
1.1 QC library(Seurat)# 多个数据读取与合并 './rawdata' '/' )function (x){'/' )[[1]][3])这里sceList是一个列表，列表里共4个元素，每个元素都是一个seurat对象，每个元素的名字就是上面的样本名。rawdata_path就是上面下载的.gz文件 所在文件夹的路径，Read10X函数只需要传入文件夹路径就可以自动读入数。
sce <- merge(sceList[[1]],sceList[-1],add.cell.ids = names(sceList),project = 'luad' )$group <- str_split(sce@meta.data$orig .ident,'_' ,simplify = T)[,2]把4个seurat对象用merge函数合并后，看一下每个细胞的基本信息，这里我根据作者给出的信息，把疾病和正常分成两组，又增加了一列，ade1是疾病，nor是正常
# 查看线粒体(MT开头)、核糖体(RPS/RPL开头)、血红细胞所占比例 '^MT' ,x=rownames(sce@assays$RNA @data),value = T)'^RP[SL]' ,x=rownames(sce@assays$RNA @data),value = T)'^HB[^(P)]' ,x=rownames(sce@assays$RNA @data),value = T)'^MT' ,col.name = 'percent_MT' )'^RP[SL]' ,col.name = 'percent_RP' )'^HB[^(P)]' ,col.name = 'percent_HB' )"nCount_RNA" ,pt.size = 0,y.max = 10000)"nFeature_RNA" ,pt.size = 0,y.max = 2500)"percent_MT" ,pt.size = 0)"percent_RP" ,pt.size = 0)"percent_HB" ,pt.size = 0,y.max = 0.1)"nCount_RNA" ,"nFeature_RNA" ,"percent_MT" ),pt.size = 0,group.by = 'orig.ident' )对于每个细胞，我又添加了三列信息：线粒体基因占比，核糖体基因占比，红细胞基因占比。这几列用来过滤数据，这里重新给sce赋值是因为，PercentageFeatureSet函数的返回值还是一个seurat对象，和原来的sce相比， 就只是在meta data里多了一列而已，所以直接重新赋值，把原来的seurat对象sce给更新掉
# 过滤细胞 # 过滤基因 $RNA @counts>0)>3,]# 过滤线粒体、核糖体、血红细胞进占比高的细胞 过滤完后22249个gene，7892个细胞
# 细胞周期评分 # S.Score较高的为S期，G2M.Score较高的为G2M期，都比较低的为G1期 $s .genes$g2m .genes'S.Score' ,'G2M.Score' ),group.by = 'orig.ident' ,pt.size = 0)'sce_qc.rds' )再看一下细胞周期，把细胞周期两个评分加到meta data里，更新sce
1.2 用harmony去除批次效应 library(harmony)'sce_qc.rds' )'LogNormalize' ,"vst" ,# 默认用variableFeature做scale 'pca' ,group.by = 'group' )从PCA降维的结果看，两个组的批次效应不大，均匀的混合在一起，但是这里还是用harmony试试
sce <- RunHarmony(sce,group.by.vars = 'orig.ident' )'harmony' )'harmony' )'harmony' )'umap' ,label = T,group.by = 'group' )'tsne' ,label = T,group.by = 'group' )选择orig.ident这一列，去除4个样本之间的批次效应。根据肘部图确定降维后的维度。
sce <- FindNeighbors(sce,reduction = 'harmony' ,dims = 1:10)# 设置不同的分辨率查看分组情况 for (i in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)){'RNA_snn_res.' )FindNeighbors函数计算细胞间的距离，FindClusters函数确定分群结果。上面的UMAP看着大概分了7坨，但实际上并不是7个分群，实际分群数量还是要看FindClusters的结果。FindClusters函数分辨率不同，分群数量会不同，具体分多少需要摸索。这里分辨率选了0.5，对应第六行蓝色，13个分群，如果选0.8的话，感觉分群数量会突然增加好多，不太合适。
sce <- FindClusters(sce,resolution = 0.5)'step2_harmony.rds' )'umap' ,group.by = 'orig.ident' )'umap' ,group.by = 'seurat_clusters' )确定好分辨率后，分群赋值给sce，在meta data里会多出来两列，这两列是一样的，表示每个细胞分在哪个cluster里，最后保存一下更新后的sce对象，再来看下分群结果
1.3 细胞注释 library(SingleR)'step2_harmony.rds' )# singleR注释 $RNA @data,$label .main,hpca.se$label .main),$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]}))'umap' ,group.by = 'singleR_label' ,label = T)'tsne' ,group.by = 'singleR_label' ,label = T)
1.4 鉴定关键细胞类型 cell_list <- split(colnames(sce), sce$singleR_label )for (i in names(cell_list)){in % cell_list[[i]]]$RNA @counts,'nFeature_RNA' ,'"percent_MT"' ))$group $cluster ,sce_temp_markers$avg_log2FC ),]$p_val_adj <0.05,'gene' ]把每种类型的细胞分别拿出来，用FindAllMarkers找出ade1和nor组之间的差异基因
1.5 单细胞通路富集 temp <- sce$singleR_label $diff <- rowMax(as.matrix(path_ways[,2:ncol(path_ways)])) - rowMin(as.matrix(path_ways[,2:ncol(path_ways)]))$diff ,decreasing = T),]'step3_sce.rds' )用pathway()提取富集出来的结果，每行是一个通路，每列是一种细胞，计算每个通路在7种细胞中表达差异的diff，根据diff排序，画出前10个通路
1.6 拟时序分析 library(Seurat)'step3_sce.rds' )# 1. 将seurat对象转化为monocle的CDS对象 $RNA @counts),'sparseMatrix' )'AnnotatedDataFrame' ,data=sce@meta.data) # 行为cell 'AnnotatedDataFrame' ,data=fData) # 行为gene 列为gene description # 2. 计算size factor和离散度 # 3.过滤低质量细胞 # 4. 找高变基因、降维分群 # 用seurat找的高变基因 'log' ,reduction_method = 'tSNE' )'singleR_label' )'~singleR_label' )# 5. 用reversed graph embedding降维 'DDRTree' )# 6. 计算细胞拟时间 这里运行到第六步，orderCells()函数大概率会报错，详细的解决办法可以参考知乎上的贴子https://zhuanlan.zhihu.com/p/589134519 ，这里说一下大概解决步骤
下载一个R文件，order_cells.R，链接在知乎贴子里，然后source调用 把orderCells()中用DDRTree计算的那部分代码单独拿出来，写成一个函数，后面直接调用，记得给出返回值。 # 6. 计算细胞拟时间 source ('order_cells.R' )'igraph' )function (cds){hash = TRUE)if (cds@dim_reduce_type == "DDRTree" ) {if (is.null(num_paths) == FALSE) {"Warning: num_paths only valid for method 'ICA' in reduceDimension()" )$Pseudotime <- cc_ordering[row.names(pData(cds)), ]$pseudo_time $pr_graph_cell_proj_tree which (V(old_mst)$name == root_cell, arr.ind = T)which (cds@auxOrderingData[["DDRTree" ]]$pr_graph_cell_proj_closest_vertex == root_cell_idx)if (length(cells_mapped_to_graph_root) == 0) {$name which (degree(minSpanningTree(cds)) == 1))in % tip_leaves][1]if (is.na(root_cell)) {$root_cell <- root_cell$Pseudotime <- cc_ordering_new_pseudotime[row.names(pData(cds)), ]$pseudo_time if (is.null(root_state) == TRUE) {"DDRTree" ]]$pr_graph_cell_proj_closest_vertex $State <- cc_ordering[closest_vertex[, 1], ]$cell_state # 可以画图了 'Pseudotime' )'State' )'singleR_label' )'singleR_label' )+facet_wrap(~singleR_label,nrow = 3)'monocle.rds' )
1.7 细胞通讯 library(CellChat)'step3_sce.rds' )# Step1.创建cellchat对象 $RNA @data'singleR_label' )# Step2.加载CellChatDB数据库 # Step3.处理表达数据 # Step4.计算通讯概率，推断细胞通讯网络 这一步很耗时 # Step5. 提取预测的细胞通讯网络为data frame 'netP' )# Step6. 在信号通路水平推断细胞通讯 # Step7. 计算加和的cell-cell通讯网络 $count ,vertex.weight = groupSize,'number of Interaction' )$weight ,vertex.weight = groupSize,'Interaction Weight' )$weight for (i in 1:nrow(mat)){# Step8. 使用层次图（Hierarchical plot），圆圈图（Circle plot）或和弦图（Chord diagram）可视化每个信号通路 $pathways 'SEMA4' )'hierarchy' )par(mfrow=c(1,2),xpd=T)'circle' )'circle' ,label.edge=T)pairLR.CXCL <- extractEnrichedLR(cellchat,signaling = pathway.show,'SEMA4A_PLXNB2' ,'circle' )'chord' ,"Chord diagram 1" )
Step 2. TCGA差异分析 library(TCGAbiolinks)# 下载LUAD表达矩阵 'TCGA-LUAD' ,'Transcriptome Profiling' ,'Gene Expression Quantification' ,'STAR - Counts' )'TCGA_LUAD_exp.Rdata' )# 下载LUAD临床信息 'TCGA-LUAD' ,'Clinical' ,'Clinical Supplement' ,'xml' )'patient' )function (x){!all(is.na(x))});table(k)'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_len.txt' )'id' ,'length' )$id ),]$id <- unlist(lapply(rownames(exp),function (x){str_split(x,'\\.' ,simplify = T)[,1]}))'id' )$id ),]$id # 计算TPM $length / 1000采用limma做差异分析，先把counts转换成TPM。不转也可以用DEseq2做差异分析。gene_length是每个基因的长度，这个文件是把gtf文件中，基因的终止位置-起始位置得到的。
# 更改行名：ENSEMBL -> SYMBOL '\\.' ,simplify = T)[,1],'ENSEMBL' )$ENSEMBL ),]$ENSEMBL '\\.' ,simplify = T)[,1],'SYMBOL' ]# 给样本分组 tumor/normal '-' ,simplify = T)[,4],function (x){as.integer(substr(x,1,2))})>=10,'normal' ,'tumor' )a是annoGene函数返回的一个ID转换的结果，是一个dataframe，多个ENSEMBL ID可能对应同一个SYMBOL ID，对symbol去重后，再映射成exp的行名。
# 过滤 # 1. 去除FFPE样本 $sample_type <- ifelse(sapply(str_split(colnames(exp),'-' ,simplify = T)[,4],function (x){substr(x,3,3)})=='A' ,'tissue' ,'ffpe' )$sample_type !='ffpe' ]$sample_type =='tissue' ,]# 2. 保留在一半以上样本中有值的gene 'TCGA_LUAD_TPM_exp.Rdata' )FFPE样本无论做测序还是蛋白组，代谢组，效果都不怎么好
# 对表达矩阵取log2 for (i in 1:ncol(exp)){# 差异筛选 $group ,levels = c('tumor' ,'normal' )))$change <- ifelse((deg$logFC >log2(2))&(deg$adj .P.Val<0.05),'up' ,$logFC <log2(0.5))&(deg$adj .P.Val<0.05),'down' ,'nochange' ))$change )
# deg比exp少了几个gene $ID ),]'TCGA_exp.Rdata' )'TCGA_DEGS.Rdata' )# volcano plot $logFC ,y=-log10(deg$adj .P.Val),color=change))+'blue' ,'gray' ,'red' ))# heatmap $logFC ,decreasing = T),]$change !='nochange' ,'ID' ]$group )$sample 'row' ,# pca plot $x ),aes(x=PC1,y=PC2,color=sample_group$group ))+'PCA Analysis' )+'Group' ))
# GO enrich 'SYMBOL' ,'BP' )'SYMBOL' ,'MF' )'SYMBOL' ,'CC' )# KEGG enrich 'SYMBOL' ,toType = 'ENTREZID' ,OrgDb = org.Hs.eg.db)$ENTREZID ,organism = 'hsa' ,keyType = 'ncbi-geneid' )
Step 3. WGCNA分析 library(WGCNA)'TCGA_DEGS.Rdata' )'TCGA_exp.Rdata' )'TCGA_LUAD_Clinical.Rdata' )# 开启多线程 # 1. 数据处理 # >>> 选取表达量前10000的gene，用全部gene的话，计算后面的邻接矩阵和TOM矩阵太慢了 $avg_exp <- rowMeans(temp)$avg_exp ,decreasing = T),]# >>> 查找并删除异常值 # 没有异常样本 # 2. 构建基因共表达网络 # >>> 2.1 选择合适的power $fitIndices ,aes(x=Power,y=-sign(slope)*SFT.R.sq))+$fitIndices $Power ,color='red' )+'red' )+"Soft Threshold (power)" )+ylab("Scale Free Topology Model Fit,signed R^2" )+"Scale independence" )+theme(plot.title = element_text(hjust = 0.5))$fitIndices ,aes(x=Power,y=mean.k.))+$fitIndices $Power ,color='red' )+"Soft Threshold (power)" )+ylab("Mean Connectivity" )+'Mean connectivity' )+theme(plot.title = element_text(hjust = 0.5))
power <- sft$powerEstimate # >>> 2.2 检验所选power是否符合无尺度网络 power=12确实也符合无尺度网络，R2较高
# >>> 2.3 计算gene在sample间的corr或者distance，构建临近矩阵。默认用dist函数计算（欧氏距离） # >>> 2.4 根据临近矩阵计算TOM（拓扑重叠矩阵） # >>> 2.5 计算gene间相异性 # 3. 对gene进行聚类并可视化 'average' )"" , sub="" , main = "Gene clustering on TOM-based dissimilarity" ,# >>> 3.2 将gene聚类树进行修剪后，把gene归到不同的模块里 'Dynamic color tree' ,'gene dendrogram and module colors' )# >>> 3.3 计算每个模块的特征基因向量(对每个模块的gene表达量降维后，只保留PC1，即MEs) $eigengenes # >>> 3.4 计算特征基因的相异度，基于相异度对原有模块进行层次聚类 'average' )"Eigengene adjacency heatmap" , # >>> 3.5 如果模块太多，可以对模块进行合并.设置剪切高度,剪切高度以下的模块进行合并 "Clustering of module eigengenes" , "" , sub = "" )'red' )
module_merge <- mergeCloseModules(exp,module_color,cutHeight = cut_height)$colors $newMEs "Dynamic Tree Cut" , "Merged dynamic" ))# 4. 计算每个模块与临床信息的相关性、pvalue。找出与临床信息相关性高的模块 ################################准备临床信息#################################### 'bcr_patient_barcode' ,'gender' ,'residual_tumor' )]$ID <- substr(traits$ID ,1,12)$ID ,clinical$bcr_patient_barcode )$gender <- as.numeric(as.factor(clinical[index,'gender' ]))$residual_tumor <- as.numeric(as.factor(clinical[index,'residual_tumor' ]))$type <- ifelse(sapply(str_split(rownames(traits),'-' ,simplify = T)[,4],function (x){as.integer(substr(x,1,2))})>=10,'normal' ,'tumor' )$type <- as.numeric(as.factor(traits$type ))################################准备临床信息#################################### clinical是前面下载的TCGA临床数据，这里随便取了两个比较关心的临床信息。clinical的每一行是一个病人的ID，而表达矩阵每一行是一个样本，一个病人身上可能有多个样本。在把临床信息和表达矩阵一一对应的时候，如果直接把exp的行名取前12位，就会产生重复值。在做的时候，先准备好样本和病人ID的对应关系，再把clinical里病人的信息给匹配过来
module_trait_cor <- cor(merged_MEs,traits,use = 'p' )"\n(" , signif(module_trait_pvalue, 1), ")" , sep = "" )"Module-trait relationships" ))模块与tumor/normal之间的相关性并不是很高。。。。。就拿粉色模块当做感兴趣的模块吧
# 5. 找出感兴趣的模块内的枢纽gene # >>> 5.1 计算gene表达量与临床性状的相关性 function (x){cor(x,traits$type )})# >>> 5.2 基于网络计算gene与临床性状的相关性 $type ,datME = merged_MEs,datExpr = exp)# >>> 5.3 计算每个gene与特征gene的相关性 # >>> 5.4 根据gs, ns, kme 找出hub gene $q .Weighted<0.01 & abs(kme$kMEmagenta )>0.4)"unsigned" ,"unsigned correlations" )看下kme，其中粉色模块里，gene的相关性并不是很高，这里筛选了>0.4的，最后找到25个hub gene。如果阈值再设高一点，hub gene就全卡掉了。。。这跟没有用所有的3w多个gene去计算有关，也跟一些参数有关，我没有一个个去调整优化，只是展示一下如何做图和分析思路
hub_gene <- intersect(hub_gene,deg$ID )'TCGA_WGCNA_hub_gene.Rdata' )最后找出hub gene和前面deg的交集部分，作为hub gene，存档
Step 4. 非负矩阵分解算法 4.1 NMF算法对样本分类 library(NMF)load('TCGA_LUAD_TPM_exp.Rdata' )'TCGA_DEGS.Rdata' )'TCGA_WGCNA_hub_gene.Rdata' )# 对表达矩阵取log2 'brunet' )注意hubgene_exp里不能有负数，不然会报错，rank为2时，波动最大
rank <- 2'brunet' )'max' )"#2874C5" ,"#EABF00" ,"#C6524A" ,"#868686" )"cluster" =group[colnames(nmf.exp)]),"1" =jco[1],"2" =jco[2],"3" =jco[3],"4" =jco[4])))NMF算法把所有样本分成两类，顺便看下最后的clueter_df
4.2 生存分析 load('TCGA_LUAD_TPM_exp.Rdata' )'TCGA_LUAD_Clinical.Rdata' )# 1. 取log（TPM） for (i in 1:ncol(exp)){# 2. 以病人为中心，去重 # 3.读取从xena下载的生存数据 'TCGA-LUAD.survival.tsv' )$ID ,sur[,'_PATIENT' ])'OS' ,'OS.time' )])$ID ,clinical$bcr_patient_barcode )'gender' ])'gender' )$gender <- as.numeric(as.factor(clin$gender ))
# 合并NMF结果 '-' ,simplify = T)[,4],function (x){as.integer(substr(x,1,2))}'normal' ,'tumor' )which (cluster_type_group=='tumor' ),]'ID' ,multiple='first' )再把NMF分类的结果添加进来，就是cluster那一列
# 7. 过滤生存数据，不要<30天的，再把天数转换成月 $OS .time>30,]$OS .time <- clin$OS .time / 30# 8. exp的列名和临床信息的行名一一对应。NMF将样本分成若干个亚型，画出亚型的生存曲线 $ID ,colnames(exp))'clin.Rdata' )4.3 免疫浸润 my_pallet <- colorRampPalette(brewer.pal(8,'Set1' ))"extdata" , "LM22.txt" , package = "CIBERSORT" )$cluster ==1,])]$cluster ==2,])]# 耗时很久 原文中免疫浸润用MCPcounter做的，这里我用CIBERSORT做，运行完的结果就是各种免疫细胞在每个样本中的含量
result1 <- TME.result1[,-c(23:25)] %>% as.data.frame() %>% 'Sample' ) %>% gather(key = cell_type,value = proportion,-Sample)'Sample' ) %>% gather(key = cell_type,value = proportion,-Sample)$cluster <- rep('cluster1' , times =nrow(result1))$cluster <- rep('cluster2' , times =nrow(result2))stat = "stratum" , discern = FALSE,aes(label = after_stat(stratum)))+"none" ,#去除刻度线和背景颜色 "black" ),#坐标轴名 "top" )把上面的透视表转换成长表，添加上cluster的信息，再把cluster1和cluster2对应的长表合并起来，就可以画图了
cell <- 'Plasma cells' $cell_type ==cell,],x='cluster' ,y='proportion' ,fill = 'cluster' ,'boxplot' ,add.params = list(fill = 'white' ,width=0.1,linetype=1),title = cell)+'t.test' ,label = 'p.format' ,aes(label = 'p.format' ),'center' ,label.y = 0.5,size=5)4.4 批量单因素Cox回归 load('TCGA_WGCNA_hub_gene.Rdata' )'TCGA_LUAD_TPM_exp.Rdata' ) 'clin.Rdata' )$ID <- substr(rownames(hubgene_exp),1,12)'ID' )把WGCNA中找到的hub gene挑出来，拿到他们在所有样本中的表达值。把表达值和生存信息合并到一起，最后的clin长这样。在合并之前要注意，删掉normal样本，我这里没有做。
res.cox <- coxph(Surv(OS.time,OS)~SAPCD2,data=clin)挑一个基因SAPCD2做单因素cox回归，最后需要的信息都藏在sumres.cox里
# 批量单因素cox回归 for (i in 6:ncol(clin)){$coefficients [,'Pr(>|z|)' ]if (p<0.05){$conf .int[,'exp(coef)' ])$conf .int[,'lower .95' ])$conf .int[,'upper .95' ])$pvalue <0.01,])'Lasso_input_gene.Rdata' )批量单因素cox回归得到每个gene的pvalue，挑出pvalue<0.01的再进行后续筛选
Step 5. 预后模型构建和验证 load('Lasso_input_gene.Rdata' )'TCGA_LUAD_TPM_exp.Rdata' )'clin.Rdata' )# 对表达矩阵取log2 # Lasso回归进一步筛选gene $ID <- substr(rownames(exp2),1,12)$ID ),]'ID' )$ID exp3是生存数据和上一步批量单因素Cox回归筛选出来的基因的表达量合并起来的结果
x <- as.matrix(exp3[,3:(ncol(exp3)-1)])$OS .time,exp3$OS )'cox' )"cox" )$lambda .min'cox' )$beta [fit2$beta [,1]!=0,1])由于之前筛选出来的gene效果不怎么好， 经过lasso过滤后，只剩2个gene了，就用这两个来做多因素cox回归
# 多因素cox回归构建预后模型 type = 'risk' ,newdata = exp4)$ID <- rownames(riskScore)$risk <- ifelse(riskScore$riskScore >median(riskScore$riskScore ),'high' ,'low' )做完多因素回归后，计算每个样本的riskscore，用riskScore把样本分成两类
# KM plot 'ID' )'hv' ,'RiskScore' ,title = 'Overall survival' ,'Cummulative survival' ,xlab='Time(Days)' )# 绘制ROC曲线 $FP [,1],$TP [,1],$times [1])for (i in 2:length(roc$times )) {$FP [,i],$TP [,i],$times [i])"gray" , "blue" , "black" )) +"dashed" ) +"1 - Specificity" ) + "Sensitivity" ) + "ROC Curve" ) + # 绘制风险图 # risk score与临床病理特征的相关性 'TCGA_LUAD_Clinical.Rdata' )'bcr_patient_barcode' ,'gender' ,'stage_event_pathologic_stage' )]'ID' ,by.y = 'bcr_patient_barcode' )$ID ),]$stage_event_pathologic_stage <- case_when(clinical_df$stage_event_pathologic_stage =='Stage IA' ~'Stage I' ,$stage_event_pathologic_stage =='Stage IB' ~'Stage I' ,$stage_event_pathologic_stage =='Stage IIA' ~'Stage II' ,$stage_event_pathologic_stage =='Stage IIB' ~'Stage II' ,$stage_event_pathologic_stage =='Stage IIIA' ~'Stage III' ,$stage_event_pathologic_stage =='Stage IIIB' ~'Stage III' ,'gender' ,y='riskScore' ,fill = 'gender' ,'boxplot' ,add.params = list(fill = 'white' ,width=0.1,linetype=1),)+'wilcox.test' ,label = 'p.format' ,aes(label = 'p.format' ),'center' ,label.y = 20,size=5)'stage_event_pathologic_stage' ,y='riskScore' ,fill = 'stage_event_pathologic_stage' ,'boxplot' ,add.params = list(fill = 'white' ,width=0.1,linetype=1),)+'wilcox.test' ,label = 'p.signif' ,aes(label = 'p.format' ),
Step 6. riskScore和临床信息的相关性 # 森林图 # 批量单因素回归 $gender <- as.numeric(as.factor(df_forest$gender ))$stage_event_pathologic_stage <- as.numeric(as.factor(df_forest$stage_event_pathologic_stage ))for (i in c('riskScore' ,'gender' ,'stage_event_pathologic_stage' )){$coefficients [,'Pr(>|z|)' ]$conf .int[,'exp(coef)' ])$conf .int[,'lower .95' ])$conf .int[,'upper .95' ])$means )$pvalue <- round(cox_df$pvalue ,3)$HR <- round(cox_df$HR ,3)'Variable' )'cm' ),graphwidth=unit(35,'mm' ),'mm' ),lwd.zero=2,lwd.ci=3,'#7AC5CD' ,lines = 'black' ,zero = 'purple' ),'solid' ,ci.vertices.height=0.05,graph.pos=2,'2' =gpar(lwd=3,columns=1:4,col='dark green' ),'5' =gpar(lwd=3,columns=1:4,col='dark green' ))# 多因素回归 $conf .int),pvalue=sum_res$coefficients [,'Pr(>|z|)' ],$conf .int[,'exp(coef)' ],lower=sum_res$conf .int[,'lower .95' ],$conf .int[,'upper .95' ],means=res$means )$pvalue <- round(df$pvalue ,3)$HR <- round(df$HR ,3)'Variable' )'cm' ),graphwidth=unit(35,'mm' ),'mm' ),lwd.zero=2,lwd.ci=3,'#7AC5CD' ,lines = 'black' ,zero = 'purple' ),'solid' ,ci.vertices.height=0.05,graph.pos=2,'2' =gpar(lwd=3,columns=1:4,col='dark green' ),'5' =gpar(lwd=3,columns=1:4,col='dark green' ))# GSEA 'sample_ID' )$sample_ID <- substr(df_exp$sample_ID ,1,12)$sample_ID ),]$sample_ID $ID ]$risk )$risk )$change <- ifelse(df_deg$logFC >log2(2) & df_deg$adj .P.Val<0.05,'up' ,$logFC <log2(0.5) & df_deg$adj .P.Val<0.05,'down' ,'nochange' ))$logFC ,decreasing = T),]'SYMBOL' ,toType = 'ENTREZID' ,OrgDb = org.Hs.eg.db)$SYMBOL ),]$SYMBOL $entrezid <- symbol_2_entrez[rownames(df_deg),'ENTREZID' ]$entrezid ),]$change =='up' ),'logFC' ]$change =='up' ),'entrezid' ]'hsa' )$change =='down' ),'logFC' ]$change =='down' ),'entrezid' ]'hsa' )先做差异分析，找到high risk vs low risk的差异基因，上调和下调的gene分别拿去做GSEA再画图
# waterfall plot 'TCGA-LUAD' ,"Simple Nucleotide Variation" ,"Masked Somatic Mutation" ,'open' )'TCGA-LUAD-SNP.Rdata' )'TCGA-LUAD-SNP.Rdata' )'median' , dashboard = TRUE, titvRaw = FALSE)
Step 7. 高低风险组间的免疫浸润 library(GSVA)#1. 获取geneSets # 下载28种免疫细胞的参考基因集 <http://cis./TISIDB/data/download/CellReports.txt> "cellReport.txt" ,header = F,sep = "\t" ,)'V1' ))for (i in colnames(geneSet)){'' ] 直接读取的geneSet长这样，有空值，需要转变成gene_list，list长度为28，其中每个元素是某个细胞所有基因的向量，名字是对应的细胞
# 载入表达矩阵,把表达矩阵分成high risk 和low risk组 'TCGA_LUAD_TPM_exp.Rdata' )'risk.Rdata' )# 对表达矩阵取log2 $risk =='high' ,'ID' ]$risk =='low' ,'ID' ]# ssGSEA 'ssgsea' ,kcdf='Gaussian' ,abs.ranking=T)'ssgsea' ,kcdf='Gaussian' ,abs.ranking=T)要注意gsva函数的输入得是matrix，不能是data.frame，不然会报错
hr_gsea_norm <- high_risk_gseafor (i in colnames(hr_gsea_norm)) {$type <- rep('high' ,times =nrow(hr_gsea_norm))for (i in colnames(lr_gsea_norm)) {$type <- rep('low' ,times =nrow(lr_gsea_norm))type ,names_to = 'cell_type' ,values_to = 'Score' )type ))+geom_boxplot()+"p.signif" ,size=3,method = 'wilcox.test' )先把high risk gsea和low risk gsea的数据分别进行min max scale，新增一列标记上high和low，然后合并起来，行是样本，列是细胞。再把透视表转成长表，就可以画图了
# 免疫检查点&TBM $ID ,]$type <- km_data$risk 'CD274' ,'PDCD1' ,'CTLA4' ,'LAG3' ,'TIGIT' ,'type' )]type ,names_to = 'ICIs' ,values_to = 'Expression' )type ))+geom_boxplot()+"p.format" ,size=3,method = 'wilcox.test' )
project参数传入的是一个字符串，一般用样本名即可，这里用str_split函数分割路径字符串得到样本名
