美文网首页生物信息学与算法R生信R语言做生信
【r<-方案】从gtf文件中抽取基因id和name

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

作者: 王诗翔 | 来源:发表于2019-01-22 18: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)
}

相关文章

网友评论

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

    本文链接:https://www.haomeiwen.com/subject/jterjqtx.html