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

相关文章

来自专栏熊二哥

深入入门系列--Data Structure--04树

终于有机会重新回头学习一下一直困扰自身多年的数据结构了,赶脚棒棒哒。一直以来,对数据结构的掌握基本局限于线性表,稍微对树有一丢丢了解,而对于图那基本上就是不懂(...

1729
来自专栏小樱的经验随笔

Codeforces 789A Anastasia and pebbles(数学,思维题)

A. Anastasia and pebbles time limit per test:1 second memory limit per test:256 ...

2798
来自专栏玉树芝兰

Review on HIT summer school

I took part in the summer school hosted by Harbin Institute of Technology. The T...

722
来自专栏SeanCheney的专栏

《大话数据结构》总结第一章 绪论第二章 算法第三章 线性表第四章 栈和队列第五章 字符串第六章 树第七章 图第八章 查找第九章 排序

第一章 绪论 什么是数据结构? 数据结构的定义:数据结构是相互之间存在一种或多种特定关系的数据元素的集合。 第二章 算法 算法的特性:有穷性、确定性、可行性、输...

2765
来自专栏计算机视觉战队

技术 | 用二进制算法加速神经网络

The original article is published on Nervana site: Accelerating Neural Networks ...

2847
来自专栏xcywt

《大话数据结构》树以及赫夫曼编码的例子

第六章 树 6.2 树的定义 树(Tree)的n个结点的有限集。当n=0时,称为空树。 任意一个非空树中: 1)有且仅有一个特定的称为根(root)的结点 2)...

3696
来自专栏小樱的经验随笔

Codeforces 810C Do you want a date?(数学,前缀和)

C. Do you want a date? time limit per test:2 seconds memory limit per test:256 m...

2515
来自专栏专注研发

poj-1005-l tanink i need a houseboat

Fred Mapper is considering purchasing some land in Louisiana to build his house ...

592
来自专栏数据结构与算法

codevs3002 石子归并 3

题目描述 Description 有n堆石子排成一列,每堆石子有一个重量w[i], 每次合并可以合并相邻的两堆石子,一次合并的代价为两堆石子的重量和w[i]+w...

2927
来自专栏WD学习记录

哈夫曼树(Huffman Tree)

给定n个权值作为n个叶子结点,构造一棵二叉树,若该树的带权路径长度达到最小,称这样的二叉树为最优二叉树,也称为哈夫曼树(Huffman Tree)。哈夫曼树是带...

934

扫码关注云+社区