ixxmu
commented
1 year ago 以Counts转TPM为例, 感受一下生信自学中阳性对照的重要性 by 果子学生信
把转录组的数据从counts转为TPM, FRPM一直是个小需求，就是没有那么重要，但是总有人需要。
恰好本次改版的TCGA转录组同时提供了counts，TPM，FRPM格式的数据，可以作为一个阳性对照来测试格式转换的流程了。
阳性对照有多重要
获取数据分析的阳性对照，是自我迭代的神器。你不再需要询问别人，这个做的对不对，那个做的对不对，你自己就可以判断。
我最怕别人问我这个代码报错了怎么办，因为代码报错千千万，能解决它的是能力，不是运气，而能力的获得是靠自己动手的。
我通常的回复是这样的(下次我就直接复制粘贴了):
第一，找到这个函数的R包，去到R包的文档那里，运行一下里面阳性的例子，确认你的 R包安装完成，以及代码书写是否正确。
第二，排查了代码的问题，那现在大概率是数据的问题，看看别人在该步骤输入的数据是什么格式，class看看属性，str看看结构，把自己的数据调整一下
第三，还解决不了，设置报错为英文，提取关键英文信息来检索，常见的问题都有答案
第四，还解决不了的就真的靠能力了，找到该函数的在github上面的原始代码，设置对应的参数，一行行运行内部的代码找到报错信息所在的位置，查找原因。
基本上前三个就能解决95%的问题了，最后一个比较考验个人能力，我习得这个技能后，报错都能自己解决了。
一个使用的例子在这里, 我找到了报错的代码修改了之后，重新打包，这样R包的所有问题都能解决了
出图神器export目前不能用了，该如何是好？
这个方法也被我用来拆解R包，掌握算法的原理：
ssgsea算法在量化免疫浸润时的运用以及原理
然后这个方法虽好，唯一的缺点大概就是耗时间，但是我比较推荐去尝试，未来我也会在答疑群给大家多演示几次。
言归正传，下面开始数据的转换。
推荐资源
首先是原理，推荐一个视频，一个网页，一本书。
进入Statquest网页，https://statquest.org/
在视频列表里面找到这个视频看一下
如果你的网络不方便，可以在哔哩哔哩检索statquest 也有好心人搬运了。
然后还要知道，这些格式之间的区别，看看这个网页
https://hbctraining.github.io/Training-modules/planning_successful_rnaseq/lessons/sample_level_QC.html
然后还有一本推荐给学员的书籍
我们使用的代码在这个章节，作者提供了测试数据，可以下载后完美运行
https://compgenomr.github.io/book/gene-expression-analysis-using-high-throughput-sequencing-technologies.html
只不过作者例子中，提供了一个基因长度，而这个关键数据，在自己的数据上是没有的，而这个数据是格式转换中最重要的部分
获取基因长度(最重要)
基因长度，指的是外显子的总和，这个可以从gtf文件中获取。
GTF文件有什么用啊？别的不谈，最起码能提lncRNA
既往的TCGA数据，使用的gtf版本是V22，xena网站用的是V23，改变后用的是V36，他们的差别在于对基因的认识，
比如哪些是非编码RNA，哪些本来是非编码RNA，现在变成了mRNA。
在这个页面，我们能看到本次使用的就是V36版本，而且提供下载
https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files
下载后解压，会变的很大，通过以下函数可以获取基因长度
library(GenomicFeatures)
## 1.读取GTF文件
txdb <- makeTxDbFromGFF("resource/gencode.v36.annotation.gtf")
## 2.提取外显子信息
exonic <- exonsBy(txdb, by="gene")
## 3.外显子长度求和
exonic.gene.sizes <- sum(width(reduce(exonic)))
## 4.结果整理
mygeneids <- data.frame(gene_id=names(exonic.gene.sizes),width=exonic.gene.sizes)
最终的结果如下，有一点要强调的是，这里的基因数目是60660，这跟我们之前下载的数据是一致的。
如果不一致怎么办，那就是没弄对，会导致后面所有的结果出现偏差。
需要注意的是，外显子求和步骤，并不是把所有外显子相加，而是经过reduce函数处理，重复的部分并没有重复计算，这是核心点。
而这里的代码基本上都是bioconductor的技能，这里挖一个坑，以后来填。
测试流程是否正确
怎么判断这样的提取过程对不对呢，我们可以来转换数据试试。
先提取一个数据来试一下, 这个数据就是我们之前下载一个数据
TCGA改版后转录组数据的下载以及整理
data <- data.table::fread("data_in_one/00063026-2a67-41de-b105-f80dca277978.rna_seq.augmented_star_gene_counts.tsv")
data <- data[-seq(1,4),c(1,2,3,4,7,8)]
把长度合并一下，便于计算和观察
data = merge(data,mygeneids,by="gene_id")
从counts到TPM
第一，矫正基因长度,这里的apply代码是可以批量处理多个列的，尽管这里只有一列。
geneLengths <- data$width
rpk <- apply( X = subset(data, select = "unstranded"),
MARGIN =2,
FUN = function(x){
x/(geneLengths/1000)
})
第二，在现有基础上，测序的深度
tpm <- apply(X = rpk,
MARGIN = 2,
FUN = function(x){
x / sum(as.numeric(x)) * 10^6
})
第三，合并数据，查看
newdata = cbind(data,round(tpm,4))
colnames(newdata)[8] = "mytpm"
现在我们确信，基因的长度是对的，转换的方法也是对的，这就叫自我迭代。
从counts到FPKM/RPKM
这个过程也是两步，但是因为他矫正测序深度，直接用的counts，或者说他是先矫正测序深度，再来矫正的基因长度
fpkm <- apply(X = subset(data, select = "unstranded"),
MARGIN = 2,
FUN = function(x) {
10^9 * x / geneLengths / sum(as.numeric(x))
})
合并后查看
newdata = cbind(data,round(fpkm,4))
colnames(newdata)[8] = "myfpkm"
几乎一样，有一点差别，但是因为TPM的成功，我们也相信这个结果是对的。
实战
说到这里，这个技能算是掌握了，下面用一个例子来展示一下实际操作。
就以之前批量下载数据中的 BRCA为例
拿去！TCGA改版后转录组数据的批量下载及分享
### 1.读取数据
data <- readRDS(file = "TCGA-BRCA.rds")
### 2.查看数据
test <- data[1:100,1:100]
### 3.获取基因长度(可以保留一份)
mygeneids <- data.frame(gene_id=names(exonic.gene.sizes),width=exonic.gene.sizes)
### 4.强制限定基因顺序，跟长度的顺序一致
rownames(data) = mygeneids$gene_id
### 转换开始
geneLengths <- mygeneids$width
rpk <- apply( X = data,
MARGIN =2,
FUN = function(x){
x/(geneLengths/1000)
})
tpm <- apply(X = rpk,
MARGIN = 2,
FUN = function(x){
x / sum(as.numeric(x)) * 10^6
})
很快就结束了。
这时候我们发现，基因长度是很重要的。而这个文件是固定的，也就是这个版本的是所有TCGA数据都可以用这个文件。
我把这个mygeneids文件，连同下载的gtf文件，以及之前的那个statquest视频一起打包了。
回复"果子转录组格式转换"自助获取。
说多了都是泪，但凡各位网络好一点，这个都不需要我自己上传。
那么既然有了gtf文件，我们就可以提取非编码RNA，看看哪些多出来了，并且之前的数据都没有采用这个新版本的gtf来比对，那么此时基因的表达量肯定有变化，会更加准确。我有一种感觉，TCGA的非编码RNA文章，流程都可以重新尝试一下了。这可能就是当前全民生信的局面下，不多的挖掘机会吧。
https://mp.weixin.qq.com/s/imJ5za173pKFJVad5j_5iw