psych igraph：共表达网络构建

WGCNA分析，简单全面的最新教程详细介绍了WGCNA构建共表达网络的原理和过程，但有时由于样本太少，效果不理想，找不到合适的Power值。实验室师弟师妹经常在讨论样本较少的时候如何筛选共表达基因（转录因子的突变体/转基因材料和野生型材料在2-3个时间点进行转录组测序），加上3个生物学重复勉强构成18个样本。这样的实验设计，已经足够算出各个时间点或者处理间的差异基因和任意基因之间的相关系数。本文介绍少样本如何构建共表达无向网络，筛选共表达基因的流程。

授权转载自SimonCat，有增删。

0 实验设计

材料	时间点	重复	数据
野生型	胁迫处理1，6，12h	3	RNA-seq SE 单端测序
突变体	胁迫处理1，6，12h	3	RNA-seq SE 单端测序

思路：（1）得到胁迫处理诱导的基因和突变体野生型间的差异基因作为候选基因集（2）两两计算候选基因的皮尔逊相关系数（3）使用igrah计算各个基因的中心度。

从RNA-seq开始构建共表达网络。

1 使用软件FastQC进行测序质量评估 (NGS基础 - FASTQ格式解释和质量评估)

$ ls *.fastq.gz | xargs -n1 -P2 fastqc $1

备注： -n1 表示一个参数传递给fastqc，-P表示线程数

2 使用trimmomatic质量清理

script.trimmo.sh，对测序数据进行质量控制

$ echo 'nohup java -jar trimmomatic-0.36.jar SE $1 $1.trim.fil.gz LEADING:20 TRAILING:20 AVGQUAL:25 SLIDINGWINDOW:10:30 MINLEN:36 > $1.trim.nohup' > script.trimmo.sh

$ ls *.fastq.gz | xargs -n1 -P2 sh script.trimmo.sh

3 使用hisat2比对 (转录组分析工具哪家强？)

使用比对软件hisat2建立索引和比对

$ hisat2-build TAIR10_Chr.all.fasta TAIR10_Chr.all
$ echo ’nohup $WD/tools/hisat2-2.0.5/hisat2 -x TAIR10_Chr.all -U $1 -S $1.sam > $1.align.stat.txt’ > script.align.sh
$ ls *.trim.fil.gz | xargs -n1 -P2 sh script.align.sh

4使用htseq-count，获得每个基因的counts数目

$ echo 'htseq-count -s no $1 TAIR10.gtf > $1.count ' > script_count.sh
$ ls *sam|xargs -n1 sh script_count.sh
# 注意用grep -v 去除掉htseq-count生成文件中的无用信息
$ csvtk join -t treat_1h_rep1.count treat_1h_rep2.count treat_1h_rep3.count control_1h_rep1.count  control_1h_rep2.count  control_1h_rep3.count > 0h_count.txt
$ csvtk join -t treat_6h_rep1.count treat_6h_rep2.count treat_6h_rep3.count control_6h_rep1.count  control_6h_rep2.count  control_6h_rep3.count > 6h_count.txt
$ csvtk join -t treat_12h_rep1.count treat_12h_rep2.count treat_12h_rep3.count control_12h_rep1.count  control_12h_rep2.count  control_12h_rep3.count > 12h_count.txt
$ csvtk join -t 0h_count.txt 6h_count.txt 12h_count.txt > total_count.txt

备注：PE数据必须先根据序列的名字或者位置进行排序，本次项目是单端，所以随意。

5 使用DESeq2 计算差异基因 （DESeq2差异基因分析和批次效应移除）

将0、6、12 h的count的table依次导入,分别计算这3个时间点的差异基因。

过滤掉低表达的基因（所有样本的和 >1），log2FC >1 with adjusted p-values<0.01 筛选差异基因

>library(DESeq2)

>table=read.table('0_H_count.txt',header=F,sep='\t',row.names =1)

> colname(table)=c('0h_c_rep1','0h_c_rep2','0h_c_rep3','0h_t_rep1','0h_t_rep2','0h_t_rep3')

>countData <- countData[rowSums(countData) > 1, ]  # 去除掉低表达的基因

>condition <- c('c','c','c','t','t','t')

>coldata <- data.frame(row.names=colnames(countData), condition)

>dds <- DESeqDataSetFromMatrix(countData = countData, colData = group, design = ~condition )

>results <- results(DESeq(dds))

>res <- na.exclude(as.data.frame(results))

> filter <- res[(abs(res$log2FoldChange)>1 & res$padj<0.01),]

> write.table(filter,'DGE_0h.txt', quote=F,sep='\t', col. names = NA)

6 使用EBSeq 对count数目进行标准化

（生信宝典注：标准化后的数据可以也直接从DESeq2获取）

获得每个基因的counts数目之后，就可以计算两两基因的皮尔逊相关系数。使用EBSeq的函数median normalization 对count数目进行标准化，再进行对数化。

> if (!requireNamespace('BiocManager', quietly = TRUE))

  install.packages('BiocManager')

BiocManager::install('EBSeq', version = '3.8')

> library(EBSeq)

> NormData <- GetNormalizedMat(counts, MedianNorm(counts))

> NormData.log <- log2(NormData + 1)

> Norm.interest <- NormData.log[rownames(filter),] # filter是步骤5获得的基因集

7 使用psych 计算两两相关性

psych在CRAN官网进行下载。使用corr.test (也可以使用WGCNA的bicor函数，如果用Python更快)函数对差异基因进行两两的相关性计算。如果用所有基因计算数据量将大得可怕。因此我们只挑选差异基因进行分析。如我们对这个2 x 3的实验组进行两两差异计算，得到一个差异基因的总表，挑选这些基因进行计算。以相关性 > 0.9 p小于0.01筛选出共表达基因

>library('psych')

>Norm.interest.corr <- corr.test( t(Norm.interest), method='pearson', ci=F)

> Norm.interest.corr$p[lower.tri( Norm.interest.corr$p,diag=TRUE)]=NA

>Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))

> Correlation <- as.data.frame(as.table(Norm.interest.corr $r))

> Cor.table <- na.exclude(cbind( Correlation, Pval.adj))[,c (1,2,3,6)]

> colnames(Cor.table) <- c('gene1','gene2','cor','p.adj')

> Cor.table.filt <- Cor.table [(abs(Cor.table[,3])>0.9 & Cor.table[,4] <0.01 ),]

8 使用igraph和Cytoscape可视化网络

使用igrah对前面筛出的互作基因Cor.table.filt进行网络分析，degree 是指节点 (这里指基因)的连接度，即一个点有多少条边相连，degree centrality是某个节点的度除以网络中所有点能构成的连接数目，能反应一个基因的中心度。介度中心性 (betweenness centrality) 是指一个节点充当其它两个节点中间人的次数“被经过”的感觉，用“被经过次数”除以总的ties，即n（n-1）/2，计算方法见下图 （来源于 易生信陈亮博士的扩增子和宏基因组授课，本期推文同期有陈亮博士的“一文学会网络分析 igraph更详细假设”）。输出的文件Node_nw_st.txt和之前获得的Cor.table.filt两个表格共同用于Cytoscape可视化了。

# By Chen Liang
node_pro <- function(igraph){  
  # 节点度
  igraph.degree <- igraph::degree(igraph)
  # 节点度中心性
  igraph.cen.degree <- round(centralization.degree(igraph)$res,3)
  # 节点介数中心性
  igraph.betweenness <- round(centralization.betweenness(igraph)$res,3)
  # 节点中心性
  igraph.closeness <- round(centralization.closeness(igraph)$res,3)
  #igraph.node.pro <- cbind(igraph.degree,igraph.closeness,igraph.betweenness,igraph.cen.degree)
  #colnames(igraph.node.pro)<-c('igraph.degree','igraph.closeness','igraph.betweenness','igraph.cen.degree')
  igraph.node.pro <- data.frame(Node=names(igraph.degree), igraph.degree,igraph.closeness,igraph.betweenness,igraph.cen.degree)
}
net_node_attr <- node_pro(net)[,1:4]
write.table(net_node_attr, file='Node_nw_st.txt', row.names=F, col.names=T,quote=F,sep='\t')

 具体可视化见：

• network3D 交互式网络生成

• Cytoscape教程1

• Cytoscape之操作界面介绍

• 新出炉的Cytoscape视频教程

• Cytoscape制作带bar图和pie图节点的网络图

• Cytoscape: MCODE增强包的网络模块化分析

• 一文学会网络分析

R统计和作图

• Graphpad，经典绘图工具初学初探

• 在R中赞扬下努力工作的你，奖励一份CheatShet

• 别人的电子书，你的电子书，都在bookdown

• R语言 - 入门环境Rstudio

• R语言 - 热图绘制 (heatmap)

• R语言 - 基础概念和矩阵操作

• R语言 - 热图简化

• R语言 - 热图美化

• R语言 - 线图绘制

• R语言 - 线图一步法

• R语言 - 箱线图（小提琴图、抖动图、区域散点图）

• R语言 - 箱线图一步法

• R语言 - 火山图

• R语言 - 富集分析泡泡图

• R语言 - 散点图绘制

• R语言 - 韦恩图

• R语言 - 柱状图

• R语言 - 图形设置中英字体

• R语言 - 非参数法生存分析

• R语言 - 绘制seq logo图

• WGCNA分析，简单全面的最新教程

• 一文看懂PCA主成分分析

• 富集分析DotPlot，可以服

• 基因共表达聚类分析和可视化

• R中1010个热图绘制方法

• 还在用PCA降维？快学学大牛最爱的t-SNE算法吧, 附Python/R代码

• 一个函数抓取代谢组学权威数据库HMDB的所有表格数据

• 文章用图的修改和排版

• network3D: 交互式桑基图

• network3D 交互式网络生成

• Seq logo 在线绘制工具——Weblogo

• 生物AI插图素材获取和拼装指导

• ggplot2高效实用指南 (可视化脚本、工具、套路、配色)

• 图像处理R包magick学习笔记

• SOM基因表达聚类分析初探

• 利用gganimate可视化全球范围R-Ladies（R社区性别多样性组织）发展情况

• 一分钟绘制磷脂双分子层：AI零基础入门和基本图形绘制

• AI科研绘图（二）：模式图的基本画法

• 你知道R中的赋值符号箭头(<-)和等号（=）的区别吗？

• R语言可视化学习笔记之ggridges包

• 利用ComplexHeatmap绘制热图(一)

• ggplot2学习笔记之图形排列

• R包reshape2，轻松实现长、宽数据表格转换

• 用R在地图上绘制网络图的三种方法

• PCA主成分分析实战和可视化 附R代码和测试数据

• iTOL快速绘制颜值最高的进化树！

• 12个ggplot2扩展包帮你实现更强大的可视化

• 编程模板-R语言脚本写作：最简单的统计与绘图，包安装、命令行参数解析、文件读取、表格和矢量图输出

• R语言统计入门课程推荐——生物科学中的数据分析Data Analysis for the Life Sciences

• 数据可视化基本套路总结

• 你知道R中的赋值符号箭头<-和等号=的区别吗？

• 使用dplyr进行数据操作30例

• 交集intersect、并集union、找不同setdiff

• R包reshape2，轻松实现长、宽数据表格转换

• 1数据类型（向量、数组、矩阵、 列表和数据框）

• 2读写数据所需的主要函数、与外部环境交互

• 3数据筛选——提取对象的子集

• 4向量、矩阵的数学运算

• 5控制结构

• 6函数及作用域

• 7认识循环函数lapply和sapply

• 8分解数据框split和查看对象str

• 9模拟—随机数、抽样、线性模型

• 1初识ggplot2绘制几何对象

• 2图层的使用—基础、加标签、注释

• 3工具箱—误差线、加权数、展示数据分布

• 4语法基础

• 5通过图层构建图像

• 6标度、轴和图例

• 7定位-分面和坐标系

• 8主题设置、存储导出

• 9绘图需要的数据整理技术

• 创建属于自己的调色板

• 28个实用绘图包，总有几个适合你

• 热图绘制

• R做线性回归

• 绘图相关系数矩阵corrplot

• 相关矩阵可视化ggcorrplot

• 绘制交互式图形recharts

• 交互式可视化CanvasXpress

• 聚类分析factoextra

• LDA分析、作图及添加置信-ggord

• 解决散点图样品标签重叠ggrepel

• 添加P值或显著性标记ggpubr

• Alpha多样性稀释曲线rarefraction curve

• 堆叠柱状图各成分连线画法：突出组间变化

• 冲击图展示组间时间序列变化ggalluvial

• 桑基图riverplot

• 微生物环境因子分析ggvegan

• 五彩进化树与热图更配ggtree

• 多元回归树分析mvpart

• 随机森林randomForest 分类Classification 回归Regression

• 加权基因共表达网络分析WGCNA

• circlize包绘制circos-plot

• R语言搭建炫酷的线上博客系统

易生信系列培训课程，扫码获取免费资料

更多阅读

画图三字经 生信视频 生信系列教程 

心得体会 TCGA数据库 Linux Python 

高通量分析 免费在线画图 测序历史 超级增强子

生信学习视频 PPT EXCEL 文章写作 ggplot2

海哥组学 可视化套路 基因组浏览器

色彩搭配 图形排版 互作网络

自学生信