Pythonバイオ? Pythonバイオ/ツール?
467   2019-07-12 (金) 17:20:37

EdgeRまわり

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

まずは、EdgeRで二群間比較

参照 二群間比較(edgeR) | RNA-Seq による比較トランスクリプトーム解析

更に元は、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)')   # 所定の手続き1
r('d <- estimateCommonDisp(d)')  # 2
r('d <- estimateTagwiseDisp(d)')  #  3
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)

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-07-12 (金) 17:20:37 (445d)