GenomicRanges 數(shù)據(jù)結(jié)構(gòu)

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

參考資料
An Introduction to the GenomicRanges Package

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末旬牲,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子搁吓,更是在濱河造成了極大的恐慌原茅,老刑警劉巖,帶你破解...
    沈念sama閱讀 218,122評論 6 505
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件堕仔,死亡現(xiàn)場離奇詭異擂橘,居然都是意外死亡,警方通過查閱死者的電腦和手機摩骨,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,070評論 3 395
  • 文/潘曉璐 我一進店門通贞,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人恼五,你說我怎么就攤上這事昌罩。” “怎么了灾馒?”我有些...
    開封第一講書人閱讀 164,491評論 0 354
  • 文/不壞的土叔 我叫張陵茎用,是天一觀的道長。 經(jīng)常有香客問我睬罗,道長轨功,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 58,636評論 1 293
  • 正文 為了忘掉前任容达,我火速辦了婚禮古涧,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘花盐。我一直安慰自己蒿褂,他們只是感情好,可當(dāng)我...
    茶點故事閱讀 67,676評論 6 392
  • 文/花漫 我一把揭開白布卒暂。 她就那樣靜靜地躺著啄栓,像睡著了一般。 火紅的嫁衣襯著肌膚如雪也祠。 梳的紋絲不亂的頭發(fā)上昙楚,一...
    開封第一講書人閱讀 51,541評論 1 305
  • 那天,我揣著相機與錄音诈嘿,去河邊找鬼堪旧。 笑死削葱,一個胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的淳梦。 我是一名探鬼主播析砸,決...
    沈念sama閱讀 40,292評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼爆袍!你這毒婦竟也來了首繁?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 39,211評論 0 276
  • 序言:老撾萬榮一對情侶失蹤陨囊,失蹤者是張志新(化名)和其女友劉穎弦疮,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體蜘醋,經(jīng)...
    沈念sama閱讀 45,655評論 1 314
  • 正文 獨居荒郊野嶺守林人離奇死亡胁塞,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,846評論 3 336
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了压语。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片啸罢。...
    茶點故事閱讀 39,965評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖胎食,靈堂內(nèi)的尸體忽然破棺而出伺糠,到底是詐尸還是另有隱情,我是刑警寧澤斥季,帶...
    沈念sama閱讀 35,684評論 5 347
  • 正文 年R本政府宣布训桶,位于F島的核電站,受9級特大地震影響酣倾,放射性物質(zhì)發(fā)生泄漏舵揭。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 41,295評論 3 329
  • 文/蒙蒙 一躁锡、第九天 我趴在偏房一處隱蔽的房頂上張望午绳。 院中可真熱鬧,春花似錦映之、人聲如沸拦焚。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,894評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽赎败。三九已至,卻和暖如春蠢甲,著一層夾襖步出監(jiān)牢的瞬間僵刮,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,012評論 1 269
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留搞糕,地道東北人勇吊。 一個月前我還...
    沈念sama閱讀 48,126評論 3 370
  • 正文 我出身青樓,卻偏偏與公主長得像窍仰,于是被迫代替她去往敵國和親汉规。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 44,914評論 2 355

推薦閱讀更多精彩內(nèi)容