搜索
查看: 254|回复: 1

生信编程直播第五题:根据GTF画基因的多个转录本结构

[复制链接]

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
发表于 2017-1-4 23:43:32 | 显示全部楼层 |阅读模式
转载自我自己的博客:画基因结构图!
可以下载各种gtf,从NCBI,ENSEMBL,UCSC,GENCODE都可以!
重点就是得到所有基因的转录本个数,以及每个转录本的外显子的坐标。
如图:


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

再比如下面这个例子:对DDX11L11这个基因来说,前两个外显子都不会翻译,直到第三个外显子才开始翻译,构成CDS。
有些转录本是没有utr的,所以该转录本的起始坐标,就是CDS的起始坐标
这个非常有用,可以更新自己的一些概念:
1,如果基因有多个转录本,基因的起始坐标,就是该基因所有转录本的第一个外显子的起始坐标的最小值,同理基因的终止坐标,就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。
2,通过这个概念,可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。(这个是我乱说的,并没有这个定义)
3,如果基因只有一个转录本,那么基因的起始终止坐标,就是转录本的起始终止坐标!
4,一个基因的一个转录本的5’utr区域可以包括多个外显子区域,前者是翻译行为,后者是转录行为
5,起始密码子和终止密码子是CDS的起止处,是基于翻译的概念
6,一个基因的多个转录本的外显子坐标不一定会排列整齐,每个转录本的剪切位点并不一定要比其它转录本一致!



R实现的代码如下:
[AppleScript] 纯文本查看 复制代码
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_min
structure$new_end=structure$end-tmp_min
tmp_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)
}













上一篇:外显子重复的基因,是如何定量的呢??
下一篇:生信编程直播第二题-hg19基因组序列的一些探究
回复

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-13 21:25:04 | 显示全部楼层
感谢群主的代码,我把这一题重复了一遍,代码太经典,我没法改,所以就不上传了。
有个问题不明白
tmp_legend表示的每条转录本的信息
if(tmp$record == 'transcript'){
    j=j+1
    tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
  }
问题是我在画图中没有发现tmp_legend在图中应该位于什么地方。标题,副标题都是明确的。所以,不知道为啥要写这个代码?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-1-22 03:30 , Processed in 0.210260 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.