TCGA的乳腺癌RNA-seq数据WGCNA分析示例
一文看懂WGCNA分析（2019更新版）
一文看懂WGCNA分析（2019更新版）发现我这个4年前的WGCNA分析教程可以排在自己最受欢迎的前10个教程里面了，而且直接以我这个授课代码出的SCI文章就有38篇了，当然不排除很多学员使用我的代码却不告知我，也不会致谢我。
不过，我这点战绩根本就算不上什么，其实这个WGCNA包已经是十多年前发表的了，仍然是广受好评及引用量一直在增加，破万也是指日可待。
大家首先可以看到3个教程：•2016-WGCNA-HCC-hub-gene.pdf 中文文章范例）•WGCNA_GBMTutorialHorvath.pdf•WGCNA_YeastTutorialHorvath.pdf其中第一个是我4年前的WGCNA分析教程最主要的参考文献，后面两个是英文教程，我相信你大概率是不会去看的，不过，我还是放在这里了。
（还是需要强调，这两个英文教程完整的展现了WGCNA的全部用法）然后你只需要简单浏览本文档，就可以在rstudio里面打开后缀是proj的文件，打开R代码，一步步跟着店铺！基本概念WGCNA其译为加权基因共表达网络分析。
该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。
适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。
一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。
基本原理从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。
第一步计算任意两个基因之间的相关系数（Person Coefficient）。
为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。
但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。
xgene：WGS，突变与癌，RNA-seq，WES
xgene：WGS，突变与癌，RNA-seq，WES⼈类全基因组测序06SNP(single nucleotide polymorphism)：有了10倍以上的覆盖深度以后，来确认SNP信息，就相当可靠了。
⼀个普通黄种⼈的基因组，与hg19这个参考基因组序列相⽐，会有350万个左右的SNP。
⼜有⼤概2万个是落在外显⼦上的，⽽⾮同义的SNP有⼤概9千个。
所谓⾮同义的SNP，就是这些SNP是会引起蛋⽩质的序列变化的。
indel：(insertion & deletion)是指⼩于50个bp以内的微⼩的插⼊、和缺失突变。
⼀个普通黄种⼈的基因组和hg19相⽐，约有50万个Indel。
其中落在外显⼦上的，⼤概在1千个左右。
那么Indel如果⼀旦落在外显⼦区域，它⼀定会引起蛋⽩质序列变化的。
如果它引起的是移码突变，那么在移码位点之后，所有氨基酸序列就和原来的序列完全不同。
如果它（基因）还能保持原来的阅读框，也会引起蛋⽩质中若⼲个氨基酸的增或者减。
SV： structure variation 染⾊体结构变异 1、染⾊体内部的位移2、染⾊体之间的位移3、⼤⽚段的缺失4、⼤⽚段的插⼊5、⼤⽚倍的加倍6、⼤⽚段的倒位 CNV ：copy number variation 拷贝数变异，是指染⾊体⽚段的拷贝数变异：包括拷贝数增加，也包括拷贝数减少。
实际上，CNV是和结构变异（也就是SV）紧密相关的。
SV 中的⼤⽚段的增加、和⼤⽚段的缺失，会直接导致CNV的变化。
突变种类与癌症04基因拷贝数异常： 例如：HER2基因，如果HER2基因的拷贝数增加到6个，或者更多，它就⽐较容易引发乳腺癌。
赫赛汀（Herceptin）这个药，可以抑制HER2蛋⽩的活性，所以赫赛洒就对于由HER2基因拷贝数异常增加引发的乳腺癌，有⾮常好的治疗作⽤。
染⾊体结构变异： 强启动⼦替换了弱启动⼦，改变了某个基因在天然条件下的表达量。
TCGA的乳腺癌RNA-seq数据WGCNA分析示例
TCGA的乳腺癌RNA-seq数据WGCNA分析示例WGCNA（WeightedCorrelationNetworkanalyi）是一个基于基因表达数据，构建基因共表达网络的方法。
WGCNA和差异基因分析（DEG）的差异在于DEG主要分析样本和样本之间的差异，而WGCNA主要分析的是基因和基因之间的关系。
WGCNA通过分析基因之间的关联关系，将基因区分为多个模块。
而最后通过这些模块和样本表型之间的关联性分析，寻找特定表型的分子特征。
网上例子千千万，但是大部分都是从文档翻译而来，要用起来还是有些费劲，要深入的可以移步这里：WGCNA##############etwd('E:/rawData/TCGA_DATA/TCGA-BRCA')ample=read.cv('ClinicalFull_matri某.t某t',ep='\\t',=1)dim(ample)#[1]1003e某pro=read.cv('Merge_matri某.t某t.cv.t某t',ep='\\t',=1)dim(e某pro)#[1]24991100数据读取完成，从上述结果可以看出100个样本，有24991个基因，这么多基因全部用来做WGCNA很显然没有必要，我们只要选择一些具有代表性的基因就够了，这里我们采取的方式是选择在100个样本中方差较大的那些基因（意味着在不同样本中变化较大）继续命令：m.var=apply(e某pro,1,var)e某pro.upper=e某pro[which(m.var>quantile(m.var,prob=eq(0,1,0.25))[4]),]##选择方差最大的前25%个基因作为后续WGCNA的输入数据集通过上述步骤拿到了6248个基因的表达谱作为WGCNA的输入数据集，进一步的我们需要看看样本之间的差异情况datE某pr=a.data.frame(t(e某pro.upper));gg=goodSampleGene(datE某pr,verboe=3);gg$allOK ampleTree=hclut(dit(datE某pr),method='average')plot(ampleTree,main='Samplecluteringtodetec toutlier',ub='',某lab='')从图中可看出大部分样本表现比较相近，而有两个离群样本，对后续的分析可能造成影响，我们需要将其去掉，共得到98个样本clut=cutreeStatic(ampleTree,cutHeight=80000,minSize=10)table(clut )#clut#01#298keepSample=(clut==1)datE某pr=datE某pr[keepSample,]nGene=ncol(datE某pr)nSample=nrow(datE某pr) ave(datE某pr,file='FPKM-01-dataInput.RData')得到最终的数据矩阵之后，我们需要确定软阈值，从代码中可以看出pickSoftThrehold 很简单，就两个参数，其他默认即可power=c(c(1:10),eq(from=12,to=20,by=2))ft=pickSoftThrehold(datE某pr,powerVector=power,verboe=5)##画图##par(mfrow=c(1,2));ce某1=0.9;plot(ft$fitIndice[,1],-ign(ft$fitIndice[,3])某ft$fitIndice[,2],某lab='SoftThrehold(power)',ylab='ScaleFreeTopologyModelFit,ignedR ^2',type='n',main=pate('Scaleindependence'));te某t(ft$fitIndice[,1],-ign(ft$fitIndice[,3])某ft$fitIndice[,2],label=power,ce某=ce某1,col='red');abline(h=0.90,col='red')plot(ft$fitIndice[,1],ft$fitIndice[,5],某lab='SoftThrehold(power)',ylab='MeanConnectivity',type='n', main=pate('Meanconnectivity'))te某t(ft$fitIndice[,1],ft$fitIndice[,5],label=power,ce某=ce某1,col='red')从图中可以看出这个软阈值选择7比较合适,选择软阈值7进行共表达模块挖掘pow=7net=blockwieModule(datE某pr,power=pow,ma某BlockSize=7000, TOMType='unigned',minModuleSize=30,reaignThrehold=0,mergeCutHeight=0.25,numericLabel=TRUE,pamRepectDendro=FALSE,aveTOM=TRUE,aveTOMFileBae='FPKM-TOM',verboe=3)table(net$color)#openagraphicwindow#izeGrWindow(12 ,9)#ConvertlabeltocolorforplottingmergedColor=label2color(net$color)#PlotthedendrogramandthemodulecolorunderneathplotDendroAndCo lor(net$dendrogram[[1]],mergedColor[net$blockGene[[1]]], groupLabel=c('Modulecolor','GS.weight'),dendroLabel=FALSE,ha ng=0.03,addGuide=TRUE,guideHang=0.05)从图中可以看出大部分基因在灰色区域，灰色部分一般认为是没有模块接受的，从这里也可以看出其实咱们选择的这些基因并不是特别好那么做到这一步了基本上共表达模块做完了，每个颜色代表一个共表达模块，统计看看各个模块下的基因个数：那么得到模块之后下一步该做啥呢，或许很多人到这就不知道如何继续分析了这里就需要咱们利用这些模块搞事情了，举个例子如果你是整合的数据（整合lnc与gene），那么同时在某个模块中的基因和lncRNA咱们可以认为是共表达的，这便是lnc-gene共表达关系的获得途径之一了，进一步你可以根据该模块的基因-lnc-基因之间的关系绘制出共表达网络今天咱们这里不讲这个，而是跟表型关联，咱们已经拿到了这98个样本的ER、PR、HER2阳性阴性信息，那么进一步的咱们可以看看哪些共表达模块跟ER、PR、HER2阴性最相关，代码如下：moduleLabelAutomatic= net$colormoduleColorAutomatic=label2color(moduleLabelAutomatic)moduleColorFemale=moduleCol orAutomaticME0=moduleEigengene(datE某pr,moduleColorFemale)$eigengeneMEFemale=orderME(ME0)ample=ample[match((datE某pr),pate0(gub('-','.',(ample)),'.01')),]#匹配98个样本数据trainDt=a.matri某(cbind(ifele(ample[,1]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,2]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,3]=='Poitive',0,1),#将阴性的样本标记为1ifele(ample[,1]=='Negative'&le[,2]=='Negative'&le[,3]= ='Negative',1,0))#将三阴性的样本标记为1)#得到一个表型的0-1矩阵modTraitCor=cor(MEFemale,trainDt,ue='p')colname(MEFemale)modTraitP=corPvalueStudent(modTraitCor,nSample)te某tMatri某=pate(ignif(modTraitCor,2),'\\n(',ignif(modTraitP,1),')',ep='')d im(te某tMatri某)=dim(modTraitCor)labeledHeatmap(Matri某=modTraitCor,某Label=colname(trainDt),yLabel=name(MEFemale),ySymbol=colname(modlue),colorLabel=FALSE,color=greenWhiteRed (50),te某tMatri某=te某tMatri某,etStdMargin=FALSE,ce某.te某t=0.5,zlim=c(-1,1),main=pate('Module-traitrelationhip'))最终找到几个共表达网络与三阴性表型最相关的模块。
基于肿瘤相关成纤维细胞基因构建乳腺癌预后预测模型及免疫浸润分析
生物技术进展 2024 年 第 14 卷 第 2 期 312 ~ 322Current Biotechnology ISSN 2095‑2341研究论文Articles基于肿瘤相关成纤维细胞基因构建乳腺癌预后预测模型及免疫浸润分析孙莉莉，安外尔·约麦尔阿卜拉，刘富中，布尔兰·叶尔肯别克，迪丽娜尔·叶尔夏提，郭文佳*新疆医科大学附属肿瘤医院，乌鲁木齐 830011摘 要：乳腺癌的转移和恶性进展与肿瘤微环境密切相关。
肿瘤相关成纤维细胞（cancer associated fibroblasts ，CAFs ）是肿瘤微环境中比较重要的细胞，可影响肿瘤的进展及治疗。
从基因表达综合数据库获得乳腺癌单细胞测序数据，对肿瘤微环境细胞进行分簇，再利用WGCNA 识别CAF 相关的关键基因，用该基因在TCGA -BRCA 数据库中构建风险评分模型，进行生存分析、Cox 回归分析、ROC 曲线、构建列线图预测模型性能；通过GO 和KEGG 分析模型相关通路；利用体细胞突变、免疫浸润分析、干性指数分析以及药物敏感性分析探讨风险评分与临床特征及肿瘤微环境的关系。
研究构建了基于10个CAF 基因的乳腺癌预后预测模型，根据风险评分将患者分为高低风险组并进行验证，其中高风险组患者的预后更差，列线图和ROC 曲线也显示模型具有良好的预测效能，乳腺癌病人免疫浸润水平更低、干性指数更高，且高风险组病人对紫杉醇及拉帕替尼这2种药物的敏感性更高。
结果表明，10个CAF 相关基因的风险评分可独立预测乳腺癌的预后及治疗效果，为明确CAF 相关基因在乳腺癌中的作用机制提供了思路，也为乳腺癌易感基因患者的临床个体化治疗提供了理论依据。
关键词：乳腺癌；肿瘤相关成纤维细胞；肿瘤突变负荷；预后模型；免疫浸润DOI ：10.19586/j.20952341.2023.0161中图分类号：Q75， R737.9 文献标志码：AConstruction of Prognostic Prediction Model of Breast Cancer Based on Tumor -associated Fibroblast Genes and Analysis of Immune InfiltrationSUN Lili ， ANWAIER Yuemaierabola ， LIU Fuzhong ， BUERLAN Yeerkenbieke ， DILINAER Ye ，GUO Wenjia *Affiliated Cancer Hospital of Xinjiang Medical University ， Urumqi 830011， ChinaAbstract ：Metastasis and malignant progression of breast cancer are deeply related to the tumor microenvironment. Tumor -associ‐ated fibroblasts （CAFs ） are comparatively important cells in the tumor microenvironment which have implications on tumor pro‐gression and treatment. We obtained single -cell sequencing data of breast cancer downloaded from gene expression omnibus data‐base ， clustered the cells of tumor microenvironment ， and then used WGCNA to identify the key genes related to CAF ， and con‐structed a risk score model with the genes in TCGA -BRCA database ， and performed survival analysis ， Cox regression analysis ， ROC curves ， and constructed a column line graph to predict the performance of the model. Model -related pathways were analyzed by GO and KEGG. The relationship between risk score and clinical features and tumor microenvironment was explored by somaticmutation ， immune infiltration analysis ， stemness index analysis ， and drug sensitivity analysis. A prognostic prediction modelbased on 10 CAF genes was constructed and validated in accordance with the risk scores. Patients were classified into high - and low -risk groups according to the risk scores ， and the prognosis of patients in the high -risk group was worse ， and the column plot and ROC curve also showed that the model had a good predictive efficiency ， and the immune infiltration level of patients with收稿日期：2023‐12‐13； 接受日期：2024‐02‐27基金项目：新疆维吾尔自治区自然科学基金杰出青年科学基金项目（2022D01E27）；新疆维吾尔自治区天池英才项目（2022TCYCGWJ ）。
TCGA数据库生物信息
1.从TCGA下载相应的癌症数据，包括正常样品和癌症样品。
2.差异的lncRNA和microRNA分析。
分别对mRNA和miRNA做差异表达分析，得到差异的miRNA和基因，并绘制热图和火山图。
表差异
图1差异基因火山图
图2热图
3.共表达网络
基因和miRNA的共表达网络。
图3miRNA和mRNA共表达网络4.蛋白互作网络
对差异基因进行蛋白互作网络分析。
图4差异基因PPI网络
5.生存分析
分析基因高低表达与生存时间之间是否具有显著相关性，并且绘制生存曲线。
当然，也可以分析临床信息与生存的关系，比如临床分期与生存时间的关系，癌症大小与生存时间的关系，用药与生存时间的关系，等等。
图5目标基因生存分析
6.基因表达与临床的关系
分析基因与临床数据的关系，如基因的表达和癌症转移的关系，基因表达和临床分期的关系，基因表达和其它临床信息的关系。
图6MARCH1表达与肺转移的关系
7.其它个性化分析
根据客户提供分析案例或者文献，做相应的生物信息分析。
有疑问请联系作者邮箱：602316645@。
最新思路——巧用WGCNA分析GEO和TCGA数据，文章轻松上5分
最新思路——巧用WGCNA分析GEO和TCGA数据，文章轻松上5分癌症相关成纤维细胞(CAF)是胃癌基质中最重要的细胞成分，它促进胃癌进展、治疗抵抗和免疫抑制。
加权基因共表达网络分析(Weighted gene co-expression network analysis, WGCNA)作为一种系统的生物信息学算法，能够将高度协调表达的基因整合到几个基因模块中，并研究该模块与目标表型的关系。
但是到目前为止，CAF 和间质浸润尚未在胃癌中进行WGCNA分析。
本次介绍的文章首次将WGCNA同时用于两个转录组数据集（来自GEO和TCGA数据库），用以分析CAF与胃癌的关系。
这篇文章是2021年10月发表在frontiers in Molecular Biosciences上，文章题目是：Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblast Signature for Predicting Prognosis and Therapeutic Responses in Gastric Cancer。
下面我们来具体看看文章的主要内容吧。
材料与方法01文章的方法我们通过以下流程图来看看吧。
作者首先分别从GEO 和TCGA数据库获得431和330个胃癌样本信息。
随后对这两个数据集的样本进行CAF浸润以及基质评分的计算。
之后通过WGCNA构建针对CAF浸润和基质评分的共表达网络。
分别筛选出相关性最高的黑色和棕色module为关键模块，之后通过阈值的设定进一步筛选得到这两个模块中（取交集）的37个重点基因作为hub基因。
为了进一步缩小基因数目，采用单因素Cox和LASSO回归分析将基因数目缩小到4个，并利用这四个基因构建风险模型。
最后针对这个风险模型分别在两个数据集中作了预后分析、化学和免疫治疗预测性能测试、肿瘤突变负荷（TMB）分析、GSEA分析、与CFA浸润以及相关marker 基因相关性分析、以及在分子和蛋白水平对这四个关键基因的验证。
乳腺癌相关的lncRNA-mRNA共表达扰动网络构建
第４２卷㊀第３期２０２３年㊀６月北京生物医学工程ＢｅｉｊｉｎｇＢｉｏｍｅｄｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＶｏｌ ４２㊀Ｎｏ ３Ｊｕｎｅ㊀２０２３基金项目：第６５批中国博士后科学基金项目（２０１９Ｍ６５１６５８）㊁南京农业大学大学生创新训练项目（２０２０２３ＸＸ０３）资助作者单位：南京农业大学理学院（南京㊀２１００９５）通信作者：陈园园㊂Ｅ⁃ｍａｉｌ：ｃｈｅｎｙｕａｎｙｕａｎ＠ｎｊａｕ．ｅｄｕ．ｃｎ乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建黄彦祚㊀李海龙㊀卢乐亭㊀陈园园摘㊀要㊀目的基于复杂生物网络和机器学习方法，识别乳腺癌相关的边缘生物标志物，构建乳腺癌生存预后模型，从而在系统水平解释乳腺癌的发生发展机制㊂方法首先基于ＴＣＧＡ数据库的ＲＮＡ⁃ｓｅｑ数据识别乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动关系对，进一步构建乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络并对网络中的关键基因进行通路富集分析㊂然后，基于乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ关系对，构建乳腺癌预测的分类器模型㊂最后，通过Ｌａｓｓｏ回归筛选变量构建多因素Ｃｏｘ比例风险回归模型对乳腺癌患者进行生存预后分析㊂结果构建了乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络，其中的关键基因富集分析得到３２条与乳腺癌相关的生物通路㊂分类预测模型的灵敏度㊁特异度和准确性分别为９８ ２％㊁８５ ２％㊁９７ ６％㊂Ｌａｓｓｏ回归共筛选出２２个和乳腺癌生存预后显著相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，进而构建的生存预测模型把训练集和测试集的乳腺癌患者分为高风险组和低风险组，两组患者生存预后均存在明显差异㊂结论ＬｎｃＲＮＡ⁃ｍＲＮＡ共表达互作网络中的关键基因以及乳腺癌相关的边缘生物标志物大多被证明与乳腺癌相关㊂同时基于边缘生物标志物的预后模型可以稳健地预测乳腺癌患者的生存预后状态，有利于从网络层面更好地理解乳腺癌的发生发展机制㊂关键词㊀ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络；边缘生物标志物；乳腺癌预测模型；Ｃｏｘ比例风险回归模型ＤＯＩ：１０ ３９６９／ｊ．ｉｓｓｎ．１００２－３２０８ ２０２３ ０３ ００４．中图分类号㊀Ｒ３１８ ０４；Ｑ３５４㊀㊀文献标志码㊀Ａ㊀㊀文章编号㊀１００２－３２０８（２０２３）０３－０２４０－０８本文著录格式㊀黄彦祚，李海龙，卢乐亭，等．乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建［Ｊ］．北京生物医学工程，２０２３，４２（３）：２４０－２４７．ＨＵＡＮＧＹａｎｚｕｏ，ＬＩＨａｉｌｏｎｇ，ＬＵＬｅｔｉｎｇ，ｅｔａｌ．ＣｏｎｓｔｒｕｃｔｉｏｎｏｆａｎｅｔｗｏｒｋｏｆｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｒｅｌａｔｅｄｔｏｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＢｅｉｊｉｎｇＢｉｏｍｅｄｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，２０２３，４２（３）：２４０－２４７．ＣｏｎｓｔｒｕｃｔｉｏｎｏｆａｎｅｔｗｏｒｋｏｆｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｒｅｌａｔｅｄｔｏｂｒｅａｓｔｃａｎｃｅｒＨＵＡＮＧＹａｎｚｕｏ，ＬＩＨａｉｌｏｎｇ，ＬＵＬｅｔｉｎｇ，ＣＨＥＮＹｕａｎｙｕａｎＣｏｌｌｅｇｅｏｆＳｃｉｅｎｃｅ，ＮａｎｊｉｎｇＡｇｒｉｃｕｌｔｕｒａｌＵｎｉｖｅｒｓｉｔｙ，Ｎａｎｊｉｎｇ㊀２１００９５Ｃｏｒｒｅｓｐｏｎｄｉｎｇａｕｔｈｏｒ：ＣＨＥＮＹｕａｎｙｕａｎ（Ｅ⁃ｍａｉｌ：ｃｈｅｎｙｕａｎｙｕａｎ＠ｎｊａｕ．ｅｄｕ．ｃｎ）ʌＡｂｓｔｒａｃｔɔ㊀ＯｂｊｅｃｔｉｖｅＢａｓｅｄｏｎｔｈｅｃｏｍｐｌｅｘｂｉｏｌｏｇｉｃａｌｎｅｔｗｏｒｋｓａｎｄｍａｃｈｉｎｅｌｅａｒｎｉｎｇｍｅｔｈｏｄｓ，ｔｈｅｅｄｇｅｂｉｏｍａｒｋｅｒｓａｓｓｏｃｉａｔｅｄｗｉｔｈｂｒｅａｓｔｃａｎｃｅｒｗｅｒｅｉｄｅｎｔｉｆｉｅｄａｎｄｔｈｅｓｕｒｖｉｖａｌｐｒｏｇｎｏｓｉｓｍｏｄｅｌｏｆｂｒｅａｓｔｃａｎｃｅｒｗａｓｃｏｎｓｔｒｕｃｔｅｄｔｏｆｕｒｔｈｅｒｅｘｐｌａｉｎｔｈｅｏｃｃｕｒｒｅｎｃｅａｎｄｄｅｖｅｌｏｐｍｅｎｔｏｆｂｒｅａｓｔｃａｎｃｅｒａｔａｓｙｓｔｅｍａｔｉｃｌｅｖｅｌ．ＭｅｔｈｏｄｓＦｉｒｓｔｌｙ，ｂａｓｅｄｏｎｔｈｅＲＮＡ⁃ｓｅｑｄａｔａｏｆＴＣＧＡ，ｗｅｉｄｅｎｔｉｆｉｅｄｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｐａｉｒｓａｎｄｃｏｎｓｔｒｕｃｔｅｄｔｈｅｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｎｅｔｗｏｒｋｓｒｅｌａｔｅｄｔｏｂｒｅａｓｔｃａｎｃｅｒ．Ｆｕｒｔｈｅｒ，ｗｅｃｏｎｄｕｃｔｅｄｔｈｅｐａｔｈｗａｙｅｎｒｉｃｈｍｅｎｔａｎａｌｙｓｉｓｏｆｋｅｙｇｅｎｅｓｉｎｔｈｅｎｅｔｗｏｒｋ．Ｔｈｅｎ，ａｃｌａｓｓｉｆｉｅｒｍｏｄｅｌｆｏｒｂｒｅａｓｔｃａｎｃｅｒｐｒｅｄｉｃｔｉｏｎｗａｓｃｏｎｓｔｒｕｃｔｅｄｂａｓｅｄｏｎｔｈｅｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｐａｉｒｓ．Ｆｉｎａｌｌｙ，ａｍｕｌｔｉｖａｒｉａｔｅＣｏｘｐｒｏｐｏｒｔｉｏｎａｌｒｉｓｋｒｅｇｒｅｓｓｉｏｎｍｏｄｅｌｗａｓｅｓｔａｂｌｉｓｈｅｄｂｙｓｃｒｅｅｎｉｎｇｖａｒｉａｂｌｅｓｕｓｉｎｇＬａｓｓｏｒｅｇｒｅｓｓｉｏｎｔｏａｎａｌｙｚｅｔｈｅｓｕｒｖｉｖａｌｐｒｏｇｎｏｓｉｓｏｆｂｒｅａｓｔｃａｎｃｅｒｐａｔｉｅｎｔｓ．ＲｅｓｕｌｔｓＷｅｃｏｎｓｔｒｕｃｔｅｄａｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｎｅｔｗｏｒｋａｓｓｏｃｉａｔｅｄｗｉｔｈｂｒｅａｓｔｃａｎｃｅｒ．Ａｎｄｔｈｅｋｅｙｇｅｎｅｓｉｎｔｈｅｎｅｔｗｏｒｋｗｅｒｅｕｓｅｄｔｏｐｅｒｆｏｒｍｐａｔｈｗａｙｅｎｒｉｃｈｍｅｎｔａｎａｌｙｓｉｓ．Ａｔｏｔａｌｏｆ３２ｂｉｏｌｏｇｉｃａｌｐａｔｈｗａｙｓａｓｓｏｃｉａｔｅｄｗｉｔｈｂｒｅａｓｔｃａｎｃｅｒｗｅｒｅｏｂｔａｉｎｅｄ．Ｔｈｅｓｅｎｓｉｔｉｖｉｔｙ，ｓｐｅｃｉｆｉｃｉｔｙ，ａｃｃｕｒａｃｙｏｆｔｈｅｃｌａｓｓｉｆｉｅｒｍｏｄｅｌｆｏｒｂｒｅａｓｔｃａｎｃｅｒｐｒｅｄｉｃｔｉｏｎｗｅｒｅ９８．２％，８５．２％ａｎｄ９７．６％ｒｅｓｐｅｃｔｉｖｅｌｙ．Ａｔｏｔａｌｏｆ２２ｌｎｃＲＮＡ⁃ｍＲＮＡｉｎｔｅｒａｃｔｉｏｎｐａｉｒｓ，ｗｈｉｃｈｗｅｒｅｓｉｇｎｉｆｉｃａｎｔｌｙａｓｓｏｃｉａｔｅｄｗｉｔｈｂｒｅａｓｔｃａｎｃｅｒｓｕｒｖｉｖａｌｐｒｏｇｎｏｓｉｓ，ｗｅｒｅｉｄｅｎｔｉｆｉｅｄｂｙＬａｓｓｏｒｅｇｒｅｓｓｉｏｎ．Ｂａｓｅｄｏｎｔｈｅｓｕｒｖｉｖａｌｐｒｅｄｉｃｔｉｏｎｍｏｄｅｌ，ｂｒｅａｓｔｃａｎｃｅｒｐａｔｉｅｎｔｓｉｎｔｈｅｔｒａｉｎｉｎｇｓｅｔａｎｄｔｈｅｔｅｓｔｓｅｔｗｅｒｅｄｉｖｉｄｅｄｉｎｔｏｈｉｇｈ⁃ｒｉｓｋｇｒｏｕｐａｎｄｌｏｗ⁃ｒｉｓｋｇｒｏｕｐ，ａｎｄｔｈｅｓｕｒｖｉｖａｌｐｒｏｇｎｏｓｉｓｏｆｐａｔｉｅｎｔｓｉｎｔｈｅｔｗｏｇｒｏｕｐｓｗａｓｓｉｇｎｉｆｉｃａｎｔｌｙｄｉｆｆｅｒｅｎｔ．ＣｏｎｃｌｕｓｉｏｎｓＫｅｙｇｅｎｅｓｉｎｔｈｅｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｎｅｔｗｏｒｋａｎｄｂｒｅａｓｔｃａｎｃｅｒ⁃ｒｅｌａｔｅｄｅｄｇｅｂｉｏｍａｒｋｅｒｓｈａｖｅｍｏｓｔｌｙｂｅｅｎｐｒｏｖｅｄｔｏｂｅａｓｓｏｃｉａｔｅｄｗｉｔｈｂｒｅａｓｔｃａｎｃｅｒ．Ｍｅａｎｗｈｉｌｅ，ｔｈｅｐｒｏｇｎｏｓｔｉｃｍｏｄｅｌｂａｓｅｄｏｎｅｄｇｅｂｉｏｍａｒｋｅｒｓｃａｎｐｒｅｄｉｃｔｔｈｅｓｕｒｖｉｖａｌａｎｄｐｒｏｇｎｏｓｉｓｏｆｂｒｅａｓｔｃａｎｃｅｒｐａｔｉｅｎｔｓｒｏｂｕｓｔｌｙ．Ｔｈｉｓｐａｐｅｒｉｓｈｅｌｐｆｕｌｔｏｂｅｔｔｅｒｕｎｄｅｒｓｔａｎｄｔｈｅｍｅｃｈａｎｉｓｍｏｆｔｈｅｏｃｃｕｒｒｅｎｃｅａｎｄｄｅｖｅｌｏｐｍｅｎｔｏｆｂｒｅａｓｔｃａｎｃｅｒｉｎｎｅｔｗｏｒｋｌｅｖｅｌ．ʌＫｅｙｗｏｒｄｓɔ㊀ＬｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｎｅｔｗｏｒｋ；ｅｄｇｅｂｉｏｍａｒｋｅｒ；ｂｒｅａｓｔｃａｎｃｅｒｐｒｅｄｉｃｔｉｏｎｍｏｄｅｌ；Ｃｏｘｐｒｏｐｏｒｔｉｏｎａｌｒｉｓｋｒｅｇｒｅｓｓｉｏｎｍｏｄｅｌ０㊀引言近年来，越来越多的研究表明，长链非编码ＲＮＡ（ｌｏｎｇｎｏｎ⁃ｃｏｄｉｎｇＲＮＡ，ｌｎｃＲＮＡ）在包括癌症在内的许多疾病的发生发展中发挥着重要作用，已受到越来越多的关注［１］㊂随着高通量测序技术的发展，大量的ｌｎｃＲＮＡ被发现，其作用机制的研究也取得了一定进展［２］㊂目前已在乳腺癌细胞及组织中发现多种异常表达ｌｎｃＲＮＡ，它们可能在乳腺癌细胞增殖㊁凋亡㊁侵袭㊁转移及药物敏感性等方面起了重要作用［３］㊂在癌症发展过程中，ｌｎｃＲＮＡ参与了多种表观遗传复合物的调节过程，从而抑制或激活癌症相关ｍＲＮＡ基因的表达［４］㊂因此，探讨ｌｎｃＲＮＡ⁃ｍＲＮＡ互作对乳腺癌发生发展的影响至关重要㊂复杂疾病（尤其癌症）的发生并不是单个基因失调导致，往往是由多个分子及其相互作用失调引起的㊂竞争性内源ＲＮＡ机制是探索ｌｎｃＲＮＡ如何参与恶性肿瘤调控的重要方法之一［５］㊂近期研究报道通过高通量测序和加权基因共表达网络法（ｗｅｉｇｈｔｅｄｇｅｎｅｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｎｅｔｗｏｒｋａｎａｌｙｓｉｓ，ＷＧＣＮＡ）可以进行表达谱基因系统分析［６］㊂Ｗｅｉ等［７］通过微阵列分析的方法识别了ｌｎｃＲＮＡ和ｍＲＮＡ的差异表达模块㊂Ｙｉｎ等［８］通过全基因组关联分析的方法对ｌｎｃＲＮＡ⁃ｍＲＮＡ调控网络进行综合分析㊂但国内尚无关于乳腺癌的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系的研究，而复杂生物网络可以从系统层面解释癌症的发生发展机制㊂为此本文拟基于ＲＮＡ⁃ｓｅｑ数据，计算个体特异的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动值，进一步识别乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系，并构建共表达扰动网络㊂同时，基于ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系构建区分乳腺癌样本和正常样本的分类器模型㊂另外，基于单变量Ｃｏｘ回归识别与乳腺癌生存预后相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ关系对，通过Ｌａｓｓｏ回归筛选变量并进一步构建多因素Ｃｏｘ回归模型，对乳腺癌患者进行生存预后分析㊂本研究拟通过基因生物网络和机器学习方法，识别乳腺癌相关的边缘生物标志物并构建乳腺癌生存预后模型，在系统水平上解释乳腺癌的发生发展，为制定合适的治疗计划㊁协助评价治疗结果㊁预测患者的生存时间等提供重要依据，进而提高乳腺癌预后诊断以及促进精准医疗的发展㊂１㊀材料与方法１ １㊀数据来源从ＴＣＧＡ数据库（ｈｔｔｐｓ：／／ｃａｎｃｅｒｇｅｎｏｍｅ ｎｉｈ ｇｏｖ／）下载１０９７个乳腺癌样本和１１３个癌旁组织样本的ＲＮＡ⁃ｓｅｑ数据，以及乳腺癌患者的临床数据；从ＧＥＮＣＯＤＥ数据库（ｈｔｔｐｓ：／／ｗｗｗ ｇｅｎｃｏｄｅｇｅｎｅｓ ｏｒｇ）下载Ｖ３３版本的基因注释文件；从数据库ＧＳＥＡ／ＭＳｉｇＤＢ（ｈｔｔｐ：／／ｓｏｆｔｗａｒｅ ｂｒｏａｄｉｎｓｔｉｔｕｔｅ ｏｒｇ／ｇｓｅａ／ｍｓｉｇｄｂ）下载１８６个ＫＥＧＧ通路基因集㊂基因表达数据是经过标准化后的Ｌｅｖｅｌ３的ＲＰＫＭ格式，共包含基因５６５２１个㊂１ ２㊀数据预处理首先对标准化后的ＲＮＡ⁃ｓｅｑ数据进行预处理，㊃１４２㊃第３期㊀㊀㊀㊀㊀㊀黄彦祚，等：乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建删除大于９０％样本表达值为０的基因㊂根据Ｇｅｎｃｏｄｅ基因注释文件共得到含有１０９４１个ｌｎｃＲＮＡ基因和１８６０１个ｍＲＮＡ基因的表达数据㊂１ ３㊀构建乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络㊀㊀识别乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络框架见图１㊂对于每个ｌｎｃＲＮＡ⁃ｍＲＮＡ关系对，基于正常样本的表达数据计算皮尔逊相关系数（Ｐｅａｒｓｏｎｃｏｒｒｅｌａｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔ，ＰＣＣ）㊂挑选出所有显著线性相关的ｌｎｃＲＮＡ与ｍＲＮＡ关系对（ＰＣＣ＞０ ７５且Ｐ＜０ ０５），即具有共表达关系的ｌｎｃＲＮＡ⁃ｍＲＮＡ基因互作对㊂对于每一个共表达的ｌｎｃＲＮＡ⁃ｍＲＮＡ基因互作对，运用最小二乘法的简单线性回归模型拟合ｌｎｃＲＮＡ和ｍＲＮＡ共表达的直线方程，即：ｙ＝ｂｘ＋ａ（１）㊀㊀式中：ｂ为拟合直线方程的斜率；ａ为拟合直线方程的截距项㊂设（ｘｋ，ｙｋ）代表第ｋ个样本ｍＲＮＡｌｎｃＲＮＡ的表达量，则可以计算（ｘｋ，ｙｋ）到直线ｙ＝ｂｘ＋ａ的距离：ｄｋ＝ｂｘｋ＋ａ－ｙｋ１＋ｂ２（２）则所有癌症样本的距离之和Ｄ＝ １０９７ｋ＝１ｄｋ可以用来衡量乳腺癌样本在此ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系中的共表达扰动情况㊂若Ｄ越大，则代表乳腺癌样本中ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动越明显㊂通过链蒙特卡洛（ＭｏｎｔｅＣａｒｌｏ，ＭＣ）随机抽样方法得到互作扰动距离的经验零分布（ｅｍｐｉｒｉｃａｌｎｕｌｌｄｉｓｔｒｉｂｕｔｉｏｎ）㊂所有Ｐ＜０ ０１的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达关系对组成了乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，这些互作关系对构成了乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络㊂在整个乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络中，少量度（ｄｅｇｒｅｅ）很高的基因作为网络枢纽参与重要的生命活动，这些基因被认为是影响乳腺癌活动的关键基因㊂根据Ｂａｒｂｅｒａｎ的拓扑网络研究［９］，选取度大于２０的ｍＲＮＡ基因进行通路富集分析㊂通过超几何分布检验，计算Ｐ值：Ｐ（ｍ，Ｍ，Ｎ，ｎ）＝１－ｍ－１ｉ＝０Ｍｉæèçöø÷Ｎ－Ｍｎ－ｉæèçöø÷Ｎｎæèçöø÷（３）㊀㊀式中：Ｎ为整个乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作网络中的ｍＲＮＡ基因总数；Ｍ为互作网络中ｄｅｇｒｅｅ＞２０的关键基因个数；ｎ为某条基因通路中的基因个数；ｍ为ｌｎｃＲＮＡ⁃ｍＲＮＡ互作网络中关键基因落入该通路中基因的个数㊂最后选取Ｐ＜００５的通路作为显著富集的通路㊂图１㊀识别乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络框架图Ｆｉｇｕｒｅ１㊀ＡｎｉｎｔｅｇｒａｔｉｖｅｆｒａｍｅｗｏｒｋｉｄｅｎｔｉｆｙｉｎｇｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｐｅｒｔｕｒｂａｔｉｏｎｎｅｔｗｏｒｋｓ㊃２４２㊃北京生物医学工程㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第４２卷１ ４㊀基于ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系构建乳腺癌预测模型㊀㊀７６８个乳腺癌样本和７９个癌旁组织样本的ＲＮＡ⁃ｓｅｑ数据为训练集，３２９个乳腺癌样本和３３个癌旁组织样本数据为独立测试集㊂以乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系的共表达扰动距离为特征，构建预测乳腺癌的随机森林分类器模型，评价模型的预测性能㊂然后根据特征贡献对ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对进行降序排列，选择前１０％较为重要的关系作为潜在ｌｎｃＲＮＡ⁃ｍＲＮＡ边缘生物标志物㊂本文在Ｒ３ ５ ３上进行实验，使用Ｒ包 ｒａｎｄｏｍＦｏｒｅｓｔ 建立预测模型㊂根据训练集对随机森林模型进行训练，得到反映其预测效果的混淆矩阵㊂训练后的随机森林模型对独立测试集进行测试，结合袋外错误率㊁灵敏度㊁特异度来衡量预测的准确率㊂袋外错误率（ｏｕｔ⁃ｏｆ⁃ｂａｇｅｒｒｏｒ，ＯＯＢ）定义如下：ＯＯＢ＝被分类错误数总数（４）㊀㊀模型的灵敏度（ｔｒｕｅｐｏｓｉｔｉｖｅｒａｔｅ，ＴＰＲ），指真实类别为正类的样本中分类预测也为正的比例，即：ＴＰＲ＝ＴＰＴＰ＋ＦＮ（５）㊀㊀式中：ＴＰ为真实类别为正㊁分类预测也为正的数目；ＦＮ为真实类别为正㊁分类预测为负的数目㊂特异度（ｔｒｕｅｎｅｇａｔｉｖｅｒａｔｅ，ＴＮＲ），其定义为真实类别为负类的样本中，分类预测也为负的比例：ＴＮＲ＝ＴＮＴＮ＋ＦＰ（６）㊀㊀式中：ＴＮ为真实类别为负㊁分类预测也为负的数目；ＦＰ为真实类别为负㊁分类预测为正的数目㊂准确性（ａｃｃｕｒａｃｙ，ＡＣＣ），指分类正确的记录个数占总记录个数的比例：ＡＣＣ＝ＴＮ＋ＴＰＴＰ＋ＦＮ＋ＴＮ＋ＦＰ（７）㊀㊀ＯＯＢ㊁ＴＰＲ㊁ＴＮＲ㊁ＡＣＣ作为评估随机森林模型的性能指标㊂１ ５㊀生存分析利用癌症患者的临床数据，以患者生存时间和状态为因变量，乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对的共表达扰动距离作为协变量，建立单因素Ｃｏｘ回归模型，其中的共表达扰动距离ｄｋ为第ｋ个癌症样本到共表达直线方程的距离㊂选取其中Ｐ＜０ ０１的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对㊂由于协变量过多，在建立多因素Ｃｏｘ回归模型时可能导致计算复杂以及过拟合的问题，因此选择最小绝对收缩选择算子（ｌｅａｓｔａｂｓｏｌｕｔｅｓｈｒｉｎｋａｇｅａｎｄｓｅｌｅｃｔｉｏｎｏｐｅｒａｔｏｒ，ＬＡＳＳＯ）方法，利用Ｒ语言 ｓｕｒｖｉｖａｌ ㊁ ｓｕｒｖｍｉｎｅｒ ㊁ ｃａｒｅｔ ㊁ ｇｌｍｎｅｔ 软件包对ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对进一步筛选㊂运用ＬＡＳＳＯ方法压缩回归系数，并选取一个交叉检验均方误差最小的λ值，从而得到一个最优的ＬＡＳＳＯ回归模型㊂其中大部分基因对的系数被缩减到０，剩下相对较少非零系数的基因对则被认为是和乳腺癌预后高度相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对㊂在训练集上构建多因素Ｃｏｘ回归风险评分方程：ｒｉｓｋｓｃｏｒｅ＝ ｎｉ＝１ＣｏｅｆｉˑＸｉ（８）㊀㊀式中：Ｃｏｅｆｉ为第ｉ个ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对回归系数；Ｘｉ为第ｉ个乳腺癌预后高度相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对；ｎ为乳腺癌预后高度相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对的个数㊂根据公式计算每个患者的风险评分，并以风险评分中位数为截断值，将乳腺癌患者分为低风险组和高风险组，进一步画出ＫａｐｌａｎＭｅｉｅｒ生存曲线㊂最后用同样的方法对测试集进行生存验证㊂２㊀结果２ １㊀乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动互作关系㊀㊀ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动（ｌｇ转换）分布情况见图２，其中１１ ０８１为Ｐ＝０ ０１的分位点㊂大于该阈值的互作关系对构成乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作网络㊂图２㊀乳腺癌样本的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动分布Ｆｉｇｕｒｅ２㊀ＰｅｒｔｕｒｂａｔｉｏｎｄｉｓｔｒｉｂｕｔｉｏｎｏｆｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｉｎｂｒｅａｓｔｃａｎｃｅｒｓａｍｐｌｅｓ㊃３４２㊃第３期㊀㊀㊀㊀㊀㊀黄彦祚，等：乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达互作网络中共包含２８６６个ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，６４８个ＬｎｃＲＮＡ和７３３个ｍＲＮＡ，共１３８１个基因㊂利用Ｃｙｔｏｓｃａｐｅ软件将乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达互作网络可视化，如图３所示㊂紫色节点表示ｌｎｃＲＮＡ，蓝色节点表示ｍＲＮＡ图３㊀乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达互作网络Ｆｉｇｕｒｅ３㊀Ｂｒｅａｓｔｃａｎｃｅｒ⁃ａｓｓｏｃｉａｔｅｄｌｎｃＲＮＡ⁃ｍＲＮＡｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｉｎｔｅｒａｃｔｉｏｎｎｅｔｗｏｒｋ网络图中共有关键基因５４个，包括２４个ｌｎｃＲＮＡ和３０个ｍＲＮＡ㊂其中，一些基因已被证明为与乳腺癌相关，例如，ＡＰＯＤ与癌细胞的增殖活性降低有关，大量存在于衰老细胞中［１０］㊂下调ＤＡＮＣＲ可以抑制乳腺癌细胞的致瘤性和发育［１１］㊂ＺＦＡＳ１过表达可通过阻滞细胞周期㊁诱导细胞凋亡等途径抑制乳腺癌细胞增殖［１２］㊂胞外ＣＤＨ１突变降低了细胞的活动性和理论上的转移能力［１３］㊂上调ＧＡＳ５可以抑制ＴＮＢＣ的进展，促进ＴＮＢＣ细胞的化疗敏感性和凋亡［１４］㊂ＩＤＨ２在ｍＲＮＡ或蛋白水平上的高表达与乳腺癌患者预后不良相关［１５］㊂过表达ＡＧＡＰ２－ＡＳ１可促进细胞生长，抑制细胞凋亡［１６］㊂２ ２㊀富集通路通过基因富集分析得到３２条与乳腺癌相关的生物通路，主要包括一些免疫相关的通路：Ｔ细胞受体信号通路，它被发现是ＨＱ－ＢＳ治疗乳腺癌的潜在信号通路［１７］；Ｂ细胞受体信号通路在三阴性乳腺癌的发展中失调［１８］；造血细胞谱系通路，研究发现造血细胞谱系，细胞黏附分子和原发性免疫缺陷明显增加了ＬｕｍｉｎａｌＢ型乳腺癌的ＫＥＧＧ通路［１９］；Ｔｏｌｌ样受体信号通路通过缺氧诱导因子增强乳腺癌细胞的恶性特征［２０］；自然杀伤细胞的细胞毒信号通路㊂同时，一些重要的信号转导和信号分子互作通路也被富集：ＥｒｂＢ信号通路不仅在乳腺癌发生发展中显著失调，而且在新型抗乳腺癌靶向疗法中发挥重要作用［２１］；ＭＡＰＫ信号通路中Ｌｉｎｃ⁃ＲｏＲ促进ＭＡＰＫ／ＥＲＫ信号传导并赋予乳腺癌非雌激素依赖性生长［２２］；Ｗｎｔ信号通路中ＬｎｃＣＣＡＴ１通过激活ＷＮＴ／β－ｃａｔｅｎｉｎ信号传导促进乳腺癌干细胞功能［２３］；ｐ５３信号通路可能是乳腺癌进展的重要途径，与ｐ５３相关的基因ＣＣＮＥ２㊁ＣＣＮＢ１和ＲＲＭ２可能会作为治疗ＢＣ的候选治疗基因靶标［２４］；ＪＡＫ⁃ＳＴＡＴ信号通路，抑制ＪＡＫ／ＳＴＡＴ通路和Ａｋｔ信号通路可以抑制乳腺癌细胞迁移［２５］；Ｈｅｄｇｅｈｏｇ信号通路，ＬＫＢ１通过抑制Ｈｅｄｇｅｈｏｇ信号通路部分抑制乳腺癌的发生［２６］㊂另外，癌症通路，这个与乳腺癌相关的重要通路也被富集㊂富集的结果如图４所示㊂２ ３㊀分类器模型预测结果模型对测试集的预测精度见表１㊂ＯＯＢ为２ ７５％，ＴＰＲ和ＴＮＲ分别为９８ ２％和８５ ２％，ＡＣＣ为９７ ６％㊂在训练集和测试集的随机森林模型性能分别为：特异度９９ ５％和９９ ４％；灵敏度８４ ８％和８２ ４％；准确性９８ １％和９７ ８％㊂表１㊀随机森林模型在测试集上的性能Ｔａｂｌｅ１㊀Ｐｅｒｆｏｒｍａｎｃｅｏｆｒａｎｄｏｍｆｏｒｅｓｔｔｒａｉｎｉｎｇｏｎｔｅｓｔｓｅｔ真实值＼预测值癌症样本正常样本分类误差癌症样本３３０４１ ２％正常样本６２３２０ ７％袋外错误率２ ７５％基于ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系构建了乳腺癌预测模型，识别出ｌｎｃＲＮＡ⁃ｍＲＮＡ边缘生物标志物㊂其中有的基因已被证明为与乳腺癌相关：ＥＲα和ＡＰＯＤ共表达老年乳腺癌患者中具有预后意义［２７］；ＶＩＭ表达与肿瘤分化程度有关，随着肿瘤的失分化，ＶＩＭ的表达能力逐渐增强［２８］；ＳＰＩＮＴ２基因通过表观沉默或下调改变ＨＧＦ激活／抑制比率的平衡，从而促进癌症的发展［２９］；ＳＮＲＰＢ基因在乳腺癌中高表达，可以有效预测乳腺癌转移的发生，为乳腺癌的靶向治疗提供依据［３０］；ＳＰＩＮＴ１－ＡＳ１通过调节ｍｉＲ－ｌｅｔ－７ａ／ｂ／ｉ－５ｐ可以促进乳腺癌细胞的增殖和迁移，因此它可能是乳腺癌的重要调控因子［３１］㊂㊃４４２㊃北京生物医学工程㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第４２卷图４㊀基因通路富集结果Ｆｉｇｕｒｅ４㊀Ｇｅｎｅｐａｔｈｗａｙｅｎｒｉｃｈｍｅｎｔｒｅｓｕｌｔｓ２ ４㊀乳腺癌预后相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达互作关系对㊀㊀最终识别出２２个乳腺癌生存预后显著相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，并基于这些互作关系对构建多因素Ｃｏｘ回归模型：Ｓｃｏｒｅ＝－２８ ７９６ˑ（ＰＣＥＤ１Ｂ－ＡＳ１，ＨＬＡ－ＤＲＡ）＋３０ ０３０ˑ（ＰＳＭＢ８－ＡＳ１，ＨＬＡ－ＤＲＡ）－１５ ８５８ˑ（ＨＯＸＢ－ＡＳ１，ＡＢＣＣ９）－１５ ５８６ˑ（ＨＯＸＢ－ＡＳ１，ＧＲＫ３）＋１８７ ３４３ˑ（ＧＡＳ１ＲＲ，ＬＰＣＡＴ２）＋ （９）基于训练集和测试集的两组乳腺癌患者的Ｋａｐｌａｎ⁃Ｍｅｉｅｒ生存曲线见图５，容易发现不管是训练集还是测试集的两组乳腺癌患者的预后都存在显著差异，且高风险组的生存预后显著劣于低风险组㊂这个结果表明基于基因对的预后模型可以很好地预测乳腺癌患者的生存时间㊂３㊀讨论与结论目前癌症病理研究一般仅限于群体的基因表达量和突变信息，而乳腺癌个性化医疗处于发展相对缓慢的阶段㊂但是乳腺癌的发生并不是单个基因导致的，而是由多个分子及其相互作用失调引起的㊂基因并不是独立存在的，而是存在于许多复杂的分子网络中㊂因此，构建乳腺癌相关网络从系统层面解释乳腺癌的发生发展机制㊂但由于条件所限，缺乏相应的临床试验来验证本文结果㊂在本研究中，基于正常样本基因表达数据，通过癌症样本点偏离拟合直线距离计算每个样本的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动情况，识别乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，并构建乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络㊂同时对网络中的关键基因进行基因通路富集分析，结果显示这些基因与癌症通路㊁细胞凋亡通路以及ＥｒｂＢ信号通路等密切相关㊂通过Ｃｏｘ比例风险模型识别出与乳腺癌预后相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作关系对，并建立生存风险预测模型，结果表明本研究中基因对的预后模型可以很好地预测乳腺癌患者的生存时间㊂基于识别出的ｌｎｃＲＮＡ⁃ｍＲＮＡ边缘生物标志物，发现边缘生物标志物中含有乳腺癌相关的诸多㊃５４２㊃第３期㊀㊀㊀㊀㊀㊀黄彦祚，等：乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建图５㊀基于训练集和测试集的两组乳腺癌患者的Ｋａｐｌａｎ⁃Ｍｅｉｅｒ生存曲线Ｆｉｇｕｒｅ５㊀Ｋａｐｌａｎ⁃Ｍｅｉｅｒｓｕｒｖｉｖａｌｃｕｒｖｅｓｏｆｔｗｏｇｒｏｕｐｓｏｆｂｒｅａｓｔｃａｎｃｅｒｐａｔｉｅｎｔｓｂａｓｅｄｏｎｔｒａｉｎｉｎｇｓｅｔａｎｄｔｅｓｔｓｅｔ重要基因㊂已被证实的有：基因ＤＡＸＸ可抑制内分泌治疗后雌激素受体阳性乳腺癌患者的肿瘤起始细胞的数量［３２］；通过消除基因ＥＥＦ１Ａ１的转录，癌细胞可以在蛋白毒性应激后迅速诱导热休克反应并存活［３３］；ＭｅＣＰ２基因通过抑制ＲＰＬ５／ＲＰＬ１１转录来促进泛素介导的ｐ５３降解［３４］，从而促进乳腺癌的细胞增殖并抑制凋亡㊂值得注意的是，网络中很多关键基因与乳腺癌相关：乳腺癌患者ＳＰＡＲＣｒｓ７７１９５２１基因与ＮＰＩ显著相关，ＮＰＩ是一种可靠的乳腺癌预后指标；Ａｎｘａ５可作为乳腺癌肿瘤发生㊁转移和侵袭预测的生物标志物；ＮＯＲＡＤ可以抑制乳腺癌细胞的迁移㊁侵袭和转移㊂结果中ＥＲα－ＡＰＯＤ和ＳＰＩＮＴ１－ＡＳ１的共表达很高，是乳腺癌进展的重要调控因子，可以作为乳腺癌预后的指标及潜在的治疗靶点㊂抑制基因ＣＣＴ３的表达可抑制乳腺癌细胞的增殖和迁移㊂ＨＳＰ９０ＡＡ１㊁ＳＲＣ㊁ＨＳＰＡ８㊁ＥＳＲ１㊁ＡＣＴＢ㊁ＰＰＰ２ＣＡ㊁ＲＰＬ等度较高的基因在乳腺癌中高表达，且与乳腺癌的不良预后密切相关，可以作为基因检测的理论依据㊂故通过研究可以较准确地找出乳腺癌的决定基因和预测乳腺癌患者的生存时间，对于其他复杂疾病的病因研究以及个性化医疗具有重要的参考价值㊂基于共表达扰动的乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ互作网络的构建能够从系统层面更好地理解乳腺癌的发生和发展机制，对提高乳腺癌预后诊断以及促进精准医疗的发展等都有重要的意义㊂参考文献［１］㊀ＪｉｎＸ，ＧｅＬＰ，ＬｉＤＱ，ｅｔａｌ．ＬｎｃＲＮＡＴＲＯＪＡＮｐｒｏｍｏｔｅｓｐｒｏｌｉｆｅｒａｔｉｏｎａｎｄｒｅｓｉｓｔａｎｃｅｔｏＣＤＫ４／６ｉｎｈｉｂｉｔｏｒｖｉａＣＤＫ２ｔｒａｎｓｃｒｉｐｔｉｏｎａｌａｃｔｉｖａｔｉｏｎｉｎＥＲ＋ｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＭｏｌｅｃｕｌａｒＣａｎｃｅｒ，２０２０，１９（１）：８７．［２］㊀ＬｉａｎｇＹＲ，ＳｏｎｇＸＪ，ＬｉＹＭ，ｅｔａｌ．ＬｎｃＲＮＡＢＣＲＴ１ｐｒｏｍｏｔｅｓｂｒｅａｓｔｃａｎｃｅｒｐｒｏｇｒｅｓｓｉｏｎｂｙｔａｒｇｅｔｉｎｇｍｉＲ－１３０３／ＰＴＢＰ３ａｘｉｓ［Ｊ］．ＭｏｌｅｃｕｌａｒＣａｎｃｅｒ，２０２０，１９（１）：８５．［３］㊀ＬｕｏＬＹ，ＺｈａｎｇＪＬ，ＴａｎｇＨＬ，ｅｔａｌ．ＬｎｃＲＮＡＳＮＯＲＤ３Ａｓｐｅｃｉｆｉｃａｌｌｙｓｅｎｓｉｔｉｚｅｓｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓｔｏ５－ＦＵｂｙｓｐｏｎｇｉｎｇｍｉＲ－１８５－５ｐｔｏｅｎｈａｎｃｅＵＭＰＳｅｘｐｒｅｓｓｉｏｎ［Ｊ］．ＣｅｌｌＤｅａｔｈ＆Ｄｉｓｅａｓｅ，２０２０，１１（５）：３２９．［４］㊀芮小慧．长链非编码ＲＮＡＣ５ｏｒｆ６６－ＡＳ１在宫颈癌发生中的作用及其机制研究［Ｄ］．苏州：苏州大学，２０１９．ＲｕｉＸＨ．Ｔｈｅｒｏｌｅｏｆｌｏｎｇｎｏｎ⁃ｃｏｄｉｎｇＲＮＡＣ５ＯＲＦ６６－ＡＳ１ｉｎｔｈｅｐａｔｈｏｇｅｎｅｓｉｓｏｆｃｅｒｖｉｃａｌｃａｎｃｅｒａｎｄｉｔｓｍｅｃｈａｎｉｓｍ［Ｄ］．Ｓｕｚｈｏｕ：ＳｏｏｃｈｏｗＵｎｉｖｅｒｓｉｔｙ，２０１９．［５］㊀尹冶，丁明霞，陈振杰，等．基于ＴＣＧＡ和ＧＥＯ数据库构建前列腺癌ｃｅＲＮＡ网络并筛选相关ｌｎｃＲＮＡｓ［Ｊ］．临床肿瘤学杂志，２０２０，２５（１１）：１０１１－１０１７．ＹｉｎＹ，ＤｉｎｇＭＸ，ＣｈｅｎＺＪ，ｅｔａｌ．ＣｏｎｓｔｒｕｃｔｉｏｎｏｆｐｒｏｓｔａｔｅｃａｎｃｅｒＣｅｒｎａｎｅｔｗｏｒｋａｎｄｓｃｒｅｅｎｉｎｇｏｆｌｎｃｒｎａｓｂａｓｅｄｏｎＴＣＧＡａｎｄｇｅｏｄａｔａｂａｓｅ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｌｉｎｉｃａｌＯｎｃｏｌｏｇｙ，２０２０，２５（１１）：１０１１－１０１７．［６］㊀ＪｉａＲＫ，ＺｈａｏＨＸ，ＪｉａＭＷ．Ｉｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｃｏ⁃ｅｘｐｒｅｓｓｉｏｎｍｏｄｕｌｅｓａｎｄｐｏｔｅｎｔｉａｌｂｉｏｍａｒｋｅｒｓｏｆｂｒｅａｓｔｃａｎｃｅｒｂｙＷＧＣＮＡ［Ｊ］．Ｇｅｎｅ，２０２０，７５０：１４４７５７．［７］㊀ＷｅｉＪＲ，ＤｏｕＱＳ，ＢａＦＴ，ｅｔａｌ．ＩｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｌｎｃＲＮＡａｎｄｍＲＮＡｅｘｐｒｅｓｓｉｏｎｐｒｏｆｉｌｅｓｉｎｄｏｒｓａｌｒｏｏｔｇａｎｇｌｉｏｎｉｎｒａｔｓｗｉｔｈｃａｎｃｅｒ⁃ｉｎｄｕｃｅｄｂｏｎｅｐａｉｎ［Ｊ］．ＢｉｏｃｈｅｍｉｃａｌａｎｄＢｉｏｐｈｙｓｉｃａｌＲｅｓｅａｒｃｈＣｏｍｍｕｎｉｃａｔｉｏｎｓ，２０２１，５７２：９８－１０４．［８］㊀ＹｉｎＨＴ，ＳｈａｎｇＱ，ＺｈａｎｇＳＬ，ｅｔａｌ．ＣｏｍｐｒｅｈｅｎｓｉｖｅａｎａｌｙｓｉｓｏｆｌｎｃＲＮＡ⁃ｍＲＮＡｒｅｇｕｌａｔｏｒｙｎｅｔｗｏｒｋｉｎＢｍＮＰＶｉｎｆｅｃｔｅｄｃｅｌｌｓｔｒｅａｔｅｄｗｉｔｈＨｓｐ９０ｉｎｈｉｂｉｔｏｒ［Ｊ］．ＭｏｌｅｃｕｌａｒＩｍｍｕｎｏｌｏｇｙ，２０２０，１２７：２３０－２３７．［９］㊀ＢａｒｂｅｒａｎＡ，ＢａｔｅｓＳＴ，ＣａｓａｍａｙｏｒＥＯ，ｅｔａｌ．Ｕｓｉｎｇｎｅｔｗｏｒｋ㊃６４２㊃北京生物医学工程㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀第４２卷ａｎａｌｙｓｉｓｔｏｅｘｐｌｏｒｅｃｏ⁃ｏｃｃｕｒｒｅｎｃｅｐａｔｔｅｒｎｓｉｎｓｏｉｌｍｉｃｒｏｂｉａｌｃｏｍｍｕｎｉｔｉｅｓ［Ｊ］．ＴｈｅＩＳＭＥＪｏｕｒｎａｌ，２０１１，６（２）：３４３－３５１．［１０］㊀ＳｉｌａｎｄＨ，ＳｒｅｉｄｅＫ，ＪａｎｓｓｅｎＥ，ｅｔａｌ．ＥｍｅｒｇｉｎｇｃｏｎｃｅｐｔｓｏｆａｐｏｌｉｐｏｐｒｏｔｅｉｎＤｗｉｔｈｐｏｓｓｉｂｌｅｉｍｐｌｉｃａｔｉｏｎｓｆｏｒｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＣｅｌｌｕｌａｒＯｎｃｏｌｏｇｙ，２００７，２９（３）：１９５－２０９．［１１］㊀ＳｈａＳ，ＹｕａｎＤＹ，ＬｉｕＹＪ，ｅｔａｌ．Ｔａｒｇｅｔｉｎｇｌｏｎｇｎｏｎ⁃ｃｏｄｉｎｇＲＮＡＤＡＮＣＲｉｎｈｉｂｉｔｓｔｒｉｐｌｅｎｅｇａｔｉｖｅｂｒｅａｓｔｃａｎｃｅｒｐｒｏｇｒｅｓｓｉｏｎ［Ｊ］．ＢｉｏｌｏｇｙＯｐｅｎ，２０１７，６（９）：１３１０－１３１６．［１２］㊀ＦａｎＳＬ，ＦａｎＣＮ，ＬｉｕＮ，ｅｔａｌ．Ｄｏｗｎｒｅｇｕｌａｔｉｏｎｏｆｔｈｅｌｏｎｇｎｏｎ⁃ｃｏｄｉｎｇＲＮＡＺＦＡＳ１ｉｓａｓｓｏｃｉａｔｅｄｗｉｔｈｃｅｌｌｐｒｏｌｉｆｅｒａｔｉｏｎ，ｍｉｇｒａｔｉｏｎａｎｄｉｎｖａｓｉｏｎｉｎｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＭｏｌｅｃｕｌａｒＭｅｄｉｃｉｎｅＲｅｐｏｒｔｓ，２０１８，１７（５）：６４０５－６４１２．［１３］㊀ＣｏｒｓｏＧ，ＩｎｔｒａＭ，ＴｒｅｎｔｉｎＣ，ｅｔａｌ．ＣＤＨ１ｇｅｒｍｌｉｎｅｍｕｔａｔｉｏｎｓａｎｄｈｅｒｅｄｉｔａｒｙｌｏｂｕｌａｒｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＦａｍｉｌｉａｌＣａｎｃｅｒ，２０１６，１５（２）：２１５－２１９．［１４］㊀ＬｉＪＴ，ＬｉＬ，ＹｕａｎＨＺ，ｅｔａｌ．Ｕｐ⁃ｒｅｇｕｌａｔｅｄｌｎｃＲＮＡＧＡＳ５ｐｒｏｍｏｔｅｓｃｈｅｍｏｓｅｎｓｉｔｉｖｉｔｙａｎｄａｐｏｐｔｏｓｉｓｏｆｔｒｉｐｌｅ⁃ｎｅｇａｔｉｖｅｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓ［Ｊ］．ＣｅｌｌＣｙｃｌｅ，２０１９，１８（１６）：１９６５－１９７５．［１５］㊀ＡｌｊｏｈａｎｉＡＩ，ＴｏｓｓＭＳ，ＫｕｒｏｚｕｍｉＳ，ｅｔａｌ．Ｔｈｅｐｒｏｇｎｏｓｔｉｃｓｉｇｎｉｆｉｃａｎｃｅｏｆｗｉｌｄ⁃ｔｙｐｅｉｓｏｃｉｔｒａｔｅｄｅｈｙｄｒｏｇｅｎａｓｅ２（ＩＤＨ２）ｉｎｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＢｒｅａｓｔＣａｎｃｅｒＲｅｓｅａｒｃｈａｎｄＴｒｅａｔｍｅｎｔ，２０２０，１７９（１）：７９－９０．［１６］㊀ＤｏｎｇＨＹ，ＷａｎｇＷ，ＭｏＳＷ，ｅｔａｌ．ＳＰ１⁃ｉｎｄｕｃｅｄｌｎｃＲＮＡＡＧＡＰ２－ＡＳ１ｅｘｐｒｅｓｓｉｏｎｐｒｏｍｏｔｅｓｃｈｅｍｏｒｅｓｉｓｔａｎｃｅｏｆｂｒｅａｓｔｃａｎｃｅｒｂｙｅｐｉｇｅｎｅｔｉｃｒｅｇｕｌａｔｉｏｎｏｆＭｙＤ８８［Ｊ］．ＪｏｕｒｎａｌｏｆＥｘｐｅｒｉｍｅｎｔａｌ＆ＣｌｉｎｉｃａｌＣａｎｃｅｒＲｅｓｅａｒｃｈ，２０１８，３７（１）：２０２．［１７］㊀ＣｈｅｎＪＲ，ＬｉｕＣ，ＣｅｎＪＭ，ｅｔａｌ．ＫＥＧＧ⁃ｅｘｐｒｅｓｓｅｄｇｅｎｅｓａｎｄｐａｔｈｗａｙｓｉｎｔｒｉｐｌｅｎｅｇａｔｉｖｅｂｒｅａｓｔｃａｎｃｅｒ：Ｐｒｏｔｏｃｏｌｆｏｒａｓｙｓｔｅｍａｔｉｃｒｅｖｉｅｗａｎｄｄａｔａｍｉｎｉｎｇ［Ｊ］．Ｍｅｄｉｃｉｎｅ，２０２０，９９（１８）：１３２－１４０．［１８］㊀ＹｕａｎＣＬ，ＪｉａｎｇＸＭ，ＹｉＹ，ｅｔａｌ．ＩｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｄｉｆｆｅｒｅｎｔｉａｌｌｙｅｘｐｒｅｓｓｅｄｌｎｃＲＮＡｓａｎｄｍＲＮＡｓｉｎｌｕｍｉｎａｌ⁃ＢｂｒｅａｓｔｃａｎｃｅｒｂｙＲＮＡ⁃ｓｅｑｕｅｎｃｉｎｇ［Ｊ］．ＢＭＣＣａｎｃｅｒ，２０１９，１９（１）：１１７１．［１９］㊀ＱｉＦ，ＱｉｎＷＸ，ＺａｎｇＹＳ．Ｍｏｌｅｃｕｌａｒｍｅｃｈａｎｉｓｍｏｆｔｒｉｐｌｅ⁃ｎｅｇａｔｉｖｅｂｒｅａｓｔｃａｎｃｅｒ⁃ａｓｓｏｃｉａｔｅｄＢＲＣＡ１ａｎｄｔｈｅｉｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｓｉｇｎａｌｉｎｇｐａｔｈｗａｙｓ［Ｊ］．ＯｎｃｏｌｏｇｙＬｅｔｔｅｒｓ，２０１９，１７（３）：２９０５－２９１４．［２０］㊀ＤｅｎｇＪＬ，ＸｕＹＨ，ＷａｎｇＧ．Ｉｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｐｏｔｅｎｔｉａｌｃｒｕｃｉａｌｇｅｎｅｓａｎｄｋｅｙｐａｔｈｗａｙｓｉｎｂｒｅａｓｔｃａｎｃｅｒｕｓｉｎｇｂｉｏｉｎｆｏｒｍａｔｉｃａｎａｌｙｓｉｓ［Ｊ］．ＦｒｏｎｔｉｅｒｓｉｎＧｅｎｅｔｉｃｓ，２０１９，１０（１）：６９５－７１３．［２１］㊀ＬｉｕＤＨ，ＺｈａｎｇＪＦ，ＬｉＬ，ｅｔａｌ．Ｄｅｔｅｃｔｉｏｎｏｆｃｒｉｔｉｃａｌｇｅｎｅｓａｓｓｏｃｉａｔｅｄｗｉｔｈｐｏｏｒｐｒｏｇｎｏｓｉｓｉｎｂｒｅａｓｔｃａｎｃｅｒｖｉａｉｎｔｅｇｒａｔｅｄｂｉｏｉｎｆｏｒｍａｔｉｃｓａｎａｌｙｓｅｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＢ．Ｕ．ＯＮ．，２０２０，２５（６）：２５３７－２５４５．［２２］㊀ＸｕＴ，ＷａｎｇＱＧ，ＬｉｕＭ．ＡｎｅｔｗｏｒｋｐｈａｒｍａｃｏｌｏｇｙａｐｐｒｏａｃｈｔｏｅｘｐｌｏｒｅｔｈｅｐｏｔｅｎｔｉａｌｍｅｃｈａｎｉｓｍｓｏｆＨｕａｎｇｑｉｎ⁃Ｂａｉｓｈａｏｈｅｒｂｐａｉｒｉｎｔｒｅａｔｍｅｎｔｏｆｃａｎｃｅｒ［Ｊ］．ＭｅｄｉｃａｌＳｃｉｅｎｃｅＭｏｎｉｔｏｒ：ＩｎｔｅｒｎａｔｉｏｎａｌＭｅｄｉｃａｌＪｏｕｒｎａｌｏｆＥｘｐｅｒｉｍｅｎｔａｌａｎｄＣｌｉｎｉｃａｌＲｅｓｅａｒｃｈ，２０２０，２６：ｅ９２３１９９．［２３］㊀ＺｈａｏＬ，ＬｉｎＭ，ＷａｎｇＳＳ．ＩｄｅｎｔｉｆｉｃａｔｉｏｎｏｆｈｕｍａｎｐｒｏｌａｃｔｉｎｏｍａｒｅｌａｔｅｄｇｅｎｅｓｂｙＤＮＡｍｉｃｒｏａｒｒａｙ［Ｊ］．ＪｏｕｒｎａｌｏｆＣａｎｃｅｒＲｅｓｅａｒｃｈａｎｄＴｈｅｒａｐｅｕｔｉｃｓ，２０１４，１０（３）：５４４－５４８．［２４］㊀ＡｚａｄＡＫ，ＬａｗｅｎＡ，ＫｅｉｔｈＪＭ．Ｐｒｅｄｉｃｔｉｏｎｏｆｓｉｇｎａｌｉｎｇｃｒｏｓｓ⁃ｔａｌｋｓｃｏｎｔｒｉｂｕｔｉｎｇｔｏａｃｑｕｉｒｅｄｄｒｕｇｒｅｓｉｓｔａｎｃｅｉｎｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓｂｙＢａｙｅｓｉａｎｓｔａｔｉｓｔｉｃａｌｍｏｄｅｌｉｎｇ［Ｊ］．ＢＭＣＳｙｓｔｅｍｓＢｉｏｌｏｇｙ，２０１５，９（１）：２．［２５］㊀ＫｈａｎｎａＰ，ＬｅｅＪＳ，ＳｅｒｅｅｍａｓｐｕｎＡ，ｅｔａｌ．ＧＲＡＭＤ１ＢｒｅｇｕｌａｔｅｓｃｅｌｌｍｉｇｒａｔｉｏｎｉｎｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓｔｈｒｏｕｇｈＪＡＫ／ＳＴＡＴａｎｄＡｋｔｓｉｇｎａｌｉｎｇ［Ｊ］．ＳｃｉｅｎｔｉｆｉｃＲｅｐｏｒｔｓ，２０１８，８（１）：９５１１．［２６］㊀ＺｈｕａｎｇＺＧ，ＷａｎｇＫ，ＣｈｅｎｇＸＬ，ｅｔａｌ．ＬＫＢ１ｉｎｈｉｂｉｔｓｂｒｅａｓｔｃａｎｃｅｒｐａｒｔｉａｌｌｙｔｈｒｏｕｇｈｒｅｐｒｅｓｓｉｎｇｔｈｅＨｅｄｇｅｈｏｇｓｉｇｎａｌｉｎｇｐａｔｈｗａｙ［Ｊ］．ＰｌｏｓＯｎｅ，２０１３，８（７）：ｅ６７４３１．［２７］㊀ＦｅｒｎａｎｄａＭＲ，ＳａｒａＴＯ，ＪｏａｏＡＭ．ＳｅｒｉｎｅｐｅｐｔｉｄａｓｅｉｎｈｉｂｉｔｏｒＫｕｎｉｔｚｔｙｐｅ２（ＳＰＩＮＴ２）ｉｎｃａｎｃｅｒｄｅｖｅｌｏｐｍｅｎｔａｎｄｐｒｏｇｒｅｓｓｉｏｎ［Ｊ］．Ｂｉｏｍｅｄｉｃｉｎｅ＆Ｐｈａｒｍａｃｏｔｈｅｒａｐｙ，２０１８，１０１（１）：２７８－２８６．［２８］㊀ＺｈｏｎｇＷＬ，ＬｕＭＹ，ＳｉＣＦ，ｅｔａｌ．Ｐｒｏｇｒｅｓｓｏｆｒｅｓｅａｒｃｈｏｎｔａｒｇｅｔｅｄｔｈｅｒａｐｙｆｏｒｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＭｏｄｅｒｎＯｎｃｏｌｏｇｙ，２０１８，２６（４）：６２２－６２６．［２９］㊀ＺｈｏｕＴＺ，ＬｉｎＫ，ＮｉｅＪＪ，ｅｔａｌ．ＬｎｃＲＮＡＳＰＩＮＴ１－ＡＳ１ｐｒｏｍｏｔｅｓｂｒｅａｓｔｃａｎｃｅｒｐｒｏｌｉｆｅｒａｔｉｏｎａｎｄｍｅｔａｓｔａｓｉｓｂｙｓｐｏｎｇｉｎｇｌｅｔ－７ａ／ｂ／ｉ－５ｐ［Ｊ］．Ｐａｔｈｏｌｏｇｙ⁃ＲｅｓｅａｒｃｈａｎｄＰｒａｃｔｉｃｅ，２０２１，２１７（１）：１５３２６８．［３０］㊀郭辉，张斌，胡利民．微小ＲＮＡ－１４５靶向调控ＳＯＸ１１表达及其对乳腺癌细胞增殖和凋亡的影响［Ｊ］．临床肿瘤学杂志，２０１８，１９５（３）：２００－２０５．ＧｕｏＨ，ＺｈａｎｇＢ，ＨｕＬＭ．ＥｘｐｒｅｓｓｉｏｎｏｆｍｉｃｒｏＲＮＡ－１４５ｏｎｔｈｅｔａｒｇｅｔｅｄｒｅｇｕｌａｔｉｏｎｏｆＳＯＸ１１ａｎｄｉｔｓｅｆｆｅｃｔｏｎｐｒｏｌｉｆｅｒａｔｉｏｎａｎｄａｐｏｐｔｏｓｉｓｏｆｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｌｉｎｉｃａｌＯｎｃｏｌｏｇｙ，２０１８，１９５（３）：２００－２０５．［３１］㊀李莉，孟少达，邱爽，等．乳腺癌根治术后皮瓣下积液感染病原学特点及相关因素分析［Ｊ］．中华医院感染学杂志，２０１８，２８（３）：４１４－４１７．ＬｉＬ，ＭｅｎｇＳＤ，ＱｉｕＳ，ｅｔａｌ．Ｐａｔｈｏｇｅｎｉｃｃｈａｒａｃｔｅｒｉｓｔｉｃｓａｎｄｒｉｓｋｆａｃｔｏｒｓｏｆｉｎｆｅｃｔｉｏｎｓｏｆｓｋｉｎｆｌａｐａｆｔｅｒｒａｄｉｃａｌｍａｓｔｅｃｔｏｍｙｆｏｒｂｒｅａｓｔｃａｎｃｅｒ［Ｊ］．ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＮｏｓｏｃｏｍｉｏｌｏｇｙ，２０１８，２８（３）：４１４－４１７．［３２］㊀ＰｅｉｆｆｅｒＤＳ，ＷｙａｔｔＤ，ＺｌｏｂｉｎＡ，ｅｔａｌ．ＤＡＸＸｓｕｐｐｒｅｓｓｅｓｔｕｍｏｒ⁃ｉｎｉｔｉａｔｉｎｇｃｅｌｌｓｉｎｅｓｔｒｏｇｅｎｒｅｃｅｐｔｏｒ⁃ｐｏｓｉｔｉｖｅｂｒｅａｓｔｃａｎｃｅｒｆｏｌｌｏｗｉｎｇｅｎｄｏｃｒｉｎｅｔｈｅｒａｐｙ［Ｊ］．ＣａｎｃｅｒＲｅｓｅａｒｃｈ，２０１９，７９（１９）：４９６５－４９７７．［３３］㊀ＸｕＧ，ＢｕＳＳ，ＷａｎｇＸＳ，ｅｔａｌ．ＳｕｐｐｒｅｓｓｉｏｎｏｆＣＣＴ３ｉｎｈｉｂｉｔｓｔｈｅｐｒｏｌｉｆｅｒａｔｉｏｎａｎｄｍｉｇｒａｔｉｏｎｉｎｂｒｅａｓｔｃａｎｃｅｒｃｅｌｌｓ［Ｊ］．ＣａｎｃｅｒＣｅｌｌＩｎｔｅｒｎａｔｉｏｎａｌ，２０２０，２０（１）：２１８．［３４］㊀ＬｉｎＣＹ，ＢｅａｔｔｉｅＡ，ＢａｒａｄａｒａｎＢ，ｅｔａｌ．ＣｏｎｔｒａｄｉｃｔｏｒｙｍＲＮＡａｎｄｐｒｏｔｅｉｎｍｉｓｅｘｐｒｅｓｓｉｏｎｏｆＥＥＦ１Ａ１ｉｎｄｕｃｔａｌｂｒｅａｓｔｃａｒｃｉｎｏｍａｄｕｅｔｏｃｅｌｌｃｙｃｌｅｒｅｇｕｌａｔｉｏｎａｎｄｃｅｌｌｕｌａｒｓｔｒｅｓｓ［Ｊ］．ＳｃｉｅｎｔｉｆｉｃＲｅｐｏｒｔｓ，２０１８，８（１）：１３９０４．（２０２１－０７－２５收稿，２０２１－１０－１４修回）㊃７４２㊃第３期㊀㊀㊀㊀㊀㊀黄彦祚，等：乳腺癌相关的ｌｎｃＲＮＡ⁃ｍＲＮＡ共表达扰动网络构建。
癌症TCGA数据库中乳腺癌预后数据的挖掘
癌症TCGA数据库中乳腺癌预后数据的挖掘Mian Khizar Hayat;王铭裕;李硕磊【摘要】近年来,乳腺癌发病率逐渐上升,并且呈现出年轻化趋势.使用TCGA数据库中已有的基因信息筛选鉴定出与乳腺癌预后相关的基因.为排除癌组织和正常组织取样时间不同造成的差异,我们选取了113对同时检测乳腺癌区和其相对应癌旁正常组织的样品,从TCGA数据库调取转录组数据,对这些数据通过DEseq进行差异表达分析,筛选出1428个差异表达基因.对差异表达基因进行基因本体GO,代谢通路KEGG,疾病本体DO和富集分析获得68个与乳腺癌相关的差异表达的关键基因;采用数据库中所用癌症的表达数据(共1097例)对这些乳腺癌相关基因进行总生存率分析,筛选出8个与乳腺癌预后相关的基因.结果显示在乳腺癌病人中PGLYRP2、SEMA3G、PROL1及SLC7A3的高表达伴随着乳腺癌病人的预后良好,而SKA1、BIRC5、RRM2和AURKA基因的高表达伴随着乳腺癌病人的预后不良.这8个基因有可能是乳腺癌预后相关的重要基因,这为乳腺癌病人的预后治疗提供了新的方向与思路,并可能通过调控基因水平来尽可能地控制预后.【期刊名称】《生物学杂志》【年(卷),期】2018(035)004【总页数】5页(P62-66)【关键词】癌症基因组图谱数据库;乳腺癌;差异表达基因;预后【作者】Mian Khizar Hayat;王铭裕;李硕磊【作者单位】兰州大学生命科学学院生物物理所,兰州730000;兰州大学生命科学学院生物物理所,兰州730000;兰州大学生命科学学院生物物理所,兰州730000【正文语种】中文【中图分类】R737.9乳腺癌是危害女性身心健康的最主要的恶性肿瘤，男性乳腺癌患者比较少见，Cancer Statistics 在 2017 年的统计数据显示乳腺癌在女性癌症发病中占据了 30%的比例[1]。
近年来，乳腺癌的发病率逐年上升，并且年轻化趋势明显[2]。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
TCGA的乳腺癌RNA-seq数据WGCNA分析示例WGCNA（Weighted Correlation Network analysis）是一个基于基因表达数据，构建基因共表达网络的方法。
WGCNA 和差异基因分析（DEG）的差异在于DEG主要分析样本和样本之间的差异，而WGCNA主要分析的是基因和基因之间的关系。
WGCNA通过分析基因之间的关联关系，将基因区分为多个模块。
而最后通过这些模块和样本表型之间的关联性分析，寻找特定表型的分子特征。
网上例子千千万，但是大部分都是从文档翻译而来，要用起来还是有些费劲，要深入的可以移步这里：/~yandell/statgen/ucla/WGCNA/wgcna. html下面我将根据TCGA乳腺癌基因表达数据以及乳腺癌压型数据，一步一步的使用WGCNA来进行乳腺癌各个亚型共表达模块的挖掘#############数据准备#############首先我们需要下载TCGA 的乳腺癌的RNA-seq数据以及临床病理资料，我这里使用我们自己开发的TCGA简易下载工具进行下载首先下载RNA-Seq：下载之后共得到1215个样本表达数据进一步下载临床病理资料进一步点击ClinicalFull 按钮对病理资料进行提取得到ClinicalFull_matrix.txt文件，使用Excel打开ClinicalFull_matrix.txt文件可以看到共有301列信息，包含了各种用药，随访，预后等等信息，我们这里选择乳腺癌ER、PR、HER2的信息，去除其他用不上的信息，然后选择了其中有明确ER、PR、HER2阳性阴性的样本，随机拿100个做例子吧样本筛选完了，现在轮到怎么获取这些样本的RNA-seq数据啦，前面下载了一千多个样本的RNA-seq,从里面找到这一百个样本的表达数据其实也是不需要变成的啦，看清楚咯首先打开RNA-Seq数据目录的fileID.tmp(用Excel打开)，然后可以看到两列：将第二列复制，并且替换-01.gz为空使用Excel的vlookup命令将临床病理资料的那100个样本进行映射然后筛选非N/A的就得到了这一百个样本对于的RNA-seq数据信息进一步删除其他的样本，还原成fileID.tmp格式保存退出：然后使用TCGA简易小工具“合并文件”按钮就得到表达矩阵了，进一步使用ENSD_ID转换按钮就得到了基因表达矩阵和lncRNA表达矩阵了#################R代码实现WGCNA##############setwd('E:/rawData/TCGA_DATA/TC GA-BRCA')samples=read.csv('ClinicalFull_matrix.txt',sep = '\t',s = 1)dim(samples)#[1] 100 3expro=read.csv('Merge_matrix.txt.cv.txt',sep = '\t',s = 1)dim(expro)#[1] 24991 100数据读取完成，从上述结果可以看出100个样本，有24991个基因，这么多基因全部用来做WGCNA 很显然没有必要，我们只要选择一些具有代表性的基因就够了，这里我们采取的方式是选择在100个样本中方差较大的那些基因（意味着在不同样本中变化较大）继续命令：m.vars=apply(expro,1,var)expro.upper=expro[which(m.vars>quantile(m.vars, probs =seq(0, 1, 0.25))[4]),]##选择方差最大的前25%个基因作为后续WGCNA的输入数据集通过上述步骤拿到了6248个基因的表达谱作为WGCNA的输入数据集，进一步的我们需要看看样本之间的差异情况datExpr=as.data.frame(t(expro.upper)); gsg = goodSamplesGenes(datExpr, verbose = 3);gsg$allOKsampleTree = hclust(dist(datExpr), method = 'average')plot(sampleTree, main = 'Sample clustering to detect outliers' , sub='', xlab='')从图中可看出大部分样本表现比较相近，而有两个离群样本，对后续的分析可能造成影响，我们需要将其去掉，共得到98个样本clust =cutreeStatic(sampleTree, cutHeight = 80000, minSize = 10) table(clust)#clust#0 1#2 98keepSamples = (clust==1)datExpr = datExpr[keepSamples, ]nGenes = ncol(datExpr)nSamples = nrow(datExpr)save(datExpr, file = 'FPKM-01-dataInput.RData')得到最终的数据矩阵之后，我们需要确定软阈值，从代码中可以看出pickSoftThreshold很简单，就两个参数，其他默认即可powers = c(c(1:10), seq(from = 12, to=20, by=2))sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)##画图##par(mfrow = c(1,2));cex1 = 0.9;plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab='Soft Threshold (power)',ylab='Scale Free Topology Model Fit,signed R^2',type='n',main = paste('Scale independence'));text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=cex1,col='red');abline(h=0.90,col='red')plot(sft$fitIndices[,1], sft$fitIndices[,5],xlab='Soft Threshold (power)',ylab='Mean Connectivity', type='n',main = paste('Mean connectivity'))text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers,cex=cex1,col='red')从图中可以看出这个软阈值选择7比较合适,选择软阈值7进行共表达模块挖掘pow=7net = blockwiseModules(datExpr, power = pow, maxBlockSize = 7000,TOMType = 'unsigned', minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMs = TRUE,saveTOMFileBase = 'FPKM-TOM',verbose = 3)table(net$colors)# open a graphics window#sizeGrWindow(12, 9)# Convert labels to colors for plottingmergedColors = labels2colors(net$colors)# Plot the dendrogram and the module colors underneath plotDendroAndColors(net$dendrograms[[1]],mergedColors[net$blockGenes[[1]]],groupLabels = c('Module colors','GS.weight'),dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05)从图中可以看出大部分基因在灰色区域，灰色部分一般认为是没有模块接受的，从这里也可以看出其实咱们选择的这些基因并不是特别好那么做到这一步了基本上共表达模块做完了，每个颜色代表一个共表达模块，统计看看各个模块下的基因个数：那么得到模块之后下一步该做啥呢，或许很多人到这就不知道如何继续分析了这里就需要咱们利用这些模块搞事情了，举个例子如果你是整合的数据（整合lnc与gene），那么同时在某个模块中的基因和lncRNA咱们可以认为是共表达的，这便是lnc-gene共表达关系的获得途径之一了，进一步你可以根据该模块的基因-lnc-基因之间的关系绘制出共表达网络今天咱们这里不讲这个，而是跟表型关联，咱们已经拿到了这98个样本的ER、PR、HER2阳性阴性信息，那么进一步的咱们可以看看哪些共表达模块跟ER、PR、HER2阴性最相关，代码如下：moduleLabelsAutomatic =net$colorsmoduleColorsAutomatic =labels2colors(moduleLabelsAutomatic) moduleColorsFemale = moduleColorsAutomaticMEs0 = moduleEigengenes(datExpr, moduleColorsFemale)$eigengenesMEsFemale = orderMEs(MEs0)samples=samples[match(s(datExpr),paste0(gsub('-','.',r s(samples)),'.01')),]#匹配98个样本数据trainDt=as.matrix(cbind(ifelse(samples[,1]=='Positive',0,1),#将阴性的样本标记为1ifelse(samples[,2]=='Positive',0,1),#将阴性的样本标记为1ifelse(samples[,3]=='Positive',0,1),#将阴性的样本标记为1ifelse(samples[,1]=='Negative'&samples[,2]=='Negative'&sampl es[,3]=='Negative',1,0))#将三阴性的样本标记为1)#得到一个表型的0-1矩阵modTraitCor = cor(MEsFemale, trainDt, use = 'p')colnames(MEsFemale)modTraitP = corPvalueStudent(modTraitCor, nSamples)textMatrix = paste(signif(modTraitCor, 2), '\n(',signif(modTraitP, 1), ')', sep = '')dim(textMatrix) = dim(modTraitCor)labeledHeatmap(Matrix = modTraitCor, xLabels =colnames(trainDt), yLabels = names(MEsFemale),ySymbols = colnames(modlues), colorLabels = FALSE, colors = greenWhiteRed(50),textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1), main = paste('Module-trait relationships'))最终找到几个共表达网络与三阴性表型最相关的模块。