以前，NCBI的genomic browser中是可以得到一个所谓的exon number的。其规则为同一个基因内所有exon按5′->3’排序，用排序序号来做为exon number。对于有重叠的exon使用同一序号，以排序先后为其加上a,b,c…的后缀。
但是这其实是一个很不靠谱的表达方式，因为随着基因的转录本的变化，有可能会新增加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)) |
运行完程序后，什么都看不到啊？请问一下这段程序运行的结果应该是怎么样的？
记数结果就在exon.numbering里保存着啊。