上一期讲到了GRanges的特点，这一期我们讲讲其具体用法。
首先我们需要引入一个新的函数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]
它也可以使用$符号，获得某一个变量的数值，如：
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文件格式有所了解的朋友，应该知道bed文件的第五列是对特定区域的一个打分（score）。既然，GRanges和bed文件格式这么相似，那GRanges是否可以对每一行添加一个score呢？答案是肯定的，我们通常使用values,elementMetadata或者mcols给GRanges赋值，这里以values方法为例
# 首先创建一个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
当然了，如果你觉得这些方法太麻烦，这里还有一个更简单的赋值方式，利用$符号：
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相比，多了正负链的信息，不过没关系，我们依然可以使用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参数，在使用该参数时，我们在取overlap时将直接忽略链的方向信息:
findOverlaps(gr, gr2,ignore.strand=TRUE)
这里如果忘记了findOverlaps的输出结果是什么，建议回顾上一篇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文件转换为GRanges对象呢？总不能自己动手一行一行打出来吧？
这里介绍一个非常实用的方法：makeGRangesFromDataFrame，用法如下：
df <- data.frame(chr = "chr1", start = 1:3, end = 4:6, score = 7:9) makeGRangesFromDataFrame(df)
### 如果像保存score信息，需要添加一个参数 keep.extra.columns
makeGRangesFromDataFrame(df, keep.extra.columns = TRUE)
我们常用的基因组都含有多个染色体，但很多情况下我们只想了解某一条染色体上的信息，那该怎么办呢？
我们有三种方法可以实现，具体如下：
### 第一种
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的区域
- dropSeqlevels(gr, "chr2")：去除chr2的区域，重新生成一个GRanges对象
- keepSeqlevels(gr, "chr1")：利用gr重新生成一个只有chr1区域的GRanges对象
最最最最最……最后，再介绍一个知识点。
许多情况下，不同数据库的染色体序号的形式各不相同，GRanges对象也提供了一个方法让命名变得统一化，我知道你们懒得看了，我也懒得再细讲了，简单给一个代码，好奇的同学可以自己了解一下：
gr <- GRanges(seqnames = "chr1", ranges = IRanges(start = 1:2, width = 2))
newStyle <- mapSeqlevels(seqlevels(gr), "NCBI")
gr <- renameSeqlevels(gr, newStyle)
讲完结束，收工……
事了拂衣去，深藏功与名
实在找不到合适的表情包了，大家自行脑补画面
网友评论