剪不断理还乱，那些年我们看过的转录关系

一文教你学会绘制转录调控网络

大家好，我是阿琛。在常见的生信分析文章中，我们常常会见到三大网络内容，分别是蛋白-蛋白互作网络（Protein-protein interaction， PPI），ceRNA（competing endogenous RNA）网络，以及转录调控网络。

今天，我们主要来看看如何绘制转录调控网络。

首先看一下概念:

所谓转录调控，是指通过改变转录速率从而改变基因表达的水平；经过转录调控过程，可以控制转录何时发生以及产生多少RNA。那么，我们来看看在生信分析中，转录调控网络是什么形式的。

0.期刊信息

1.文献解析

在这里，我找了一篇今年11月9号发表在Cancer Cell International杂志上的文章，题为“Identification of an immune gene signature for predicting the prognosis of patients with uterine corpus endometrial carcinoma.”。

接下来，我们来看看在这篇文章中，作者是如何进行转录调控分析的。根据实验方法部分的描述，作者使用limma包进行了差异分析，根据P<0.05和｜logFC｜>=1的标准，得到了显著差异表达基因（Differential expression genes， DEGs）。接着，分别从IMMPORT数据库和Cistrome数据库中获取了免疫基因和转录因子列表；将两个基因列表分别与DEGs取交集，得到各自的数据集。随后，作者检测了差异转录因子和差异免疫相关基因表达之间的共表达关系；设置｜相关系数｜>0.4和P <0.001的标准，得到共表达的基因信息，最后使用Cytoscape（version 3.7.1）将共表达网络信息可视化。

看到这里，我们可以发现，在生信分析中，所谓的转录调控网络实际上是基于转录因子和基因之间的共表达关系；通过两个基因之间的表达相关性，来预测可能存在正调控或者负调控关系的转录因子和基因之间的关系。

下面，我们一起来看一下如何绘制转录调控网络。

2.差异转录因子的表达

2.1 差异基因提取

根据实验方法，我们首先需要进行差异表达分析，获得DEGs列表。这里，我们已经完成了差异表达分析；并且，根据P<0.05和｜logFC｜>1.5的标准，得到了显著差异表达基因。直接使用load()函数将差异分析结果和基因表达矩阵读取进来。

rm(list=ls()) #清空环境变量options(stringsAsFactors = F)load("exp.Rdata")

在.RData文件中，一共包含两个数据框，可以看到差异基因列表中一共包含了760个DEGs。

table(diff$change)

查看一下基因的变化结果，可以看到，在差异基因中，其中343个基因的表达显著上调，417个基因显著下调。

tpms <- tpms[rownames(diff),]rownames(tpms) <- diff$Gene

接着，根据差异基因情况，提取差异基因的TPM值，并将数据框的行名从Ensemble id转换成Gene symbol。

tpms[1:3,1:3]

查看一下tpms文件的前3行前3列。

这样，差异基因的表达文件就准备好了；接下来，我们需要提取其中的转录因子。

2.2 提取差异转录因子表达

首先，我们需要去Cistrome数据库获取转录因子信息。Cistrome数据库（http://www.cistrome.org/）是一个较为全面且公开的人类和小鼠ChIP-seq及开放染色质信息的数据库；从数量上来看，Cistrome数据库可以说是收录chip类型最多的数据库之一。

此外，在该数据库中，基于对对TCGA表达图谱和公共ChIP-seq图谱的综合分析，提供了肿瘤中预测转录因子（TF）靶标和增强子谱的全面资源。点击Visit site进入Cistrome Cancer板块内容中。

进入Cistrome Cancer板块后，首先可以看到一段介绍的视频，随后是两大板块，分别是Cancer Transcription Factor Targets和Cancer Enhancer Prediction。在此，我们选择第一个转录因子板块，点击进入其中。

在转录因子板块中，我们可以看到，左边是肿瘤相关的转录因子列表，右边是相应的一些可视化分析内容。在此，我们直接将转录因子名称进行保存即可，一共得到318个肿瘤相关转录因子。

接下来，将转录因子列表读取到R中，获取差异表达的转录因子信息。

rt <- read.table("TF.txt", sep = "\t", header = T)TF_diff <- intersect(rt$gene, rownames(tpms))TF1 <- tpms[TF_diff,]Gene1 <- tpms[!(rownames(tpms) %in% TF_diff),]

可以看到，通过intersect()函数，最终得到9个差异表达的转录因子信息，分别是E2F1，EGR1，EPO，HOXC11，HOXC9，KAT2B，MAF，SOX17，以及SOX9基因。

3.相关性检验

3.1 设置相关性阈值

corFilter=0.5 #相关系数过滤标准pvalueFilter=0.001 #p值过滤标准在这里，我们将相关性分析的阈值设定为｜cor｜>0.5和P<0.001。

3.2 相关性分析

接下来，我们逐个提取转录因子表达矩阵TF1和差异基因Gene1中的基因，两两进行相关性分析。

outTab=data.frame()for(i in row.names(TF1)){for(j in row.names(Gene1)){ x=as.numeric(TF1[i,]) y=as.numeric(Gene1[j,]) corT=cor.test(x,y) cor=corT$estimate pvalue=corT$p.valueif((cor > corFilter) & (pvalue < pvalueFilter)){ outTab=rbind(outTab,cbind(TF=i,Gene=j,cor,pvalue,Regulation="postive")) }if((cor < -corFilter) & (pvalue < pvalueFilter)){ outTab=rbind(outTab,cbind(TF=i,Gene=j,cor,pvalue,Regulation="negative")) } }}

根据相关性系数cor>0.5和P<0.001，设置为正相关性；以及cor< -0.5和P<0.001，设置为负相关性。最终，一共得到了54条相互调控关系网络。

write.table(file="TFs.corResult.txt",outTab,sep="\t",quote=F,row.names=F)

最后，我们将相关性分析结果输出到工作目录中。

3.3 注释文件的准备

另外，我们还需要准备一份注释文件，作为Cytoscape软件的输入文件。

tf <- unique(outTab$TF)gene <- unique(outTab$Gene)type <- data.frame(gene = c(tf, gene), group = c(rep("TF",times=length(tf)), rep("Gene",times=length(gene))))write.table(type, file = "type.txt", sep = "\t", quote = F, row.names = F)

4.转录调控网络的制作

然后，我们将相关性结果导入Cytoscape软件（version 3.8.0）中，进行可视化展示。

打开Cytoscape软件，选择File---Import---Network from File，选择文件“TFs.corResult.txt”；