[[Pythonバイオ]] [[Pythonバイオ/ツール]]~
&counter();   &lastmod();~

*EdgeRまわり [#l8b20a9d]
**まずは、EdgeRで二群間比較 [#ee83bcbb]
参照 [[二群間比較(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]]

このページによると、
 edgeR のバージョンが異なると、解析結果も少し異なってくる。このページでは 3.10.0 を利用してる。
ということだが、どうだろうか?

このページに従って、自前データで試す。
 > source("http://bioconductor.org/biocLite.R")
 > biocLite("edgeR")
 # 依存ライブラリが多いので、今まで使っていないと結構時間がかかる
 > library(edgeR)

 > count <- read.table("../all.counts.txt", sep="\t", header = T, row.names = 1)
 > head(count)
                   Chr Start   End Strand Length SRR453566.sorted.bam
 gene_0001 NC_001133.9  1807  2169      -    363                    0
 gene_0002 NC_001133.9  2480  2707      +    228                    0
 gene_0003 NC_001133.9  7235  9016      -   1782                    0
 > count <- count[,c('SRR453566.sorted.bam','SRR453567.sorted.bam','SRR453568.sorted.bam','SRR453569.sorted.bam','SRR453570.sorted.bam', 'SRR453571.sorted.bam')]
 > colnames(count) <- c("batch_1","batch_2","batch_3", "chemostat_1", "chemostat_2", "chemostat_3")
 > head(count)
           batch_1 batch_2 batch_3 chemostat_1 chemostat_2 chemostat_3
 gene_0001       0       2       6           0           0           1
 gene_0002       0       0       0           0           0           0
 gene_0003       0       0       0           0           0           0

 > count <- as.matrix(count)
 > dim(count)   # 6420   6
 > group <- factor(c("b", "b", "b", "c", "c", "c"))

<exact test>
 # まず、countデータとgroupを1つの変数に入れる(ラベル付け)
 > d <- DGEList(counts = count, group = group)
 # 発現数の非常に少ないgeneを除去  [[edgeRUsersGuide:https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf]] p55
 > keep <- rowSums(cpm(d)>1) >= 2    # CPM(count per million)が1以上が2つ以上
 > d <- d[keep, , keep.lib.sizes=FALSE]
 > dim(d)   # 5899   6

 # 正規化
 > d <- calcNormFactors(d)
 # 分散の推測
 > d <- estimateCommonDisp(d)
 > d <- estimateTagwiseDisp(d)
 
 # で、exact testの検定
 > result <- exactTest(d)
 > dim(result)    # 5899    3
 
 # 結果表示のためには
 > topTags(result)
 Comparison of groups:  c-b 
               logFC    logCPM PValue FDR
 gene_2989 11.273497  8.543906      0   0
 gene_4667  9.683678  9.715407      0   0
 gene_3497  8.659467 11.419660      0   0
 gene_3421  8.501696  7.911698      0   0
 gene_4882  7.545279 10.715972      0   0
 gene_0615  7.399665  8.409708      0   0
 gene_1031  7.233161  7.473834      0   0
 gene_5986  7.169033  6.784983      0   0
 gene_2035  6.893378  8.849378      0   0
 gene_0017  6.841963 10.445186      0   0
 # logFCとは:log fold-change、wild type (WT) と knock out (KO) のシグナル値を比較するとき、log2(KO)/log2(WT)
 # locCPMとは:log counts per million
 # FDR:  False Discovery Rate
 
 # 結果をファイルに格納する
 > table <- as.data.frame(topTags(result, n = nrow(count)))
 > write.table(table, file = "result.txt", col.names = T, row.names = T, sep = "\t")
 # ファイルの、途中のエントリーの例
         "logFC"     "logCPM"        "PValue"        "FDR"
 "gene_2901"     0.167058457314462       6.8615588901153 0.0493154259214582      0.0667381733220193
 "gene_0340"     0.289892217524922       5.33859450515354        0.0494251531593589      0.0668713253410684
 "gene_0519"     0.189496458586853       6.08722498147699        0.0494834466115673      0.0669348432840255
 "gene_0954"     0.256044030124387       5.65989391337004        0.0496916298717389      0.0672010372795479
 "gene_3371"     -0.205143944867366      9.08724668881671        0.0498370683360759      0.0673822750663561
 "gene_0287"     -0.255406659792859      5.95521336812612        0.0501243894523602      0.0677552184645904
 "gene_5879"     0.208805500563976       5.93528971411533        0.0503683786485931      0.0680694308472052
 "gene_5595"     -0.172878902469866      6.9102371649123 0.0503961888612053      0.0680914150463238
 "gene_2980"     -0.174306726898056      6.36273018338751        0.0505773479860026      0.0683201796940756
 "gene_3461"     0.190387033748589       6.37468042292314        0.0505886667068524      0.0683201796940756

[[exactTest From edgeR v3.14.0:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/exactTest]]


[[DGEList:https://www.rdocumentation.org/packages/edgeR/versions/3.10.5/topics/DGEList-class]]

 # 結果を図示する
 > is.DEG <- as.logical(table$FDR < 0.01)
 > DEG.names <- rownames(table)[is.DEG]
 > plotSmear(result, de.tags = DEG.names)
&ref(PlotSmear.png);
かなりはみ出している?

<尤度比検定>

 > design <- model.matrix(~ group)
 > design
   (Intercept) groupc
 1           1      0
 2           1      0
 3           1      0
 4           1      1
 5           1      1
 6           1      1
 attr(,"assign")
 [1] 0 1
 attr(,"contrasts")
 attr(,"contrasts")$group
 [1] "contr.treatment"
 
 # 解析する
 # まずもう一度データをcountから作り直し
 > d <- DGEList(counts = count, group = group)
 # 前と同様に、発現数の非常に少ないgeneを除去
 > keep <- rowSums(cpm(d)>1) >= 2    # CPM(count per million)が1以上が2つ以上
 > d <- d[keep, , keep.lib.sizes=FALSE]
 # ここから、分布(分散)の推定 GLMを使う
 > d <- calcNormFactors(d)
 > d <- estimateGLMCommonDisp(d, design)
 > d <- estimateGLMTrendedDisp(d, design)
 > d <- estimateGLMTagwiseDisp(d, design)
 # 分布にフィット
 > fit <- glmFit(d, design)
 # 尤度比検定
 このとき、full model と reduced model を比べるので、この情報を関数に与える。
 glmFit を利用して作成した fit には full model のデータ含まれている。従って、あとは
 reduced model の情報だけを与えればよい。ここで、full model と reduced model の
 デザイン行列を見比べてみると、reduced model は full model のデザイン行列の 2 列を
 すべて 0 にしたものに等しいことがわかる。そこで、coef = 2 と指定して、reduced model
 の情報を与える。
 > lrt <- glmLRT(fit, coef = 2)
 > dim(lrt)  # 5899   4
 # 結果は
 > topTags(lrt)
 Coefficient:  groupc 
               logFC    logCPM       LR PValue FDR
 gene_2989 11.272397  8.543271 2260.090      0   0
 gene_3497  8.658186 11.419157 1934.825      0   0
 gene_3421  8.500384  7.911092 2271.734      0   03859
 gene_4882  7.543818 10.715354 2673.692      0   0
 gene_0615  7.398431  8.409080 1897.589      0   0
 
 # ファイルに出力
 > lrttable <- as.data.frame(topTags(lrt, n = nrow(count)))
 > write.table(lrttable, file="lrtresult.txt", col.names=T, row.names=T, sep="\t")
                logFC           logCPM       LR        PValue              FDR
 "gene_3295"     0.247645680731094       5.40510860572214        3.87945550718096        0.0488804073717157      0.065772245229414
 "gene_4943"     -0.194295943127542      7.18848127443337        3.8631344790163 0.0493580361534112      0.0663997845539276
 "gene_3461"     0.190360832039291       6.37466110822439        3.85464691865404        0.0496083672460962      0.0667213311410673
 "gene_6181"     -0.204622554470669      5.77230370218388        3.85415389282258        0.0496229496429794      0.0667257305548064
 "gene_1983"     0.856380812164874       1.74233409044916        3.8422391885909 0.0499767355704516      0.0671861356267307
 "gene_5577"     -0.263126831966537      5.60034347444897        3.84011138664067        0.050040196820857       0.0672561223618673
 "gene_3658"     -0.194852891673291      6.72394023500093        3.82819449899496        0.0503971930016531      0.0677205105960711
 "gene_2310"     -0.197601718589681      6.24609836497735        3.82407278015195        0.0505212936731495      0.0678718085579387
 "gene_2980"     -0.17430474127026       6.36273796901708        3.81688233743943        0.0507385643479778      0.0681481764773954
 "gene_0447"     0.325651368079434       4.40033732555139        3.81349708835367        0.0508411968215156      0.0682704803209925
 
増減の全体を見ると ([[EdgeUsersGuide:https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf]] p58)
 >summary(decideTests(lrt))    # 1% FDR
        groupc
 Down     2072
 NotSig   1620
 Up       2207

MD図を描くことができる  ([[EdgeUsersGuide:https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf]] p58)
 > plotMD(lrt)
 > abline(h=c(-1,1), col="blue")

&ref(MDchart.png);

[[glmFit/glmLRT:https://www.rdocumentation.org/packages/edgeR/versions/3.14.0/topics/glmFit]]


<GO analysis>

 > biocLite("GO.db")
 > biocLite("org.Hs.eg.db")
 > go <- goana(lrt)
  goana.default(de = DEGenes, universe = universe, ...) でエラー: 
   No annotated genes found in universe

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS