人基因组注释包提取转录本cDNA长度

bioconductor包TxDb.Hsapiens.UCSC.hg38.knownGene提取转录本cDNA长度。
 

# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")

library("TxDb.Hsapiens.UCSC.hg38.knownGene") 
library(dplyr)



## 转录本cDNA长度
if(F){
  # 拿到所有基因的所有转录本,CDS,5UTR,3UTR等长度信息
  t_l=transcriptLengths(txdb)
  t_l=transcriptLengths(txdb,with.cds_len=TRUE,
                        with.utr5_len=TRUE, 
                        with.utr3_len=TRUE)
  head(t_l)
  t_l=na.omit(t_l)
  
  #增加一列，mRNA长度
  t_l <- mutate(t_l,mRNA_len = cds_len + utr5_len + utr3_len)
      
  # 按基因名，mRNA长度排序
  t_l=t_l[order(t_l$gene_id,t_l$mRNA_len,decreasing = T),]
  str(t_l)
  ## 取最长的mRNA代表转录本加工后的长度
  # 已经排序，第一位最长，去除后面冗余
  t_l_max=t_l[!duplicated(t_l$gene_id),] 
  head(t_l_max)
  t_l_max=t_l_max[,c(3,9)]
  t_l_max=t_l_max[order(t_l_max$gene_id),]
  head(t_l_max)
  ## 取cds平均长度位基因长度
  library(dplyr)
  
  t_l_mean <- t_l %>% 
    group_by(gene_id) %>% 
    summarise(
      avg_len = mean(cds_len, na.rm = TRUE)) 
  head(t_l_mean)
  
  ## 取mRNA长度中位数作为基因长度
  t_l_median <- t_l %>% 
    group_by(gene_id) %>% 
    summarise(
       median_len = median(cds_len, na.rm = TRUE)) 
  
  head(t_l_median)
  
}

head(t_l_max)

## id转基因symbol
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")  
s2g=toTable(org.Hs.egSYMBOL)
head(s2g)
t_l_max=merge(t_l_max,s2g,by='gene_id')
head(t_l_max)
## 查看某一基因的cDNA长度
t_l_max[which(t_l_max$symbol=='BRAF'),]