生物信息学技能面试题(第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_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)
}