count格式的数据转换(count to FPKM,count to TPM) 【GEO数据库】

news2024/9/23 15:28:29

在正式分析之前，对于数据的处理是至关重要的，这种重要性是体现在很多方面，其中有一点是要求分析者采用正确的数据类型。

对于芯片数据，原始数据进行log2处理之后可以进行很多常见的分析，比如差异分析、热图、箱线图、PCA分析、生存分析、模型构建，聚类分析和相关性分析等。

对于转录组数据，在上述的常见分析中只有差异分析时需要采用count值，其他的分析是需要采用log2后的cpm，tpm，fpkm，rpkm数据。

count数据转换为cpm数据非常简单
# exprSet是count表达矩阵
# 一句代码搞定
exprSet = log2(edgeR::cpm(exprSet)+1)

比较难的就是count数据转换为tpm数据，因此搬运了常规的流程和R包的方法，做个对比

首先要去获取基因长度文件，因为后续需要用这个数据去矫正基因长度。

网址：https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files

1.Count值

对给定的基因组参考区域，计算比对上的read数，又称为raw count（RC）。在RNAseq数据中，raw reads count一般是指mapped到基因外显子区域的reads数目。

2.RPKM/FPKM

RPKM

RPKM(reads per kilobase per million),对测序深度和长度做了一个矫正。用于单端测序数据。

RPKM = (gene_read_count / (total_count*gene_length) * 10^6 * 10^3

基因count除以总read count，除以gene长度。然后避免数值太小，根据测序深度大小乘10^6, 基因长度就乘以一个10^3。

FPKM

RPKM(fragments per kilobase per million) 用于双端数据，一个fragment是一对reads。

故FPKM = RPKM / 2

FPKM(Fragment Per Kilobase of transcript, per Million mapped reads）：每千碱基片段每百万映射读取的 reads 数)，是针对双端测序的一个normalization方法。通常来讲，当paired reads同时匹配到一个位置，记为fragment(注：即便是双端测序，RPKM也不完全是FPKM的2倍)。

RPKM: Reads Per Kilobase of exon model per Million mapped reads (每千个碱基的转录每百万映射读取的reads)

RPKM/FPKM方法：10^3标准化了基因长度的影响，10^6标准化了测序深度的影响。FPKM方法与RPKM类似，主要针对双末端RNA-seq实验的转录本定量。在双末端RNA-seq实验中，有左右两个对应的read来自相同的DNA片段。在进行双末端read进行比对时，来自同一DNA片段的高质量的一对或单个read可以定位到参考序列上。为避免混淆或多次计数，统计一对或单个read比对上的参考序列片段（Fragment），来计算FPKM，计算方法同RPKM。

RPKM与FPKM的区别：RPKM值适用于单末端RNA-seq实验数据，FPKM适用于双末端RNA-seq测序数据。

3.TPM (Transcript per million)

TPM

TPM(transcript per million), 基因FPKM 占总的FPKM的比例， TPM衡量基因在样本中的相对表达量，对同一个基因不同样本，其FPKM可能相同，但相对所有基因而言，其在不同样本中可能所占比例不同。

TPM = gene_fpkm / total_fpkm * 10^6

TPM（Transcripts Per Million）是一种常用的基因表达量归一化方法，它将基因的表达量调整为每百万条转录本的数量。TPM 值考虑了基因的长度和测序深度，通过将每个基因的 Counts 值除以其长度，并进行适当的归一化，将基因的表达量转换为每百万转录本数，以便进行样本间的比较和分析。TPM 值消除了样本间测序深度的差异和基因长度的影响。

TPM的计算方法也同RPKM/FPKM类似，首先使用式2计算每个基因的表达值，去除基因长度的影响。随后计算每个基因的表达量的百分比，最后再乘以10^6，TPM可以看作是RPKM/FPKM值的百分比。

直接说事情，我有一个基因A，它在这个样本的转录组数据中被测序而且mapping到基因组了 5000个的reads，而这个基因A长度是10K，我们总测序文库是50M，所以这个基因A的RPKM值是 5000除以10，再除以50，为10. 就是把基因的reads数量根据基因长度和样本测序文库来normalization 。那么它的TPM值是多少呢？这个时候这些信息已经不够了，需要知道该样本其它基因的RPKM值是多少，加上该样本有3个基因，另外两个基因的RPKM值是5和35，那么我们的基因A的RPKM值为10需要换算成TPM值就是 1,000,000 *10/(5+10+35)=200,000，看起来是不是有点大呀，其实主要是因为我们假设的基因太少了，一般个体里面都有两万多个基因的，总和会大大的增加，这样TPM值跟RPKM值差别不会这么恐怖的。

TPM与RPKM/FPKM的区别：从计算公式来说，唯一的不同是计算操作的顺序，TPM是先去除了基因长度的影响，而RPKM/FPKM是先去除测序深度的影响，具体可看这篇博文，有计算步骤的详细说明；TPM实际上改进了RPKM/FPKM方法在跨样品间定量的不准确性。

TPM的使用范围与RPKM/FPKM相同。

#设置工作目录左边/，右边\\

setwd("D:\\R_Script")

#安装↓这个包#install.packages("tidyverse")
library(tidyverse)

#读gtf文件，计算所有外显子的长度

gtf <- read_tsv("human.gtf", comment="#", col_names=c('chr','source','type','start','end','score','strand','phase','attributes')) %>% filter(type=='exon') %>% mutate(len = end - start + 1) %>% select(start, end, attributes,len)

#计算基因的非冗余外显子的长度，获得有效基因长度

gtf$attributes %>% str_extract(., "gene_id \"[\\w|\\.]+") %>% str_remove(., "gene_id \"") -> gtf$gene_id

gtf %>% select(start, end, gene_id, len) %>% distinct(start,end,gene_id, .keep_all = T) %>% select(gene_id,len) %>% group_by(gene_id) %>% summarise(est_len=sum(len)) -> gtf

#储存基因长度文件备用

write.csv(gtf, "gene_length.csv", row.names = TRUE)

#读取文件：1.自己的count.txt矩阵；2.上面生成的基因长度

rt <-read.table("GSEXXXXXX.txt", row.names = 1,header = TRUE,sep="\t")

str(rt)

eff_length <- read.csv("gene_length.csv", row.names = 1, header = T)

rownames(eff_length)<-eff_length$gene_id

rownames(eff_length) <- do.call(rbind,strsplit(as.character(eff_length$gene_id),'\\.'))[,1]eff_length[1:3,]

#读取基因长度文件

gen<- intersect(rownames(rt), rownames(eff_length))

rt<- rt[gen,]

eff_length<- eff_length[gen, ]

#定义函数

countToFpkm <- function(counts, effLen) { N <- sum(counts) exp(log(counts) + log(1e9) - log(effLen) - log(N)) }##count转换为FPKM值

fpkms <- as.data.frame(apply(rt, 2, countToFpkm, effLen = eff_length$est_len))

write.table(fpkms, "data_fpkms.txt", sep="\t", quote=F, row.names=T)