2015-03-25 56 views
2

我想提取染色體的基因間座標。我做的代碼塊,但因爲我是新來的這些包,我不知道我是否按照正確的邏輯在這裏:查找基因間區域

library(IRanges) 
library(GenomicFeatures) 
library(TxDb.Hsapiens.UCSC.hg19.knownGene) 
txdb = transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, by = "gene",use.names=TRUE) 

#For example, only I am interested in intergenic coordinates in chromosome 1 
seqlevels(txdb, force=TRUE) = c("chr1") 

#Creates IRanges 

ir = IRanges(start=unlist(start(txdb)), end=unlist(end(txdb)), names=names(txdb)) 

# Reduce overlapping gene positions to continuous range 
ir.red = reduce(ir) # Collapses ranges of overlapping genes to continuous range. 

#Identify overlaps among the initial and reduced range data sets 
overlap = findOverlaps(ir, ir.red) 

我必須採取的「+」和關懷「 - 」鏈?

回答

2

最簡單的方法是先從基因()存取,減少()的,並採取了補充:

library(TxDb.Hsapiens.UCSC.hg19.knownGene) 
genic <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene) 
genic <- reduce(genic, ignore.strand=T) 
intergenic <- gaps(genic) 
intergenic <- intergenic[strand(intergenic) == "*"] #This is important!!! 

最後一步是非常重要的,否則你會得到一個額外的2項每個染色體(每個+和 - )。

0

雖然不是您的問題的直接答案,但值得一提的是,您可以直接從UCSC表瀏覽器http://genome.ucsc.edu/cgi-bin/hgTables獲取這些數據。表格瀏覽器允許您設置兩個軌道或它們的補充之間的交集,甚至可以將它們返回到表瀏覽器中的自定義軌道以進行進一步的操作,從而允許針對多個軌道和表的任意複雜的嵌套查詢。步驟爲您的目的如下...

1)選擇您的基因組和感興趣的大會。 2)選擇「基因和基因預測」作爲組,「UCSC基因」作爲軌道,「已知基因」作爲表格。

3)如果你想有一個染色體的單染色體或一部分,對「區域」行框中輸入內容,否則將其設置爲「基因組」

4)單擊「創建」上5)進入下一頁後,選擇group =「Mapping and Sequencing」,track =「Assembly」,table =「Assembly(gold)」,點擊單選按鈕「Base-pair- (UCSC基因和組裝)的智能相交(AND)「,勾選」在基礎對明智的相交/聯合之前補充UCSC基因「框,然後單擊」提交「。

6)回到主表瀏覽器屏幕上,輸入任何你想要的輸出選項,並點擊「獲取輸出」來檢索基因間區域。

快樂狩獵!