WGCNA升级：CCNA+WGCNA筛选更可靠的候选基因

生信技能树

发布于 2025-02-19 13:31:45

11300

代码可运行

文章被收录于专栏：生信技能树

运行总次数：0

代码可运行

今天给大家介绍一个WGCNA的新玩法，即一致性共表达网络分析（Consensus co-expression network analysis，CCNA），该方法于2022年9月份发表在 Journal of Advanced Research 杂志（IF=11.4）上，文献标题为：《Consensus co-expression network analysis identifies AdZAT5 regulating pectin degradation in ripening kiwifruit》。

发表团队来自浙江大学园艺植物综合生物学浙江省重点实验室，Xueren Yin(殷学仁) 教授团队，他的介绍见：https://www.aminer.cn/profile/xueren-yin/542c39addabfae1bbfd2380a。

文章背景（为什么做这个？）

如何从高通量测序数据中筛选出更可靠的候选基因用于后续的实验验证一直是一个比较难的事情。为了提高从转录组数据（来自乙烯处理和非处理果实）中筛选基因的可靠性，作者结合了一致性共表达网络分析（CCNA）和WCGNA方法，得到 85 consensus clusters，并与差异分析结果取交集，鉴定出6个负责果胶降解的细胞壁基因和4个转录因子，后续的一系列实验验证了AdZAT5对AdPL5和Ad-Gal5的调控作用。

一致性共表达网络分析（CCNA）有助于以更高的效率和准确性发现未知的调控因子！

Fig. 1. Transcriptome analysis and cluster-trait associations

候选基因筛选：

分析步骤

作者将他的分析代码放在了github上面：https://github.com/czawora/WGCNA_consensus_clustering，代码来自 Shahan et al., 2018 [35]，并进行了修改：

1、首先，低表达基因过滤：过滤所有样本中FPKM<1的基因。
2、subsampled 1000 times：对过滤后的基因随机抽取 80% genes，抽样1000次，用于WGCNA分析，分析参数为 randomized parameters (power transformation [1,2,4,8,12,16], minModuleSize [40, 60, 90, 120, 150, 180, 210], merge on eigengenes [true/false]). 对应的脚本是make_subsamp_wgcna.py 和 subsamp_wgcna.R。
3、定义矩阵A：the number of times each gene pair clustered together were defined as matrix A
4、定义矩阵B：the number of times each gene pair subsampled together were defined as matrix B，对应的脚本是 seq_adding_mat.R 和 seq_indicator_mat.R
5、定义邻接矩阵C：The adjacency matrix C, representing connection strength between gene pairs, was calculated by matrix A dividing matrix B with scripts merge_mats.R.
6、WGCNA分析：the matrix C as an input to replace gene pair connections of standard WGCNA and correlated with physiological trait data.
7、PCA分析与绘图：PCA分析R包为FactoMineR和factoextra，ggpubr 和 ggcor 包用来挑选候选基因

拿到表达矩阵

因为作者并没有给我们提供表达矩阵，这里先来做一下上游分析，拿到表达矩阵！作者的数据为：收集了乙烯处理组和空气对照组在1天储存和4天储存时的十二个样本。数据编号：GenBank SRR6885590–SRR6885601，查询得到项目编号 PRJNA445209。

1、物种是 ‘Hong Yang’

学名：Actinidia deliciosa 人家就真的叫美味！

美味猕猴桃（学名：Actinidiadeliciosa (A.Chev.) C.F.Liang & A.R.Ferguson）是猕猴桃科、猕猴桃属大型落叶藤本；花枝多数较长，达15-20厘米，被黄褐色长硬毛，毛落后仍可见到硬毛残迹。叶倒阔卵形至倒卵形，长9-11厘米，宽8-10厘米，顶端常具突尖，叶柄被黄褐色长硬毛。花较大，直径3.5厘米左右；子房被刷毛状糙毛。…

中文名： 美味猕猴桃
拉丁学名： Actinidia deliciosa (A.Chev.) C.F.Liang & A.R.Ferguson
别名： 硬毛猕猴桃

下载参考基因组，下载地址：http://kiwifruitgenome.org/ftp/A_chinensis/Hongyang/v3.0/

# 我的下载目录：/nas1/zhangj/database/genome/Hongyang/v3.0/
axel -n 100 http://kiwifruitgenome.org/ftp/A_chinensis/Hongyang/v3.0/Hongyang_cds_v3.0.fa.gz
axel -n 100 http://kiwifruitgenome.org/ftp/A_chinensis/Hongyang/v3.0/Hongyang_genome_v3.0_update.fa.gz
axel -n 100 https://kiwifruitgenome.org/ftp/A_chinensis/Hongyang/v3.0/Hongyang_mRNA_v3.0.fa.gz
axel -n 100 http://kiwifruitgenome.org/ftp/A_chinensis/Hongyang/v3.0/Hongyang_v3.0_update.gff3.gz

解压：

gunzip *
ls

# Hongyang_cds_v3.0.fa  Hongyang_genome_v3.0_update.fa  Hongyang_mRNA_v3.0.fa  Hongyang_v3.0_update.gff3

构建hisat2的索引：

nohup hisat2-build -p 30 -f Hongyang_genome_v3.0_update.fa Hongyang_genome 1>index.log  2>&1 &

2、数据下载

使用 kingfisher 一行下载，kingfisher详细介绍见帖子：一行代码下载原始数据—Kingfisher

# 我的下载目录 /nas1/zhangj/project/PRJNA445209/data/rawdata
nohup kingfisher get -p PRJNA445209 -m ena-ascp ena-ftp prefetch aws-http 1>down_PRJNA445209.log 2>&1 &

这个数据应该上传的是cleandata，别问为啥知道，问就是一眼看出来了文件的命名方式是以前工作过的公司的风格：在ENA下载的文件表格 filereport_read_run_PRJNA445209_tsv.txt，我们可以使用fastqc简单看一眼：

nohup fastqc -t 6 -o ./ SRR*.fastq.gz 1>qc.log  2>&1 &
multiqc *.zip -o ./ -n qc_fastqc

3、数据比对：

# 目录/nas1/zhangj/project/PRJNA445209/Mapping/
mkdir Mapping/
cd Mapping/
## ----2.比对
# 进入比对文件夹
cd /nas1/zhangj/project/PRJNA445209/Mapping

## 已经构建好的索引，版本为 Hongyang/v3.0
## 定义索引前缀index
index=/nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome

## 使用echo生成 hisat.sh
ls ../data/rawdata/*_1.fastq.gz | while read id
do
  base_name="${id##*/}" # 去掉路径，只保留文件名
  name="${base_name%_1.fastq.gz}"  # 去掉后缀之前的部分
  echo "hisat2 -p 5 -x ${index} -1 ../data/rawdata/${name}_1.fastq.gz -2 ../data/rawdata/${name}_2.fastq.gz 2>${name}.log |samtools sort -@ 3 -o ${name}.Hisat_aln.sorted.bam - "
done >hisat.sh

hisat.sh内容如下：

hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885590_1.fastq.gz -2 ../data/rawdata/SRR6885590_2.fastq.gz 2>SRR6885590.log |samtools sort -@ 3 -o SRR6885590.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885591_1.fastq.gz -2 ../data/rawdata/SRR6885591_2.fastq.gz 2>SRR6885591.log |samtools sort -@ 3 -o SRR6885591.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885592_1.fastq.gz -2 ../data/rawdata/SRR6885592_2.fastq.gz 2>SRR6885592.log |samtools sort -@ 3 -o SRR6885592.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885593_1.fastq.gz -2 ../data/rawdata/SRR6885593_2.fastq.gz 2>SRR6885593.log |samtools sort -@ 3 -o SRR6885593.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885594_1.fastq.gz -2 ../data/rawdata/SRR6885594_2.fastq.gz 2>SRR6885594.log |samtools sort -@ 3 -o SRR6885594.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885595_1.fastq.gz -2 ../data/rawdata/SRR6885595_2.fastq.gz 2>SRR6885595.log |samtools sort -@ 3 -o SRR6885595.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885596_1.fastq.gz -2 ../data/rawdata/SRR6885596_2.fastq.gz 2>SRR6885596.log |samtools sort -@ 3 -o SRR6885596.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885597_1.fastq.gz -2 ../data/rawdata/SRR6885597_2.fastq.gz 2>SRR6885597.log |samtools sort -@ 3 -o SRR6885597.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885598_1.fastq.gz -2 ../data/rawdata/SRR6885598_2.fastq.gz 2>SRR6885598.log |samtools sort -@ 3 -o SRR6885598.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885599_1.fastq.gz -2 ../data/rawdata/SRR6885599_2.fastq.gz 2>SRR6885599.log |samtools sort -@ 3 -o SRR6885599.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885600_1.fastq.gz -2 ../data/rawdata/SRR6885600_2.fastq.gz 2>SRR6885600.log |samtools sort -@ 3 -o SRR6885600.Hisat_aln.sorted.bam -
hisat2 -p 5 -x /nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_genome -1 ../data/rawdata/SRR6885601_1.fastq.gz -2 ../data/rawdata/SRR6885601_2.fastq.gz 2>SRR6885601.log |samtools sort -@ 3 -o SRR6885601.Hisat_aln.sorted.bam -

并行投递：

# 运行，一次投递6个样本
# -CPU：一次投递几个任务
nohup ParaFly -c hisat.sh -CPU 6 1>hisat.log  2>&1 &

第二天看看：

每个样本的总比对率，没有那么差：

4、featureCounts 定量

这里需要用到gtf格式的注释文件，我们下载的只有gff3，需要转一下：

conda activate rna
conda install gffread -y # 安装软件

# 转换
gffread Hongyang_v3.0_update.gff3 -T -o Hongyang_v3.0.gtf

万事具备：

mkdir /nas1/zhangj/project/PRJNA445209/Expression/
cd /nas1/zhangj/project/PRJNA445209/Expression/

## 定义输入输出文件夹
gtf=/nas1/zhangj/database/genome/Hongyang/v3.0/Hongyang_v3.0.gtf

# featureCounts对bam文件进行计数 v2.0.3以上版本
nohup featureCounts -T 30 -p -t exon -g gene_id -a $gtf -o all.id.txt ../Mapping/*.sorted.bam 1>count.log 2>&1 &

# 对定量结果质控
multiqc ./all.id.txt.summary -n qc_count

# 得到表达矩阵
# 处理表头，
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#../Mapping/##g' |sed 's#.Hisat_aln.sorted.bam##g' > raw_counts.txt