R 包 GenomicRanges 是 Bioconductor 項目承載基因組位置信息的基礎(chǔ)包租副,為諸如 BSgenome,rtracklayer,VariantAnnotation 及其他 R 包提供了基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)倒槐。
GenomicRanges 包含了 3 類對象:GRanges,GPos,GRangesList. 本文主要討論 GRanges 結(jié)構(gòu)和操作函數(shù)。
GRanges 對象是基因組范圍(Genomic Ranges)的集合附井,每個基因組范圍有一個在基因組的起點和終點位置讨越。GRanges 可用于存儲諸如轉(zhuǎn)錄本、外顯子等基因組特征永毅。
下面代碼用 genes
函數(shù)提取 hg38 基因。
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
hg38Gene <- genes(txdb)
用 class
可以查看它屬于 GRanges 對象沼死。
> class(hg38Gene)
[1] "GRanges"
attr(,"package")
[1] "GenomicRanges"
> show(hg38Gene)
GRanges object with 25750 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58345178-58362751 - | 1
10 chr8 18391282-18401218 + | 10
100 chr20 44619522-44652233 - | 100
1000 chr18 27950966-28177130 - | 1000
100009613 chr11 70072434-70075433 - | 100009613
... ... ... ... . ...
9991 chr9 112217716-112333667 - | 9991
9992 chr21 34364024-34371389 + | 9992
9993 chr22 19036282-19122454 - | 9993
9994 chr6 89829894-89874436 + | 9994
9997 chr22 50523568-50526461 - | 9997
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
用 show
函數(shù)顯示 GRanges 結(jié)構(gòu)着逐,可見用豎線分為 2 部分,左邊為基因組座標(biāo)信息意蛀,包含 seqnames,ranges 等耸别;右邊為 metadata 存儲相關(guān)聯(lián)的信息,具備非常高的可拓展性县钥。
左邊信息根據(jù)列名秀姐,GenomicRanges 定義了相應(yīng)函數(shù)提取詳細(xì)信息。
> seqnames(hg38Gene)
factor-Rle of length 25750 with 19826 runs
Lengths: 1 1 1 1 1 1 ... 1 1 1 1 1
Values : chr19 chr8 chr20 chr18 chr11 chr10 ... chr9 chr21 chr22 chr6 chr22
Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
> ranges(hg38Gene)
IRanges object with 25750 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
1 58345178 58362751 17574
10 18391282 18401218 9937
100 44619522 44652233 32712
1000 27950966 28177130 226165
100009613 70072434 70075433 3000
... ... ... ...
9991 112217716 112333667 115952
9992 34364024 34371389 7366
9993 19036282 19122454 86173
9994 89829894 89874436 44543
9997 50523568 50526461 2894
用 granges
函數(shù)提取左邊信息若贮。
> granges(hg38Gene)
GRanges object with 25750 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
1 chr19 58345178-58362751 -
10 chr8 18391282-18401218 +
100 chr20 44619522-44652233 -
1000 chr18 27950966-28177130 -
100009613 chr11 70072434-70075433 -
... ... ... ...
9991 chr9 112217716-112333667 -
9992 chr21 34364024-34371389 +
9993 chr22 19036282-19122454 -
9994 chr6 89829894-89874436 +
9997 chr22 50523568-50526461 -
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
用 mcols
函數(shù)提取 metadata 信息省有。
> mcols(hg38Gene)
DataFrame with 25750 rows and 1 column
gene_id
<character>
1 1
10 10
100 100
1000 1000
100009613 100009613
... ...
9991 9991
9992 9992
9993 9993
9994 9994
9997 9997
也可以修改 metadata, 比如:
> mcols(hg38Gene)$test <- "something"
> show(hg38Gene)
GRanges object with 25750 ranges and 2 metadata columns:
seqnames ranges strand | gene_id test
<Rle> <IRanges> <Rle> | <character> <character>
1 chr19 58345178-58362751 - | 1 something
10 chr8 18391282-18401218 + | 10 something
100 chr20 44619522-44652233 - | 100 something
1000 chr18 27950966-28177130 - | 1000 something
100009613 chr11 70072434-70075433 - | 100009613 something
... ... ... ... . ... ...
9991 chr9 112217716-112333667 - | 9991 something
9992 chr21 34364024-34371389 + | 9992 something
9993 chr22 19036282-19122454 - | 9993 something
9994 chr6 89829894-89874436 + | 9994 something
9997 chr22 50523568-50526461 - | 9997 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
GRanges 可以取子集,類似于向量進行切片操作谴麦。
> hg38Gene[1:3]
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | gene_id test
<Rle> <IRanges> <Rle> | <character> <character>
1 chr19 58345178-58362751 - | 1 something
10 chr8 18391282-18401218 + | 10 something
100 chr20 44619522-44652233 - | 100 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
# 指定需要的 metadata 列
> hg38Gene[1:3, "gene_id"]
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58345178-58362751 - | 1
10 chr8 18391282-18401218 + | 10
100 chr20 44619522-44652233 - | 100
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
而用 split
函數(shù)可以將 GRanges 分割為多個 GRanges 對象蠢沿,從而組織成 GRangesList 對象。
> sp <- split(hg38Gene, rep_len(1:2, length.out = 25750))
> show(sp)
GRangesList object of length 2:
$`1`
GRanges object with 12875 ranges and 2 metadata columns:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58345178-58362751 - | 1
100 chr20 44619522-44652233 - | 100
100009613 chr11 70072434-70075433 - | 100009613
100009676 chr3 101676475-101679217 + | 100009676
10002 chr15 71792638-71818259 + | 10002
... ... ... ... . ...
9987 chr4 82422564-82430408 - | 9987
9989 chr18 9546791-9615240 - | 9989
9990 chr15 34229784-34338060 - | 9990
9992 chr21 34364024-34371389 + | 9992
9994 chr6 89829894-89874436 + | 9994
test
<character>
1 something
100 something
100009613 something
100009676 something
10002 something
... ...
9987 something
9989 something
9990 something
9992 something
9994 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
$`2`
GRanges object with 12875 ranges and 2 metadata columns:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
10 chr8 18391282-18401218 + | 10
1000 chr18 27950966-28177130 - | 1000
100009667 chr10 68010205-68010862 - | 100009667
10001 chr14 70581257-70641204 - | 10001
100033411 chr16 75693893-75701461 - | 100033411
... ... ... ... . ...
9988 chr7 87152361-87196337 + | 9988
999 chr16 68737292-68835541 + | 999
9991 chr9 112217716-112333667 - | 9991
9993 chr22 19036282-19122454 - | 9993
9997 chr22 50523568-50526461 - | 9997
test
<character>
10 something
1000 something
100009667 something
10001 something
100033411 something
... ...
9988 something
999 something
9991 something
9993 something
9997 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
多個 GRanges 對象可以用 c
或者 append
函數(shù)合并匾效。
> c(sp[[1]], sp[[2]])
GRanges object with 25750 ranges and 2 metadata columns:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58345178-58362751 - | 1
100 chr20 44619522-44652233 - | 100
100009613 chr11 70072434-70075433 - | 100009613
100009676 chr3 101676475-101679217 + | 100009676
10002 chr15 71792638-71818259 + | 10002
... ... ... ... . ...
9988 chr7 87152361-87196337 + | 9988
999 chr16 68737292-68835541 + | 999
9991 chr9 112217716-112333667 - | 9991
9993 chr22 19036282-19122454 - | 9993
9997 chr22 50523568-50526461 - | 9997
test
<character>
1 something
100 something
100009613 something
100009676 something
10002 something
... ...
9988 something
999 something
9991 something
9993 something
9997 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
對基因區(qū)域記錄的修改分為 intra-range methods,inter-range methods,between-range methods 三類舷蟀。intra-range methods 單獨修改每一個區(qū)域記錄,有 flank,resize,shift 等函數(shù)。比如用 shift
函數(shù)移動每個區(qū)域 100bp.
> shift(hg38Gene, 100)
GRanges object with 25750 ranges and 2 metadata columns:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58345278-58362851 - | 1
10 chr8 18391382-18401318 + | 10
100 chr20 44619622-44652333 - | 100
1000 chr18 27951066-28177230 - | 1000
100009613 chr11 70072534-70075533 - | 100009613
... ... ... ... . ...
9991 chr9 112217816-112333767 - | 9991
9992 chr21 34364124-34371489 + | 9992
9993 chr22 19036382-19122554 - | 9993
9994 chr6 89829994-89874536 + | 9994
9997 chr22 50523668-50526561 - | 9997
test
<character>
1 something
10 something
100 something
1000 something
100009613 something
... ...
9991 something
9992 something
9993 something
9994 something
9997 something
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
Inter-range methods 涉及多個區(qū)域記錄的比較野宜,包含 reduce,gaps,coverage 等函數(shù)碗殷。如 coverage
函數(shù)統(tǒng)計區(qū)域間重疊情況。
> coverage(hg38Gene)
RleList of length 595
$chr1
integer-Rle of length 248956422 with 5195 runs
Lengths: 11868 2535 6 15161 ... 47051 13751 36476
Values : 0 1 2 1 ... 0 1 0
$chr2
integer-Rle of length 242193529 with 3397 runs
Lengths: 38813 8057 171265 46731 ... 5123 34215 424713
Values : 0 1 0 1 ... 0 1 0
$chr3
integer-Rle of length 198295559 with 2954 runs
Lengths: 11744 13105 84857 87056 ... 914 4486 251839
Values : 0 1 0 1 ... 2 1 0
$chr4
integer-Rle of length 190214555 with 2071 runs
Lengths: 53285 34923 36292 ... 1278990 3363 146691
Values : 0 1 0 ... 0 1 0
$chr5
integer-Rle of length 181538259 with 2377 runs
Lengths: 92150 97822 1522 4840 ... 94960 995 169997
Values : 0 1 0 1 ... 0 1 0
...
<590 more elements>
Between-range methods 是多個 GRanges 對象間操作速缨,包含 union,intersect,setdiff 等函數(shù)。如 intersect 取 2 個 GRanges 交集代乃。
> hg38Gene2 <- hg38Gene[3:5]
> intersect(hg38Gene, hg38Gene2)
GRanges object with 3 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr11 70072434-70075433 -
[2] chr18 27950966-28177130 -
[3] chr20 44619522-44652233 -
-------
seqinfo: 595 sequences (1 circular) from hg38 genome