Hatena::Groupbioconductor

Bioconductorノート このページをアンテナに追加 RSSフィード

 | 

2007-06-18Bioconductorチュートリアル @ NIG

ケーススタディCELファイルから興味ある遺伝子リストの作成  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート を含むブックマーク はてなブックマーク -  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート のブックマークコメント

目次

つぎのfromcels.Rを参考にして、CELファイルから興味ある遺伝子リストの作成についてのケーススタディを説明します。

▼のついたソースコード片をRコンソールにコピーペーストしてすすめます。

ストーリー

  1. データ:ある条件下で発現変動のある遺伝子をねらった Affymetrix GeneChip による遺伝子発現データ
  2. シグナル/ノイズ比を改善するための前処理
  3. 発現変動遺伝子の順位付けのための適切な統計量の選択
  4. 候補遺伝子選択のためのフィルタの決定
  5. 関連するアノテーションの付加
  6. 出口:ウェブブラウザで閲覧できる発現変動のある遺伝子リスト

データ読み込み

library("affy")
library("SpikeInSubset")
data(spikein95)

このケーススタディでは spikein95 という AffyBatch クラスデータを利用します。このデータスパイクイン実験データなので、発現変動している遺伝子あきらかな状態になっています。


前処理

▼このケーススタディ用にデータのラベルを操作します。

pd <- data.frame(population = c(1,1,1,2,2,2),
                 replicate = c(1,2,3,1,2,3))
rownames(pd) <- sampleNames(spikein95)
vl <- list(population = "1 is control, 2 is treatment",
           replicate = "arbitrary numbering") 
phenoData(spikein95) <- new("phenoData", pData = pd, varLabels = vl)

RMAを適用します。

eset <- rma(spikein95)

△実行結果

Background correcting
Normalizing
Calculating Expression

RMAバックグラウンド補正、正規化、プローブレベルデータの要約をおこないます。ここで eset は exprSet クラスデータです。


データ次元を確認

e <- exprs(eset)
dim(e)

△実行結果

[1] 12626     6

行が遺伝子、列がアレイを示しています。

▼それぞれのアレイがどの集団(population)に属しているかを示す

pData(eset)

△実行結果

                 population replicate
1521a99hpp_av06           1         1
1532a99hpp_av04           1         2
2353a99hpp_av08           1         3
1521b99hpp_av06           2         1
1532b99hpp_av04           2         2
2353b99hpp_av08r          2         3

▼集団毎のインデクス・ベクトルを作成

Index1 <- which(eset$population == 1)
Index2 <- which(eset$population == 2)

遺伝子の順位付け

▼各遺伝子の発現変動の対数変換された平均値を計算する

d <- rowMeans(e[, Index2]) - rowMeans(e[, Index1])

▼発現量の対数変換された平均値を計算する

a <- rowMeans(e)

▼発現変動が二倍以上ある遺伝子の数を計算する

sum(abs(d) > 1)

△実行結果

[1] 1

一つしか無かった。


遺伝子の順位づけのための統計検定量の比較

発現変動が興味深い遺伝子を選択するために、いくつかの基準を比較する。

MAプロットデータをみる

plot(a, d, ylim = c(-1, 1), main = "A) MA-plot", pch = ".")

△実行結果

f:id:nakao_mitsuteru:20070616221318p:image


▼t統計量を計算する

library("genefilter")
tt <- rowttests(e, factor(eset$population))

▼t検定の結果をVolcanoプロットでみる

lod <- -log10(tt$p.value)
plot(d, lod, cex =.25, main = "B) Volcano plot for $t$-test")
abline(h = 2)

△実行結果

f:id:nakao_mitsuteru:20070616221317p:image

Volcanoプロットは、(対数オッズ化した)p値と平均の差を同時にみるためのものです。

▼t検定の結果をVolcanoプロットで拡大してみる

o1 <- order(abs(d), decreasing = TRUE)[1:25]
o2 <- order(abs(tt$statistic), decreasing = TRUE)[1:25]
o <- union(o1, o2)
plot(d[-o], lod[-o], cex = .25, xlim = c(-1, 1), ylim = range(lod), main = "C) Close up of B)")
points(d[o1], lod[o1], pch = 18, col = "blue")
points(d[o2], lod[o2], pch = 1, col = "red")

△実行結果

f:id:nakao_mitsuteru:20070616221316p:image

対数スケールの倍数変化の平均によって上位25位までに順位付けられた遺伝子は青いひし形で示し,p値の小さい順に選択した25遺伝子は赤い円で示している。

倍数変化によって選択される遺伝子とt統計量によって選択される遺伝子には比較的大きな違いがあることが読み取れる。


改良されたt統計量のひとつ(moderated t検定)を試す。

limmaを実行する

library("limma")
design <- model.matrix(~factor(eset$population))
fit <- lmFit(eset, design)
ebayes <- eBayes(fit)

▼moderated t検定の結果をVolcanoプロットでみる

lod <- -log10(ebayes$p.value[,2])
mtstat<- ebayes$t[,2]
o1 <- order(abs(d), decreasing = TRUE)[1:25]
o2 <- order(abs(mtstat), decreasing = TRUE)[1:25]
o <- union(o1, o2)
plot(d[-o], lod[-o], cex = .25, xlim = c(-1,1), ylim = c(0,4), main = "D) Volcano plot for moderated $t$-test")
points(d[o1], lod[o1], pch = 18, col = "blue")
points(d[o2], lod[o2], pch = 1, col = "red")

△実行結果

f:id:nakao_mitsuteru:20070616221314p:image


カットオフの選択

順位付けられた遺伝子をどこまで利用するかをきめるカットオフについて考えます。

遺伝子数が多い(12,626個)ので、0.01の有意水準でも 126個の偽陽性が期待されてしまいます。

▼p値が 0.01 以下の遺伝子数を計算

sum(tt$p.value <= 0.01)

△実行結果

[1] 46

比較

データスパイクイン実験なので、発現変動遺伝子の数は既知です。それを取り出して比較します。

data(spikein95)
spikedin <- colnames(pData(spikein95))
spikedIndex <- match(spikedin, geneNames(eset))

スパイクインの16遺伝子がそれぞれの統計量でどのような順位になっているかを示します。

d.rank <- sort(rank(-abs(d))[spikedIndex])
t.rank <- sort(rank(-abs(tt$statistic))[spikedIndex])
mt.rank <- sort(rank(-abs(mtstat))[spikedIndex])
ranks <- cbind(mt.rank, d.rank, t.rank)
rownames(ranks) <- NULL
ranks

△実行結果

      mt.rank d.rank t.rank
 [1,]       1      1      1
 [2,]       2      2      3
 [3,]       3      3      5
 [4,]       6      4      8
 [5,]       7      5     13
 [6,]       8      6     17
 [7,]       9      9     19
 [8,]      10     11     27
 [9,]      11     12     28
[10,]      12     14     29
[11,]      16     16     45
[12,]      25     53     70
[13,]      28     86     71
[14,]      48    226     93
[15,]      77    331    390
[16,]     465    689    900

比較的、modarated t検定が良い結果でした。


関連するアノテーションの付加

統計検定量の付加

▼modareted t統計量をつかって、上位10遺伝子をとりだします。

tab <- topTable(ebayes, coef = 2, adjust = "fdr", n = 10)

▼tabの上位5つを表示します。

tab[1:5,]

△実行結果

           ID  logFC AveExpr      t  P.Value adj.P.Val    B
779   1708_at -7.061    7.95 -73.53 7.82e-17  9.87e-13 8.65
6252 36202_at  0.853    9.37   9.98 4.94e-07  3.12e-03 4.59
6362 36311_at  0.832    8.56   8.36 3.02e-06  1.27e-02 3.57
3284 33264_at  0.712    4.92   7.43 9.67e-06  2.71e-02 2.84
2674 32660_at  0.655    8.68   7.36 1.07e-05  2.71e-02 2.77

▼Affymetrix IDを取り出します。

genenames <- as.character(tab$ID)
PubMedアブストラクトの付加

▼eset が必要とするアノテーションパッケージを調べる。

annotation(eset)

△実行結果

[1] "hgu95a"

▼必要なパッケージを読み込み、関連した論文PubMed アブストラクトを取得する。

library("hgu95av2")
library("XML")
library("annotate")
absts <- pm.getabst(genenames, "hgu95av2")

△実行結果

Read 8255 items
Read 4181 items
Read 4277 items
Read 1329 items
Read 1157 items
Read 8379 items
Read 22358 items
Read 3816 items
Read 11971 items
Read 3653 items

アブストラクトの一部をみる。

absts[[1]][[4]]

△実行結果

An object of class 'pubMedAbst':
Title: Absence of excitotoxicity-induced apoptosis in
     the hippocampus of mice lacking the Jnk3 gene.
PMID: 9349820
Authors: DD Yang, CY Kuan, AJ Whitmarsh, M Rincón,
     TS Zheng, RJ Davis, P Rakic, RA Flavell
Journal: Nature
Date: Oct 1997

▼関連した論文タイトルを表示する。

titl <- sapply(absts[[2]], articleTitle)
strwrap(titl, simplify = FALSE)

△実行結果

[[1]]
[1] "Isolation and characterization of cDNA clones for an"
[2] "inhibitor protein of cAMP-dependent protein kinase." 

[[2]]
[1] "Inhibition of protein kinase-A by overexpression of"
[2] "the cloned human protein kinase inhibitor."         

[[3]]
[1] "Structure of a peptide inhibitor bound to the"
[2] "catalytic subunit of cyclic adenosine"        
[3] "monophosphate-dependent protein kinase."      

[[4]]
[1] "Tyrosine kinase catalyzed phosphorylation and"
[2] "inactivation of the inhibitor protein of the" 
[3] "cAMP-dependent protein kinase."               

[[5]]
[1] "Identification of an inhibitory region of the"      
[2] "heat-stable protein inhibitor of the cAMP-dependent"
[3] "protein kinase."                                    

[[6]]
[1] "Primary-structure requirements for inhibition by the"
[2] "heat-stable inhibitor of the cAMP-dependent protein" 
[3] "kinase."                                             

[[7]]
[1] "The expression and intracellular distribution of the"
[2] "heat-stable protein kinase inhibitor is cell cycle"  
[3] "regulated."                                          

[[8]]
[1] "Glutamic acid 203 of the cAMP-dependent protein" 
[2] "kinase catalytic subunit participates in the"    
[3] "inhibition by two isoforms of the protein kinase"
[4] "inhibitor."                                      

[[9]]
[1] "Evidence for the importance of hydrophobic residues"
[2] "in the interactions between the cAMP-dependent"     
[3] "protein kinase catalytic subunit and the protein"   
[4] "kinase inhibitors."                                 

[[10]]
[1] "HIV Rev uses a conserved cellular protein export"    
[2] "pathway for the nucleocytoplasmic transport of viral"
[3] "RNAs."                                               

[[11]]
[1] "PrKX is a novel catalytic subunit of the"      
[2] "cAMP-dependent protein kinase regulated by the"
[3] "regulatory subunit type I."                    

[[12]]
[1] "Protein inhibitor of neuronal nitric oxide synthase"
[2] "interacts with protein kinase A inhibitors."        

[[13]]
[1] "Generation and initial analysis of more than 15,000"
[2] "full-length human and mouse cDNA sequences."        

[[14]]
[1] "Mutants of protein kinase A that mimic the" 
[2] "ATP-binding site of protein kinase B (AKT)."

[[15]]
[1] "Development of human protein reference database as an"
[2] "initial platform for approaching systems biology in"  
[3] "humans."                                              

[[16]]
[1] "Human protein reference database as a discovery"
[2] "resource for proteomics."                       

[[17]]
[1] "Human protein reference database--2006 update."

▼R のなかで PubMed アブストラクトキーワード検索をおこなう。ここではキーワードが Protein もしくは protein です。

pro.res <- sapply(absts, function(x) pm.abstGrep("[Pp]rotein", x))
pro.res[[2]]

△実行結果

 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[10]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE

PubMed アブストラクトを閲覧するために HTML レポートを作成する。

pmAbst2HTML(absts[[2]], filename = "pm.html")

△実行結果

pm.html

f:id:nakao_mitsuteru:20070616231151p:image


レポートの生成

▼genenames 遺伝子リストに対応する LocusLink ID遺伝子名(Symbol)を取り出します。

ll <- getLL(genenames, "hgu95av2")
sym <- getSYMBOL(genenames, "hgu95av2")

▼つぎに、LocusLink ID (gene ID) を付与した HTML 形式のレポートを作成します。HTML 形式のレポートウェブブラウザで閲覧できるので便利です。

tab <- data.frame(sym, tab[, -1])
htmlpage(ll, filename = "report.html", 
                title = "HTML report",
           othernames = tab,
           table.head = c("Locus ID", colnames(tab)),
         table.center = TRUE)

△実行結果

report.html

f:id:nakao_mitsuteru:20070616231150p:image

▼report.html は tab の内容にリンクをつけたものです。

tab

△実行結果

            sym  logFC AveExpr      t  P.Value adj.P.Val     B
1708_at  MAPK10 -7.061    7.95 -73.53 7.82e-17  9.87e-13  8.65
36202_at   PKIA  0.853    9.37   9.98 4.94e-07  3.12e-03  4.59
36311_at  PDE1A  0.832    8.56   8.36 3.02e-06  1.27e-02  3.57
33264_at ENOSF1  0.712    4.92   7.43 9.67e-06  2.71e-02  2.84
32660_at   LBA1  0.655    8.68   7.36 1.07e-05  2.71e-02  2.77
38734_at    PLN  0.747    6.26   7.19 1.35e-05  2.83e-02  2.62
1024_at  CYP1A1  0.843    9.70   6.73 2.50e-05  4.40e-02  2.20
36085_at  GNAT1  0.645   12.19   6.65 2.79e-05  4.40e-02  2.12
33818_at    VCP  0.532   12.29   6.45 3.70e-05  5.19e-02  1.92
39058_at    ABR  0.609    7.53   6.28 4.77e-05  5.69e-02  1.74

KEGGやGOへのリンクを付与するには annaffy パッケージを使用します。

library("KEGG")
library("GO")
library("annaffy")
atab <- aafTableAnn(genenames, "hgu95av2", aaf.handler())
saveHTML(atab, file = "report2.html")

△実行結果

report2.html

f:id:nakao_mitsuteru:20070616231149p:image

report2.html は atab の内容をテーブル状に表示しています。atab は annaffy パッケージの aafTable クラスオブジェクトです。

▼カラム名を表示

colnames(atab['1708_at',])

△実行結果

 [1] "Probe"               "Symbol"             
 [3] "Description"         "Function"           
 [5] "Chromosome"          "Chromosome Location"
 [7] "GenBank"             "LocusLink"          
 [9] "Cytoband"            "UniGene"            
[11] "PubMed"              "Gene Ontology"      
[13] "Pathway" 

▼1708_at プローブ遺伝子シンボルを表示

atab['1708_at',]$Symbol

△実行結果

[[1]]
[1] "MAPK10"
attr(,"class")
[1] "aafSymbol"

attr(,"class")
[1] "aafList"

まとめ

CELファイルから、有意に発現変動する遺伝子を列挙する方法をみました。

発現変動のさまざまな基準が利用できます。

どの方法が良いかはデータ依存しています。

データパッケージがあると、さまざまな解析が簡単におこなえます。

PubMedアブストラクトタイトルはAffymetrix IDには簡単に対応づきます。

LocusLink や GO や KEGG などのメタデータを簡単に付与できます。

行列HTML のテーブルページにするメソッドがあります。




もくじにもどる

unjgoigpbcunjgoigpbc2013/07/27 00:33mxuftcjpdpoevdups, <a href="http://www.qfdljtypod.com/">lmcioritxd</a> , [url=http://www.gzzikngyus.com/]bdfczpsmri[/url], http://www.zftisrmfht.com/ lmcioritxd

nwnmrjzdklnwnmrjzdkl2013/07/30 09:34chgdqcjpdpoevdups, <a href="http://www.dtpjtpmcwz.com/">nwnjygkzzb</a> , [url=http://www.kdfgrbwcvr.com/]bypzwjsiji[/url], http://www.igffaiytve.com/ nwnjygkzzb

cdmzgyfdmrcdmzgyfdmr2013/08/27 20:48jjpbfcjpdpoevdups, <a href="http://www.sykbutykif.com/">ukuuydefwr</a> , [url=http://www.uhpsolyxko.com/]fzsrybsnxm[/url], http://www.ltbqizswiw.com/ ukuuydefwr

rzneokixsarzneokixsa2013/12/19 06:07ojooxcjpdpoevdups, <a href="http://www.mbhcuaxduh.com/">kcoetjdwiz</a>

lfenndqlkxlfenndqlkx2013/12/21 07:18qqkzkcjpdpoevdups, http://www.ezysqwcddb.com/ ahraopcipd

rcnnhmqvoercnnhmqvoe2013/12/23 12:07bxpfacjpdpoevdups, <a href="http://www.szrshwiogp.com/">cxeugadcxs</a> , [url=http://www.zzekrrgrud.com/]dxqnzipezs[/url], http://www.dydxvalriy.com/ cxeugadcxs

gbzpsabhzdgbzpsabhzd2014/01/12 06:32svcyrcjpdpoevdups, <a href="http://www.fmvxfxjkfd.com/">tnaexwiklw</a> , [url=http://www.gtxtbpmhhq.com/]rvsokjanyr[/url], http://www.hymjvzsvzt.com/ tnaexwiklw

JasonCoofeJasonCoofe2017/01/25 04:28печать журналов http://wkrolik.com.ua/products/otkrytki

GuestAcirmGuestAcirm2017/09/06 10:27guest test post
<a href="http://googlee.te/">bbcode</a>
<a href="http://googlee.te/">html</a>
http://googlee.te/ simple

 |