生物信息学技能面试题(第5题)-根据GTF画基因的多个转录本结构

可以下载各种gtf,从NCBI,ENSEMBL,UCSC,GENCODE都可以!(记住,你下载什么样的gtf就需要修改成什么样的代码!!!)本文来源于我的个人博客:

画基因结构图! http://www.bio-info-trainee.com/1404.html

重点就是得到所有基因的转录本个数,以及每个转录本的外显子的坐标。 如图:

比如对这个ANXA1基因来说,非常多的转录本,但是基因的起始终止坐标,是所有转录本起始终止坐标的极大值和极小值!同时,它是一个闭合基因,因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的,虽然多个转录本会共享某些外显子,但是也存在某些外显子比同区域其它外显子长的现象!

再比如下面这个例子:对DDX11L11这个基因来说,前两个外显子都不会翻译,直到第三个外显子才开始翻译,构成CDS。

有些转录本是没有utr的,所以该转录本的起始坐标,就是CDS的起始坐标 这个非常有用,可以更新自己的一些概念:

1. 如果基因有多个转录本,基因的起始坐标,就是该基因所有转录本的第一个外显子的起始坐标的最小值,同理基因的终止坐标,就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。 2. 通过这个概念,可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。(这个是我乱说的,并没有这个定义) 3. 如果基因只有一个转录本,那么基因的起始终止坐标,就是转录本的起始终止坐标! 4. 一个基因的一个转录本的5’utr区域可以包括多个外显子区域,前者是翻译行为,后者是转录行为 ‍5. 起始密码子和终止密码子是CDS的起止处,是基于翻译的概念

6‍. ‍一个基因的多个转录本的外显子坐标不一定会排列整齐,每个转录本的剪切位点并不一定要比其它转录本一致!

R实现的代码如下:

rm(list=ls())## [url=http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gz]http://www.broadinstitute.org/ca ... n_goodContig.gtf.gz[/url]setwd('tmp')gtf <- read.table('gencode.v7.annotation_goodContig.gtf.gz',stringsAsFactors = F,                  header = F,comment.char = "#",sep = '\t'                  )table(gtf[,2])gtf <- gtf[gtf[,2] =='HAVANA',]gtf <- gtf[grepl('protein_coding',gtf[,9]),]lapply(gtf[1:10,9], function(x){  y=strsplit(x,';') })gtf$gene <- lapply(gtf[,9], function(x){  y <- strsplit(x,';')[[1]][5]  strsplit(y,'\\s')[[1]][3]  })draw_gene = 'TP53'structure = gtf[gtf$gene==draw_gene,]colnames(structure) =c(  'chr','db','record','start','end','tmp1','tmp2','tmp3','tmp4','gene')gene_start <- min(c(structure$start,structure$end))gene_end <- max(c(structure$start,structure$end))tmp_min=min(c(structure$start,structure$end))structure$new_start=structure$start-tmp_minstructure$new_end=structure$end-tmp_mintmp_max=max(c(structure$new_start,structure$new_end))num_transcripts=nrow(structure[structure$record=='transcript',])tmp_color=rainbow(num_transcripts)x=1:tmp_max;y=rep(num_transcripts,length(x))#x=10000:17000;y=rep(num_transcripts,length(x))plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))title(main = draw_gene,sub = paste("chr",structure$chr,":",gene_start,"-",gene_end,sep=""))j=0;tmp_legend=c()for (i in 1:nrow(structure)){  tmp=structure[i,]  if(tmp$record == 'transcript'){    j=j+1    tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))  }  if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)}

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-02-20

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏hotqin888的专栏

request+goquery+mahonia实现自动抓取网页数据

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/hotqin888/article/det...

544
来自专栏王小雷

SAS进阶《深入分析SAS》之数据汇总和展现

SAS进阶《深入分析SAS》之数据汇总和展现 1. 通过Print过程制作报表 proc print <data=数据集>; run; 选项: obs=修改观测...

15610
来自专栏潇涧技术专栏

When Math meets Android Animation (3)

当数学遇上动画:讲述ValueAnimator、TypeEvaluator和TimeInterpolator之间的恩恩怨怨(3)

782
来自专栏阿炬.NET

【人在江湖飘,哪有不带刀】神器Jumony

2686
来自专栏菩提树下的杨过

Flash/Flex学习笔记(22):滤镜学习

Silverlight中称之为“效果(Effect)”的东东,在Flash里叫“滤镜(Filter)",而且Flash里内置的滤镜要比Silverlight丰富...

1789
来自专栏生信技能树

(13)Hg19基因组的一些分析-生信菜鸟团博客2周年精选文章集

搞生信研究的,大部分数据都是针对于人类的,那么人类的参考基因组就不得不知了! 与hg19的突变相关的一些数据解释。 Hg19基因组的分析 R的bioconduc...

3056
来自专栏菩提树下的杨过

“AS3.0高级动画编程”学习:第三章等角投影(上)

什么是等角投影(isometric)? 刚接触这个概念时,我也很茫然,百度+google了N天后,找到了一些文章: [转载]等角(斜45度)游戏与数学 [转载]...

1975
来自专栏大数据文摘

你的数据科学python编程能力过关吗?看看这40道题你能得几分

1243
来自专栏大数据智能实战

tensorflow.models.rnn.rnn_cell.linear在tensorflow1.0版本之后找不到(附tensorflow1.0 API新变化)

由于版本更新关系,从原来的tensorflow低版本到升级到tensorflow1.0以上时,发现有很多API函数变化是很正常的事情,大多碰到的如: 如其中tf...

2437
来自专栏R语言___生物信息

基因功能富集分析-R语言

BgRation:所有background基因中与该Term相关的基因数与所有background基因的比值

56110

扫描关注云+社区