R语言实现基因表达量单位转换：从Counts到FPKM、RPKM、TPM和CPM

基本概念解释

• Counts：直接从测序数据得到的读数计数。

• FPKM (Fragments Per Kilobase of transcript per Million mapped reads)：每百万映射读数中每千碱基转录本的片段数，用于RNA-seq数据标准化。

• RPKM (Reads Per Kilobase of transcript per Million mapped reads)：与FPKM类似，但用于单端测序数据。

• TPM (Transcripts Per Million)：每百万转录本中的转录本数，解决了FPKM/RPKM中先标准化后比较的问题。

• CPM (Counts Per Million)：每百万读数中的读数计数，用于简化的数据标准化。

从GTF文件提取基因长度信息

首先，我们需从GENCODE下载所需的GTF文件。在本例中，我使用的示例数据源自GSE224127的小鼠RNA-seq数据。由于文章中采用的参考基因组为mm10，我们因此从GENCODE下载mm10版本的GTF文件，以便在后续的标准化过程中使用，计算有效基因长度。

gtf_path <- "data/gencode.vM25.basic.annotation.gtf.gz"
gtf <-
  as.data.frame(x = rtracklayer::import(con = gtf_path))
table(gtf$type)

          gene     transcript           exon            CDS    start_codon 
         55401          81540         527731         428717          45617 
    stop_codon            UTR Selenocysteine 
         45382         117716             59 

exon <- gtf[
  gtf$type == "exon",
  c("start", "end", "gene_name")
]

exon_by_gene_name <- split(exon, exon$gene_name)

cl <- parallel::makeCluster(0.5 * parallel::detectCores())
gene_length <-
  parallel::parLapply(
    cl = cl,
    X = exon_by_gene_name,
    fun = function(x) {
      tmp <-
        apply(x, 1, function(y) {
          y[1]:y[2]
        })
      length(unique(unlist(tmp)))
    }
  )

gene_length <- data.frame(
  gene_name = names(gene_length),
  length = as.numeric(gene_length)
)
write.table(
  x = gene_length,
  file = "data/gencode.vM25.basic.gene_length.tsv",
  sep = "\t",
  row.names = FALSE,
  quote = FALSE
)

gene_length <-
  read.table(
    file = "data/gencode.vM25.basic.gene_length.tsv",
    header = TRUE,
    row.names = 1
  )

读取counts表达矩阵

counts <-
  read.table(
    file = "data/GSE224127_2_counts_CLP_timecourse.txt.gz",
    row.names = 1,
    header = TRUE,
    check.names = FALSE
  )
counts <- counts[rowSums(x = counts) > 0, ]
rownames(x = counts) <-
  gsub(
    pattern = "\\.",
    replacement = "-",
    rownames(x = counts)
  )
intersect_features <- intersect(
  x = rownames(x = counts),
  y = rownames(x = gene_length)
)

counts <- counts[intersect_features, ]
gene_length <-
  as.matrix(x = gene_length[intersect_features, , drop = FALSE])

Counts到CPM的转换

y <- edgeR::DGEList(counts = counts)
cpm <- edgeR::cpm(y = y, log=FALSE)

Counts到FPKM和RPKM的转换

total_counts <- colSums(x = counts)
gene_length_kb <- as.vector(x = gene_length / 1000)

fpkm <- sweep(counts, 2, total_counts, FUN = "/") * 1e6 / gene_length_kb

Counts到TPM的转换

counts_per_kb <- sweep(counts, 1, gene_length_kb, FUN = "/")
tpm <- sweep(counts_per_kb, 2, colSums(x = counts_per_kb), FUN = "/") * 1e6