关键问题答疑：WGCNA的输入矩阵到底是什么格式

虽然我们全国巡讲课程并不讲解WGCAN内容，因为时间的确有限，短短的3天要传授给大家R语言，linux还有RNA-seq数据分析实战，希望给大家打造好的基础成为合格的生信工程师，但是我们公众号有数不胜数的高级分析教程，比如WGCNA，有了基础的大家看教程就容易很多。今天收到生信技能树201908北京站学员提问，问题描述是：

  • 请问用tcga做wgcna分析，原始数据输入tpm和fpkm格式都行吗？

  • 如果下的raw_count有r包转换吗？

这样的问题我其实被问过好多次了，因为这次是学员提问，虽然已经过了一个月的答疑期，但是情谊还在，所以就系统性的回复一下。

首先TCGA目前的确是以count格式的矩阵下载为主

至于能不能找到RPKM这样的矩阵，肯定是可以的，但是我教大家的主要是count值，因为对RNA-seq数据的差异分析以这个count为input，大家可以看我B站的教学视频

然后问题就是，用tcga做wgcna分析，是不是原始数据输入一定要是tpm和fpkm格式？

(PS, 类似的基因表达量的归一化还有很多，详细见：https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html)

那么问题就是，用tcga做wgcna分析，是不是原始数据输入一定要是tpm和fpkm格式？

其实呢，我最开始的教程，的确是fpkm，所以大家会以为必须要这样的输入格式，详细教程见：一文看懂WGCNA 分析(2019更新版)

实际上，WGCNA首先会对全部基因的表达量计算两两之间的相关性，这个时候，只需要基因的表达量是适合计算相关性的即可，如果是 原始 counts值，可以直接转为 log(cpm+1) 的格式 ，更为重要的其实是挑选多少个基因进入后续的wgcna流程。

以及我们的基因被WGCNA算法分成了不同模块后，哪些是有生物学意义的，跟表型相关性。

接着什么样的程序一定要tpm和fpkm格式呢？

类似tpm和fpkm的基因表达量的归一化还有很多，详细见：https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html  如果是需要对基因表达量进行排序，这个时候，基因长度就有影响，所以需要使用tpm和fpkm，比如：http://xcell.ucsf.edu/

The expression matrix should be a matrix with genes in rows and samples in columns. The rownames should be gene symbols. If the data contains non-unique gene symbols, rows with same gene symbols will be averaged. xCell uses the expression levels ranking and not the actual values, thus normalization does not have an effect, however normalizing to gene length (RPKM/FPKM/TPM/RSEM) is required.

最后如果下的raw_count有r包转换为tpm和fpkm

其实我GitHub有代码的，而且我还提出了3种方法，全部代码如下:

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)
## 载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
# 注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。

library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
?TxDb
## 下面是定义基因长度为 非冗余exon长度之和
if(F){
  exon_txdb=exons(txdb)
  genes_txdb=genes(txdb)
  genes_txdb
  ?GRanges
  o = findOverlaps(exon_txdb,genes_txdb)
  o
  ## exon - 1 : chr1 4807893-4807982
  ## 1        6523
  #  genes_txdb[6523]  # chr1 4807893-4846735 , 18777
  t1=exon_txdb[queryHits(o)]
  t2=genes_txdb[subjectHits(o)]
  t1=as.data.frame(t1)
  t1$geneid=mcols(t2)[,1]
  # 如果觉得速度不够，就参考R语言实现并行计算
  # http://www.bio-info-trainee.com/956.html
  #lapply : 遍历列表向量内的每个元素，并且使用指定函数来对其元素进行处理。返回列表向量。
  #函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组;
  #它的返回值是一个列表，代表分组变量每个水平的观测。
  g_l = lapply(split(t1,t1$geneid),function(x){
    # x=split(t1,t1$geneid)[[1]]
    head(x)
    tmp=apply(x,1,function(y){
      y[2]:y[3]
    })
    length(unique(unlist(tmp)))
    # sum(x[,4])
  })
  head(g_l)
  g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))

save(g_l,file = 'step7-g_l.Rdata')
}
load(file = 'step7-g_l.Rdata')
## 下面是定义基因长度为 最长转录本长度
if(F){

t_l=transcriptLengths(txdb)
  head(t_l)
  t_l=na.omit(t_l)
  head(t_l)
  t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
  head(t_l)
  str(t_l)
  t_l=t_l[!duplicated(t_l$gene_id),]
  head(t_l)
  g_l=t_l[,c(3,5)]

}

head(g_l)
library(org.Mm.eg.db)
s2g=toTable(org.Mm.egSYMBOL)
head(s2g)
g_l=merge(g_l,s2g,by='gene_id') #把g_l,s2g两个数据框以'gene_id'为连接进行拼接

#merge函数可以实现对两个数据表进行匹配和拼接的功能。

# 参考counts2rpkm，定义基因长度为非冗余CDS之和
# http://www.bio-info-trainee.com/3298.html  
a[1:4,1:4]
ng=intersect(rownames(a),g_l$symbol) #取a数据框的行名与g_l数据框的symbol列的交集
#intersect()取交集

# 有了counts矩阵和对应的基因长度信息，就很容易进行各种计算了：
exprSet=a[ng,]
lengths=g_l[match(ng,g_l$symbol),2]
head(lengths)
head(rownames(exprSet))
# http://www.biotrainee.com/thread-1791-1-1.html
exprSet[1:4,1:4]
total_count<- colSums(exprSet)
head(total_count)
head(lengths)
total_count[4]
lengths[1]
1*10^9/(1122*121297)
rpkm <- t(do.call( rbind,
                   lapply(1:length(total_count),
                          function(i){
  10^9*exprSet[,i]/lengths/total_count[i]
}) ))

rpkm[1:4,1:4]
# 下面可以比较一下 自己根据counts值算出来的RPKM和作者提供的RPKM区别。
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',
             header = T ,sep = '\t')
# 每次都要检测数据
a[1:4,1:4]
rpkm_paper=a[ng,] 
rpkm_paper[1:4,1:4]

rpkm[1:4,1:4]

上面的代码有点复杂，如果R语言水平不够，不建议去理解了。其它知识点代码是：https://github.com/jmzeng1314/scRNA_smart_seq2

1

南京场（正在进行）

10.12-10.14

2

南宁场（马上开始）

10.26-10.28

课程内容

1

生信R语言入门

2

GEO数据库挖掘

5

生信-Linux基础

6

转录组课题设计与流程分析
赞 (0)

相关推荐