小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看

搜索
分享

Bioinformatic Data Skills 學(xué)習(xí)專題(6) Genomic Range之二

 微笑如酒 2018-06-24

概念回溯

  • 為什么要用ranges?

  • 圖像畫思維：雖然在考慮結(jié)構(gòu)變異的時(shí)候我們已經(jīng)將線性的基因組縮短到了插入、刪除、轉(zhuǎn)位和拷貝數(shù)目變異等類別，但是更好的方法就是把一維序列當(dāng)作多位圖像來(lái)看

6. Run Length Encoding and Views

除了可以表示特定位置的核酸序列之外，一個(gè)區(qū)間也可以有其他含義，每個(gè)堿基位置上可以有不同的數(shù)值含義，比如：

  • Coverage覆蓋度，即堿基上reads的總數(shù)

  • Conservation tracks，用于表示某堿基在不同物種間的的進(jìn)化保守性（可以用phastCons計(jì)算）

  • 以人群為單位，計(jì)算某特定位置堿基的多樣性（哈哈，沒(méi)錯(cuò)就是SNP） 當(dāng)然還有很多其他的信息，但是作者在這張主要科普coverage的計(jì)算，從序列數(shù)據(jù)中獲得區(qū)間信息，以及一個(gè)工具view。

6.1 Run-length encoding and coverage()

由于coverage數(shù)據(jù)通常都比較連續(xù)延綿（run），IRanges就用了一個(gè)名為 run-length encoding的算法把值一樣的位點(diǎn)重疊起來(lái)算作一個(gè)'run'（有沒(méi)有想到一種數(shù)據(jù)格式？答案見(jiàn)文末，也可以參考我們生信菜鳥團(tuán)公眾號(hào)的“數(shù)據(jù)格式專題”欄目）。這樣就可以把分辨率從單堿基縮小到若干堿基乃至一個(gè)很大的區(qū)域，這就大大節(jié)省了存儲(chǔ)空間。打個(gè)比方就是把一條線很有規(guī)律地折疊起來(lái)，數(shù)值一樣的點(diǎn)就算一個(gè)run(每個(gè)run長(zhǎng)度可能不一樣，但堿基上的覆蓋度必然一樣)，像一個(gè)個(gè)山峰：

6.1.1

首先用 Rle將包含單堿基覆蓋度信息序列轉(zhuǎn)換為一個(gè)run的信息

  1. > x <->as.integer(c(4, 4, 4, 3, 3, 2, 1, 1, 1, 1, 1, 0, 0, 0,

  2. 0, 0, 0, 0, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4))

  3. > xrle <->Rle(x)

  4. > xrle

  5. integer-Rle of length 28 with 7 runs

  6. Lengths: 3 2 1 5 7 3 7

  7. Values : 4 3 2 1 0 1 4

由于xrle本質(zhì)上也是IRange的子類，所以同樣地，我們也可以對(duì)壓縮好的數(shù)據(jù) xrle進(jìn)行IRanges的常規(guī)操作：

  1. > as.vector(xrle)

  2. [1] 4 4 4 3 3 2 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 4 4 4 4 4 4 4

  3. > xrle/2  ## 注意runlength不變，值減少為原來(lái)的一半

  4. numeric-Rle of length 28 with 7 runs

  5. Lengths: 3 2 1 5 7 3 7

  6. Values : 2 1.5 1 0.5 0 0.5 2

  7. > xrle[xrle > 3]

  8. > runLength(xrle) ## 獲得每個(gè)run的長(zhǎng)度  

  9. [1] 3 2 1 5 7 3 7

  10. > runValue(xrle) ## 獲得每個(gè)run里堿基的覆蓋度

  11. [1] 4 3 2 1 0 1 4

6.1.2

rle能直接壓縮一條序列上的信息，但是如何得到一段序列上若干reads的覆蓋度信息呢？

另一種更接地氣的問(wèn)法就是“如何從原始測(cè)序數(shù)據(jù)獲得一段特定序列上run-length格式的coverage數(shù)據(jù)？”。答案如下，關(guān)鍵是 coverage這個(gè)函數(shù)。同時(shí)IRanges也提供了統(tǒng)計(jì)一段區(qū)域上run信息的方法，參考第三段代碼。

  1. # 1. 造reads，獲得含有coverage信息的序列

  2. > set.seed(0)

  3. > rngs <->IRanges(start=sample(seq_len(60), 10), width=7) ## 模擬70條長(zhǎng)度為7的序列

  4. > names(rngs)[9] <->'A' # label one range for examples later

  5. > rngs_cov <- coverage(rngs)="">## coverage提取reads覆蓋度信息

  6. > rngs_cov

  7. integer-Rle of length 63 with 18 runs

  8. Lengths: 11 4 3 3 1 6 4 2 5 2 7 2 3 3 1 3 2 1

  9. Values : 0 1 2 1 2 1 0 1 2 1 0 1 2 3 4 3 2 1

  10. # 這是一個(gè)長(zhǎng)度為63bp的序列，由18個(gè)run組成，參考上圖，可以發(fā)現(xiàn)有18條橫線

  12. # 2. 同樣可以取子集

  13. > rngs_cov > 3 # where is coverage greater than 3?

  14. logical-Rle of length 63 with 3 runs

  15. Lengths: 56 1 6

  16. Values : FALSE TRUE FALSE

  17. > rngs_cov[as.vector(rngs_cov) > 3] # extract the depths that are greater than 3

  18. integer-Rle of length 1 with 1 run

  19. Lengths: 1

  20. Values : 4

  22. # 3. 取出某個(gè)區(qū)域的run信息

  24. > rngs_cov[rngs['A']]  #step1中命名為A一個(gè)reads

  25. integer-Rle of length 7 with 2 runs

  26. Lengths: 5 2

  27. Values : 2 1

  29. # 4. 看看這個(gè)reads周圍的特征，比如豐度之類

  30. > mean(rngs_cov[rngs['A']])  

  31. [1] 1.714286

小結(jié)-Coverage節(jié)

很多生物學(xué)的問(wèn)題追根溯源可以歸結(jié)為遺傳和表觀遺傳問(wèn)題，但他們都離不開二級(jí)結(jié)構(gòu)的特征分析，我們可以問(wèn)一些問(wèn)題，比如發(fā)生遺傳突變的位置是不是少了某些蛋白質(zhì)的結(jié)合？是不是突變的序列導(dǎo)致的？那么這些序列有什么特征呢？GC content/nucleotide diversity等等都是值得追問(wèn)的；除了探索為止序列（測(cè)序結(jié)果），我們還可以看實(shí)驗(yàn)變量帶來(lái)的已知序列上的特征改變，比如重復(fù)元件、蛋白編碼區(qū)、低重組區(qū)域上的甲基化水平、組蛋白修飾水平或者染色質(zhì)開放程度的改變。關(guān)鍵還是要探究變量可能帶來(lái)的影響。 一旦我們把改變量化到區(qū)間或者序列上，分辨率降低，計(jì)算因此也變得相對(duì)高效。GenomicRanges也提供相類似的功能。

6.2 Going from run-length encoded sequences to ranges with slice()

Rle可以把區(qū)域信息pileup起來(lái)，pileup的數(shù)據(jù)需要做一定的過(guò)濾，比如堿基的coverge至少要為某個(gè)值，這就需要用 slice函數(shù)去處理獲得的run信息。 slice很形象，可以想象在峰圖上畫一條水平線，保留在水平線以上的區(qū)域，因此獲得的信息多是連續(xù)的run。被slice過(guò)的rle稱為 view，屬于一個(gè)序列的一部分，特征是存在一個(gè)個(gè)峰。

  1. > min_cov2 <- slice(rngs_cov,="" lower="">2)

  2. > min_cov2

  3. Views on a 63-length Rle subject

  4. views:

  5. start end width

  6. [1] 16 18 3 [2 2 2]

  7. [2] 22 22 1 [2]

  8. [3] 35 39 5 [2 2 2 2 2]

  9. [4] 51 62 12 [2 2 2 3 3 3 4 3 3 3 2 2]

6.3 Advanced IRanges: Views

我們使用 slice獲得 rngs_cov中run的子集或者并集，得到了一個(gè)覆蓋度大于一定值的 view。某種意義上， view是的 rngs_cov子類，這個(gè)子類繼承了父類的方法，但是在 view中，可用的方法前面都要加一個(gè)view，比如 viewMaxs, viewMeans, viewApply。

  1. > viewMeans(min_cov2) # 統(tǒng)計(jì)coverage>2連續(xù)區(qū)間上覆蓋reads的平均值

  2. [1] 2.000000 2.000000 2.000000 2.666667

  3. > viewMaxs(min_cov2) # 統(tǒng)計(jì)coverage>2連續(xù)區(qū)間上覆蓋reads數(shù)的最大值

  4. [1] 2 2 2 4

  5. > viewApply(min_cov2, median)

  6. [1] 2 2 2 3

除了根據(jù)覆蓋度可以篩選特定區(qū)域，我們也可以把一個(gè)區(qū)域分成若干個(gè)window/bin，統(tǒng)計(jì)每個(gè)window信息

  1. > length(rngs_cov)

  2. [1] 63

  3. > bwidth <->5L # 設(shè)定bin寬度為5bp

  4. > end <- bwidth="" *="" floor(length(rngs_cov)="" bwidth)=""># 獲得序列最后一個(gè)位置信息

  5. > windows <->IRanges(start=seq(1, end, bwidth), width=bwidth) # 根據(jù)start, end和binWidth制定windows

  6. > head(windows)

  7. IRanges of length 6

  8. start end width

  9. [1] 1 5 5

  10. [2] 6 10 5

  11. [3] 11 15 5

  12. [4] 16 20 5

  13. [5] 21 25 5

  14. [6] 26 30 5

  15. > cov_by_wnd <->Views(rngs_cov, windows) # 用Views獲得rngs_cov上我們指定的windows上的信息

  16. > head(cov_by_wnd)

  17. Views on a 63-length Rle subject

  18. views:

  19. start end width

  20. [1] 1 5 5 [0 0 0 0 0]

  21. [2] 6 10 5 [0 0 0 0 0]

  22. [3] 11 15 5 [0 1 1 1 1]

  23. [4] 16 20 5 [2 2 2 1 1]

  24. [5] 21 25 5 [1 2 1 1 1]

  25. [6] 26 30 5 [1 1 1 0 0]

  26. > viewMeans(cov_by_wnd) #就可以得到均等區(qū)間上的平均信號(hào)啦！

  27. [1] 0.0 0.0 0.8 1.6 1.2 0.6 0.8 1.8 0.2 0.4 2.4 3.2

7. Storing Genomic Ranges with GenomicRanges

GenomicRanges在 IRanges基礎(chǔ)上加了一類 GRanges對(duì)象用于儲(chǔ)存序列信息，它在 IRanges基礎(chǔ)上加入了其他兩類信息用來(lái)指明基因組信息：序列名稱比如染色體號(hào)以及正負(fù)鏈信息。 小編再次提示，假如不明白什么是genomic range，請(qǐng)閱讀我們生信菜鳥團(tuán)公眾號(hào)里“生物數(shù)據(jù)格式專題”。

在指明序列的基礎(chǔ)上，我們需要用 metadata columns指明GRanges上的除了序列信息之外的額外信息，比如gc含量，覆蓋度。所有metadata數(shù)據(jù)都存儲(chǔ)在 DataFrame里，區(qū)別于R里的data.frame，支持更多了列類型，比如：

  • run-length encoded vectors

  • identifiers and names: 轉(zhuǎn)錄本、snp信息、外顯子名稱等等

  • annotation data： GC 含量、重復(fù)序列信息、保守性打分等等

  • 假如是比對(duì)的reads信息，還可以存儲(chǔ)Q30和gap數(shù)啊，真是妙用!

'the union of genomic location with any type of data is what makes GRanges so powerful'

  1. > library(GenomicRanges)

  2. > gr <->GRanges(seqname=c('chr1', 'chr1', 'chr2', 'chr3'),

  3. + ranges=IRanges(start=5:8, width=10),

  4. + strand=c('+', '-', '-', '+'))

  5. > gr

  6. GRanges object with 4 ranges and 0 metadata columns:

  7.      seqnames    ranges strand

  8.         Rle> IRanges>  Rle>

  9.  [1]     chr1   [5, 14]      +

  10.  [2]     chr1   [6, 15]      -

  11.  [3]     chr2   [7, 16]      -

  12.  [4]     chr3   [8, 17]      +

  13.  -------

  14.  seqinfo: 3 sequences from an unspecified genome; no seqlengths

  1. #2. 再加點(diǎn)料，比如gc數(shù)目, 注意gc是metadata哦，所以有`|`分割

  2. > gr <->GRanges(seqname=c('chr1', 'chr1', 'chr2', 'chr3'),

  3.  ranges=IRanges(start=5:8, width=10),

  4.  strand=c('+', '-', '-', '+'), gc=round(runif(4), 3))

  5. > gr

  6. GRanges object with 4 ranges and 1 metadata column:

  7.      seqnames    ranges strand |        gc

  8.         Rle> IRanges>  Rle> |

  9.  [1]     chr1   [5, 14]      + |     0.238

  10.  [2]     chr1   [6, 15]      - |      0.45

  11.  [3]     chr2   [7, 16]      - |     0.236

  12.  [4]     chr3   [8, 17]      + |      0.49

  13.  -------

  14.  seqinfo: 3 sequences from an unspecified genome; no seqlengths

  16. #3. 如何解決“seqinfo: 3 sequences from an unspecified genome; no seqlengths”？

  18. ## strategy1

  19. seqlens <- c(chr1="">152, chr2=432, chr3=903)

  20. gr <->GRanges(seqname=c('chr1', 'chr1', 'chr2', 'chr3'),

  21. ranges=IRanges(start=5:8, width=10),

  22. strand=c('+', '-', '-', '+'),

  23. gc=round(runif(4), 3),

  24. seqlengths=seqlens)

  26. ## strategy2

  27. > seqlengths(gr) <- seqlens=""># another way to do the same as above

  29. #4.基本操作

  30. start(gr)

  31. end(gr)

  32. width(gr)

  33. seqnames(gr)

  34. strand(gr)

  35. ranges(gr)

  36. length(gr) #4

  37. names(gr) <->1:length(gr)]

  38. start(gr) > 7

  39. table(seqnames(gr)) # 有多少條seq，比如染色體

  41. #5. metadata相關(guān)操作,可以計(jì)算子集的各種metadata的統(tǒng)計(jì)信息

  42. mcols(gr) #metadata的列信息

  43. mcols(gr)$gc

  44. > mcols(gr)$gc #等價(jià)

  45. [1] 0.897 0.266 0.372 0.573

  46. > gr$gc

  47. [1] 0.897 0.266 0.372 0.573

  48. > mcols(gr[seqnames(gr) == 'chr1'])$gc

  49. [1] 0.897 0.266

  50. > mean(mcols(gr[seqnames(gr) == 'chr1'])$gc)

  51. [1] 0.5815

8. Grouping Data with GRangesList

我們可以使用 GRangesList或 c()把 GRanges對(duì)象整合到列表中。 seqnames(), start(), end(), width(), ranges(), strand()也都可以用。

  1. > gr1 <->GRanges(c('chr1', 'chr2'), IRanges(start=c(32, 95), width=c(24, 123)))

  2. > gr2 <->GRanges(c('chr8', 'chr2'), IRanges(start=c(27, 12), width=c(42, 34)))

  3. > grl <->GRangesList(gr1, gr2)

  4. > doubled_grl <- c(grl,="">

  5. > length(doubled_grl)

  7. # 操作grl： 根據(jù)染色體信息split grl

  8. > chrs <->'chr3', 'chr1', 'chr2', 'chr2', 'chr3', 'chr1')

  9. > gr <->GRanges(chrs, IRanges(sample(1:100, 6, replace=TRUE),

  10. width=sample(3:30, 6, replace=TRUE)))

  11. > head(gr)

  12. GRanges with 6 ranges and 0 metadata columns:

  13. seqnames ranges strand

  14. Rle> IRanges> Rle>

  15. [1] chr3 [90, 93] *

  16. [2] chr1 [27, 34] *

  17. [3] chr2 [38, 44] *

  18. [4] chr2 [58, 79] *

  19. [5] chr3 [91, 103] *

  20. [6] chr1 [21, 44] *

  21. ---

  22. seqlengths:

  23. chr3 chr1 chr2

  24. NA NA NA

  25. > gr_split <- split(gr,="">

  26. > gr_split[[1]]

  27. GRanges with 4 ranges and 0 metadata columns:

  28. seqnames ranges strand

  29. Rle> IRanges> Rle>

  30. [1] chr3 [90, 93] *

  31. [2] chr3 [91, 103] *

  32. [3] chr3 [90, 105] *

  33. [4] chr3 [95, 117] *

  34. ---

  35. seqlengths:

  36. chr3 chr1 chr2

  37. NA NA NA

  38. > names(gr_split)

  39. [1] 'chr3' 'chr1' 'chr2

  42. # lapply操作 grl

  43. > lapply(gr_split, function(x) order(width(x)))

  44. $chr3

  45. [1] 1 2 3 4

  46. $chr1

  47. [1] 1 2

  48. $chr2

  49. [1] 1 4 2 3

  50. > sapply(gr_split, function(x) min(start(x)))

  51. chr3 chr1 chr2

  52. 90 21 38

  53. > sapply(gr_split, length)

  54. chr3 chr1 chr2

  55. 4 2 4

  56. > elementLengths(gr_split)

  57. chr3 chr1 chr2

  58. 4 2 4

9. Working with Annotation Data: GenomicFeatures and rtracklayer

這節(jié)主要講了兩個(gè)包，他們的核心就是GenomicRanges的各種擴(kuò)展。

  • 一個(gè)是 GenomicFeatures，其工作基礎(chǔ)是菜鳥團(tuán)前期介紹過(guò)的一個(gè)注釋數(shù)據(jù)庫(kù) TranscriptDb，有各種transcript的信息。

  • 另一個(gè)包 rtracklayer能夠?qū)牒蛯?dǎo)出許多常見(jiàn)的文件類型，BED,WIG,GTF,還能查詢和導(dǎo)航UCSC genome Browser rtracklayer包含測(cè)序所用BED文件。

emmm，其實(shí)，這兩個(gè)包我們之前都已經(jīng)有專門提到過(guò)啦，詳見(jiàn)前期推送的兩個(gè)相關(guān)專題： R專題和Bioconductor注釋專題，小菜鳥們務(wù)必參考一下

Practices

6～9 小節(jié)主要講了兩個(gè)range對(duì)象和注釋相關(guān)的包，我們對(duì)Range已經(jīng)相對(duì)了解，知道它可以以Rle這種壓縮數(shù)據(jù)的方法存儲(chǔ)，接下來(lái)就是一些簡(jiǎn)單的使用。Buffalo在github上共享了一個(gè)gtf文件（Musmusculus.GRCm38.75chr1.gtf.gz），書里就以這個(gè)數(shù)據(jù)為例講了一些簡(jiǎn)單的實(shí)戰(zhàn)。

無(wú)論如何，先讀取數(shù)據(jù)吧

  1. > library(rtracklayer);library(GenomicRanges)

  2. > mm_gtf <->import('Mus_musculus.GRCm38.75_chr1.gtf.gz') #rtracklayer功能請(qǐng)看之前的推送～

  3. > colnames(mcols(mm_gtf)) # metadata columns read in

  4. [1] 'source' 'type' 'score'

  5. [4] 'phase' 'gene_id' 'gene_name'

  6. [7] 'gene_source' 'gene_biotype' 'transcript_id'

  7. [10] 'transcript_name' 'transcript_source' 'tag'

  8. [13] 'exon_number' 'exon_id' 'ccds_id'

  9. [16] 'protein_id'

隨機(jī)取一些基因做測(cè)試對(duì)象，以及輸出方法示例

  1. > set.seed(0)

  2. > pseudogene_i <- which(mm_gtf$gene_biotype="=">'pseudogene' & mm_gtf$type == 'gene')

  3. > pseudogene_sample <- sample(pseudogene_i,="">5)

  4. > export(mm_gtf[pseudogene_sample], con='five_random_pseudogene.gtf', format='GTF') #輸出GTF

  5. > bed_data <->

  6. > mcols(bed_data) <- null=""># clear out metadata columns

  7. > export(bed_data, con='five_random_pseudogene.bed', format='BED') #輸出BED




下期終于要引來(lái)Ranges的終章：練習(xí)篇！


自己看了一些，感覺(jué)有些還是需要多拿數(shù)據(jù)練習(xí)，邊學(xué)邊用才比較能理解其含義，大家不妨自己下些數(shù)據(jù)來(lái)練練手。

    本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn)。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購(gòu)買等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊一鍵舉報(bào)。
    轉(zhuǎn)藏 分享 獻(xiàn)花（0）

    0條評(píng)論

    發(fā)表

    請(qǐng)遵守用戶 評(píng)論公約

    類似文章 更多