专栏首页生信技能树关键问题答疑:WGCNA的输入矩阵到底是什么格式

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

  • 请问用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

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-14

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 《GEO数据挖掘课程》配套练习题

    我把3年前的收费视频课程:3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月,免费到B站:

    生信技能树
  • 有趣的基因命名

    还有,如果你看到HS.开头的基因,它是unigene的ID了,已经不再是symbol啦。

    生信技能树
  • 使用Mfuzz包做时间序列分析

    既然是讲解时间序列分析,那么就不得不提一下Mfuzz包了,恰好生信技能树创始人jimmy的200篇生物信息学文献阅读活动分享过的一篇文章就有这个,作者主要使用了...

    生信技能树
  • 【Rust 日报】2020-03-25 regex crate 的计划

    测试驱动的 Rust 学习项目,适合有其他语言编程经验的 Rust 新手. 在这个项目中,你可以通过一系列测试驱动的练习以及阅读材料来学习如何构建一个 JIRA...

    MikeLoveRust
  • 实战 | 一行代码让你的电脑可以看图说话

    Image Caption 任务是一个需要综合计算机视觉和自然语言处理的任务,需要使用计算机建立某种映射方式,将处于视觉模态当中的数据映射到文本模态当中,即让视...

    潘永斌
  • 机器学习 TOP 10 必读论文 | 资源

    用户1737318
  • 干货!机器学习 必读论文 TOP 10

    Medium上的机器学习深度爱好者必关注的账号Mybridge照例对11月发表的学术论文进行了排名,整理出了10篇必读论文,建议收藏深读。 1. Alpha Z...

    企鹅号小编
  • 可视化算法VxOrd论文研读

    摘要 本文介绍了一种适合挖掘超大型数据库的聚类和排序ordination算法,包括微阵列表达式研究microarray expression studies产生...

    ZONGLYN
  • 泛微 e-cology OA 远程代码执行漏洞复现

    Poc已在github公开,由于环境搭建较为复杂,所以我在空间搜索引擎中找了国外的网站进行复现

    用户1631416
  • 黑客法律知识课堂

    大部分学网络安全的人,都有个致命的缺点,就是不懂这方面的法律知识,装逼挂个黑页,等下就被警察叔叔发了通告,身边也不少例子,今天就帮大家补一下这方面的知识

    周俊辉

扫码关注云+社区

领取腾讯云代金券