ID转换

作者: 小胡同学ime | 来源:发表于2022-01-24 16:50 被阅读0次

    前面处理步骤之后,行名时官方格式,转换成对应的genesymble


    image.png
    image.png

    操作步骤可见下文件夹


    image.png

    title: "基因ID转换"
    author: "Sun Xiaojie"
    output: rmarkdown::html_vignette
    editor_options:
    chunk_output_type: console


    knitr::opts_chunk$set(
      collapse = TRUE,
      comment = "#>"
    )
    knitr::opts_chunk$set(fig.width = 6,fig.height = 6,collapse = TRUE)
    knitr::opts_chunk$set(message = FALSE)
    options(rmarkdown.html_vignette.check_title = FALSE)
    

    01.需求

    TCGA的RNA-seq数据使用的geneid是ensembl id,两个常见的需求:

    1.差异分析结果中每个ensembl id对应的symbol和类型(mRNA/lncRNA或其它)

    2.将行名从ensembl id 转换为symbol

    02.思路

    1.找到TCGA数据对应的参考基因组注释版本。

    2.下载该版本的参考基因组注释文件,提取ensembl id 与symbol的对应关系及每个基因的gene type信息。

    3.可以将symbol和gene type 用merge添加到差异分析结果中,也可以在差异分析前先转换矩阵的行名。

    03.动起来

    1.找参考基因组版本

    gdc首页的support

    image

    about the data - GDC Reference Files

    image

    可以看到,使用的参考基因组版本是genecode的v22。(版本很多,这个是14年的版本了)

    image

    2.找区分类型的列

    在gtf文件里并不是直接分出了lncRNA,需要找gtf文件里对biotype的说明,不看不知道,一看发现这是一个很长的表格。

    image

    其中对lncRNA的说明是:

    Generic long non-coding RNA biotype that replaced the following biotypes: 3prime_overlapping_ncRNA, antisense, bidirectional_promoter_lncRNA, lincRNA, macro_lncRNA, non_coding, processed_transcript, sense_intronic and sense_overlapping.

    所以需要将genetype里这些类型对应的行挑出来,就是lncRNA了。
    然后与表达矩阵行名进行匹配替换,就可以分别得到mRNA和lncRNA的矩阵了。

    options(stringsAsFactors = F)
    if(!file.exists("gtf_gene.Rdata")){
      #step1:读取并探索gtf文件----
      #BiocManager::install("rtracklayer")
      library(rtracklayer)
      gtf = rtracklayer::import("gencode.v22.annotation.gtf")
      class(gtf)
      gtf = as.data.frame(gtf);dim(gtf)
      colnames(gtf)
      table(gtf$type)
      #step2:先筛选出gene对应的行
      gtf_gene = gtf[gtf$type=="gene",]
      save(gtf_gene,file = "gtf_gene.Rdata")
    }
    load("gtf_gene.Rdata")
    load("TCGA-CHOL_DEG.Rdata")
    deg = DESeq2_DEG
    table(rownames(deg) %in% gtf_gene$gene_id)  #测试一下deg(作为模版必须齐全)是否包含所有gene_id,下面用merge的话就可以不全部包含,match必须包含
    
    an = gtf_gene[,c("gene_name","gene_id","gene_type")]
    deg = merge(deg,an,by.x = "row.names",by.y = "gene_id")
    
    # mRNA和lncRNA总共有多少个?
    
    lnc = c("3prime_overlapping_ncRNA", "antisense", "bidirectional_promoter_lncRNA", "lincRNA", "macro_lncRNA", "non_coding", "processed_transcript", "sense_intronic" , "sense_overlapping")
    
    k1 = gtf_gene$gene_type %in% lnc;table(k1)
    k2 = gtf_gene$gene_type == "protein_coding";table(k2)
    
    # deg中有多少mRNA和lncRNA?
    
    k3 = deg$gene_type %in% lnc;table(k3)
    k4 = deg$gene_type =="protein_coding";table(k4)
    
    # 差异的 mRNA和lncRNA 各有多少
    k5 = deg$change !="NOT"
    table(k3&k5)
    table(k4&k5)
    

    表达矩阵的行名id转换

    rm(list = ls())
    load("TCGA-CHOL_gdc.Rdata")
    load("gtf_gene.Rdata")
    an = gtf_gene[,c("gene_name","gene_id","gene_type")]
    exp = exp[rownames(exp) %in% an$gene_id,]  #exp里面存在an里所有基因ID挑选出来
    an = an[match(rownames(exp),an$gene_id),] #match要求是两向量内容相同顺序可以不一样,若内容不一样则要求模版基因都要在比对数据里面
    identical(an$gene_id,rownames(exp))   #此时exp和an和行名完全一致
    
    k = !duplicated(an$gene_name);table(k)   #去除相同的行名,R语言是不允许行名重复的
    
    an = an[k,]
    exp = exp[k,]
    
    rownames(exp) = an$gene_name
    
    # 最终得到的结果
    exp[1:2,1:2]
    
    save(exp,file = paste0(cancer_type,"_symbol_exp.Rdata"))
    

    相关文章

      网友评论

          本文标题:ID转换

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