1. 程式人生 > 其它 >生物資訊學技能面試題(第5題)-根據GTF畫基因的多個轉錄本結構

生物資訊學技能面試題(第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)}