上一期講到了GRanges的特點蜀撑，這一期我們講講其具體用法看铆。
首先我們需要引入一個新的函數(shù)DataFrame徽鼎，它與data.frame類似，但是它的列變量可以是任何類型弹惦，比如IRanges否淤。舉個例子：
library(GenomicRanges)
ir <- IRanges(start = 1:2, width = 3)
df1 <- DataFrame(iranges = ir)
df1
## iranges
## <IRanges>
## 1 [1, 3]
## 2 [2, 4]
它也可以使用$符號，獲得某一個變量的數(shù)值棠隐，如：
df1$iranges
## IRanges of length 2
## start end width
## [1] 1 3 3
## [2] 2 4 3
我們使用data.frame對IRanges也進行一次處理石抡，這樣我們將更加清楚DataFrame與data.frame之間的差別：
df2 <- data.frame(iranges = ir)
df2
## iranges.start iranges.end iranges.width
## 1 1 3 3
## 2 2 4 3
對bed文件格式有所了解的朋友，應(yīng)該知道bed文件的第五列是對特定區(qū)域的一個打分（score）助泽。既然汁雷，GRanges和bed文件格式這么相似，那GRanges是否可以對每一行添加一個score呢报咳？答案是肯定的侠讯，我們通常使用values,elementMetadata或者mcols給GRanges賦值，這里以values方法為例
# 首先創(chuàng)建一個GRanges對象
gr <- GRanges(seqnames = "chr1", strand = c("+", "-", "+"),
ranges = IRanges(start = c(1,3,5), width = 3))
values(gr) <- DataFrame(score = c(0.1, 0.5, 0.3))
gr
## GRanges object with 3 ranges and 1 metadata column:
## seqnames ranges strand | score
## <Rle> <IRanges> <Rle> | <numeric>
## [1] chr1 [1, 3] + | 0.1
## [2] chr1 [3, 5] - | 0.5
## [3] chr1 [5, 7] + | 0.3
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
當(dāng)然了暑刃，如果你覺得這些方法太麻煩厢漩，這里還有一個更簡單的賦值方式，利用$符號：
gr$score
## [1] 0.1 0.5 0.3
gr$score2 = gr$score * 0.2
gr
## GRanges object with 3 ranges and 2 metadata columns: ## seqnames ranges strand | score score2
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 [1, 3] + | 0.1 0.02
## [2] chr1 [3, 5] - | 0.5 0.1
## [3] chr1 [5, 7] + | 0.3 0.06
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
上一篇bioconductor教程中岩臣，我們使用findOverlap來尋找不同的IRanges對象之間的overlap溜嗜。GRanges與IRanges相比，多了正負鏈的信息架谎，不過沒關(guān)系炸宵，我們依然可以使用findOverlap方法查找不同的GRanges對象的overlap。用法如下：
# 大家可以自己重新定義一個新的gr2對象
gr2 <- GRanges(seqnames = c("chr1", "chr2", "chr1"), strand = "*", ranges = IRanges(start = c(1, 3, 5), width = 3))
findOverlaps(gr, gr2)
## Hits object with 4 hits and 0 metadata columns:
## queryHits subjectHits
## <integer> <integer>
## [1] 1 1
## [2] 2 1
## [3] 2 3
## [4] 3 3
## -------
## queryLength: 3
## subjectLength: 3
這里值得一提的是谷扣，如果一個GRanges對象的鏈方向信息使用*號表示土全。這種情況捎琐，我們在尋找overlap時，它既能與另一個GRanges對象的正鏈重疊裹匙，也能與負鏈重疊瑞凑。此外，findOverlaps方法概页，還提供了一個ignore.strand參數(shù)籽御，在使用該參數(shù)時，我們在取overlap時將直接忽略鏈的方向信息:
findOverlaps(gr, gr2,ignore.strand=TRUE)
這里如果忘記了findOverlaps的輸出結(jié)果是什么惰匙，建議回顧上一篇bioconductor教程
還有一個方法subsetByOverlaps技掏，它能夠直接輸出我們想要得到的重疊的ranges：
subsetByOverlaps(gr, gr2)
## GRanges object with 3 ranges and 2 metadata columns:
## seqnames ranges strand | score score2
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 [1, 3] + | 0.1 0.02
## [2] chr1 [3, 5] - | 0.5 0.1
## [3] chr1 [5, 7] + | 0.3 0.06
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
講到這里，愛思考的你可能會問：我該如何將bed文件轉(zhuǎn)換為GRanges對象呢项鬼？總不能自己動手一行一行打出來吧哑梳？
這里介紹一個非常實用的方法：makeGRangesFromDataFrame，用法如下：
df <- data.frame(chr = "chr1", start = 1:3, end = 4:6, score = 7:9) makeGRangesFromDataFrame(df)
### 如果像保存score信息秃臣，需要添加一個參數(shù) keep.extra.columns
makeGRangesFromDataFrame(df, keep.extra.columns = TRUE)
我們常用的基因組都含有多個染色體涧衙，但很多情況下我們只想了解某一條染色體上的信息哪工，那該怎么辦呢奥此？
我們有三種方法可以實現(xiàn)，具體如下：
### 第一種
gr <- GRanges(seqnames = c("chr1", "chr2"), ranges = IRanges(start = 1:2, end = 4:5))
seqlevels(gr, force=TRUE) <- "chr1"
gr
### 第二種
gr <- GRanges(seqnames = c("chr1", "chr2"), ranges = IRanges(start = 1:2, end = 4:5))
dropSeqlevels(gr, "chr2")
### 第三種
keepSeqlevels(gr, "chr1")
其實雁比，我們從這三個方法的名稱就可以簡單知道其具體含義：
- seqlevels(gr, force=TRUE) <- "chr1" ：對gr對象本身進行操作稚虎，只留下chr1的區(qū)域
- dropSeqlevels(gr, "chr2")：去除chr2的區(qū)域，重新生成一個GRanges對象
- keepSeqlevels(gr, "chr1")：利用gr重新生成一個只有chr1區(qū)域的GRanges對象
最最最最最……最后偎捎，再介紹一個知識點蠢终。
許多情況下，不同數(shù)據(jù)庫的染色體序號的形式各不相同茴她，GRanges對象也提供了一個方法讓命名變得統(tǒng)一化寻拂，我知道你們懶得看了，我也懶得再細講了丈牢，簡單給一個代碼祭钉，好奇的同學(xué)可以自己了解一下：
gr <- GRanges(seqnames = "chr1", ranges = IRanges(start = 1:2, width = 2))
newStyle <- mapSeqlevels(seqlevels(gr), "NCBI")
gr <- renameSeqlevels(gr, newStyle)
講完結(jié)束，收工……
事了拂衣去己沛，深藏功與名
實在找不到合適的表情包了慌核，大家自行腦補畫面
