[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[Pythonバイオ]] [[Pythonバイオ/ツール]]~
&counter(); &lastmod();~
*発現解析 [#q2eb5d88]
[[発現量解析 | RNA-Seq を利用した発現変動遺伝子の検出:https://bi.biopapyrus.jp/rnaseq/analysis/]] BioPapyrusサイトでのまとめ
**正規化の説明 [#qae69fd7]
-[[RNA-Seq | 遺伝子発現量解析:https://bi.biopapyrus.jp/rnaseq/]] ⇒ [[FPKM / RPKM | RNA-Seq リードカウントデータに転写物の長さなどを補正した発現量:https://bi.biopapyrus.jp/rnaseq/analysis/normalizaiton/fpkm.html]] の中にRを利用して FPKM を計算する方法
-[[次世代シーケンサーでの遺伝子発現量解析 | PictBio:https://www.pictbio.com/tips/2554.html]]
-[[機能ゲノム学（第6回）:http://www.iu.a.u-tokyo.ac.jp/~kadota/20110929_kadota.pdf]] 門田先生による各手法の比較実験(2011/09/29)~
RPM(Reads per million mapped reads) ⇒ RPKM(Reads per kilobase of exon per million mapped reads) / TMM正規化法 / 門田法
-[[機能ゲノム学（第6回）:https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2011_illumina_rna-seq_session3.pdf]](2011/11/17)
-[[Question: What is the reason why we usually use normalized values from RNA-Seq (FPKM, RPKM, etc.) ?:https://www.biostars.org/p/270537/]] (BioStars)~
R(F)PKM/TPM values are used to normalize read counts by library size (total number of reads you have in a given RNAseq experiment) and the length of the feature (gene/transcript). But remember that commonly used software for differential expression analysis (DESEQ2/EdgeR) are using raw counts instead of normalized values (they do their internal normalization steps).
-[[マイクロアレイより少し複雑な、RNA-Seqのデータ解析手順 | Subio:https://www.subioplatform.com/ja/info_technical/293/an-rna-seq-data-analysis-procedure-a-bit-more-complicated-than-microarrays]]
-[[An integrative method to normalize RNA-Seq data:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4067528/]](2014/6/14)~
Since RNA-Seq emergence, a number of normalization methods have been developed to address one or two of the different biases [1-12,14]. Our aim was to develop an integrated method able to correct all these sources of bias.
-----------------------------------
~
~
**Ballgown [#u8416342]
[[R上でのインストール:https://github.com/alyssafrazee/ballgown#installation]]
source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")
使おう� �[[データのロード:https://github.com/alyssafrazee/ballgown#loading-data-into-r]]
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:https://github.com/alyssafrazee/ballgown#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:https://github.com/alyssafrazee/ballgown#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')
終了行:
[[Pythonバイオ]] [[Pythonバイオ/ツール]]~
&counter(); &lastmod();~
*発現解析 [#q2eb5d88]
[[発現量解析 | RNA-Seq を利用した発現変動遺伝子の検出:https://bi.biopapyrus.jp/rnaseq/analysis/]] BioPapyrusサイトでのまとめ
**正規化の説明 [#qae69fd7]
-[[RNA-Seq | 遺伝子発現量解析:https://bi.biopapyrus.jp/rnaseq/]] ⇒ [[FPKM / RPKM | RNA-Seq リードカウントデータに転写物の長さなどを補正した発現量:https://bi.biopapyrus.jp/rnaseq/analysis/normalizaiton/fpkm.html]] の中にRを利用して FPKM を計算する方法
-[[次世代シーケンサーでの遺伝子発現量解析 | PictBio:https://www.pictbio.com/tips/2554.html]]
-[[機能ゲノム学（第6回）:http://www.iu.a.u-tokyo.ac.jp/~kadota/20110929_kadota.pdf]] 門田先生による各手法の比較実験(2011/09/29)~
RPM(Reads per million mapped reads) ⇒ RPKM(Reads per kilobase of exon per million mapped reads) / TMM正規化法 / 門田法
-[[機能ゲノム学（第6回）:https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2011_illumina_rna-seq_session3.pdf]](2011/11/17)
-[[Question: What is the reason why we usually use normalized values from RNA-Seq (FPKM, RPKM, etc.) ?:https://www.biostars.org/p/270537/]] (BioStars)~
R(F)PKM/TPM values are used to normalize read counts by library size (total number of reads you have in a given RNAseq experiment) and the length of the feature (gene/transcript). But remember that commonly used software for differential expression analysis (DESEQ2/EdgeR) are using raw counts instead of normalized values (they do their internal normalization steps).
-[[マイクロアレイより少し複雑な、RNA-Seqのデータ解析手順 | Subio:https://www.subioplatform.com/ja/info_technical/293/an-rna-seq-data-analysis-procedure-a-bit-more-complicated-than-microarrays]]
-[[An integrative method to normalize RNA-Seq data:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4067528/]](2014/6/14)~
Since RNA-Seq emergence, a number of normalization methods have been developed to address one or two of the different biases [1-12,14]. Our aim was to develop an integrated method able to correct all these sources of bias.
-----------------------------------
~
~
**Ballgown [#u8416342]
[[R上でのインストール:https://github.com/alyssafrazee/ballgown#installation]]
source("http://bioconductor.org/biocLite.R")
biocLite("ballgown")
使おう� �[[データのロード:https://github.com/alyssafrazee/ballgown#loading-data-into-r]]
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:https://github.com/alyssafrazee/ballgown#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:https://github.com/alyssafrazee/ballgown#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')
ページ名:
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
