美文网首页生物信息学与算法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