「R」从gtf文件中抽取基因id和name

王诗翔呀

发布于 2020-07-03 04:10:46

4.9K00

代码可运行

文章被收录于专栏：优雅R优雅R

运行总次数：0

代码可运行

参考文章http://www.bioinfo-scrounger.com/archives/342计算FPKM值，发现计算完每个基因下所有外显子的总长度后，记录的都是ENSEMBL gene id，而我需要的是symbol。奇怪的是GenomicFeatures既然把GTF文件读取进去了还抽取基因id了，但它就是不提供抽gene symbol的功能。

尝试使用clusterProfiler包装的转换器进行转换，发现基因丢了一半，这可不行。谷歌了一波没有发现满意的答案，有个refGenome包好像可以做，但读取文件半天卡死了，特别奇怪。最后还是自己动手，完成了6万个gene feature的转换。

整个提取操作包装为函数了，输入可以是文件名或已经导入的gtf文件数据框（最好还是文件吧）。由data.table包支持，速度杠杠的！

get_map = function(input) {
  if (is.character(input)) {
    if(!file.exists(input)) stop("Bad input file.")
    message("Treat input as file")
    input = data.table::fread(input, header = FALSE)
  } else {
    data.table::setDT(input)
  }
  
  input = input[input[[3]] == "gene", ]
  
  pattern_id = ".*gene_id \"(ENSG[0-9]+)\";.*"
  pattern_name = ".*gene_name \"([^;]+)\";.*"
  
  gene_id = sub(pattern_id, "\\1", input[[9]])
  gene_name = sub(pattern_name, "\\1", input[[9]])
  
  data.frame(gene_id = gene_id,
             gene_name = gene_name,
             stringsAsFactors = FALSE)
}