劉小澤寫于19.4.2
主要看看IRanges的基本用法
上次寫了關(guān)于RangeData的介紹：http://www.itdecent.cn/p/9150213743d7
總的來說，主要會(huì)利用5個(gè)R包，來處理關(guān)于基因組區(qū)間（Genomic Range）的信息
IRanges、GenomicRanges：顧名思義，就是處理基因組區(qū)間信息 【GenomicRanges可以將染色體坐標(biāo)范圍、序列名稱和正負(fù)鏈信息結(jié)合起來。許多Bioconductor包都高度依賴IRanges和GenomicRanges提供的底層數(shù)據(jù)結(jié)構(gòu)】
-
GenomicFeatures：處理基因組上的gene model或其他序列特征（gene feature）信息(比如：genes、exons、UTRs、transcripts)
A Gene Model is defined as any description of a gene product from a variety of sources including computational prediction, mRNA sequencing, or genetic characterization.
The gene feature is meant to approximately cover the region of nucleic acid considered by workers in the field to be the gene. Biostrings、BSgenome：處理序列，比如提取子集【Biostrings包含了序列比對(duì)、模式匹配等基本序列分析函數(shù)；BSgenome針對(duì)有注釋的全基因組數(shù)據(jù)進(jìn)行操作】
rtracklayer：讀取像BED、GTF、WIG這樣的文件 【rtracklayer包將數(shù)據(jù)導(dǎo)入U(xiǎn)SCS基因組瀏覽器(http://genome-asia.ucsc.edu/)進(jìn)行瀏覽、操作、輸出】
Before diving into working with genomic ranges, we’re going to get our feet wet with generic ranges (i.e., ranges that represent a contiguous subsequence of elements over any type of sequence)
從簡單的基因區(qū)間入手，就像批量跑流程一樣，由小及大，逐漸培養(yǎng)區(qū)間思維"range thinking"
利用IRanges存儲(chǔ)基因區(qū)間信息
IRanges包是GenomicRanges的一個(gè)依賴包，它會(huì)創(chuàng)建一個(gè)同名的對(duì)象，有兩個(gè)基本的組成：起始和終止位點(diǎn)
# 例如：創(chuàng)建一個(gè)從4-13的區(qū)間
> rng <- IRanges(start = 4,end = 13)
> rng
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 4 13 10
IRanges和GenomicRanges的坐標(biāo)都是1-based系統(tǒng)，也可以給定任意的起始或終止坐標(biāo)，并給出區(qū)間長度來創(chuàng)建
# 給定一個(gè)起始位置
> IRanges(start=4, width=3)
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 4 6 3
# 給定一個(gè)終止位置
> IRanges(end=5, width=5)
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 5 5
除了指定一個(gè)起始、終止位點(diǎn)，還可以使用向量創(chuàng)建多個(gè)區(qū)間：
> x <- IRanges(start=c(4, 7, 2, 20), end=c(13, 7, 5, 23))
> x
IRanges object with 4 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 4 13 10
[2] 7 7 1
[3] 2 5 4
[4] 20 23 4
既然是一個(gè)range對(duì)象，那么其中的每個(gè)range都可以指定一個(gè)名稱，使用names 函數(shù)
> names(x) <- paste0("gene",1:4)
> x
IRanges object with 4 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
gene1 4 17 14
gene2 7 11 5
gene3 2 9 8
gene4 20 27 8
對(duì)于一個(gè)對(duì)象，最直接的查看它的內(nèi)容就是通過str()[就像使用levels()查看因子型變量內(nèi)容一樣]
> str(x)
Formal class 'IRanges' [package "IRanges"] with 6 slots
..@ start : int [1:4] 4 7 2 20
..@ width : int [1:4] 14 5 8 8
..@ NAMES : chr [1:4] "gene1" "gene2" "gene3" "gene4"
..@ elementType : chr "ANY"
..@ elementMetadata: NULL
..@ metadata : list()
有了起始、終止、區(qū)間，就可以抽出其中的任意部分，或者對(duì)其進(jìn)行操作[操作過程和數(shù)據(jù)框很像]：
# 找到起始位點(diǎn)
> start(x)
[1] 4 7 2 20
# 找到終止位點(diǎn)
> end(x)
[1] 13 7 5 23
# # 找到區(qū)間長度
> width(x)
[1] 10 1 4 4
# 將所有的終止位點(diǎn)增加4
> end(x) <- end(x) + 4
> x
IRanges object with 4 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
gene1 4 17 14
gene2 7 11 5
gene3 2 9 8
gene4 20 27 8
# 使用range()可以得到總區(qū)間
> range(x)
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 2 27 26
# 找到起始位點(diǎn)小于5的
> x[start(x) < 5]
IRanges object with 2 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
gene1 4 17 14
gene3 2 9 8
# 找到長度大于8的基因
> x[width(x) > 8]
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
gene1 4 17 14
# 查看gene3的長度信息
> x['gene3']
IRanges object with 1 range and 0 metadata columns:
start end width
<integer> <integer> <integer>
gene3 2 9 8
# 還可以輕松組合(merge)
> a <- IRanges(start=7, width=4)
> b <- IRanges(start=2, end=5)
> c(a,b)
IRanges object with 2 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 7 10 4
[2] 2 5 4
上面可以看到，IRanges和數(shù)據(jù)框也差不了多少，而且生成的方式好像比數(shù)據(jù)框還要復(fù)雜一點(diǎn)，"存在即合理"，既然一開始構(gòu)建稍微麻煩一點(diǎn)，那么肯定他能做的事情也更多
IRanges結(jié)合基本運(yùn)算、轉(zhuǎn)換
-
改變序列兩側(cè)范圍：比如想在一個(gè)編碼區(qū)域上、下游加一段序列
# 現(xiàn)有的編碼區(qū)假設(shè)是這樣 > x <- IRanges(start=c(40, 80), end=c(67, 114)) > x IRanges object with 2 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 40 67 28 [2] 80 114 35 # 上、下游同時(shí)增加4bp > x+4L IRanges object with 2 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 36 71 36 [2] 76 118 43 # 看一下變化：start和end都同時(shí)向自己的方向增加了4bp -
改變中間：有時(shí)我們主要想保留序列中間的某一段區(qū)域，而不是增加/縮減兩端序列
# 原始序列范圍 > y <- IRanges(start=c(4, 6, 10, 12), width=13) > y IRanges object with 4 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 4 16 13 [2] 6 18 13 [3] 10 22 13 [4] 12 24 13 # 現(xiàn)在只想保留y的所有基因的5-13bp區(qū)域(因?yàn)槊總€(gè)基因的起始、終止都不同，如果基因起始位點(diǎn)在5之前的，截取后就從5開始；如果起始位置在5之后的，截取后就從現(xiàn)在的位置開始；終止位置也是如此) > restrict(y,5,13) IRanges object with 4 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 5 13 9 [2] 6 13 8 [3] 10 13 4 [4] 12 13 2 -
獲得兩端部分(flank)，比如想得到左側(cè)：轉(zhuǎn)錄起始位點(diǎn)( transition start site，TSS)；右側(cè)：轉(zhuǎn)錄終止位點(diǎn)(transcription termination site， TTS)
> x IRanges object with 2 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 40 67 28 [2] 80 114 35 # 現(xiàn)在要獲得左側(cè)7bp的起始位點(diǎn)(flank默認(rèn)是計(jì)算上游) > flank(x,width = 7) IRanges object with 2 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 33 39 7 [2] 73 79 7 # 如果要計(jì)算下游的終止位點(diǎn)(就要設(shè)置start=FALSE) -
如果ranges對(duì)象中包含多個(gè)基因，并且它們的區(qū)間存在重疊，比如基因1是10-15bp，gene2是15-18bp，那么我們想看總體覆蓋的區(qū)間，也就是10-18bp。這時(shí)使用
reduce函數(shù)，將重疊的區(qū)域進(jìn)行縮減#模擬數(shù)據(jù)(隨機(jī)產(chǎn)生20個(gè)起始位點(diǎn)，然后區(qū)間長度為5) > set.seed(12) > alns <- IRanges(start=sample(seq_len(50), 20), width=5) > head(alns, 4) IRanges object with 4 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 4 8 5 [2] 41 45 5 [3] 46 50 5 [4] 13 17 5 # 將重疊區(qū)域縮減，最后的總覆蓋是1-26和28-54 > reduce(alns) IRanges object with 2 ranges and 0 metadata columns: start end width <integer> <integer> <integer> [1] 1 26 26 [2] 28 54 27 -
上面已經(jīng)將多個(gè)區(qū)域縮減到了幾個(gè)大區(qū)域，怎么找這幾個(gè)大區(qū)域的位置呢？就好像通過基因區(qū)域找到基因間區(qū)一樣，利用
gap函數(shù)得到（從前一個(gè)大區(qū)域的結(jié)尾到后一個(gè)大區(qū)域的開頭，這算一個(gè)gap）> gaps(alns) IRanges object with 1 range and 0 metadata columns: start end width <integer> <integer> <integer> [1] 27 27 1 # 因?yàn)樯厦嬷坏玫絻山M大區(qū)域(可以想象成兩個(gè)基因區(qū)域)，然后找基因間區(qū)，就是27bp這個(gè)位置了 -
可以尋找兩個(gè)IRange對(duì)象的交集、補(bǔ)集、并集
> a <- IRanges(start=4, end=13) > b <- IRanges(start=12, end=17) # 求a、b交集 > intersect(a, b) IRanges object with 1 range and 0 metadata columns: start end width <integer> <integer> <integer> [1] 12 13 2 # 求a對(duì)b的補(bǔ)集(就是a中有b中沒有的部分)；顛倒參數(shù)就是求b對(duì) # a補(bǔ)集 > setdiff(a,b) IRanges object with 1 range and 0 metadata columns: start end width <integer> <integer> <integer> [1] 4 11 8 # 并集是union
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
