之前已经简单学习了测序raw counts的FPKM标准化，接下来就一个实际例子，利用R语言手动计算原始数据的FPKM值；并与官方给出的FPKM文件对比验证。

1、准备数据

• GSE81861上传的CRC_NM_epithelial的count与fpkm文件
  https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81861
  1-1
• 原始count矩阵

nm_count <- read.csv("GSE81861_CRC_NM_epithelial_cells_COUNT.csv/GSE81861_CRC_NM_epithelial_cells_COUNT.csv")
row.names(nm_count) <- nm_count[,1]
nm_count <- nm_count[,-1]
dim(nm_count)
nm_count <- as.matrix(nm_count)
#原始的基因名太长，我们仅留ensemble ID
rownames(nm_count) <- sapply(strsplit(rownames(nm_count),"_"),"[",3)
rownames(nm_count) <- sapply(strsplit(rownames(nm_count),"[.]"),"[",1)
#注意R好像不能直接以点号为分隔符，需要在两边加上英文中括号，才能识别
dim(nm_count) #57241个基因
nm_count[1:4,1:4]

1-2

计算FPKM的三要素：原始counts矩阵，样本总reads数，基因长度。其中样本总reads数直接使用colSums()函数即可。基因长度是需要加载特定的生物dataset package进行计算。因为基因长度有不同的定义标准，所以也有不同的计算思路，有意思的是结果可能也有很大的不同。

2、基因长度

2.1 TxDb.Hsapiens.UCSC.hg19.knownGene 包

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
if (!requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE))
  BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)
##The GRanges class is a container for the genomic locations and their associated annotations.

定义1：非冗余外显子长度之和

o <- findOverlaps(exon_txdb,genes_txdb)
o
genes_txdb[581] #chr1 11874-14409
exon_txdb[1] #chr1 11874-12227
exon_txdb[2] #chr1 12595-12721
exon_txdb[5] #chr1 13221-14409

如上可知道，581的gene包含了ID1,2,3,4,5的exon。所以基因长度之和就是求所有外显子的长度和。但注意要考虑外显子间可能存在交叉，所以不能直接把长度相加。
思路：假设外显子1序列3-10，2序列3-5。计算总长度即uniq所有的具体的数即可。3,4,5,6,7,8,9,10,3,4,5 unique下还是8个

2-1

  t1=exon_txdb[queryHits(o)]
  t2=genes_txdb[subjectHits(o)]
  t1=as.data.frame(t1)
  t1$geneid=mcols(t2)[,1]
# 得到exon_id与geneid的对应关系

2-2

g_l.1 <- lapply(split(t1,t1$geneid),function(x){
    #按gene id拆分表格
    head(x)
    tmp=apply(x,1,function(y){
      y[2]:y[3]
    }) #根据每一个gene所有exon的区间，生成区间内的整数，返回的为list。
    length(unique(unlist(tmp)))
    #计算共有多少种整数，即为最终的非冗余exon长度之和
  })
head(g_l.1) #为一个list
g_l.1=data.frame(gene_id=names(g_l.1),length=as.numeric(g_l.1))
dim(g_l.1)
head(g_l.1)

#为基因ID增添ENSEMBLE ID
library(org.Hs.eg.db)
s2g=toTable(org.Hs.egENSEMBL)
head(s2g)
g_l.1=merge(g_l.1,s2g,by='gene_id') 
#把g_l,s2g两个数据框以'gene_id'为连接进行拼接
head(g_l.1)

2-3

定义2：最长转录本

t_l=transcriptLengths(txdb)
head(t_l)
t_l=na.omit(t_l)
#先按基因ID，再按转录本长度从大到小排序
t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
head(t_l);dim(t_l)

2-4

#根据gene_id去重，选择第一个，也就是最长的那个
t_l=t_l[!duplicated(t_l$gene_id),]
head(t_l);dim(t_l)
g_l=t_l[,c(3,5)]
g_l.2=merge(g_l.2,s2g,by='gene_id') 
head(g_l.2)

2-5

2.2 biomaRt 包

https://bioconductor.org/packages/release/bioc/vignettes/biomaRt/inst/doc/biomaRt.html

if (!requireNamespace("biomaRt", quietly = TRUE))
  BiocManager::install("biomaRt")
ensembl <- useMart("ensembl") #connect to a specified BioMart database
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
#use the hsapiens(人类) dataset.或者直接如下设置
#ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
test <- getBM(attributes=c('ensembl_gene_id', 'start_position',
                           'end_position','ensembl_transcript_id',
                           'transcript_length'),
              mart = ensembl)

如图，分别找到基因ensemble ID、基因起止位点、转录本ID，转录本长度

2-6

test <- test[order(test$ensembl_gene_id,test$transcript_length,
                   decreasing = T),]
g_l.3 <- test[!duplicated(test$ensembl_gene_id),]
g_l.3 <- g_l[,c(1,5)]
head(g_l.3)

比较上述三种方法计算的基因长度

dim(g_l.1);dim(g_l.2);dim(g_l.3)
#[1] 23729     3
#[1] 25159     3
#[1] 67130     2
#TxDb.Hsapiens.UCSC.hg19.knownGene包的两种方法结果都只有2w+
#biomaRt包的结果有6w+的基因数量

g_l.1 <- g_l.1[,-1]
colnames(g_l.1) <- c("g_l.1","ensembl_id")
head(g_l.1)
g_l.2 <- g_l.2[,-1]
colnames(g_l.2) <- c("g_l.2","ensembl_id")
head(g_l.2)
colnames(g_l.3) <- c("ensembl_id","g_l.3")
head(g_l.3)

g_l_all <- merge(g_l.1, g_l.2, by="ensembl_id")
g_l_all <- merge(g_l_all, g_l.3, by="ensembl_id")
head(g_l_all,10)
summary(g_l_all) 

• 如下结果，三种方法总体分布大致相同，但是存在部分基因长度差别还是挺大的。此外数量方面两个包得到的结果也是有较大差异。最终计算FPKM值只有第三中方法比较理想（与官方给的FPKM文件相近），可能与数量差异存在一定关联。

  
  2-7
  2-8

3、批量计算FPKM

• 如上所说，我们选取第三种方法得出的基因长度进行演示

g_l <- g_l.3
dim(nm_count) #57241
dim(g_l) #67130
#选取二者交集基因
ng=intersect(rownames(nm_count),g_l$ensembl_gene_id) 
length(ng) #最后保留50813个
lengths=g_l[match(ng,g_l$ensembl_gene_id),2]
names(lengths) <- g_l[match(ng,g_l$ensembl_gene_id),1]
head(lengths)
#最终得到排好序的gene length数值向量，并进行命名
#ENSG00000000003 ENSG00000000005 ENSG00000000419 
#           3796            1205            1161 
#ENSG00000000457 ENSG00000000460 ENSG00000000938 
#           6308            4355            2637

#根据最终选取的counts矩阵计算每个样本的文库大小
nm_count <- nm_count[names(lengths),]
total_count <- colSums(nm_count)
head(total_count)
#RHC4104__stemTA__.2749FE RHC6087__stemTA__.2749FE RHL2880__stemTA__.2749FE 
#                593518.0                 702180.9                 103072.8 
#RHC5949__stemTA__.2749FE RHC4975__stemTA__.2749FE RHC4099__stemTA__.2749FE 
#               1636163.9                 964027.4                 991415.4

• 最终得到nm_count、lengths、·total_count，根据公式批量计算即可。
  3-1

nm_fpkm <- t(do.call( rbind,
                   lapply(1:length(total_count),
                          function(i){
                            10^9*nm_count[,i]/lengths/total_count[i]
                            #lengths向量自动遍历
                          }) ))
nm_fpkm[1:4,1:4]
#最终得到的fpkm矩阵

3-2

• 比较下官方给的fpkm矩阵

nm_fpkm_paper <- read.csv("GSE81861_CRC_NM_epithelial_cells_FPKM.csv/GSE81861_CRC_NM_epithelial_cells_FPKM.csv")
row.names(nm_fpkm_paper) <- nm_fpkm_paper[,1]
nm_fpkm_paper <- nm_fpkm_paper[,-1]
dim(nm_fpkm_paper)
nm_fpkm_paper <- as.matrix(nm_fpkm_paper)
rownames(nm_fpkm_paper) <- sapply(strsplit(rownames(nm_fpkm_paper),"_"),"[",3)
rownames(nm_fpkm_paper) <- sapply(strsplit(rownames(nm_fpkm_paper),"[.]"),"[",1)
colnames(nm_fpkm_paper) <- c(1:160)

• 如下图，总体来说还是蛮接近的，也验证了计算公式。

  
  3-3