FPKM-UQ（Fragments Per Kilobase Million - Upper Quartile）是一种用于量化基因表达的方法，主要用于RNA-seq数据的分析。
FPKM（Fragments Per Kilobase of transcript per Million mapped reads）这是一种常用的标准化方法，用于比较不同样本中基因的表达水平。FPKM考虑了每个基因的长度和测序深度，使得不同基因和不同样本之间的表达量可以直接比较。
然后，FPKM-UQ是FPKM的一个变体，它在计算FPKM之后，还进行了一个额外的标准化步骤。在这个步骤中，所有基因的FPKM值被除以所有基因FPKM值的上四分位数（Upper Quartile）。这样可以进一步减小高表达基因对总体表达水平估计的影响，使得结果更加稳健。上四分位数是指将一个数据集分为四等份后，处于第三个区间的值。
FPKM-UQ是一种更加稳健的基因表达量化方法，特别适合于处理有大量高表达基因或表达差异大的数据。
# 计算FPKM
computeFPKM <- function(counts, lengths) {
mappedReads = sum(counts)
geneLengthsKb = lengths / 1000
fpkm = (counts / geneLengthsKb) / (mappedReads / 1e6)
return(fpkm)
}
# 计算FPKM-UQ
computeFPKM_UQ <- function(fpkm) {
upperQuartile = quantile(fpkm, 0.75)
fpkm_uq = fpkm / upperQuartile
return(fpkm_uq)
}
# 假设我们有一些基因表达计数和对应的基因长度
counts = c(10, 20, 30, 40, 50)
lengths = c(1000, 2000, 3000, 4000, 5000)
# 计算FPKM
fpkm = computeFPKM(counts, lengths)
# 计算FPKM-UQ
fpkm_uq = computeFPKM_UQ(fpkm)
print(fpkm_uq)
说明：虽然下面有两种计算基因长度的方法，但是目前看来都没有上面这个好
转换公式
countToTpm <- function(counts, effLen)
{
rate <- log(counts) - log(effLen)
denom <- log(sum(exp(rate)))
exp(rate - denom + log(1e6))
}
### FPKM = [counts / (gene length in Kb * total counts in millions)]
countToFpkm <- function(counts, effLen)
{
N <- sum(counts)
exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
countToEffCounts <- function(counts, len, effLen)
{
counts * (len / effLen)
}
这里的efflen就是基因长度或者说总外显子长度
如果你要计算的话，切记用apply对每一列，也就是每个样本进行计算,不要直接把矩阵丢进去
apply(dataexp,2,function)
获取外显子的长度方法1
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb<-TxDb.Hsapiens.UCSC.hg38.knownGene
## 定义基因的长度为非冗余exon长度之和
if (F) {
exon_txdb=exons(txdb) #取出人类基因组中的全部外显子
genes_txdb=genes(txdb)
o=findOverlaps(exon_txdb,genes_txdb) # 找出外显子和基因之间的重叠部分
o
t1=exon_txdb[queryHits(o)]# 将与gene重叠的exon提取出来
t2=genes_txdb[subjectHits(o)]
t1=as.data.frame(t1)
t1$geneid=mcols(t2)[,1]
#lapply:历遍列表向量中的每个元素，并且使用指定的函数来对其元素进行处理。返回向量列表。
# 函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组；返回值为列表
g_l=lapply(split(t1,t1$geneid), function(x){
head(x)
tmp=apply(x, 1,function(y){
y[2]:y[3]
})
length(unique(unlist(tmp)))
})
head(g_l)
g_l=data.frame(gene_id=names(g_l), length=as.numeric(g_l))
save(g_l,file ='hg38_g_l.RData')
}
load('hg38_g_l.RData')
## 下面是定义基因长度为最长转录本长度
if(F) {
t_l=transcriptLengths(txdb)
head(t_l)
t_l=na.omit(t_l)
t_l=t_l[order(t_l$gene_id, t_l$tx_len, decreasing = T),]
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.Hs.eg.db)
s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
gl=merge(g_l,s2g,by='gene_id') # merge函数可以实现对两个数据框根据共同的列名来进行连接。
获取外显子的长度方法2
Htseq Count To Fpkm | KeepNotes blog (bioinfo-scrounger.com)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("hg38.gtf",format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
如何优雅的统计基因外显子长度 - 云+社区 - 腾讯云 (tencent.com)
能不能反过来算
很多的时候我们从数据库中得到的不是原始数据count，而是fpkm或者是log2(fpkm+1)，那么我们就想知道能不能反过来进行计算？
这里说一下结论，之前简单推了一下（不知道有没有错），发现是没有办法反过去算的，比如你只知道FPKM，想逆过来算count，按照上面的公式，你会发现，你得先知道Count的内容，比如总的count的数量，才能逆过来计算Count,否则没法计算。这就是死循环。
