在正式分析之前，对于数据的处理是至关重要的，这种重要性是体现在很多方面，其中有一点是要求分析者采用正确的数据类型。
对于芯片数据，原始数据进行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)
##FPKM转TPM
fpkmToTpm <- function(fpkm){ exp(log(fpkm) - log(sum(fpkm)) + log(1e6))}tpms <- apply(fpkms,2,fpkmToTpm)
#！！！检查，所有样本的总和必须一样，10的6次方，才代表转化成功
colSums(tpms)
range(tpms)
#保存tpm文件
write.table(tpms, "GSEXXXXXX_tpms.txt", sep="\t", quote=F, row.names=T)
#log化
boxplot(tpms,las=2)
tpms <- log2(tpms+1)
boxplot(tpms,las=2)
#保存log化文件write.table(tpms, "GSEXXXXXX_tpms_log2.txt", sep="\t", quote=F, row.names=T)
GEO中有些数据集直接上传的count文件，我们目前较多拿来分析的都是FPKM或者TPM格式
扫一扫
8246
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
