生物信息学技能面试题(第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 条评论
登录 后参与评论

相关文章

来自专栏码匠的流水账

聊聊hikari连接池的isAllowPoolSuspension

本文主要研究一下hikari连接池的isAllowPoolSuspension属性

1502
来自专栏余生开发

echarts太阳分布图-饼图来回穿梭

var dom = document.getElementById("container");

1202
来自专栏Hadoop数据仓库

Oracle sqlldr 如何导入一个日期列

1. LOAD DATA INFILE * INTO TABLE test FIELDS TERMINATED BY X'9' TRAILING NULLCO...

1816
来自专栏一个会写诗的程序员的博客

SQLNonTransientConnectionException: CLIENT_PLUGIN_AUTH is required

java.sql.SQLNonTransientConnectionException: CLIENT_PLUGIN_AUTH is required at ...

2722
来自专栏lonelydawn的前端猿区

js多级联动示例(省份和城市二级联动)

本文描述的省份和城市二级联动下拉菜单,在js部分设置省份和城市信息,当然也可以将数据保存在数据库中,当需要使用的时候在读取。这里用的是js原生代码 在html文...

3495
来自专栏Pulsar-V

Save Camera Document

#pragma once #include "HCCamera.h" #include <time.h> #include <cstdio> #incl...

2828
来自专栏一个会写诗的程序员的博客

java.sql.BatchUpdateException: Can not issue SELECT via executeUpdate() or executeLargeUpdate().

java.sql.BatchUpdateException: Can not issue SELECT via executeUpdate() or execu...

893
来自专栏曾大稳的博客

MediaCodec硬编码pcm2aac

MediaCodec是Android(api>=16)提供的一个多媒体硬解编码库,能实现音视频的编解码。

1382
来自专栏Flutter入门

Android平台硬件编码总结一

1482
来自专栏一个会写诗的程序员的博客

java.base.jmod

/Library/Java/JavaVirtualMachines/jdk-9.jdk/Contents/Home/jmods$ jmod list java....

1112

扫码关注云+社区