前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >关键问题答疑:WGCNA的输入矩阵到底是什么格式

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

作者头像
生信技能树
发布2019-10-15 16:28:21
3.3K0
发布2019-10-15 16:28:21
举报
文章被收录于专栏:生信技能树
  • 请问用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种方法,全部代码如下:

代码语言:javascript
复制
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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-10-14,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 首先TCGA目前的确是以count格式的矩阵下载为主
  • 然后问题就是,用tcga做wgcna分析,是不是原始数据输入一定要是tpm和fpkm格式?
  • 接着什么样的程序一定要tpm和fpkm格式呢?
  • 最后如果下的raw_count有r包转换为tpm和fpkm
相关产品与服务
GPU 云服务器
GPU 云服务器(Cloud GPU Service,GPU)是提供 GPU 算力的弹性计算服务,具有超强的并行计算能力,作为 IaaS 层的尖兵利器,服务于生成式AI,自动驾驶,深度学习训练、科学计算、图形图像处理、视频编解码等场景。腾讯云随时提供触手可得的算力,有效缓解您的计算压力,提升业务效率与竞争力。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档