Pythonバイオ? Pythonバイオ/ツール?
1420 2019-07-12 (金) 17:20:37
Pythonバイオ/ツール/RNAseq��-X EdgeRまわりとほぼ同じことを、PythonからRへ放って行う。
参照 二群間比較（edgeR） | RNA-Seq による比較トランスクリプトーム解析
更に元は、edgeRUsersGuide.pdf
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
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)
