豆豆寫于19.5.28
跟著Bioconductor書再看一遍他們是如何介紹GenomicRanges的，不同的作者不同的見解，也有不同的教學(xué)方法。多學(xué)多看，這一部分內(nèi)容來自Bioconductor 2018的第四章
Bioconductor沒想象的那么簡單(part1)
Bioconductor沒想象的那么簡單(part2)
Bioconductor沒想象的那么簡單(part3)
Bioconductor沒想象的那么簡單-(part4)
Bioconductor沒想象的那么簡單-(part5)
Bioconductor沒想象的那么簡單(part6)
Bioconductor沒想象的那么簡單-part7
Bioconductor沒想象的那么簡單-part8
前言
什么是Ranges基礎(chǔ)？
它就是處理組學(xué)數(shù)據(jù)使用的，包含了數(shù)據(jù)結(jié)構(gòu)和算法。它包含了標(biāo)準(zhǔn)的基因組數(shù)據(jù)'容器“(像GRanges和SummarizedExperiment)，優(yōu)化的數(shù)據(jù)重現(xiàn)方式(Rle)，還有快速的算法(像計(jì)算overlap、找近鄰位點(diǎn)、統(tǒng)計(jì)區(qū)間信息)?？梢钥吹剑@是一個(gè)很大的整體，其中我們 之前接觸到的GRanges只是其中一小部分。
為什么要學(xué)Ranges基礎(chǔ)？
因?yàn)閿?shù)百個(gè)R包都是建立在Ranges基礎(chǔ)上的，它可以方便整合多個(gè)包、多種數(shù)據(jù)類型以構(gòu)建一個(gè)流程。代碼拓展性強(qiáng)(可以從流程腳本向R包腳本轉(zhuǎn)換)。
入門——GenomicRanges
核心的數(shù)據(jù)結(jié)構(gòu)就是GenomicRanges，它包含了許多代表基因組位置的信息，一般是起始和終止信息(當(dāng)然這個(gè)位置信息的重合是少不了的，這會給統(tǒng)計(jì)一些基礎(chǔ)信息帶來困擾，這個(gè)問題之后再說)。除了這個(gè)，還可以包含結(jié)合位點(diǎn)、序列比對、轉(zhuǎn)錄本信息等。
實(shí)戰(zhàn)
生成
最簡單的生成GRanges方法應(yīng)該就是數(shù)據(jù)框轉(zhuǎn)換了
library(GenomicRanges)
df <- data.frame(
seqnames = rep(c('chr1', 'chr2', 'chr1', 'chr3'), c(1, 3, 2, 4)), start = c(101, 105, 125, 132, 134, 152, 153, 160, 166, 170),
end = c(104, 120, 133, 132, 155, 154, 159, 166, 171, 190),
strand = rep(strand(c('-', '+', '*', '+', '-')), c(1, 2, 2, 3, 2)), score = 1:10,
GC = seq(1, 0, length=10),
row.names = head(letters, 10))
View(df)
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)
這樣的結(jié)果就是建立了包含10個(gè)ranges的GRanges對象
之前介紹過這里，以|為分隔，左邊是基因組坐標(biāo)信息(序列名、范圍、鏈)，右邊是其他的注釋信息(這里是score和GC含量，當(dāng)然它可以容納更多類型的數(shù)據(jù))
也可以從標(biāo)準(zhǔn)文件格式(如：BED、GFF、BigWig)導(dǎo)入為GRanges，利用rtracklayer包，part7：https://www.jianshu.com/p/37448531b533有過介紹
另外，利用GenomicAlignments包可以利用BAM文件構(gòu)建GAlignments對象，這個(gè)對象和GRanges很像，當(dāng)然也可以轉(zhuǎn)換成GRanges
查看
對左側(cè)基本信息操作
# 先看三大項(xiàng)(名稱、區(qū)間、鏈)
seqnames(gr)
ranges(gr)
strand(gr)
# 另外，可以同時(shí)看這三項(xiàng)，并排除其余的metadata(也就是只返回豎線左側(cè)的內(nèi)容)
granges(gr)
# 如果要分別查看起始、終止、區(qū)間位置(也就是把上面ranges函數(shù)的結(jié)果拆開)
start(gr)
end(gr)
width(gr)
對右側(cè)注釋信息操作
mcols(gr) #將metadata返回為數(shù)據(jù)框，當(dāng)然可以直接從結(jié)果中選取某一列
mcols(gr)$score
#上面結(jié)果等同于 score(gr)
總體操作
# 看行名
names(gr)
# 看多少行
length(gr)
取子集
這個(gè)對象取子集的方法和向量的操作類似，可以直接gr[2:3]
取出來是這樣的：
然后我們可以指定后面取出哪個(gè)metadata：
> gr[2:3, 'GC']
GRanges object with 2 ranges and 1 metadata column:
seqnames ranges strand | GC
<Rle> <IRanges> <Rle> | <numeric>
b chr2 105-120 + | 0.888888888888889
c chr2 125-133 + | 0.777777777777778
當(dāng)然這只是一種最基礎(chǔ)的辦法，更方便的辦法如：subset()
> subset(gr, strand == '+' & score > 5, select = score)
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | score
<Rle> <IRanges> <Rle> | <integer>
f chr1 152-154 + | 6
g chr3 153-159 + | 7
h chr3 160-166 + | 8
對GRanges對象的修改
比如可以復(fù)制并替換一行
grMod <- gr
grMod[2] <- gr[1]
# 結(jié)果新的grMod對象的第二行就和第一行一樣
還可以對行進(jìn)行重復(fù)、反轉(zhuǎn)、選特定行
rep(gr[2], times = 3) # 重復(fù)第二行3次
rev(gr) #將原來的所有行反轉(zhuǎn)
head(gr,n=2) # 從頭取
tail(gr,n=2) # 從尾端取
window(gr, start=2,end=4) # 從中間連續(xù)取
gr[IRanges(start=c(2,7),end=c(3,9))] # 從中間不連續(xù)取(也就是取出2~3行 + 7~9行)
拆分與合并
拆分用split()，結(jié)果返回一個(gè)列表GRangesList 【很好記：對象下面是列表，列表下面是數(shù)據(jù)框，數(shù)據(jù)框下面是向量】。只不過這里的列表中包含的還是小對象(而不是一般意義上的數(shù)據(jù)框)
首先最直觀的是按行拆分
sp <- split(gr, rep(1:2, each=5))
sp
# 這個(gè)意思就是說：將原來的對象按rep(1:2, each=5)這個(gè)格式拆，而這個(gè)格式是1 1 1 1 1 2 2 2 2 2一共10個(gè)數(shù)，因此這里也就明白，只要我們告訴split函數(shù)我們要拆分的結(jié)果是怎么組合的，就可以想怎么拆就怎么拆
另外除了按行，還可以按注釋信息metadata拆分
split(gr, ~ strand)
# 這個(gè)結(jié)果就分成了正鏈、負(fù)鏈和未知鏈三個(gè)
合并用c()、append()
# 因?yàn)橹安鸱值慕Y(jié)果是一個(gè)列表，那么選擇其中的元素還是要用[[]]
c(sp[[1]], sp[[2]])
另外stack()函數(shù)可以將GRangesList 列表中的小對象重新組裝到一起，并可以為每個(gè)小對象添加一個(gè)分組信息
stack(sp, index.var='group')
整合統(tǒng)計(jì)
利用aggregate函數(shù)，其中要指定一個(gè)公式y ~ x (x表示按什么分組，y表示分組后每組的值)
有兩種方法進(jìn)行統(tǒng)計(jì)：
# 第一種
aggregate(gr, score ~ strand, mean)
# 第二種
aggregate(gr, ~ strand, n_score = lengths(score), mean_score = mean(score))
# 注意到其中統(tǒng)計(jì)長度用的是lengths而不是length，因?yàn)檫@里是先按照strand進(jìn)行分組，然后再統(tǒng)計(jì)score長度與均值信息。因此score的信息是類似于一個(gè)列表，其中有正鏈的score信息、負(fù)鏈的score信息以及未知鏈的score信息，如果我們對其中任一個(gè)score信息進(jìn)行長度統(tǒng)計(jì)，是可以用length()的，但是放在一起統(tǒng)計(jì)，就要用lengths()
