Pythonバイオ/ツール/RNAseq④-Y EdgeRまわりをPythonからドライブの変更点

追加された行はこの色です。
削除された行はこの色です。
Pythonバイオ/ツール/RNAseq④-Y EdgeRまわりをPythonからドライブへ行く。
Pythonバイオ/ツール/RNAseq④-Y EdgeRまわりをPythonからドライブの差分を削除
[[Pythonバイオ]]　[[Pythonバイオ/ツール]]~
&counter();　　　&lastmod();~

*EdgeRまわり [#p0dc7d75]
[[Pythonバイオ/ツール/RNAseq④-X EdgeRまわり]]とほぼ同じことを、PythonからRへ放って行う。

**まずは、EdgeRで二群間比較 [#z436b5c9]
参照　[[二群間比較（edgeR） | RNA-Seq による比較トランスクリプトーム解析:https://bi.biopapyrus.jp/rnaseq/analysis/de-analysis/2g-edger.html]]

更に元は、[[edgeRUsersGuide.pdf:https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf]]

 import pandas as pd
 import rpy2
 from rpy2.robjects.packages import importr
 print(rpy2.__version__)
 from rpy2.robjects import pandas2ri
 pandas2ri.activate()
 from rpy2.robjects import r, StrVector, FactorVector
 
 base = importr('base')
 base.source("http://www.bioconductor.org/biocLite.R")
 biocinstaller = importr("BiocInstaller")  # 最初だけ
 biocinstaller.biocLite("edgeR")  # 最初だけ
 
 edgeR = importr("edgeR")
 
 count = pd.read_csv('../all.counts.txt', sep='\t', skiprows=1)  # 初めはGeneidを付けておく
 count = count[['Geneid', 'SRR453566.sorted.bam','SRR453567.sorted.bam','SRR453568.sorted.bam','SRR453569.sorted.bam','SRR453570.sorted.bam', 'SRR453571.sorted.bam']]
 count.columns = ['Geneid', "batch_1","batch_2","batch_3", "chemostat_1", "chemostat_2", "chemostat_3"]
 
 r.assign("count", pandas2ri.py2ri(count))  # countをRへアップする
 r('genelist <- count[,c("Geneid")]')    # R上でGeneidを覚えておく
 r('count <- count[, colnames(count) != "Geneid"]')  # R上でGeneid列を消す
 
 group = FactorVector(StrVector(['b', 'b', 'b', 'c', 'c', 'c']))   # groupを作っておく
 r.assign("group", group)
 
 r('d <- DGEList(counts = count, group = group, genes = genelist)')  # countとgroupからDEGlistを作る、genesにGeneid(～genelist)を入れる
 r('d <- calcNormFactors(d)')   # 所定の手続き１
 r('d <- estimateCommonDisp(d)')  # ２
 r('d <- estimateTagwiseDisp(d)')  #  ３
 r('result <- exactTest(d)')
 
 r('result0 <- result$table')  # R側でtableだけにする。全部持ってくると自動変換時に pandasクラスではなくrpy2クラスにしかならない
 r('result0ix <- result$genes')  # result$table にはGeneidが入っていない
 r0 = r['result0']
 print('r0\n', r0.head())
 r0ix = r['result0ix']
 r0.index = r0ix['genes']   # r0にindexとしてGeneidを付ける
 print('result0\n', r0.head())
 
 print('pvalues\n', (r0.sort_values(['PValue', 'logFC'], ascending=[True, False])).head(20))
 ## sort_valuesで複数キーを指定できる。更にascending条件が違ってもよい
 
 r('r2 <- topTags(result, n=nrow(count))')   # R上でtopTagsを通す。
 r('r22 <- r2$table')    # 結果のtable部分だけを抽出する
 r2 = r['r22']           # Python側へ持ってくる
 print('r2\n', r2.head())



 ################################
 # 完全に最初から
 count = pd.read_csv('../all.counts.txt', sep='\t', skiprows=1)
 count = count[['Geneid', 'SRR453566.sorted.bam','SRR453567.sorted.bam','SRR453568.sorted.bam','SRR453569.sorted.bam','SRR453570.sorted.bam', 'SRR453571.sorted.bam']]
 count.columns = ['Geneid', "batch_1","batch_2","batch_3", "chemostat_1", "chemostat_2", "chemostat_3"]
 r.assign("count", pandas2ri.py2ri(count))
 r('genelist <- count[,c("Geneid")]')
 r('count <- count[, colnames(count) != "Geneid"]')
 group = FactorVector(StrVector(['b', 'b', 'b', 'c', 'c', 'c']))
 r.assign("group", group)
 
 r('d <- DGEList(counts = count, group = group, genes=genelist)')
 r('row.names(d$count) <- genelist')
 r('row.names(d$counts) <- genelist')
 r('design <- model.matrix(~ group)')
 r('d <- estimateGLMCommonDisp(d, design)')
 r('d <- estimateGLMTrendedDisp(d, design)')
 r('d <- estimateGLMTagwiseDisp(d, design)')
 r('fit <- glmFit(d, design)')
 r('lrt <- glmLRT(fit, coef = 2)')   # LRTを出力
 
 r('lrt0 <- topTags(lrt, n = nrow(count))')   # topTagを通す。
 # print('lrt0\n', r['lrt0'])
 r('lrtgenes <- lrt0$table[, "genes"]')
 lrtgenes = r['lrtgenes']
 r('lrttable <- lrt0$table')    # table部分だけを抽出する
 lrttable = r['lrttable']    # Pythonの側へ持ってくる
 print(type(lrttable))
 print('lrttable\n', lrttable.head(20))
 lrttable.to_csv("lrttable.csv", sep='\t')
 
 #decideTestsでGeneidが出力に渡されない。
 
 r('dt <- decideTests(lrt)')
 dt = pd.DataFrame(r['dt'])
 dt.index = lrtgenes
 print('dt\n', dt)

-[[ExactTests:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/exactTest]]
-[[topTags function | R Documentation:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/topTags]]
-[[DGELRT-class function | R Documentation:https://www.rdocumentation.org/packages/edgeR/versions/3.10.5/topics/DGELRT-class]]
-[[DGEList function | R Documentation:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/DGEList]]
-[[glmFit function | R Documentation:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/glmFit]]
-[[decideTests function | R Documentation: https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/decideTests]]
-[[decideTests function | R Documentation:https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/decideTests]]
-[[TestResults-class function | R Documentation:https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/TestResults-class]]

-[[dimnames function | R Documantation:https://www.rdocumentation.org/packages/base/versions/3.6.1/topics/dimnames]]

-[[Rpy2の使い方の備忘録 - Qiita:https://qiita.com/MTNakata/items/04a045f4beb0fa9e3a4d]]
-[[rpy2 / R interface &#8212; pandas 0.22.0 documentation:https://pandas.pydata.org/pandas-docs/version/0.22/r_interface.html]]
-[[Python と R の違い (データフレーム編) &#8211; Python でデータサイエンス:https://pythondatascience.plavox.info/python%E3%81%A8r%E3%81%AE%E9%81%95%E3%81%84/%E3%83%87%E3%83%BC%E3%82%BF%E3%83%95%E3%83%AC%E3%83%BC%E3%83%A0]]
-[[RのデータフレームをPythonのpandasデータフレームにする方法 - Qiita:https://qiita.com/tatekaiten/items/e29ac74660b249fdd42f]]
-[[Converting DataFrames to and from R - python:http://www.jamesloach.com/python/pandas2ri.html]]
-[[pythonから統計解析ソフトRを利用する方法・第１回:http://cup.sakura.ne.jp/pyrcmd/pyrcmd01.htm]]
-[[pythonから統計解析ソフトRを利用する方法・第２回:http://cup.sakura.ne.jp/pyrcmd/pyrcmd02.htm]]
-[[Python rpy2 で pandas の DataFrame を R の data.frame に変換する - StatsFragments:http://sinhrks.hatenablog.com/entry/2014/10/16/224948]]



-[[R-Source データ型とデータ構造:http://cse.naro.affrc.go.jp/takezawa/r-tips/r/25.html]]
-[[R_列名で列を削除する - RとLinuxと...:http://rmecab.jp/wiki/index.php?R_%CE%F3%CC%BE%A4%C7%CE%F3%A4%F2%BA%EF%BD%FC%A4%B9%A4%EB]]
-[[R: Get and Set Row Names for Data Frames:https://stat.ethz.ch/R-manual/R-devel/library/base/html/row.names.html]]
-[[列名を指定してデータを抽出する【R】 - QIITA:https://qiita.com/hitsujisuke/items/d71ee40daa0786ae1680]]

-[[Pandasで列データをソートするsort_values関数の使い方 - DeepAge:https://deepage.net/features/pandas-sort-values.html]]

-[[p値と有意水準 | ブログ | 統計WEB:https://bellcurve.jp/statistics/blog/14004.html]]
-[[有意 - Wikipedia:https://ja.wikipedia.org/wiki/%E6%9C%89%E6%84%8F]]
-[[尤度比検定 - Wikipedia:https://ja.wikipedia.org/wiki/%E5%B0%A4%E5%BA%A6%E6%AF%94%E6%A4%9C%E5%AE%9A]]
-[[正規母集団（分散既知）の平均に関する検定 | 有意に無意味な話:http://starpentagon.net/analytics/z_test_via_lrt/]]
Pythonバイオ/ツール/RNAseq④-Y EdgeRまわりをPythonからドライブ の変更点

Pythonバイオ/ツール/RNAseq④-Y EdgeRまわりをPythonからドライブの変更点
