Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Summary for Computational Genomics with R chapter 6

Funakoshi
December 05, 2020

Summary for Computational Genomics with R chapter 6

Funakoshi

December 05, 2020
Tweet

Other Decks in Technology

Transcript

  1. GRangesオブジェクト_ゲノム範囲データの入れ物 作り方 library(GenomicRanges) gr=GRanges(seqnames=c(“chr1”,“chr2”,“chr2”), ① ranges=IRanges(start=c(50,150,200), ➁ end=c(100,200,300)), strand=c(“+”,“-”,“-”) ➂

    ) 出力結果 ## GRanges object with 3 ranges and 0 metadata columns: ## seqnames ranges strand ## <Rle> <IRanges> <Rle> ## [1] chr1 50-100 + ## [2] chr2 150-200 - ## [3] chr2 200-300 - ## ------- ## seqinfo: 2 sequences from an unspecified genome; no seqlengths GRangesに入力する情報 ①ゲノム範囲の名前 ➁位置 ➂2本鎖のどちらの配列か ※メタデータも追加可能 https://bio-megane.com/
  2. タスク1:異なる二つのデータを重ね合わせてみる Chip-seqとは? 着目するタンパク質の結合配列を調べる →転写因子SP1に着目 (ref: https://www.rhelixa.com/seq/chip-seq/) 重ね合わせるデータ① 重ね合わせるデータ➁ CpG islandsとは?

    5’-CG-3’と続いた配列。Cのメチル化状態と遺伝子 発現の関係性が報告されている (ref: https://en.wikipedia.org/wiki/CpG_site) 転写因子が結合するDNAはプロモーター配列 のみとは限らない。ノイズが大きい可能性。 確度の高いデータへ!
  3. CpG配列データの準備 コード # データの読み込み filePath=system.file("extdata", "cpgi.hg19.chr21.bed", package="compGenomRData") cpgi.df = read.table(filePath,

    header = FALSE, stringsAsFactors=FALSE) # データの前処理(余分なデータを含めない) cpgi.df =cpgi.df [grep("_",cpgi.df[,1],invert=TRUE),] # データフレームから、GRangeオブジェクトへ変換 cpgi.gr=GRanges(seqnames=cpgi.df[,1], ranges=IRanges(start=cpgi.df[,2], end=cpgi.df[,3])) 出力結果 > cpgi.gr GRanges object with 205 ranges and 0 metadata columns: seqnames ranges strand <Rle> <IRanges> <Rle> [1] chr21 9825442-9826296 * [2] chr21 9909011-9909218 * [3] chr21 9968264-9968620 * [4] chr21 10989913-10991413 * [5] chr21 14409412-14410501 * ... ... ... ... [201] chr21 47918497-47918728 * [202] chr21 48018542-48018791 * [203] chr21 48055199-48056060 * [204] chr21 48068517-48068808 * [205] chr21 48081241-48081849 * -------
  4. Chip-seqデータの準備 コード # データの読み込み library(genomation) filePathPeaks=system.file("extdata", "wgEncodeHaibTfbsGm12878Sp1Pcr1xPkRep1.broadPeak.gz", package="compGenomRData") # Grangesオブジェクトへの変換(データフレームを経由しなくて良い)

    pk1.gr=readBroadPeak(filePathPeaks) 補足 データの読み込みは、webサイトから直接読み込むことも可能 (例)rtracklayerパッケージ CpG配列データが205個だったのに比べるとデータ数がかなり多い理由は、 Chip-seqの配列データは、全染色体由来だから。CpGと同じChr21に限れば620個。 使用データの作成方法 > pk1.gr GRanges object with 26121 ranges and 5 metadata columns: seqnames ranges strand | name score signalValue <Rle> <IRanges> <Rle> | <character> <integer> <numeric> [1] chr1 9991-10480 * | peak1 143 464.20 [2] chr1 11021-12230 * | peak2 347 1123.56 [3] chr1 13883-14300 * | peak3 68 221.82 [4] chr1 22721-23011 * | peak4 33 109.62 [5] chr1 23480-23786 * | peak5 40 129.51 ... ... ... ... . ... ... ... [26117] chrX 154611046-154611315 * | peak26117 55 179.08 [26118] chrX 154688675-154688944 * | peak26118 55 179.08 [26119] chrX 154689871-154690497 * | peak26119 91 296.94 [26120] chrX 155242821-155243301 * | peak26120 83 270.31 [26121] chrX 155256146-155256507 * | peak26121 47 154.07
  5. Chip-seqデータとCpG配列データの重ね合わせ コード # get the peaks that overlap with CpG

    islands subsetByOverlaps(pk1.gr,cpgi.gr) 出力結果 ## GRanges object with 44 ranges and 5 metadata columns: ## seqnames ranges strand | name score signalValue ## <Rle> <IRanges> <Rle> | <character> <integer> <numeric> ## [1] chr21 9825360-9826582 * | peak14562 56 183.11 ## [2] chr21 9968469-9968984 * | peak14593 947 3064.92 ## [3] chr21 15755368-15755956 * | peak14828 90 291.90 ## [4] chr21 19191579-19192525 * | peak14840 290 940.03 ## [5] chr21 26979619-26980048 * | peak14854 32 104.67 ## ... ... ... ... . ... ... ... ## [40] chr21 46237464-46237809 * | peak15034 32 106.36 ## [41] chr21 46707702-46708084 * | peak15037 67 217.02 ## [42] chr21 46961552-46961875 * | peak15039 38 124.31 ## [43] chr21 47743587-47744125 * | peak15050 353 1141.58 ## [44] chr21 47878412-47878891 * | peak15052 104 338.78 ※この結果は第1引数のデータの中で、第2引数のデータと重複した物しか出力されない。 CpGの結果がどれが重複したか分かりにくい。
  6. Chip-seqデータとCpG配列データの重ね合わせ コード # Chip-seqとCpGの個々の重複を出力 findOverlaps(pk1.gr,cpgi.gr) 出力結果 ## Hits object with

    45 hits and 0 metadata columns: ## queryHits subjectHits ## <integer> <integer> ## [1] 14562 1 ## [2] 14593 3 ## [3] 14828 8 ## [4] 14840 13 ## [5] 14854 16 ## ... ... ... ## [41] 15034 155 ## [42] 15037 166 ## [43] 15039 176 ## [44] 15050 192 ## [45] 15052 200 ## ------- ## queryLength: 26121 / subjectLength: 205 補足 ①結果は各データのインデックスで出力される ➁Chip-seqの結果で、CpGと重複している個数を調べたところ、 subsetByOverlaps関数では44個だったのに、 findOverlaps関数では45個なのはなぜか? ➁に対する答え findOverlaps関数と異なり、subsetByOverlaps関数は ユニーク(固有)なデータしか返さないから。
  7. 最近傍のChip-seqデータとCpG配列データの距離 コード # 2種類のデータで、どのペアが一番近いか、とその距離(最近傍距離)の計算 dists=distanceToNearest(pk1.gr,cpgi.gr,select="arbitrary") dists 出力結果 ## Hits object

    with 620 hits and 1 metadata column: ## queryHits subjectHits | distance ## <integer> <integer> | <integer> ## [1] 14440 1 | 384188 ## [2] 14441 1 | 382968 ## [3] 14442 1 | 381052 ## [4] 14443 1 | 379311 ## [5] 14444 1 | 376978 ## ... ... ... . ... ## [616] 15055 205 | 26212 ## [617] 15056 205 | 27402 ## [618] 15057 205 | 30468 ## [619] 15058 205 | 31611 ## [620] 15059 205 | 34090 ## ------- ## queryLength: 26121 / subjectLength: 205 コード # 最近傍距離の結果をヒストグラムで表示 dist2plot=mcols(dists)[,1] hist(log10(dist2plot),xlab="log10(dist to nearest TSS)",main="Distances") 出力結果
  8. テキスト中では距離をそのまま対数変換している #対数変換前のChip-seqとCpGデータの距離データ SP1_CpGDist=dists@elementMetadata@listData[["distance"]] > SP1_CpGDist [1] 384188 382968 381052 379311

    376978 374944 374145 [8] 373317 372647 371220 370004 369582 367584 366264 [15] 365610 358845 339002 338427 337502 336888 257422 [22] 253151 250232 243192 240688 240311 239926 238821 …… [582] 20641 25538 49421 50280 43945 32914 0 [589] 555 5001 0 0 13507 2372 0 [596] 11905 34066 0 7656 0 7713 72171 [603] 99759 164844 59426 2383 1065 1323 18214 [610] 23869 0 17058 0 9697 22640 26212 [617] 27402 30468 31611 34090 #対数変換後のChip-seqとCpGデータの距離データ log10(SP1_CpGDist) [1] 5.584544 5.583162 5.580984 5.578995 5.576316 5.573966 [7] 5.573040 5.572078 5.571298 5.569631 5.568206 5.567711 [13] 5.565357 5.563794 5.563018 5.554907 5.530202 5.529465 [19] 5.528276 5.527486 5.410646 5.403380 5.398343 5.385949 …… [583] 4.407187 4.693912 4.701395 4.642909 4.517381 -Inf [589] 2.744293 3.699057 -Inf -Inf 4.130559 3.375115 [595] -Inf 4.075729 4.532321 -Inf 3.884002 -Inf [601] 3.887223 4.858363 4.998952 5.217073 4.773976 3.377124 [607] 3.027350 3.121560 4.260405 4.377834 -Inf 4.231928 [613] -Inf 3.986637 4.354876 4.418500 4.437782 4.483844 [619] 4.499838 4.532627 対数変換により、Chip-seqとCpGデータで一致している配列は 距離が本来0のところを無限にすり替えられてしまった。
  9. ①プロモーター配列の定義作業1 転写開始地点を含むbedファイルの読み込み コード # データを読み込んでGrangesオブジェクトに変換 filePathRefseq=system.file("extdata", "refseq.hg19.chr21.bed", package="compGenomRData") ref.df =

    read.table(filePathRefseq, header = FALSE, stringsAsFactors=FALSE) ref.gr=GRanges(seqnames=ref.df[,1], ranges=IRanges(start=ref.df[,2], end=ref.df[,3]), strand=ref.df[,6],name=ref.df[,4]) 補足 NR: RNA配列 NM: mRNA配列 出力結果 >ref.gr GRanges object with 571 ranges and 1 metadata column: seqnames ranges strand | name <Rle> <IRanges> <Rle> | <character> [1] chr21 41384342-42219039 - | NR_073202 [2] chr21 41384342-42219039 - | NM_001389 [3] chr21 41384342-42219039 - | NM_001271534 [4] chr21 17442841-17982094 + | NR_027790 [5] chr21 17566698-17982094 + | NR_027791 ... ... ... ... . ... [567] chr21 48055506-48075276 + | NM_001242865 [568] chr21 48055506-48085155 + | NM_001242864 [569] chr21 48018530-48025035 - | NM_006272 [570] chr21 48055506-48085155 + | NM_001535 [571] chr21 48055506-48085155 + | NM_206962
  10. ①プロモーター配列の定義作業2 転写開始地点を設定 コード # tss (transcription start site)用のオブジェクト作成 tss.gr=ref.gr #

    strandの+と-で転写開始点を設定 end(tss.gr[strand(tss.gr)=="+",]) =start(tss.gr[strand (tss.gr)=="+",]) start(tss.gr[strand(tss.gr)=="-",])=end(tss.gr[strand (tss.gr)=="-",]) # トランスクリプションバリアントなどの重複を削除 tss.gr=tss.gr[!duplicated(tss.gr),] 出力結果 >tss.gr GRanges object with 378 ranges and 1 metadata column: seqnames ranges strand | name <Rle> <IRanges> <Rle> | <character> [1] chr21 42219039 - | NR_073202 [2] chr21 17442841 + | NR_027790 [3] chr21 17566698 + | NR_027791 [4] chr21 30396937 + | NM_001001992 [5] chr21 27543138 - | NM_001204303 ... ... ... ... . ... [374] chr21 47705236 - | NM_003906 [375] chr21 47743785 - | NM_058180 [376] chr21 47882383 + | NR_046400 [377] chr21 48055506 + | NM_001242866 [378] chr21 48025035 - | NM_006272 ------- seqinfo: 1 sequence from an unspecified genome; no seqlengths https://bio-megane.com/ http://www.seibutsushi.net/blog/2008/08/535.html ・Irangesが範囲から1点に変わっている ・データ数が571→378に減っている
  11. ①プロモーター配列の定義作業3 転写開始地点前後1 kbpをプロモーター配列と定義 コード #転写開始地点の前後1 kbpをプロモーター配列と設定 promoter.gr=tss.gr start(promoter.gr)=start(promoter.gr)-1000 end(promoter.gr) =end(promoter.gr)+1000

    promoter.gr=promoter.gr[seqnames(promoter.gr) =="chr21"] 出力結果 >promoter.gr GRanges object with 378 ranges and 1 metadata column: seqnames ranges strand | name <Rle> <IRanges> <Rle> | <character> [1] chr21 42218039-42220039 - | NR_073202 [2] chr21 17441841-17443841 + | NR_027790 [3] chr21 17565698-17567698 + | NR_027791 [4] chr21 30395937-30397937 + | NM_001001992 [5] chr21 27542138-27544138 - | NM_001204303 ... ... ... ... . ... [374] chr21 47704236-47706236 - | NM_003906 [375] chr21 47742785-47744785 - | NM_058180 [376] chr21 47881383-47883383 + | NR_046400 [377] chr21 48054506-48056506 + | NM_001242866 [378] chr21 48024035-48026035 - | NM_006272 出力結果のIrangesにて、2001 bp の幅があることが確認できる
  12. ➁限定した領域だけのデータの読み出し bigwigファイルから、Chip-seqデータの読み込み コード library(rtracklayer) # Chip-seqデータの読み込み bwFile=system.file("extdata", "wgEncodeHaibTfbsA549.chr21.bw", package="compGenomRData") #

    プロモーター領域に合致するChip-seqデータの取り出し bw.gr=import(bwFile, which=promoter.gr) bw.gr 出力結果 ## GRanges object with 9205 ranges and 1 metadata column: ## seqnames ranges strand| score ## <Rle> <IRanges> <Rle>| <numeric> ## [1] chr21 9825456-9825457 * | 1 ## [2] chr21 9825458-9825464 * | 2 ## [3] chr21 9825465-9825466 * | 4 ## [4] chr21 9825467-9825470 * | 5 ## [5] chr21 9825471 * | 6 ## ... ... ... ... . ... ## [9201] chr21 48055809-48055856 * | 2 ## [9202] chr21 48055857-48055858 * | 1 ## [9203] chr21 48055872-48055921 * | 1 ## [9204] chr21 48055944-48055993 * | 1 ## [9205] chr21 48056069-48056118 * | 1 ※Scoreは各Chip-seqのリード数を表している
  13. マッピングされた配列の可視化と要約のための Viewオブジェクトの作成 出力結果 ## RleViewsList object of length 1: ##

    $chr21 ## Views on a 48129895-length Rle subject ## ## views: ## start end width ## [1] 42218039 42220039 2001 [2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [2] 17441841 17443841 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [3] 17565698 17567698 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [4] 30395937 30397937 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [5] 27542138 27544138 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 1 ...] ## [6] 27511708 27513708 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [7] 32930290 32932290 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [8] 27542446 27544446 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [9] 28338439 28340439 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## ... ... ... ... ... ## [370] 47517032 47519032 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [371] 47648157 47650157 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [372] 47603373 47605373 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [373] 47647738 47649738 2001 [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 ...] ## [374] 47704236 47706236 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [375] 47742785 47744785 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [376] 47881383 47883383 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [377] 48054506 48056506 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [378] 48024035 48026035 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] コード # bigwigフォーマットをRleListに変換 cov.bw=import(bwFile, which=promoter.gr,as = "RleList") #Viewオブジェクトに変換 myViews=Views(cov.bw,as(promoter.gr, "IRangesList")) myViews 補足 ・RlelistはRle vectorをリストにしたもの ・インデックスは、プロモーター領域のもの ・[]で囲まれた2001個の数字は、プロモーター領域 の各塩基と一致するChip-seqデータの数
  14. Viewオブジェクトを用いたChi-seqデータの可視化 作成したViewオブジェクトに含まれる全データ ## RleViewsList object of length 1: ## $chr21

    ## Views on a 48129895-length Rle subject ## ## views: ## start end width ## [1] 42218039 42220039 2001 [2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [2] 17441841 17443841 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [3] 17565698 17567698 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [4] 30395937 30397937 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [5] 27542138 27544138 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 1 ...] ## [6] 27511708 27513708 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [7] 32930290 32932290 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [8] 27542446 27544446 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [9] 28338439 28340439 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## ... ... ... ... ... ## [370] 47517032 47519032 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [371] 47648157 47650157 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [372] 47603373 47605373 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [373] 47647738 47649738 2001 [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 ...] ## [374] 47704236 47706236 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [375] 47742785 47744785 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [376] 47881383 47883383 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [377] 48054506 48056506 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [378] 48024035 48026035 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] コード # 5番目のプロモーター領域の可視化 plot(myViews[[1]][[5]],type="l") 出力結果
  15. Viewオブジェクトを用いたChi-seqデータの要約 作成したViewオブジェクトに含まれる全データ ## RleViewsList object of length 1: ## $chr21

    ## Views on a 48129895-length Rle subject ## ## views: ## start end width ## [1] 42218039 42220039 2001 [2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [2] 17441841 17443841 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [3] 17565698 17567698 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [4] 30395937 30397937 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [5] 27542138 27544138 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 1 ...] ## [6] 27511708 27513708 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [7] 32930290 32932290 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [8] 27542446 27544446 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [9] 28338439 28340439 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## ... ... ... ... ... ## [370] 47517032 47519032 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [371] 47648157 47650157 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [372] 47603373 47605373 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [373] 47647738 47649738 2001 [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 ...] ## [374] 47704236 47706236 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [375] 47742785 47744785 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [376] 47881383 47883383 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] ## [377] 48054506 48056506 2001 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...] ## [378] 48024035 48026035 2001 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...] コード # Viewオブジェクトの平均値の取得 head( viewMeans(myViews[[1]])) 出力結果 ## [1] 0.2258871 0.3498251 1.2243878 0.4997501 2.0904548 0.6996502 コード # Viewオブジェクトの最大値の取得 head(viewMaxs(myViews[[1]])) 出力結果 ## [1] 2 4 12 4 21 6