[[Pythonバイオ]] [[Pythonバイオ/ツール]]~ &counter(); &lastmod();~ *発現量の解析 [#q2eb5d88] *発現解析 [#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]")