本研究对两大公共数据库（GEO和ArrayExpress）进行了系统检索，共纳入了3个与肺癌相关转录组高通量测序（RNAsequencing，RNA-seq）数据和两个肺癌相关的TCGARNA-seq数据（LUAD和LUSC），
根据目前主流研究推荐重新搭建了RNA-seq数据分析流程，对3个GEO数据中的2个原始数据重新进行标准化流程分析，得到了转录组基因计数文件；对于2个肺癌相关的TCGA数据，由于没有获取原始测序文件的权限，因此直接利用GDC的API下载了TCGA提供的LUAD和LUSC转录组表达计数（counts）表达矩阵。
随后对五个数据集进行了合并，利用DESeq2和edgeR进行差异表达分析，进而利用limma程序包去除批次效应（batch effect）并利用DESeq2程序包中vst功能进行正态化转换，获得可用于后续WGCNA和机器学习的基因表达矩阵。
使用加权基因共表达网络分析（Weighted Gene Co-expression NetworkAnalysis，WGCNA）方法，对1327例NSCLC组织和231例癌旁正常对照的基因转录组表达谱构建基因共表达网络、划分基因模块并寻找与NSCLC密切相关的基因模块，
对与NSCLC密切相关的模块进行基因本体（Gene Ontology，GO）和KEGG通路分析。将基因差异表达分析得出的结果与WGCNA结果进行联合分析，可以得到一批与NSCLC密切相关的差异表达基因，
随后从去除批次效应并进行正态化转换的转录组表达数据中获取这些基因的表达数据，利用十折交叉验证结合机器学习的方法，构建NSCLC预测模型，在验证组对预测模型效果进行评价。
DESeq2和edgeR程序包差异表达基因分析结果显示，当差异表达基因定义为|log2FC| > 1且校正P < 0.01时，共有2956个基因在NSCLC中高表达，其中2124个基因为蛋白编码基因（mRNA），254个基因为lncRNA，578个基因为其他类型基因；共有1790个基因低表达，其中1565个基因为mRNA，96个为lncRNA，129个基因为其他类型基因。WGCNA网络中共划分了39个基因模块，其中2个模块与非小细胞肺癌呈强相关（宝石绿模块R2= 0.60，蓝色模块R2= -0.79，均有P < 0.001），其中宝石绿
模块与NSCLC最为密切。
对宝石绿模块中基因的GO分析结果显示，这些基因为核染色体、染色体、中心体、微管组织中心、细胞骨架、微管、微管细胞骨架等组分，DNA结合、转录调控、结合ATP等生物学功能，参与增殖、细胞骨架和微观组织、有丝分裂细胞周期、核分裂、姐妹染色体分离、DNA代谢过程、DNA复制、DNA修复以及细胞DNA损伤刺激反应等生物学过程；
KEGG通路分析显示宝石绿模块基因主要富集在细胞周期、卵母细胞减数分裂、细胞衰老等信号通路，模块中差异表达基因主要参与细胞周期、卵母细胞减数分裂、孕酮介导的卵母细胞成熟、细胞衰老、P53信号通路、同源重组等信号通路。这进一步揭示了NSCLC的分子学机制。
WGCNA分析结果联合差异表达基因分析结果显示，与NSCLC最密切的宝石绿模块中，共有988个差异表达基因。
利用十折交叉验证结合机器学习方法对1558例研究对象的988个(1558*988)基因表达矩阵分析结果显示，构建的多个NSCLC预测模型具有很好的分辨能力，这些模型在验证组中也表现良好，其中SVM、XGBoost、C5.0、PLS、AdaBoost和gbm等算法构建的模型在验证组数据中预测准确率可高达0.98以上；尽管 JRip、PART、和rpart算法构建的半透明模型在验证组中准确率也较高，但是特异度较低，综合比较，选取SVM和XGBoost这类黑盒子算法模型作为最终NSCLC预测模型。本研究成功构建了多个准确度在0.98以上的NSCLC预测模型。
GEO Gene Expression Omnibus 基因表达综合数据库（美国）
GO Gene Ontology 基因本体
GS Gene Significance 基因显著性
KEGG Kyoto Encyclopedia of Genes and Genomes 京都基因和基因组百科全书
MM Module Membership 模块隶属度
MS Module Significance 模块显著性
ROC Receiver Operating Characteristic Curve 受试者工作曲线
TCGA The Cancer Genome Atlas 美国癌症基因图谱计划
WGCNA Weighted Gene Co-Expression Analysis 加权基因共表达网络分
引言的写法：背景—>WGCNA—>机器学习—>总结
如果训练集不足，获得的模型预测效果不佳。如果训练集过多，又会降低任务处理效率，同时也容易导致过拟合。交叉验（crossvalidation，CV）方法应运而生，既可以避免拟合不佳或过拟合的情况出现，也可以在保证预测结果的情况下提高处理效率[18,19]。交叉验证又称为循环估计，它可以将数据样本切割成较小的子集，并将其中一部分作为训练集来学习或者训练模型，其他子集作为验证集或测试集对模型效果进行评价。常见的交叉验证方法有 k 折交叉验证（k-fold corss-validation）、留一验证（leave one out crossvalitaion，LOOCV）。在本研究中，我们采用文献推荐的十折交叉验证（10 foldcross validation）来进行预测模型构建[19, 20]
本研究中建立的标准化数据分析流程概述如下：首先利用prefetch 软件下载相关数据并转为 fastq 文件，然后利用 HISAT2 将 fastq 文件为 sam 文件，用 samtools 将 sam 文件转为 bam 文件，随后利用 featureCounts 软件将 bam 文件转为 counts 文件。（数据一步步转化） 数据下载和处理代码详见附录 A
counts 数据用于后续差异表达基因分析，利用edgeR 标准化处理之后的 counts 数据用于 WGCNA 分析和 caret 构建肺癌预测模型。
GEO 和 ArrayExpress 是目前最常用的两个存储基因测序数据的公共数据库
芯片快车（ArrayExpress，https://www.ebi.ac.uk/arrayexpress/）是一个由多家科学杂志推荐用于存储基因芯片和测序平台的功能基因组学数据的一个数据库。
TCGA 由美国国家癌症研究所和人类基因组研究所创建，收录美国和加拿大等地区共11000 多个病人癌症组织和与之匹配的正常组织，包含 36 种癌症基因组测序数据，数据类型有单核苷酸变异数据、 RNA-Seq 数据及甲基化数据等[8]
{宝贝的开题}2000 年美国国立研究院生物技术信息中心建立基因表达数据库（GEO），其收录了世界各地实验室提交的 1871121 个样本实验数据及 3848 种研究类型的基因表达谱数据，数据类型多为微阵列数据，尤其是基因芯片数据； 2005 年美国国家癌症研究所和国家人类基因组研究所创建癌症基因数据库（TCGA），其数据来源于美国和加拿大地区 11000 多个病人的肿瘤和匹配正常组织，共包括 34 种癌症基因组测序数据，数据类型包括 RNA-Seq 数据及甲基化数据等，
{https://shengxin.ren/article/454}GEO主要包含各种芯片数据，也有少部分测序数据，与TCGA的差别在于TCGA只包含人的数据，而GEO是多物种的，GEO上有各种平台的数据，而TCGA只有测序数据，芯片数据的数据量较小，而TCGA的测序数据数据量较大。
Anders 等人对多种差异表达分析软件比较研究结果表明，传统基因芯片数据一般推荐在 R 软件中用limma 包进行差异表达基因分析，而转录组测序数据采用 DEseq2 和 edgeR 软件包进行差异表达分析效果优于其他同类软件[29]。这两个软件包在使用计数数据进行差异表达分析时使用类似的策略，edgeR 软件相对 DEseq2 来说更加保守，筛选到的差异表达基因范围更小。然而，他们还是有一些重要的区别。首先，两者默认的数据正态化方法不一样：edgeR 使用截断的平均 M 值，而 DESeq2 通过创造一个虚拟的库，并计算每个样本的相对 log 表达水平。在实际使用中，DESeq2 更为保守，而 edgeR 对离群值更为敏感。有比较研究表明这两种差异表达分析方法效果相当[30]
dds2<-DESeq(dds) #对原始的dds进行标准化
KM曲线和cox回归 ROC曲线的关系
由于转录组基因表达计数文件为离散型变量，容易受建库质量、测序深度等多种因素影响，而且来自不同研究的数据之间可能存在批次效应，因此在进行下游分析之前，需要先去除批次效应，然后对数据进行正态化转换。本研究中用 R软件中 limma 包中集成的removeBatchEffect 功能去除批次效应，用 DESeq 包中集成的 vst（varianceStabilizingTransformation）功能对去除批次效应后的数据进行正态化转换，用于后续 WGCNA 分析和 NSCLC 预测模型构建。
2.2.3.2 网络构建的前提
a 去除离群样本
b 确定软阈值 β
根据 Zhang 等人的描述，基因共表达网络应当符合无尺度特征，即服从幂律分布[14]。构建基因共表达网络需要选择软阈值 β，将共表达相似度提高到该阈值来计算邻接度（adjacency）。Zhang 等人建议通过接近无尺度拓扑的标准来选择软阈值，即可以使无尺度拓扑模型的 R2接近 0.9 并达到平台期，同时平均连接度（mean connectivity）也不至于太低（需要显著大于零）。
c 加权基因共表达网络构建
本研究中采用一步法来构建基因共表达网络、确定基因模块。基因共表达矩阵由基因两两之间 Pearson 相关系数的绝对值构成。通过选择合适的软阈值，可以通过指数函数变换的形式将相似矩阵中的相关系数进行连续变换以得到邻接函数，在此基础上计算拓扑重叠矩阵。基因模块的划分基于模块之间的连接度，因此需要将拓扑连接重合度转化为相异度。对于拓扑重叠度较高（相异度低于0.25）的模块，采用自适应性的动态剪枝算法进行合并，然后重新计算基因模块。
d 基因模块与临床信息的关联分析
通过将基因模块与临床信息进行关联分析，有助于发现与研究对象性状相关的基因模块，这些模块中的基因很有可能对于该性状具有重要意义。对于分组表型（如样本类型等），可以先计算基因显著性（GeneSignificance，GS），再将模块的显著性（ModuleSignificance，MS）定义为模块所包含所有基因的显著性平均值。一般而言，MS 值越高，说明这个模块与样本类型的关联度越高。模块中的基因在该基因模块内的模块隶属度（ModuleMembership，MM）为该基因与本模块性状特征基因的相关系数，可以用于筛选模块内的重要基因。
e 寻找枢纽基因
模块中连接度较高的基因也被称为枢纽基因（hub gene），枢纽基因在模块中可能扮演重要角色[13]。枢纽基因具有一定的保守性，它们可以作为减少其他基因突变影响的遗传缓冲器[32]。已有研究表明，枢纽基因在基因共表达网络中处于核心位置，然而却不总是有显著的生物学意义[33]。我们将与 NSCLC 密切相关的模块中前 30 个枢纽基因导出，并利用 VisAnt 软件绘制基因-基因相互作用网络图[34]
2.2.4 基因本体富集分析和 KEGG 通路富集分析
2.2.4.1 基因本体富集分析
2.2.4.2 KEGG 通路富集分析
本研究利用 R 软件中的 clusterProfiler 包进行 GO 富集分析和 KEGG 功能富集分析，对本研究中筛选的与 NSCLC 密切相关的基因模块进行注释，探索模块基因的功能、位置及所参与的信号通路
本研究使用 R 软件中的 caret （ClassificationAndREgressionTraining）软件包[38]，构建 NSCLC 预测模型。caret 软件包是一款专门用于预测模型和评估的 R软件包[39]。
一般通过以下指标来评价预测模型预测效果：真阳性率（truepositive，TP）：又称敏感度（sensitivity），即实际患病而按该模型预测结果可以正确的被判为患病的比例，它可以反映预测模型发现患病患者的能力；
通常用准确度和符合率来评价模型性能：准确度（accuracy），又叫准确性，用于反映预测模型预测结果与真实情况的接近程度；符合率，又称一致率，指预测模型预测结果符合实际情况的人数占测试/验证数据集总人数的比例。Kappa值适用于计数资料的一致性评测，如果 Kappa 值大于 0.75，说明一致性极好，如果 Kappa 值=1，说明预测结果完全符合实际情况。
这五个纳入的数据经过整理后最终形成的表达数据包含 1558 个样本和53728 个基因，其中包含 735 例肺腺癌、568 例肺鳞癌、24 例大细胞肺癌和 231例癌旁正常组织。表 3.1 中概括了这五组实验数据的细节。
DESeq2 和 edgeR 包在进行转录组差异表达分析时可以直接利用原始基因表达计数数据，因此不需要对数据进行转换。
本研究利用 limma 程序包去除批次效应，然后用DESeq2集成的vst功能对数据进行正态化转换。
利用 org.Hs.eg.db 程序包和GENECODE 提供的注释文件对筛选出来的差异表达基因进行注释[64, 65]
当筛选条件设定为|log2FoldChange|>1、FDR<0.01 时，DESeq2 筛选出 5085 个差异表达基因，其中上调基因为 3092 个，下调基因为 1993 个；edgeR 筛选出来 6175 个差异表达基因，其中上调基因为 4093 个，下调基因为 2082 个。对两个软件包筛选出来的差异表达基因求交集，两个软件同时筛选出来的差异表达基因为 4746个，其中上调基因为 2956 个，下调基因为 1790 个。 分析代码详见附录 B
{前面的}{edgeR 软件相对 DEseq2 来说更加保守，筛选到的差异表达基因范围更小。然而，他们还是有一些重要的区别。首先，两者默认的数据正态化方法不一样：edgeR 使用截断的平均 M 值，而 DESeq2 通过创造一个虚拟的库，并计算每个样本的相对 log 表达水平。在实际使用中，DESeq2 更为保守（少），而 edgeR 对离群值更为敏感。有比较研究表明这两种差异表达分析方法效果相当[30] }#宝贝的代码
dds<- DESeqDataSetFromMatrix(countData = exprSet,colData = coldata,design = ~condition) ##进行差异分析 dds2<-DESeq(dds) #对原始的dds进行标准化 resultsNames(dds2) #查看结果名称 res<-results(dds2) #用results函数提取结果，并赋值给res变量
https://www.jianshu.com/p/b7e55bacbede
#构建dds dds <- DESeqDataSetFromMatrix(countData=indata, colData = state, design= ~ condition ) #归一化，因为样本量超过了30，因此用vst vsd <- vst(dds, blind = FALSE)
datTraits 是样本特征数据 ，datExpr 是样本表达数据 （这两种数据什么关系） 来自下列数据
femData = read.csv("LiverFemale3600.csv"); datExpr = datExpr0[keepSamples, ] traitData = read.csv("ClinicalTraits.csv"); datTraits = allTraits[traitRows, -1];
先删除后转至改成frame,然后去列名 names () 列：基因 行：样本
WGCNA 关键代码
gsg = goodSamplesGenes(datExpr0, verbose = 3); gsg$allOK #是否全部通过检查 gsg$goodGenes #通过检查的基因 返回数据是 bool 型数组，通过检查的为 true，未通过的为 false gsg$goodSamples #通过检查的样本 返回数据是 bool 型数组，通过检查的为 true，未通过的为 false
sampleTree = hclust(dist(datExpr0), method = "average");clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
附录 B：差异表达基因分析 R 代码
原文：https://www.cnblogs.com/mohuishou-love/p/10707046.html