生物信息学与算法R生信R语言做生信

【r<-方案】从gtf文件中抽取基因id和name

2019-01-22  本文已影响15人  王诗翔

参考文章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)
}

上一篇下一篇

猜你喜欢

热点阅读