TCGA系列学习笔记(4)差异分析结果注释

TCGA系列学习笔记(4)差异分析结果注释

内容目录

前言载入数据关于ID的结构处理gtf文件下载注释文件用R处理gtf文件后记

前言

知识笔记来源:

简书总址:https://www.jianshu.com/u/c93ac360691a

fasta、fastq、vcf、bam、bed、gtf–多种生信格式的R语言读取:https://www.jianshu.com/p/18aa866a8c20

关于不同ID之间的介绍:https://www.jianshu.com/p/13c0a04fd507

同样,吃水不忘挖井人,还是要远程云感谢jimmy老师和美小洁老师的代码教程。

这里其实很简单,但是考虑到,这里介绍了多种linux文件在R里的处理方法,所以,还是整理下好了。

载入数据

if (T) {
  rm(list = ls())
  options(stringsAsFactors = F)
  load("step1_gdc.Rdata")
  load("DEG.Rdata")
}

# 基本检查数据及数据转换
> head(DESeq2_DEG,3)
baseMean log2FoldChange     lfcSE      stat       pvalue         padj change
ENSG00000105607.11 2295.7385      -3.312555 0.1836103 -18.04123 9.246210e-73 2.413168e-68    NOT
ENSG00000052802.11 7540.8294      -3.770359 0.2186743 -17.24189 1.287541e-66 1.680177e-62    NOT
ENSG00000080709.13  468.7602      -5.608751 0.3340184 -16.79174 2.804891e-63 2.440161e-59   DOWN
> head(edgeR_DEG,3)
logFC   logCPM       LR       PValue          FDR change
ENSG00000080709.13 -5.584815 2.877132 343.7979 9.500838e-77 8.906011e-73   DOWN
ENSG00000042781.11 -6.435770 3.109970 343.6720 1.011989e-76 8.906011e-73   DOWN
ENSG00000105607.11 -3.279389 5.180216 320.8182 9.609167e-72 5.637698e-68    NOT
> head(limma_voom_DEG,3)
logFC    AveExpr         t      P.Value    adj.P.Val        B change
ENSG00000042781.11 -6.535457 -0.2487210 -24.49253 2.734489e-28 2.406487e-24 53.76319   DOWN
ENSG00000266964.4  -6.884717  0.1108023 -22.99443 4.173883e-27 1.685428e-23 51.28812   DOWN
ENSG00000251424.1  -8.250455 -4.0615968 -24.57640 2.357502e-28 2.406487e-24 51.27016   DOWN

可以看到这里的差异分析结果中,我们只知道基因的ID是ENSGxxxxxx。但是我们需要的是基因symbol,所以就需要我们进行ID转换。

关于ID转换,其实之前就介绍过很多方法,有兴趣的自己去找找技能树里关于ID转换的工具吧,包括jimmy老师特意写了几个ID转换的R包供我们使用

这里要介绍的是最可靠的注释手段——利用gtf进行注释。

关于ID的结构

结构如下:

[species prefix][feature type prefix][a unique eleven digit number]

关于对应关系可以在官网查看:http://asia.ensembl.org/info/genome/stable_ids/prefixes.html

常用的物种编号:

species prefix 物种
ENS Homo sapiens (Human)
ENSMUS Mus musculus (Mouse)
ENSDAR Danio rerio (Zebrafish)

关于测序类型:

第4个字母 含义
E exon
FM Ensembl protein family
G gene
GT gene tree
P protein
R regulatory feature
T transcript

处理gtf文件

下载注释文件

去到GDC官网,点击support

然后进入到reference file里:

选择RNA-seq用到的注释信息:

点击下载gencode.v22.annotation.gtf.gz即可:

https://api.gdc.cancer.gov/data/25aa497c-e615-4cb7-8751-71f744f9691f

大概也就39MB,也不大,因为这是纯文本文件,所以压缩比很高,解压后大概1个多GB。不过由于很多R包在读入文件时可以不需要解压,所以我们就不要解压了,直接用函数读入就好了。

用R处理gtf文件

主要是使用一个R包——rtracklayer

> library(rtracklayer)
# x是个Granges对象
> x = rtracklayer::import("gencode.v22.annotation.gtf.gz")
# 转化为数据框,非常好用
> x2 = as.data.frame(x)
> tj = as.data.frame(table(x2$type))
> tj
            Var1    Freq
1           gene   60483
2     transcript  198442
3           exon 1172082
4            CDS  699443
5    start_codon   82228
6     stop_codon   74337
7            UTR  276542
8 Selenocysteine     114

这里直接读入时的数据是个Granges对象,这种对象我之前也介绍过:

Ranges对象和GRanges对象的介绍

虽然也很好用,但是还是不如直接一个as.data.frame转为我们最熟悉的数据框结构!

后面的操作就很简单了,但是需要你对gtf文件非常了解,同时还要懂一丢丢R语言,所以,不懂的自己去补习知识吧。

关于gtf文件中每种基因类型的介绍需要自己去官网认真阅读说明书

https://www.gencodegenes.org/pages/biotypes.html

如果是lncRNA,我们需要拿到:

如果是mRNA(包括ORF),我们只需要拿到:

不过这里倒是让我涨了个奇怪的知识——mRNA是包括ORF的,而不包括ORF的mRNA则成了lncRNA了

> lnc_bype = c("3prime_overlapping_ncRNA", "antisense", "bidirectional_promoter_lncRNA", "lincRNA", "macro_lncRNA", "non_coding", "processed_transcript", "sense_intronic" , "sense_overlapping")
> table(x2$gene_type %in% lnc_bype)

FALSE  TRUE
45657 14826
> table(x2$gene_type == “protein_coding”)

FALSE  TRUE
40669 19814

所以有14826个lncRNA,有19814个mRNA。

接下来把所有的都收集并保存起来即可:

lnc_anno = x2[x2$gene_type %in% lnc_bype,c("gene_name","gene_id","gene_type")]
mRNA_anno = x2[x2$gene_type == "protein_coding",c("gene_name","gene_id","gene_type")]
save(lnc_anno,mRNA_anno,file = "gtf.anno.Rdata")

让我们看看得到的数据长什么样:

mRNA:

lncRNA:

处理到这一步,后面的就没什么好讲的了,直接用mergematch这几个函数就可以完成几乎你想要的一切操作了。

后记

这些应该都算是基础知识了,大家可以多研究下最常用的R语言函数,相信也会很快弄明白这些东西的~

生物信息学

使用seaborn绘制热图

2020-10-7 21:45:06

生物信息学

TCGA系列学习笔记(6)多因素生存分析

2020-10-7 22:22:50

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索