怀疑了不该怀疑的人

生信技能树

发布于 2023-02-27 21:40:25

4020

发布于 2023-02-27 21:40:25

前面介绍了：一行命令将count转为CPM/TPM/FPKM ，是一个Python软件，也是读取gtf文件，根据id来自动计算每个基因的长度后进行相对应的公式的转换：

rnanorm sample.count.tsv --annotation gencode.vM25.annotation.gtf \
    --cpm-output sample.count.cpm.tsv
    --tpm-output sample.count.tpm.tsv \
    --fpkm-output sample.count.fpkm.tsv \

位置参数为基因count文件
--annotation 基因组注释GTF文件
--cpm-output CPM输出文件
--tpm-output TPM输出文件
--fpkm-output FPKM输出文件

但是我们公众号读者可能是更熟悉R代码，恰好看到了我们的马拉松授课王牌讲师小洁也整理了她的转换方法，借花献佛送给大家哈。

0.起因

想要重新整理一下count转换为tpm的代码。隐约记得之前有反馈说计算出来的结果和哪里不一样，所以我想检查一下。

1.基因长度

之前写过：基因长度并不是end-start，有4种计算方式，其中非冗余外显子长度之和是更推荐的。

2.非冗余外显子长度之和的计算方法

找到了两种方法，曾老板的代码是我之前一直在用的。

2.1 曾老板的代码

下载genecodev36版本的gtf文件，即新版TCGA数据使用的参考基因组注释文件

if(!require(rtracklayer))BiocManager::install("rtracklayer")
library(rtracklayer)
# 读取文件是不需要解压的。
gtf = rtracklayer::import("gencode.v36.annotation.gtf.gz")
gtf = as.data.frame(gtf)
exon = gtf[gtf$type=="exon",
           c("start","end","gene_id")]
length(unique(exon$gene_id)) # 有多少个基因
## [1] 60660
f = "gle.Rdata"
if(!file.exists(f)){
  gle = lapply(split(exon,exon$gene_id),function(x){
    tmp=apply(x,1,function(y){
      y[1]:y[2]
    })
    length(unique(unlist(tmp)))
  })
  gle=data.frame(gene_id=names(gle),
                 length=as.numeric(gle))
  save(gle,file = f)
}
load(f)

head(gle)
##              gene_id length
## 1 ENSG00000000003.15   4536
## 2  ENSG00000000005.6   1476
## 3 ENSG00000000419.13   1207
## 4 ENSG00000000457.14   6883
## 5 ENSG00000000460.17   5970
## 6 ENSG00000000938.13   3382

2.2一个是网上搜到的方法

使用R包GenomicFeatures

f2 = "gfe.Rdata"
if(!file.exists(f2)){
    # First, import the GTF-file that you have also used as input for htseq-count
  if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
  library(GenomicFeatures)
  txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",format="gtf")
  # then collect the exons per gene id
  exons.list.per.gene <- exonsBy(txdb,by="gene")
  # then for each gene, reduce all the exons to a set of non overlapping exons, calculate their lengths (widths) and sum then
  exonic.gene.sizes <- sum(width(GenomicRanges::reduce(exons.list.per.gene)))
  
  gfe <- data.frame(gene_id=names(exonic.gene.sizes),
                    length=exonic.gene.sizes)
  save(gfe,file = f2)
}
load(f2)
head(gfe)
##                               gene_id length
## ENSG00000000003.15 ENSG00000000003.15   4536
## ENSG00000000005.6   ENSG00000000005.6   1476
## ENSG00000000419.13 ENSG00000000419.13   1207
## ENSG00000000457.14 ENSG00000000457.14   6883
## ENSG00000000460.17 ENSG00000000460.17   5970
## ENSG00000000938.13 ENSG00000000938.13   3382

看两个结果的前几行，是完全一致的。

但是全部六万个基因比较起来，居然有6个不一样的。。。

k = gle$length!=gfe$length;table(k)
## k
## FALSE  TRUE 
## 60654     6
m = cbind(gle[k,],gfe[k,])
m
##                 gene_id length           gene_id length
## 12605 ENSG00000169488.6   2432 ENSG00000169488.6   1216
## 14181 ENSG00000176925.8   2272 ENSG00000176925.8   1136
## 14862 ENSG00000180433.5   2344 ENSG00000180433.5   1172
## 17035 ENSG00000196248.5   2242 ENSG00000196248.5   1121
## 19788 ENSG00000204246.3   2216 ENSG00000204246.3   1108
## 53388 ENSG00000275553.4   2436 ENSG00000275553.4   1218

有趣，刚好是两倍的关系啊。那么问题来了，谁做的对呢？

错的那个，怎么就万分之一的基因错了呢？

3.对答案

TCGA提供了tpm的，这个很权威，不太可能会出错。

CHOL_rnaseq.Rdata里的数据是TCGAbiolinks下载整理好的结果，里面包含了count和tpm。

load("CHOL_rnaseq.Rdata") 
library(SummarizedExperiment)
exp_count = assay(dat_chol)
exp_tpm = assay(dat_chol,4)
countToTpm <- function(counts, effLen)
{
  rate <- log(counts) - log(effLen)
  denom <- log(sum(exp(rate)))
  exp(rate - denom + log(1e6))
}

tpm1 <- apply(exp_count,2,countToTpm,gle$length)
tpm2 <- apply(exp_count,2,countToTpm,gfe$length)
# 把答案不一样的六个基因提取出来。
tpmaa = exp_tpm[k,]
which(colSums(tpmaa)!=0)

## TCGA-4G-AAZO-01A-12R-A41I-07 TCGA-YR-A95A-01A-12R-A41I-07 
##                           25                           33

re = data.frame(ans = exp_tpm[k,33],
                le = tpm1[k,33],
                fe = tpm2[k,33])
re

##                      ans         le         fe
## ENSG00000169488.6 0.0000 0.00000000 0.00000000
## ENSG00000176925.8 0.0497 0.02484599 0.04969198
## ENSG00000180433.5 0.0482 0.02408280 0.04816560
## ENSG00000196248.5 0.0000 0.00000000 0.00000000
## ENSG00000204246.3 0.0000 0.00000000 0.00000000
## ENSG00000275553.4 0.0000 0.00000000 0.00000000

在这个结果表格里，第一列是答案，第二列是曾老板代码的计算结果，第三列是R包的计算结果，很明显是第三列和答案一样啊，第二列不一样。。。

好了，现在我得出了一个不成熟的结论：曾老板的代码有bug。

正当我打算发个邮件告诉他的时候，转念一想，不对啊。曾老板的代码什么时候有bug?我是不是冤枉他了。

不行，发现了问题得完整解决，我还是继续深究一下好了。把那六个出错的基因提取出来

miniexon = exon[exon$gene_id %in% rownames(re),]
miniexon
##             start       end           gene_id
## 190671  158754720 158755891 ENSG00000180433.5
## 190678  158754720 158755891 ENSG00000180433.5
## 1403912 104535749 104536856 ENSG00000204246.3
## 1403919 104535749 104536856 ENSG00000204246.3
## 1588707   4821321   4822456 ENSG00000176925.8
## 1588714   4821321   4822456 ENSG00000176925.8
## 1740534 123976661 123977781 ENSG00000196248.5
## 1740541 123976661 123977781 ENSG00000196248.5
## 1972167  19975444  19976659 ENSG00000169488.6
## 1972174  19975444  19976659 ENSG00000169488.6
## 2540391  46968695  46969912 ENSG00000275553.4
## 2540396  46968695  46969912 ENSG00000275553.4

好家伙，gtf文件里居然有重复的行啊！

所以其实不是代码有问题，是输入文件有问题！问题已经得到解决，只要用distinct函数给整个exon数据框去个重复，重新计算就行了。

4. 属于R语言高级玩家的快乐

必须要研究一下曾老板高难度的R语言代码，看看是个什么原理，为什么碰上这种问题数据会算出错误结果。

找两个基因出来，一个是计算正确的，一个是计算错误（有重复行的）的，解析一下循环里的代码就可以了。

出了新手村的R语言玩家可以先学习一下apply和lapply这两个函数，查帮助文档即可。

g = c("ENSG00000000003.15","ENSG00000169488.6")
miniexon = exon[exon$gene_id %in% g,]
a = split(miniexon,miniexon$gene_id)
lapply(a, head)

## $ENSG00000000003.15
##             start       end            gene_id
## 2964638 100636608 100636806 ENSG00000000003.15
## 2964641 100635558 100635746 ENSG00000000003.15
## 2964643 100635178 100635252 ENSG00000000003.15
## 2964645 100633931 100634029 ENSG00000000003.15
## 2964647 100633405 100633539 ENSG00000000003.15
## 2964649 100632485 100632568 ENSG00000000003.15
## 
## $ENSG00000169488.6
##            start      end           gene_id
## 1972167 19975444 19976659 ENSG00000169488.6
## 1972174 19975444 19976659 ENSG00000169488.6
#把正常的基因代入进lapply循环里的代码：
tmp=apply(a[[1]],1,function(y){
    y[1]:y[2]
})
lapply(tmp[1:3], head)

## $`2964638`
## [1] 100636608 100636609 100636610 100636611 100636612
## [6] 100636613
## 
## $`2964641`
## [1] 100635558 100635559 100635560 100635561 100635562
## [6] 100635563
## 
## $`2964643`
## [1] 100635178 100635179 100635180 100635181 100635182
## [6] 100635183

#结果是一个列表，每个元素是一串连续的整数组成的向量，代表每个外显子的位置。
#因为重复序列不重复算长度，所以用下面的代码拆列表，unique即可
length(unique(unlist(tmp)))

## [1] 4536

把不正常的基因代进去：

tmp=apply(a[[2]],1,function(y){
    y[1]:y[2]
})
head(tmp)

##       1972167  1972174
## [1,] 19975444 19975444
## [2,] 19975445 19975445
## [3,] 19975446 19975446
## [4,] 19975447 19975447
## [5,] 19975448 19975448
## [6,] 19975449 19975449

#因为不正常基因的两行输入数据是一样的，:生成的而向量长度相同，导致结果不是列表了，是个只有两列，且两列完全一致的矩阵啊！
#所以下面的代码就不能再拆列表成为向量，然后去重复了。
length(unique(unlist(tmp)))

## [1] 2432

如果不给exon数据框去重，要从代码的角度解决它，那就把上面的代码换成：

length(unique(as.integer(unlist(tmp))))

## [1] 1216

这样，不管是正常的基因还是不正常的基因，就都能搞定了。

代码不能适应问题数据是正常的，因为写代码的时候没碰上这种问题数据啊！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2023-01-09，如有侵权请联系 cloudcommunity@tencent.com 删除

黑石物理服务器

编程算法

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

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

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

黑石物理服务器

编程算法

登录后参与评论

0 条评论

热度