Bioconductor做生信分析入門介紹（上）

這是我聽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新進展：

  • 完成了質量非常好的基因組測序

  • 完善了很多基因組的注釋

我們可以說Bioconductor已經(jīng)是做生信分析最火的工具了，沒有之一

Bioconductor里面的干貨

主要在三個部分：

1.Bioconductor的Learning部分

主要包括一些課程和課程代碼，還有一些會議的介紹，最干的干貨都在這里了，可以說只要有谷歌翻譯和你的耐心，你就可以學會任何這里面有的包。

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個，咋又多了兩個？）
Common work flows部分

我們常見的分析流程它都有。大家遇到不懂的經(jīng)常會各種問，比如：什么什么怎么分析呀？其實我想說的是，不是所有搞生信的都像Jimmy老師那樣是一個全棧式、全能型人才，大部分人都是有自己所擅長和專攻的版塊，各個版塊之間軟件和R包的用法差別很大，我也不能根據(jù)我的經(jīng)驗去解決你遇到的問題。最好的解決辦法就是這里找，比如，你想做最近非常火的simpleSingleCell分析，你就點進去，之后找到HTML點進去，然后你就按照人家的流程一步一步的走就行了，完全不用改代碼，到最后就會跟人家得到一模一樣的結果，然后再把數(shù)據(jù)換成自己的就可以。比問任何人都好使~
這里面最重要的是你要理解人家代碼的意圖，這樣換自己數(shù)據(jù)的時候才知道怎么換。至于理解代碼，那就看自己的R語言的基本功啦。

2.Bioconductor的工具包說明文檔

一個好用的工具包一定有一個好用的說明文檔。
如果你要學習網(wǎng)上沒有現(xiàn)成中文教程的Bioconductor的工具包，那閱讀包說明文檔可以幫助你快速學會它的用法，這里有一點需要切記：我們使用包只是當它是一個工具讓它幫我們解決實際問題，我們只需要它可以跑起來就行了，不需要精通它，所以千萬不要在閱讀說明文檔的時候太耗神，因為說明文檔一般都很詳細告訴你各種高階的用法，我們只是用哪個功能就看哪個功能的用法就好，其他的就當沒看見。

image.png

與上張圖在同一個網(wǎng)頁上

DESeq2是目前做差異表達分析最常用最常用的一個R工具包。我們看看人家文檔寫的特別特別的詳細，你只需要一個谷歌翻譯，和一個復制功能，把人家給的代碼復制到自己的R里面完全能夠重現(xiàn)人家的結果，而且人家這個文檔大概兩個月就更新一次，不用擔心代碼過時不能用了。
image.png

3.Bioconductor的在線提問內容

在Bioconductor的Learning部分的Support site版塊，點開后注冊一下。比如你有DESeq2的問題，直接輸入search一下就可以了，它會告訴你前人遇到過哪些問題，這個論壇上有很多世界級的大神~
在你使用包時，很多你遇到的問題大家已經(jīng)遇到過了，這時我們可以通過這一板塊解決我們的問題。很多問題都是包的作者在回答，所以只要你在上面正確的提問基本上都能得到權威的解答。
所以以后遇到關于Bioconductor里包的問題，先在這里找下能不能解決，再去咨詢別人。

接下來使用Bioconductor里的包，做一些圖，對這些包有一個簡單的認知

  • GenomicRanges：這個包就是把基因組開始到結束的結構儲存下來，包括的基礎信息就是在染色體在基因組上開始的位置，結束的位置，在正鏈還是負鏈，最主要的就是這4個信息，它可以進行基因組的定位。這個包最大的優(yōu)點就是可以做非常復雜的操作：去找overlap算coverage等等非常方便。這個包是我們現(xiàn)在做組學分析的基礎，基本上你只要做組學分析用Bioconductor就離不開這個包，很多包都是在這個基礎上構建出來的

  • Biostrings：這個包主要是處理基因組的一些序列信息。這個序列信息可以快速的幫我們計算GC含量，最早這個包設計的時候是方便探針相關的操作，現(xiàn)在我們芯片數(shù)據(jù)用的少了漸漸的這個功能就淡化了。現(xiàn)在主要用于：正負鏈、找反向互補鏈、統(tǒng)計GC含量、統(tǒng)計各個堿基的含量、三連字母的含量.....這些都是一行命令可以解決的。

  • BSgenome：說白了它不算是一個包，它是一個標準的框架：我怎樣把基因組的信息存成一個BSgenome對象? 存成這樣一個對象之后主要包括兩個信息：1. 各條染色體的序列信息和注釋信息 2.整個基因組的SNP（單核苷酸多態(tài)性）信息

目前已經(jīng)有存在的BSgenome對象不到一百個，主要都是模式生物和模式作物的。非模式生物和作物的BSgenome怎么獲得呢？有兩種方法：

1.自己構建

2.Biomart 這個包里提供的內容我們可以直接下載，這個包里提供了目前所有已經(jīng)測序的生物的基因組信息

如果想讀取每一個染色體的長度和序列信息同時計算GC含量，就需要結合剛才的Biostrings包來完成。

  • GenomicFeatures：它也是一個開放性的框架。里面記錄一些基因注釋信息，最后會生成一個對象叫，txdb對象

基因注釋的核心用大白話說就是，比如，在1號染色體上，在1-100bp這個位置上是CDS還是exon。

exon：外顯子； CDS：編碼一段蛋白產物的序列

一會兒介紹：UCSC網(wǎng)站上的GenomicFeatures，以及如何通過自己的gff或者gtf文件生成GenomicFeatures文件。

  • rtracklayer：這個包簡單一句話，他是讀文件用的。比如我們常見的bed格式，bigwig格式，bigbed格式，壓縮后的gff格式和gtf格式....怎樣把這些格式讀到R里面并且高效的存儲呢？rtracklayer可以幫助我們，里面的函數(shù)都是已經(jīng)封裝好了，直接讀，我們就能拿到我們想要的東西。它還會把很多內容展成GenomicRanges的對象方便我們的操作。我們配合GenomicFeatures包使用可以把我們自己的gff文件讀進來然后進行相關的操作。

  • Gviz：幫助我們畫圖的。使用IGV看圖優(yōu)點是不用編程導入到里面就能看，缺點是必須一個一個看。這個包就可以幫我們批量生成我們想要的圖。

我們知道Bioconductor在創(chuàng)建包的時候最主要是做差異表達分析的，我們不涉及到這些，因為那個是一些 specific 的包，我們今天從最根本的包講起，在有了今天這些包的基礎上，再學那些包就會容易很多。因此我們從最根本的包出發(fā)，然后涉及到基因組的注釋、基因的注釋、轉錄本的注釋、exon的注釋、snp的注釋包括怎樣從公共數(shù)據(jù)庫上下載公共數(shù)據(jù)，帶著大家把這些知識點融合到一起，去解決一些具體的生物信息學問題，這樣的話希望能夠給大家?guī)硪恍┓e極的影響，也希望能夠起到拋磚引玉的作用

如何安裝Bioconductor中的包？

以我們一會兒要用的GenomicRanges包為例，去Bioconductor上搜GenomicRanges點進這個包的詳情頁面，找到安裝代碼復制我們RStudio里直接運行就好啦~（我去年安裝的時候還不是這行代碼呢，所以問別人不如自己知道在哪里學會，教程和別人腦子里的存貨都會過時，但是官網(wǎng)資料會隨時更新）

image.png

我們先舉一個簡單的例子理解下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必須要定義的內容有：seq、ranges、strand。注：strand如果不定義的話默認是正負鏈都可以。實際應用過程中會復雜一些，不但要設置基礎信息而且還會加上復雜的信息，比如：

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）訪問其中的一列。
大家應該也發(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

輸出結果會告訴你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"

運行結果：

> 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賦值為更復雜的信息

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))

展示運行結果

> 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就變得稍微復雜了一點了，最前面那一列（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()我們用它的結束減開始，這里width()函數(shù)可以辦到的。

> width(gr3)
 [1] 11 11 11 11 11 11 11 11 11 11

這個函數(shù)這么簡單為什么還要單獨拿出來說呢？大家有沒有想到它能發(fā)揮什么樣的作用？
可以統(tǒng)計全基因組所有基因的長度的分布。這里先有個印象如果現(xiàn)在不理解沒關系
同理，length()函數(shù)可以統(tǒng)計行數(shù)

> length(gr3)
[1] 10

經(jīng)過了上面基礎的鋪墊，接下來講點實用的，也就是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)重復，所以需要去重處理：
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里面的信息應該能看懂了吧，這是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ī)思維應該這樣取：對于正鏈我們在它start位置減10，這樣就拿到了上游10bp的區(qū)域；對于負鏈我們在end的位置加10，這樣就拿到了上游10bp的區(qū)域。因為正負鏈是相對存在的嘛~
但是我們這里使用GRanges()會有一個參數(shù)，非常方便的取某個range：flank參數(shù)
flank(gr5,width = 10) 一句命令搞定，這樣就取到了上游10bp的區(qū)間
問：這步操作有什么用？ 就是人家設計這個功能是為了什么？
答：promter 啟動子
啟動子一般被認為是基因上游的2000k以內(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ū)域一般被認為是轉錄起始位點（TSS）上游2K，到轉錄起始位點下游100bp-200bp的這段區(qū)域，并不是說TSS以前才是，所以說如果不設置這樣一個函數(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()

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()一下就可以把有共同轉錄區(qū)的都去掉，最后剩下來的都是specific的內容，在研究Alternative splicing（可變剪切） 的時候這個操作是非常有用的。

思路：gtf -> GRanges -> disjoin -> remove -> overlap region -> specific region

下面我們再介紹一個GRanges里非常好用的一個功能——overlap

overlap應用

我們經(jīng)常遇到這樣的問題，我有 region1和 region2，我想找region1里面哪段有region2；region2里哪段有region1 ，這就是一個overlap的問題。再具體一點比如我有一堆reads，是個bam文件，我想算它map上了多少reads，這也會典型的是一個overlap的問題。overlap的問題都可以轉成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的結果，這個結果是什么意思呢？
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))
發(fā)現(xiàn)絕大多數(shù)的exon長度大概在250多bp

看下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ù)一定要記?。喝说幕蚪Mexon的平均長度是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ū)域。

總結一下，我們目前要記住的數(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做生信分析入門介紹（下）

最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容