Pythonバイオ? Pythonバイオ/ツール?
233   2019-04-20 (土) 14:22:38

発現解析

発現量解析 | RNA-Seq を利用した発現変動遺伝子の検出 BioPapyrusサイトでのまとめ

正規化の説明




Ballgown

R上でのインストール

source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")

使おう データのロード

library(ballgown)

別の資料によると以下も導入(詳細未確認)

library(RSkittleBrewer)
library(genefilter)
library(dplyr)

ballgownを使う。

データ構造が次のような形になっているという前提

extdata/
   sample01/
       e2t.ctab
       e_data.ctab
       i2t.ctab
       i_data.ctab
       t_data.ctab
   sample02/
       e2t.ctab
       e_data.ctab
       i2t.ctab
       i_data.ctab
       t_data.ctab
   ...
   sample20/
       e2t.ctab
       e_data.ctab
       i2t.ctab
       i_data.ctab
       t_data.ctab

ここではstringtieを処理した時に、

$ stringtie -e -B -p 16 -G s288c_e.gff -o ballgown/SRR453566/SRR453566.gtf SRR453566.sorted.bam
$ stringtie -e -B -p 16 -G s288c_e.gff -o ballgown/SRR453567/SRR453567.gtf SRR453567.sorted.bam

としたので、下のようになっている。マニュアルAccessing assembly data参照。

ballgown/
   SRR453566/
       SRR453566.gtf
       e2t.ctab
       e_data.ctab
       i2t.ctab
       i_data.ctab
       t_data.ctab
   SRR453567/
       SRR453567.gtf
       e2t.ctab
       e_data.ctab
       i2t.ctab
       i_data.ctab
       t_data.ctab

これをballgownに食わせる。

bg = ballgown(dataDir="~/src/RNAseq-Saccha/Saccha/ballgown", samplePattern='SRR', meas='all')

dataDirはデータの置いてあるディレクトリ(ballgown)、samplePatternは中のディレクトリ中のサンプルの共通接頭辞。ここではSRRxxxxxなのでSRRにした。measは不明。

この処理が終わると、bgが使えるようになる。bg中のstructureには、Exon, intron, and transcript structuresがある。それぞれを取出してみるには、

structure(bg)$exon
GRanges object with 6801 ranges and 2 metadata columns:
            seqnames      ranges strand |        id transcripts
               <Rle>   <IRanges>  <Rle> | <integer> <character>
     [1] NC_001133.9   1807-2169      - |         1           1
     [2] NC_001133.9   2480-2707      + |         2           2
     [3] NC_001133.9   7235-9016      - |         3           3
     [4] NC_001133.9 11565-11951      - |         4           4
     [5] NC_001133.9 12046-12426      + |         5           5
     ...         ...         ...    ... .       ...         ...
  [6797] NC_001224.1 78089-78162      - |      6797        6441
  [6798] NC_001224.1 78533-78608      + |      6798        6442
  [6799] NC_001224.1 79213-80022      + |      6799        6443
  [6800] NC_001224.1 85035-85112      + |      6800        6444
  [6801] NC_001224.1 85295-85777      + |      6801        6445
  -------
  seqinfo: 17 sequences from an unspecified genome; no seqlengths
structure(bg)$trans
GRangesList object of length 6445:
$1 
GRanges object with 1 range and 2 metadata columns:
         seqnames    ranges strand |        id transcripts
            <Rle> <IRanges>  <Rle> | <integer> <character>
  [1] NC_001133.9 1807-2169      - |         1           1

$2 
GRanges object with 1 range and 2 metadata columns:
         seqnames    ranges strand | id transcripts
  [1] NC_001133.9 2480-2707      + |  2           2

$3 
GRanges object with 1 range and 2 metadata columns:
         seqnames    ranges strand | id transcripts
  [1] NC_001133.9 7235-9016      - |  3           3

...
<6442 more elements>
-------
seqinfo: 17 sequences from an unspecified genome; no seqlengths

など。

次に、exprスロットを取出す。t/e/i/gに対して、texpr, eexpr, iexpr, gexprが対応し、 それぞれに取り出したいものをtexpr(bg, 'FPKM')のように指定する。

具体的には、次のようなものが取り出せる。

transcript_fpkm = texpr(bg, 'FPKM')
transcript_cov = texpr(bg, 'cov')
whole_tx_table = texpr(bg, 'all')
exon_mcov = eexpr(bg, 'mcov')
junction_rcount = iexpr(bg)
whole_intron_table = iexpr(bg, 'all')
gene_expression = gexpr(bg)

たとえば、transcript_fpkmは

>transcript_fpkm
     FPKM.SRR453566 FPKM.SRR453567
1          0.251292       1.106585
2          0.000000       0.000000
3          0.000000       0.000000
4          0.071457       0.042556
5          1.218477       2.062575
6          0.000000       0.000000
7          0.000000       0.000000
8          0.000000       0.000000
9          1.410877       1.360651
10         9.473732       8.717489
11        18.461433      12.996510
12       197.994659     219.502182
13        95.170197     101.183212
14        35.876900      39.537365
15        23.731741      23.713844
...
 [ reached getOption("max.print") -- 5945 行を無視しました ] 

となり、gene_expressionは

> gene_expression
          FPKM.SRR453566 FPKM.SRR453567
gene_0001       0.251292       1.106585
gene_0002       0.000000       0.000000
gene_0003       0.000000       0.000000
gene_0004       0.071457       0.042556
gene_0005       1.218477       2.062575
gene_0006       0.000000       0.000000
gene_0007       0.000000       0.000000
gene_0008       0.000000       0.000000
gene_0009       1.410877       1.360651
gene_0010       9.473732       8.717489
gene_0011      18.461433      12.996510
gene_0012     197.994659     219.502182
gene_0013      95.170197     101.183212
gene_0014      35.876900      39.537365
gene_0015      23.731741      23.713844
...
 [ reached getOption("max.print") -- 5945 行を無視しました ] 

のような結果が得られる。

indexスロットは、もう少し勉強必要。indexes(bg)に対して、indexes(bg)$e2t, indexes(bg)$i2t, indexes(bg)$t2gなどが可能。(テーブルを見ているだけ?)

> indexes(bg)$e2t
     e_id t_id
1       1    1
2       2    2
3       3    3
4       4    4
5       5    5
6       6    6
7       7    7
8       8    8
9       9    9
10     10   10
11     11   11
12     12   12
13     13   13
14     14   14
15     15   15
> indexes(bg)$i2t
    i_id t_id
1      1   41
2      2   66
3      3   68
4      4   71
5      5   86
6      6  104
7      7  117
8      8  124
9      9  129
10    10  154
11    11  155
12    12  163
13    13  173
14    14  188
15    15  189
> indexes(bg)$t2g
     t_id      g_id
1       1 gene_0001
2       2 gene_0002
3       3 gene_0003
4       4 gene_0004
5       5 gene_0005
6       6 gene_0006
7       7 gene_0007
8       8 gene_0008
9       9 gene_0009
10     10 gene_0010
11     11 gene_0011
12     12 gene_0012
13     13 gene_0013
14     14 gene_0014
15     15 gene_0015

あと、マニュアルによると、phenotype情報を見るコンポーネントpDataがある。 pDataの情報は手で作る必要があり、いろいろと(記述順番とかの)設定制限があるらしい。

描画については、マニュアル Plotting transcript structuresによると、次の通り。

plotTranscriptsを使うと、描画される。

plotTranscripts(gene='XLOC_000454', gown=bg, samples='sample12', 
   meas='FPKM', colorby='transcript', 
   main='transcripts from gene XLOC_000454: sample 12, FPKM')

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-04-20 (土) 14:22:38 (127d)