- 参考资料：
https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html
https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GRanges_and_GRangesList_slides.pdf
一、GenomicRanges包主要作用
- 通过创造一个
GRanges对象，storing a set of genomic ranges，来描述一段段序列的基本特征。 - 主要由以下两部分构成：
1、basic info (genomic coordinates)
（1）包含a chromosome name(染色体名字), a start, an end, and a strand.四列描述；
（2）主要基于IRanges这个基础包提供序列的头尾信息
library(IRanges)
test <- testIRanges(16:21, 20)
#返回start、end、width三列信息
#可分别用三个函数去取对应的信息
start(test)
#[1] 16 17 18 19 20 21
end(test)
#[1] 20 20 20 20 20 20
width(test)
#[1] 5 4 3 2 1 0
#其实GenomicRanges相当部分的函数来自IRanges包
如下图，可知start and end are both 1-based positions relative to the 5’ end of the plus strand of the chromosome。
The width of the range is the number of genomic positions included in it. So width = end - start + 1.
IRanges
2、metadata columns (annotation)
- 用来存储序列的来源，名字，分类等个性化特征信息
- 可有可无，一般在
GRanges对象中，以|相隔
二、GRanges对象及其基本操作
library(GenomicRanges)
gr1 <- GRanges(seqnames=Rle(c("ch1", "chMT"), c(2, 4)),
ranges=IRanges(16:21, 20),
strand=rep(c("+", "-", "*"), 2))
gr1
gr1
1、seqnames（所在染色体）
seqnames(gr1)
seqlevels(gr1)
#若想修改染色体名字，直接修改level即可
seqlevels(gr1) <- c("a","b")
as.vector(seqnames(gr1))
seqlevels(gr1) <- c("ch1","chMT")
as.vector(seqnames(gr1))
此外还可以添加染色体的长度信息于对象
seqlengths(gr1)
seqlengths(gr1) <- c(50000, 800)
seqlengths(gr1)
2、start end操作同前
3、strand 链的方向
描述该序列在DNA的哪一条链上，正链(5'→3')还是负链(3'→5')
正负链
strand(gr1)
注意正负链在后面的进阶操作时影响很大，要记住的是start and end 的位置都是相对于所在链的5'端所说的。
4、序列命名
names(gr1)
# NULL
names(gr1) <- LETTERS[1:6]
names(gr1)
# [1] "A" "B" "C" "D" "E" "F"
5、meta column
- 如前所述，
GRanges对象包含两大部分，其中第二部分为序列的meta信息，可有可无。 - 在上例中就没有，可利用
mcols函数进行查找、添加
gr1
mcols(gr1)
mcols(gr1) <- DataFrame(score=11:16, GC=seq(1, 0, length=6))
gr1 # | 分隔
meta column
二、GRanges对象之增删改查
相关操作基本可以类比data.frame的操作
gr1[c("F", "A")] #根据序列名取子集
gr1[strand(gr1) == "+"] #筛选正链的子集
gr1 <- gr1[-5] #删除第五列
gr2 <- GRanges(seqnames="ch2",
ranges=IRanges(start=c(2:1,2), width=6),
score=15:13,
GC=seq(0, 0.4, length=3))
gr12 <- c(gr1, gr2) #合并两个对象
gr12[length(gr12)] == gr12
duplicated(gr12)
unique(gr12)
sort(gr12)
三、GRanges对象之range-based operations
gr1
gr1
-
shift序列平移
shift(gr1, 50)
shift(gr1, 50)
-
resize修改长度
resize(gr1, 12)
#首先序列的长度都变为12
#其次无论正负链，都是左边的位点不同，改变右边位点的位置
resize(gr1, 12)
-
flank取序列上下游
#默认取序列上游部分，即5'端
flank(gr1, 3)
#通过设置start = F参数，可用于取下游部分
flank(gr1, 3)
-
range概括序列跨越的最大区域，注意正负链
gr3 <- shift(gr1, c(35000, rep(0, 3), 100))
range(gr3)
gr3，range(gr3)
-
reduce概括序列覆盖的区域，注意正负链
#注意对比与上面range的区别
reduce(gr3)
reduce(gr3)
-
gaps,与reduce相反，结果为后者的补集
gaps(gr3)
gaps(gr3)
-
disjoin不重复的最小片段，有可能要切下
disjoin(gr3)
disjoin(gr3)
四、GRanges对象之覆盖深度
cvg1
seqinfo(gr1)
-
coverage以染色体为界，统计所有碱基的覆盖情况
cvg1 <- coverage(gr1)
cvg1
cvg1
mean(cvg12) #染色体碱基平均覆盖深度
max(cvg12) #染色体碱基最大覆盖深度
-
slice查询指定覆盖度以上的碱基
sl1 <- slice(cvg1, lower=1)
sl1
#最低覆盖1次的，中括号里的为实际覆盖的次数
sl1
存疑：好像不考虑正负链的差异，直接根据start与end来考虑覆盖深度，我觉得应该不可以这样计算吧
四、GRanges对象之overlap
- 准备示例数据
#安装相关包
#BiocManager::install("pasillaBamSubset")
#BiocManager::install("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
library(pasillaBamSubset)
library(GenomicAlignments)
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
#示例测序reads数据
untreated1_chr4()
reads <- readGAlignments(untreated1_chr4())
reads <- as(reads, "GRanges")
reads[1:4]
#示例参考基因组数据
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
dm3_genes <- genes(txdb)
- 统计overlap情况
hits <- findOverlaps(reads, dm3_genes) #注意正负链
hits
#默认只要有交集，就算overlap；可通过
#可以通过调整参数，设置判断标准
?findOverlaps
hits
#举例查看下
reads[6298]
reads[6296]
dm3_genes[11499]
举例
四、GRanges对象之GRangesList
1、将不同GRanges对象组成list
grl <- GRangesList(gr3, gr2)
#最一开始的基础函数也同样适用于list，只是结果也是list的格式
names(grl) <- c("TX1", "TX2")
mcols(grl)$geneid <- c("GENE1", "GENE2")
2、将某一GRanges对象按照染色体，拆为GRangesList
split(gr12, seqnames(gr12)) #按seq
以此类推，用
split函数还可以按照GRanges对象内其它元素进行拆解
