如何得到exon numbering

~ admin

以前，NCBI的genomic browser中是可以得到一个所谓的exon number的。其规则为同一个基因内所有exon按5′->3’排序，用排序序号来做为exon number。对于有重叠的exon使用同一序号，以排序先后为其加上a,b,c…的后缀。
art9-cc-figure-1(fix)
但是这其实是一个很不靠谱的表达方式，因为随着基因的转录本的变化，有可能会新增加exon。所以一些文章中使用exon number来标注alternative splicing的方式已经渐渐被淘汰。所以我们现在很难在ncbi或者ucsc的genomic browser中再找到exon number了。

但是因为使用exon number标记仍然是一种习惯，很多时候还不得不在文章中使用到。

下面代码示例如何在bioconductor中得到正确的exon number.

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(Gviz)
library(org.Hs.eg.db)
exons <- GeneRegionTrack(TxDb.Hsapiens.UCSC.hg19.knownGene)
exons <- exons@range
exons <- exons[!is.na(exons$gene)]
exons$symbol <- as.character(mget(exons$gene, org.Hs.egSYMBOL, ifnotfound=NA))
genelist <- c("RALGAPA1", "RPS6KB1", "SYNCRIP", "KIF1B", "CSPP1", "TPR", "MYO6", "DENND4A", "TRRAP", "CHD2", "PRKDC", "MBNL1", "KIAA1109", "MEF2A", "MYCBP2")
exon.numbering <- list()
for(gene in genelist){
	exons.sel <- exons[exons$symbol==gene]
	exons.sel <- exons.sel[exons.sel$feature=="CDS"]
	exons.sel <- unique(exons.sel)
	exons.sel.rd <- reduce(exons.sel, with.revmap=T, min.gapwidth=0)
	exons.sel.rd$ExonNumbering <- ifelse(strand(exons.sel.rd)=="+", 1:length(exons.sel.rd), length(exons.sel.rd):1) 
	exons.sel.rev <- mapply(function(.id, .num){
		.ele <- exons.sel[.id]; .ele$ExonNumbering <- as.character(.num)
		if(length(.ele)>1) .ele$ExonNumbering <- paste(.ele$ExonNumbering, letters[1:length(.ele)], sep="")
		.ele
		}, exons.sel.rd$revmap, exons.sel.rd$ExonNumbering)
	exons.sel.rev <- unlist(GRangesList(exons.sel.rev))
	exon.numbering[[gene]] <- exons.sel.rev
}

新代码：

library(EnsDb.Hsapiens.v75)
exons <- exonsBy(EnsDb.Hsapiens.v75, by="gene")
genelist <- c("ENSG00000012048", "ENSG00000089280", "ENSG00000111247", "ENSG00000144554", "ENSG00000165392")
exons <- exons[genelist]
exon.numbering <- list()
for(i in 1:length(exons)){
	exons.sel <- exons[[i]]
	exons.sel <- unique(exons.sel)
	exons.sel.rd <- reduce(exons.sel, with.revmap=T, min.gapwidth=0)
	exons.sel.rd$ExonNumbering <- ifelse(strand(exons.sel.rd)=="+", 1:length(exons.sel.rd), length(exons.sel.rd):1) 
	exons.sel.rev <- mapply(function(.id, .num){
		.ele <- exons.sel[.id]; .ele$ExonNumbering <- as.character(.num)
		if(length(.ele)>1) .ele$ExonNumbering <- paste(.ele$ExonNumbering, letters[1:length(.ele)], sep="")
		.ele
		}, exons.sel.rd$revmap, exons.sel.rd$ExonNumbering)
	exons.sel.rev <- unlist(GRangesList(exons.sel.rev))
	exon.numbering[[names(exons)[i]]] <- exons.sel.rev
}
 exons1 <- mapply(function(query, subject) {
         ol <- findOverlaps(query, subject, type="equal")
         query <- query[queryHits(ol)]
         query$ExonNumbering <- subject[subjectHits(ol)]$ExonNumbering; query}, 
    exons, exon.numbering, SIMPLIFY=FALSE)
exons2 <- unlist(GRangesList(exons1))

发布者：admin

2 thoughts on “如何得到exon numbering”

  1. 运行完程序后，什么都看不到啊？请问一下这段程序运行的结果应该是怎么样的？

    回复

发表评论

邮箱地址不会被公开。 必填项已用*标注