這是我聽B站鯪魚不會飛視頻（R與Bioconductor的入門課）里的筆記哦~
Bioconductor是什么示弓？
Bioconductor是一個專門做生信R包的平臺次企，可以把它看成一個R工具包管理組織；里面發(fā)布了各種生信分析的R包。
Bioconductor建立時的兩個宗旨：
reducing barriers 降低科研人員做生信分析的門檻婚肆，提供更方便的工具
reproducibility 實現(xiàn)已經(jīng)發(fā)表文章的重現(xiàn)性
2015年在Nature Methods上發(fā)表了Bioconductor新進展：
完成了質(zhì)量非常好的基因組測序
完善了很多基因組的注釋
我們可以說Bioconductor已經(jīng)是做生信分析最火的工具了，沒有之一
Bioconductor里面的干貨
主要在三個部分：
1.Bioconductor的Learning部分
主要包括一些課程和課程代碼坐慰，還有一些會議的介紹较性，最干的干貨都在這里了，可以說只要有谷歌翻譯和你的耐心，你就可以學會任何這里面有的包赞咙。
特別注意：Bioconductor上Learning部分的course责循，特別特別有用。就是說Bioconductor上有自己的課程攀操，這些課程不但給代碼而且給PPT有的時候還會錄視頻院仿。而且可以說是最官方的資料了因為很多時候都是這些包的作者親自告訴你該怎么用。這里面有單細胞測序的崔赌、甲基化數(shù)據(jù)分析的意蛀、CNV、SNP健芭、RNA-Seq分析的....等等....全都有县钥，非常全，它會告訴你用什么包畫什么圖怎么解決問題慈迈，代碼非常詳細若贮，最最最最重要的是這里面所有的course完全免費！Ｑ髁簟谴麦！
Package vignettes部分是各個Package的“驚鴻一瞥” 它最精華、最特色的部分展示在這里Common work flows部分展示比如：RNA-seq數(shù)據(jù)分析伸头、DNA的SNP分析....我們常用的這些流程匾效，Bioconductor推薦的流程放在這里。這里可以說是Bioconductor干貨中的干貨了恤磷，我們花點時間介紹下
現(xiàn)在它有23個
workflow（去年看的時候才21個面哼，咋又多了兩個？）
我們常見的分析流程它都有扫步。大家遇到不懂的經(jīng)常會各種問魔策，比如：什么什么怎么分析呀？其實我想說的是河胎，不是所有搞生信的都像Jimmy老師那樣是一個全棧式闯袒、全能型人才，大部分人都是有自己所擅長和專攻的版塊游岳，各個版塊之間軟件和R包的用法差別很大政敢，我也不能根據(jù)我的經(jīng)驗去解決你遇到的問題。最好的解決辦法就是這里找胚迫，比如喷户，你想做最近非常火的simpleSingleCell分析晌区，你就點進去摩骨，之后找到HTML點進去，然后你就按照人家的流程一步一步的走就行了朗若，完全不用改代碼恼五，到最后就會跟人家得到一模一樣的結(jié)果，然后再把數(shù)據(jù)換成自己的就可以哭懈。比問任何人都好使~
這里面最重要的是你要理解人家代碼的意圖灾馒，這樣換自己數(shù)據(jù)的時候才知道怎么換。至于理解代碼遣总，那就看自己的R語言的基本功啦睬罗。
2.Bioconductor的工具包說明文檔
一個好用的工具包一定有一個好用的說明文檔。
如果你要學習網(wǎng)上沒有現(xiàn)成中文教程的Bioconductor的工具包旭斥，那閱讀包說明文檔可以幫助你快速學會它的用法容达，這里有一點需要切記：我們使用包只是當它是一個工具讓它幫我們解決實際問題，我們只需要它可以跑起來就行了垂券，不需要精通它花盐，所以千萬不要在閱讀說明文檔的時候太耗神，因為說明文檔一般都很詳細告訴你各種高階的用法菇爪，我們只是用哪個功能就看哪個功能的用法就好算芯，其他的就當沒看見。
DESeq2是目前做差異表達分析最常用最常用的一個R工具包凳宙。我們看看人家文檔寫的特別特別的詳細熙揍，你只需要一個谷歌翻譯，和一個復(fù)制功能氏涩，把人家給的代碼復(fù)制到自己的R里面完全能夠重現(xiàn)人家的結(jié)果届囚，而且人家這個文檔大概兩個月就更新一次，不用擔心代碼過時不能用了削葱。
3.Bioconductor的在線提問內(nèi)容
在Bioconductor的Learning部分的Support site版塊奖亚，點開后注冊一下。比如你有DESeq2的問題析砸，直接輸入search一下就可以了昔字，它會告訴你前人遇到過哪些問題，這個論壇上有很多世界級的大神~
在你使用包時首繁，很多你遇到的問題大家已經(jīng)遇到過了作郭，這時我們可以通過這一板塊解決我們的問題。很多問題都是包的作者在回答弦疮，所以只要你在上面正確的提問基本上都能得到權(quán)威的解答夹攒。
所以以后遇到關(guān)于Bioconductor里包的問題，先在這里找下能不能解決胁塞，再去咨詢別人咏尝。
接下來使用Bioconductor里的包压语，做一些圖，對這些包有一個簡單的認知
GenomicRanges：這個包就是把基因組開始到結(jié)束的結(jié)構(gòu)儲存下來编检，包括的基礎(chǔ)信息就是在染色體在基因組上開始的位置胎食，結(jié)束的位置，在正鏈還是負鏈允懂，最主要的就是這4個信息厕怜，它可以進行基因組的定位。這個包最大的優(yōu)點就是可以做非常復(fù)雜的操作：去找overlap算coverage等等非常方便蕾总。這個包是我們現(xiàn)在做組學分析的基礎(chǔ)粥航，基本上你只要做組學分析用Bioconductor就離不開這個包，很多包都是在這個基礎(chǔ)上構(gòu)建出來的
Biostrings：這個包主要是處理基因組的一些序列信息生百。這個序列信息可以快速的幫我們計算GC含量递雀，最早這個包設(shè)計的時候是方便探針相關(guān)的操作，現(xiàn)在我們芯片數(shù)據(jù)用的少了漸漸的這個功能就淡化了蚀浆。現(xiàn)在主要用于：正負鏈映之、找反向互補鏈、統(tǒng)計GC含量蜡坊、統(tǒng)計各個堿基的含量杠输、三連字母的含量.....這些都是一行命令可以解決的。
BSgenome：說白了它不算是一個包秕衙，它是一個標準的框架：我怎樣把基因組的信息存成一個BSgenome對象? 存成這樣一個對象之后主要包括兩個信息：1. 各條染色體的序列信息和注釋信息 2.整個基因組的SNP（單核苷酸多態(tài)性）信息
目前已經(jīng)有存在的BSgenome對象不到一百個蠢甲，主要都是模式生物和模式作物的。非模式生物和作物的BSgenome怎么獲得呢据忘？有兩種方法：
1.自己構(gòu)建
2.Biomart 這個包里提供的內(nèi)容我們可以直接下載鹦牛，這個包里提供了目前所有已經(jīng)測序的生物的基因組信息
如果想讀取每一個染色體的長度和序列信息同時計算GC含量，就需要結(jié)合剛才的Biostrings包來完成勇吊。
- GenomicFeatures：它也是一個開放性的框架曼追。里面記錄一些基因注釋信息，最后會生成一個對象叫汉规，txdb對象
基因注釋的核心用大白話說就是礼殊，比如，在1號染色體上针史，在1-100bp這個位置上是CDS還是exon晶伦。
exon：外顯子； CDS：編碼一段蛋白產(chǎn)物的序列
一會兒介紹：UCSC網(wǎng)站上的GenomicFeatures啄枕，以及如何通過自己的gff或者gtf文件生成GenomicFeatures文件婚陪。
rtracklayer：這個包簡單一句話，他是讀文件用的频祝。比如我們常見的bed格式泌参，bigwig格式脆淹，bigbed格式，壓縮后的gff格式和gtf格式....怎樣把這些格式讀到R里面并且高效的存儲呢沽一？rtracklayer可以幫助我們未辆，里面的函數(shù)都是已經(jīng)封裝好了，直接讀锯玛，我們就能拿到我們想要的東西。它還會把很多內(nèi)容展成GenomicRanges的對象方便我們的操作兼蜈。我們配合GenomicFeatures包使用可以把我們自己的gff文件讀進來然后進行相關(guān)的操作攘残。
Gviz：幫助我們畫圖的。使用IGV看圖優(yōu)點是不用編程導入到里面就能看为狸，缺點是必須一個一個看歼郭。這個包就可以幫我們批量生成我們想要的圖。
我們知道Bioconductor在創(chuàng)建包的時候最主要是做差異表達分析的辐棒，我們不涉及到這些病曾，因為那個是一些 specific 的包，我們今天從最根本的包講起漾根，在有了今天這些包的基礎(chǔ)上泰涂，再學那些包就會容易很多。因此我們從最根本的包出發(fā)辐怕，然后涉及到基因組的注釋逼蒙、基因的注釋、轉(zhuǎn)錄本的注釋寄疏、exon的注釋是牢、snp的注釋包括怎樣從公共數(shù)據(jù)庫上下載公共數(shù)據(jù)，帶著大家把這些知識點融合到一起陕截，去解決一些具體的生物信息學問題驳棱，這樣的話希望能夠給大家?guī)硪恍┓e極的影響，也希望能夠起到拋磚引玉的作用
如何安裝Bioconductor中的包农曲？
以我們一會兒要用的GenomicRanges包為例社搅，去Bioconductor上搜GenomicRanges點進這個包的詳情頁面，找到安裝代碼復(fù)制我們RStudio里直接運行就好啦~（我去年安裝的時候還不是這行代碼呢乳规，所以問別人不如自己知道在哪里學會罚渐，教程和別人腦子里的存貨都會過時，但是官網(wǎng)資料會隨時更新）
我們先舉一個簡單的例子理解下GRanges函數(shù)
> gr1 = GRanges(seq="chr1",ranges = IRanges(start = 1,end = 10))
> gr1
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 1-10 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
GRanges里有一個對象驯妄，是1號染色體（chr1）從1到10（1-10）星號 * 表示正負鏈都可以荷并。
一個GRanges必須要定義的內(nèi)容有：seq、ranges青扔、strand源织。注：strand如果不定義的話默認是正負鏈都可以翩伪。實際應(yīng)用過程中會復(fù)雜一些，不但要設(shè)置基礎(chǔ)信息而且還會加上復(fù)雜的信息谈息，比如：
gr2 <- GRanges(seqnames = c("chr1","chr2","chr2","chr2","chr1","chr1","chr3","chr3","chr3","chr3"),
ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
strand = c("-","+","+","*","*","+","+","+","-","-"),
score = 1:10,
GC = seq(1, 0, length=10))
gr2
在RStudio運行完成后缘屹，gr2就變?yōu)榱耍?/p>
GRanges object with 10 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
a chr1 101-111 - | 1 1
b chr2 102-112 + | 2 0.888888888888889
c chr2 103-113 + | 3 0.777777777777778
d chr2 104-114 * | 4 0.666666666666667
e chr1 105-115 * | 5 0.555555555555556
f chr1 106-116 + | 6 0.444444444444444
g chr3 107-117 + | 7 0.333333333333333
h chr3 108-118 + | 8 0.222222222222222
i chr3 109-119 - | 9 0.111111111111111
j chr3 110-120 - | 10 0
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
中間會有一個豎線 " | " 前面是基本信息，后面是你添加的信息侠仇。
用gr2$GC（GC content）訪問其中的一列轻姿。
大家應(yīng)該也發(fā)現(xiàn)了，這么創(chuàng)建信息有點冗余逻炊，比如我們有1個Chr1,100個Chr2互亮，1個Chr3且所有的信息都是正鏈， 要是按這么個辦法創(chuàng)建下去余素，額....
所以豹休，就誕生出了一個函數(shù)RLe
比如我們有這樣一個向量a
a = c(1,2,3,3,3,4,4,4,5,5,5,5)
如何快速記錄這樣的信息呢？就用到了RLe函數(shù)
> a = c(1,2,3,3,3,4,4,4,5,5,5,5)
> Rle(a)
numeric-Rle of length 12 with 5 runs
Lengths: 1 1 3 3 4
Values : 1 2 3 4 5
輸出結(jié)果會告訴你1有一個桨吊；2有一個威根；3有三個；4有三個视乐；5有四個
學會了這個洛搀，我們就用RLe創(chuàng)建一個GRanges object
gr3 <- GRanges(seqnames = Rle(c("chr1", "chr2", "chr3", "chr1"), c(1, 2, 3, 2)),
ranges = IRanges(start = 100:107,width = 10),
strand = "+")
### 理解一下 Rle 函數(shù)做了什么 ###
> Rle(c("chr1", "chr2", "chr3", "chr1"), c(1, 2, 3, 2))
character-Rle of length 8 with 4 runs
Lengths: 1 2 3 2
Values : "chr1" "chr2" "chr3" "chr1"
運行結(jié)果：
> gr3
GRanges object with 8 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 100-109 +
[2] chr2 101-110 +
[3] chr2 102-111 +
[4] chr3 103-112 +
[5] chr3 104-113 +
[6] chr3 105-114 +
[7] chr1 106-115 +
[8] chr1 107-116 +
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
這樣是不是方便很多
以下代碼理解下對GRanges object的基本操作，以及as.vector函數(shù)
> seqnames(gr3)
factor-Rle of length 8 with 4 runs
Lengths: 1 2 3 2
Values : chr1 chr2 chr3 chr1
Levels(3): chr1 chr2 chr3
# seqnames()函數(shù)可以告訴我們有1個chr1佑淀；2個chr2姥卢；3個chr3；又接著2個chr1
> as.vector(seqnames(gr3))
[1] "chr1" "chr2" "chr2" "chr3" "chr3" "chr3" "chr1" "chr1"
# 這種形式是不是看著順眼多啦~
通過as.vector函數(shù)渣聚，還可以讓Rle(a)變回去
> a
[1] 1 2 3 3 3 4 4 4 5 5 5 5
> Rle(a)
numeric-Rle of length 12 with 5 runs
Lengths: 1 1 3 3 4
Values : 1 2 3 4 5
> as.vector(Rle(a))
[1] 1 2 3 3 3 4 4 4 5 5 5 5
如果想要獲得start:
> start(gr3)
[1] 100 101 102 103 104 105 106 107
同理
> end(gr3)
[1] 109 110 111 112 113 114 115 116
> strand(gr3)
factor-Rle of length 8 with 1 run
Lengths: 8
Values : +
Levels(3): + - *
為了進一步學習独榴，我們把gr3賦值為更復(fù)雜的信息
gr3 <- GRanges(seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
ranges = IRanges(start = 101:110, end = 111:120, names = head(letters, 10)),
strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
score = 1:10,
GC = seq(1, 0, length=10))
展示運行結(jié)果
> gr3
GRanges object with 10 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
a chr1 101-111 - | 1 1
b chr2 102-112 + | 2 0.888888888888889
c chr2 103-113 + | 3 0.777777777777778
d chr2 104-114 * | 4 0.666666666666667
e chr1 105-115 * | 5 0.555555555555556
f chr1 106-116 + | 6 0.444444444444444
g chr3 107-117 + | 7 0.333333333333333
h chr3 108-118 + | 8 0.222222222222222
i chr3 109-119 - | 9 0.111111111111111
j chr3 110-120 - | 10 0
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
此時，gr3就變得稍微復(fù)雜了一點了奕枝，最前面那一列（a,b,c....j）是names棺榔；第二列（ chr1, chr2....chr3）是seqnames 他們兩個不一樣，需要注意一下隘道。
我們用names(gr3)可以獲取names那一列症歇，也就是第一列：
> names(gr3)
[1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j"
而且names這一列是可以指定的：
names(gr3) <- 1:10
我們來看下指定names之后的gr3：
> gr3
GRanges object with 10 ranges and 2 metadata columns:
seqnames ranges strand | score
<Rle> <IRanges> <Rle> | <integer>
1 chr1 101-111 - | 1
2 chr2 102-112 + | 2
3 chr2 103-113 + | 3
4 chr2 104-114 * | 4
5 chr1 105-115 * | 5
6 chr1 106-116 + | 6
7 chr3 107-117 + | 7
8 chr3 108-118 + | 8
9 chr3 109-119 - | 9
10 chr3 110-120 - | 10
GC
<numeric>
1 1
2 0.888888888888889
3 0.777777777777778
4 0.666666666666667
5 0.555555555555556
6 0.444444444444444
7 0.333333333333333
8 0.222222222222222
9 0.111111111111111
10 0
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
用 gr3加$ 可以取附加列（豎線 " | " 之后的都是附加列）
> gr3$score
[1] 1 2 3 4 5 6 7 8 9 10
> gr3$GC
[1] 1.0000000 0.8888889 0.7777778 0.6666667 0.5555556 0.4444444
[7] 0.3333333 0.2222222 0.1111111 0.0000000
使用函數(shù)mcols()可以把所有的附加列取出來（也就是只要豎線 " | " 之后的部分）
> mcols(gr3)
DataFrame with 10 rows and 2 columns
score GC
<integer> <numeric>
a 1 1
b 2 0.888888888888889
c 3 0.777777777777778
d 4 0.666666666666667
e 5 0.555555555555556
f 6 0.444444444444444
g 7 0.333333333333333
h 8 0.222222222222222
i 9 0.111111111111111
j 10 0
我們已經(jīng)使用函數(shù)mcols()把附加部分取出來了，想取附加部分的某一列也很方便
> mcols(gr3)[,1] #取附加部分的第一列
[1] 1 2 3 4 5 6 7 8 9 10
> mcols(gr3)[,2] #取附加部分的第二列
[1] 1.0000000 0.8888889 0.7777778 0.6666667 0.5555556 0.4444444
[7] 0.3333333 0.2222222 0.1111111 0.0000000
我想知道我的gr3每一個對象到底多寬谭梗？如果是data.frame()我們用它的結(jié)束減開始忘晤，這里width()函數(shù)可以辦到的。
> width(gr3)
[1] 11 11 11 11 11 11 11 11 11 11
這個函數(shù)這么簡單為什么還要單獨拿出來說呢激捏？大家有沒有想到它能發(fā)揮什么樣的作用设塔？
可以統(tǒng)計全基因組所有基因的長度的分布。這里先有個印象如果現(xiàn)在不理解沒關(guān)系
同理远舅，length()函數(shù)可以統(tǒng)計行數(shù)
> length(gr3)
[1] 10
經(jīng)過了上面基礎(chǔ)的鋪墊闰蛔，接下來講點實用的痕钢，也就是GRanges()函數(shù)的特點——GRanges()在處理基因組數(shù)據(jù)的時候，它比data.frame()要快序六！
快在哪里呢任连？
GRanges()函數(shù)可以很方便的過濾數(shù)據(jù)
比如：我們想要獲得score小于5的：gr3[gr3$score < 5]
> gr3[gr3$score < 5]
GRanges object with 4 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
1 chr1 101-111 - | 1 1
2 chr2 102-112 + | 2 0.888888888888889
3 chr2 103-113 + | 3 0.777777777777778
4 chr2 104-114 * | 4 0.666666666666667
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
再比如：我們想要獲得score小于5并且GC count大于0.7：
gr3[gr3$score < 5 & gr3$GC > 0.7]
> gr3[gr3$score < 5 & gr3$GC > 0.7]
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
1 chr1 101-111 - | 1 1
2 chr2 102-112 + | 2 0.888888888888889
3 chr2 103-113 + | 3 0.777777777777778
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
再比如：我們把所有不分正負鏈的篩出來：gr3[strand(gr3) == "*"]
同理穴翩，把所有分正負鏈的取出來：gr3[strand(gr3) != "*"]
> gr3[strand(gr3) == "*"]
GRanges object with 2 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
4 chr2 104-114 * | 4 0.666666666666667
5 chr1 105-115 * | 5 0.555555555555556
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
假如這里面存有全基因組所有的數(shù)據(jù)萎河，我只想分析1號染色體的數(shù)據(jù)，咋辦荸型？：
gr3[seqnames(gr3) == "chr1"]
> gr3[seqnames(gr3) == "chr1"]
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
1 chr1 101-111 - | 1 1
5 chr1 105-115 * | 5 0.555555555555556
6 chr1 106-116 + | 6 0.444444444444444
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
sort(gr3) 就可以按照染色體排序了繁涂，同一條染色體先排正鏈拱她，再排負鏈，最后排不分正負鏈的爆土。
> sort(gr3)
GRanges object with 10 ranges and 2 metadata columns:
seqnames ranges strand | score GC
<Rle> <IRanges> <Rle> | <integer> <numeric>
6 chr1 106-116 + | 6 0.444444444444444
1 chr1 101-111 - | 1 1
5 chr1 105-115 * | 5 0.555555555555556
2 chr2 102-112 + | 2 0.888888888888889
3 chr2 103-113 + | 3 0.777777777777778
4 chr2 104-114 * | 4 0.666666666666667
7 chr3 107-117 + | 7 0.333333333333333
8 chr3 108-118 + | 8 0.222222222222222
9 chr3 109-119 - | 9 0.111111111111111
10 chr3 110-120 - | 10 0
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
按照GC含量排序 gr3[order(gr3$GC)] 默認是從小到大排序；
gr3[order(gr3$GC,decreasing = T)] 這樣就是從大到小排序诸蚕。
所謂多列排序就是先排一列步势，再排一列。
合并兩個GRanges 對象：merge GRange obj
用 'c()' 這個函數(shù)就可以把兩個GRanges 對象合并在一起背犯，效率非常高坏瘩，我們可以看到整個GRanges包效率非常之高。
比如漠魏，我們的需求是倔矾，將gr2和gr3合并在一起，將合并后形成的新對象賦值給gr4柱锹，操作方法：
gr4 = c(gr2,gr3)
兩個對象合并后難免會出現(xiàn)重復(fù)哪自，所以需要去重處理：
unique(gr4)
介紹針對GRanges的四個函數(shù)：flank、shift禁熏、reduce壤巷、disjoin
首先準備一個gr5
gr5 = GRanges(seqnames = "chr2",
ranges = IRanges(start = c(6,8,12,14,21,22,23),width = c(11,4,2,5,7,7,7)),
strand = c("-","-","+","+","+","+","+"))
現(xiàn)在gr5里面的信息應(yīng)該能看懂了吧，這是2號染色體瞧毙，起始位置分別是 6,8,12,14,21,22,23胧华，寬度分別是 11,4,2,5,7,7,7 ，前兩段是負鏈后面都是正鏈宙彪。
> gr5
GRanges object with 7 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 6-16 -
[2] chr2 8-11 -
[3] chr2 12-13 +
[4] chr2 14-18 +
[5] chr2 21-27 +
[6] chr2 22-28 +
[7] chr2 23-29 +
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
flank()函數(shù)功能
上面矩动，每一行都是2號染色體上的基因。
怎樣取基因上游10bp的region（upsteam 10bp）释漆？
想取這段區(qū)域悲没，常規(guī)思維應(yīng)該這樣取：對于正鏈我們在它start位置減10男图，這樣就拿到了上游10bp的區(qū)域檀训；對于負鏈我們在end的位置加10柑潦，這樣就拿到了上游10bp的區(qū)域。因為正負鏈是相對存在的嘛~
但是我們這里使用GRanges()會有一個參數(shù)峻凫，非常方便的取某個range：flank參數(shù)
flank(gr5,width = 10) 一句命令搞定渗鬼，這樣就取到了上游10bp的區(qū)間
問：這步操作有什么用？ 就是人家設(shè)計這個功能是為了什么荧琼？
答：promter 啟動子
啟動子一般被認為是基因上游的2000k以內(nèi)(promter譬胎，gene upstream 2000k)，所以說當我們把一個基因賦值給gr對象以后（gr.obj）我們可以直接通過flank()取出它的啟動子區(qū)域命锄。
思路如下：
gene -> gr.obj -> flank(gr.obj,width = 2000)
是不是這樣很方便堰乔？不止這樣，還有一個更方便的
promoters(gr5,upstream = 2000,downstream = 10) 這樣直接取出來的就是啟動子區(qū)域
啟動子區(qū)域一般被認為是轉(zhuǎn)錄起始位點（TSS）上游2K脐恩，到轉(zhuǎn)錄起始位點下游100bp-200bp的這段區(qū)域镐侯，并不是說TSS以前才是，所以說如果不設(shè)置這樣一個函數(shù)的話驶冒，得分兩步：先取前面苟翻，再取后面。有了promoters函數(shù)骗污，就整合到一起了非常方便崇猫。
下游怎么取需忿？ flank(gr5,width = 10,start = F) 取到了下游10bp的區(qū)間
> flank(gr5,width = 10,start = F)
GRanges object with 7 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 -4-5 -
[2] chr2 -2-7 -
[3] chr2 14-23 +
[4] chr2 19-28 +
[5] chr2 28-37 +
[6] chr2 29-38 +
[7] chr2 30-39 +
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
有負數(shù)怎么辦诅炉？用一步操作把負數(shù)變成1即可：
gr6 = flank(gr5,width = 10,start = F)
start(gr6[start(gr6) < 1]) = 1
> gr6
GRanges object with 7 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 1-5 -
[2] chr2 1-7 -
[3] chr2 14-23 +
[4] chr2 19-28 +
[5] chr2 28-37 +
[6] chr2 29-38 +
[7] chr2 30-39 +
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
shift()函數(shù)功能
如果要對一個區(qū)域進行整體的平移，就用到shift
我們將這里的gr5整體的向正向平移：
shift(gr5,shift = 10)
這樣無論是正鏈還是負鏈屋厘，這個操作都會像正向平移涕烧；
shift(gr5,shift = -10) 向染色體的上游平移
在講reduce()、reduce()之前汗洒，先來張圖理解下shift() VS reduce() VS disjoin()
如圖澈魄，第一行，藍色條就是我們的gr5仲翎；
第二行shift(gr5,5)就是把整體平移了5痹扇；
第三行reduce(gr5)就是把所有重疊（overlap）區(qū)域合并（merge）到一起；
第四行disjoin(gr5)就是把有overlap的地方切割開
問題來了溯香，reduce()能干什么鲫构？disjoin()能干什么？
reduce()函數(shù)功能
計算整個exon的長度玫坛，exon上有重疊的區(qū)域结笨，如果把兩個exon的長度直接相加的話就會算重，這時候就用到reduce()
比如，我們找到一個gtf文件炕吸，把它讀到R里面賦值給GRanges()函數(shù)伐憾，再用reduce()即可算出全部的exon的長度。
思路：gtf -> GRanges -> reduce -> total exon length
disjoin()函數(shù)功能
比如說我有一個gtf文件把它讀到R里面賦值給GRanges()函數(shù)赫模，把它disjoin()一下就可以把有共同轉(zhuǎn)錄區(qū)的都去掉树肃，最后剩下來的都是specific的內(nèi)容，在研究Alternative splicing（可變剪切） 的時候這個操作是非常有用的瀑罗。
思路：gtf -> GRanges -> disjoin -> remove -> overlap region -> specific region
下面我們再介紹一個GRanges里非常好用的一個功能——overlap
overlap應(yīng)用
我們經(jīng)常遇到這樣的問題胸嘴，我有 region1和 region2，我想找region1里面哪段有region2斩祭；region2里哪段有region1 劣像，這就是一個overlap的問題。再具體一點比如我有一堆reads摧玫，是個bam文件耳奕，我想算它map上了多少reads，這也會典型的是一個overlap的問題诬像。overlap的問題都可以轉(zhuǎn)成GRanges的操作屋群。
*我們的chip-seq分析中，經(jīng)常會問到：這個信號在promter是否有富集颅停？
第一步我們得先拿到promter（啟動子）區(qū)域谓晌；
第二步我們得拿到chip-seq的region掠拳；
第三步我們找它們倆有沒有overlap癞揉；
第四步比較它們之間的overlap和基因組的overlap之間有沒有顯著的富集。
所以說這也是一個overlap的問題溺欧。
GRanges如何做overlap喊熟？我們演示下：
我們有兩個GRanges，分別是：gr6 和 gr7
gr6 = GRanges(seqnames = "chr2",
ranges = IRanges(start = c(6,8,12,14,21,22,23),width = c(11,4,2,5,7,7,7)),
strand = "*")
gr7 = GRanges(seqnames = "chr2",
ranges = IRanges(start = c(6,15),width = 10),
strand = "*")
它們分別是：
> gr6
GRanges object with 7 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 6-16 *
[2] chr2 8-11 *
[3] chr2 12-13 *
[4] chr2 14-18 *
[5] chr2 21-27 *
[6] chr2 22-28 *
[7] chr2 23-29 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> gr7
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 6-15 *
[2] chr2 15-24 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
find overlap 操作: findOverlaps(gr6,gr7)
over.gr6 = findOverlaps(gr6,gr7)
over.gr6
overlap的結(jié)果姐刁，這個結(jié)果是什么意思呢芥牌？
gr6的1和gr7的1有overlap；gr6的1和gr7的2有overlap聂使；gr6的2和gr7的1有overlap......
> over.gr6
Hits object with 9 hits and 0 metadata columns:
queryHits subjectHits
<integer> <integer>
[1] 1 1
[2] 1 2
[3] 2 1
[4] 3 1
[5] 4 1
[6] 4 2
[7] 5 2
[8] 6 2
[9] 7 2
-------
queryLength: 7 / subjectLength: 2
那有沒有一個簡單的辦法壁拉，我只要gr6和gr7有overlap的：
gr6[gr6 %over% gr7]
> gr6[gr6 %over% gr7]
GRanges object with 7 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr2 6-16 *
[2] chr2 8-11 *
[3] chr2 12-13 *
[4] chr2 14-18 *
[5] chr2 21-27 *
[6] chr2 22-28 *
[7] chr2 23-29 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
如何統(tǒng)計human全部exon的總長度？
這個要用到GRanges里的reduce()
首先下載含有基因組信息的包柏靶，再加載：
biocLite("EnsDb.Hsapiens.v75")
library(EnsDb.Hsapiens.v75)
ensembl.hg19 = EnsDb.Hsapiens.v75
ensembl.hg19.exon = exons(ensembl.hg19)
# 取這個包里所有的exons的region弃理，取exon會稍微慢一點，因為exon的序列非常的多屎蜓。
看一下exon的長度：
> length(ensembl.hg19.exon)
[1] 745593
來看下exon長度的分布：width(ensembl.hg19.exon)
將長度分布畫成直方圖痘昌，更直觀的顯示：
hist(width(ensembl.hg19.exon))
發(fā)現(xiàn)exon長度有長有短，長的太長了，短的都小到忽略了辆苔，所以我們修改下代碼算灸，讓圖有更好的觀看價值：
hist(log2(width(ensembl.hg19.exon) + 1))
看下exon的region信息：
> ensembl.hg19.exon
GRanges object with 745593 ranges and 1 metadata column:
seqnames ranges strand | exon_id
<Rle> <IRanges> <Rle> | <character>
ENSE00002234944 1 11869-12227 + | ENSE00002234944
ENSE00002234632 1 11872-12227 + | ENSE00002234632
ENSE00002269724 1 11874-12227 + | ENSE00002269724
ENSE00001948541 1 12010-12057 + | ENSE00001948541
ENSE00001671638 1 12179-12227 + | ENSE00001671638
... ... ... ... . ...
ENSE00001741452 Y 28774418-28774584 - | ENSE00001741452
ENSE00001681574 Y 28776794-28776896 - | ENSE00001681574
ENSE00001638296 Y 28779492-28779578 - | ENSE00001638296
ENSE00001797328 Y 28780670-28780799 - | ENSE00001797328
ENSE00001794473 Y 59001391-59001635 + | ENSE00001794473
-------
seqinfo: 273 sequences from GRCh37 genome
算下exon的平均長度：
> mean(width(ensembl.hg19.exon))
[1] 300.8443
這個數(shù)一定要記住：人的基因組exon的平均長度是300
求中位數(shù)：median(width(ensembl.hg19.exon)) 中位數(shù)是141
如果我們直接：sum(width(ensembl.hg19.exon))算全部exon的總長度驻啤，肯定不對菲驴，因為它里面有很多overlap的情況。
所以我們需要用 reduce()
sum(width(reduce(ensembl.hg19.exon)))
> sum(width(ensembl.hg19.exon))
[1] 224307403
> sum(width(reduce(ensembl.hg19.exon)))
[1] 134663597
我們對比上面的數(shù)值街佑，發(fā)現(xiàn)reduce()之后數(shù)量減少了一般谢翎，說明確實去除了overlap。按照嚴格意義來說其實這個數(shù)也不算完全正確沐旨，因為我們沒有區(qū)分正負鏈森逮，不過這么算也有一定的道理，正負鏈分開算嘛磁携。如果正負鏈分開算的話大概有100多兆的區(qū)域是exon的區(qū)域褒侧。
總結(jié)一下，我們目前要記住的數(shù)：
1.exon的平均長度是300
2.中位數(shù)是141
3.如果區(qū)分正負鏈的話exon的總長度大概是134
參考資料：
https://www.bilibili.com/video/av49363776?from=search&seid=10658324414632164518
https://www.bilibili.com/video/av26731585?from=search&seid=10658324414632164518
https://www.bilibili.com/video/av25643438?from=search&seid=10658324414632164518
https://www.bilibili.com/video/av37568990?from=search&seid=10658324414632164518
https://www.bilibili.com/video/av24355734
https://www.bilibili.com/video/av29255401?from=search&seid=10658324414632164518
還有Bioconductor做生信分析入門介紹（下）
