给定一些基因，从中找出联系，并能很好的显示给人们看，是一个生物信息学中非常常见的问题。其中pathway分析与gene ontology(GO)分析是两个很重要方面。基于种种原因，人们对GO分析往往比较注重，而随之而来的分析工具也多。而对PATH的分析，却因数据库的原因，手段寥寥。为此，我收集了一些关于pathway数据库方面的资料，总结如下：
1）代谢数据库有哪些？
这个问题很难全面[......]
Tags: bioconductor, 生物信息学, 编程
我们在进行GO分析时，时常需要得到某一特定的GO term下所有涉及到的基因，但是R当中并没有明确的工具。下面的脚本片断可以解决这个问题。
> library("goProfiles") > GOTermsList <- as.GOTerms.list(as.character(genelist.entrez.id),"Entrez","org.Hs.eg.db",onto="BP") > ancestorsList<-getAncestorsLst(GOTermsList,onto="BP") > yapply <- function(X,FUN, ...) { + index <- seq(length.out=length(X)) + namesX <- names(X) + if(is.null(namesX)) + namesX <- rep(NA,length(X)) + + FUN <- match.fun(FUN) + fnames <- names(formals(FUN)) + if( ! "INDEX" %in% fnames ){ + formals(FUN) <- append( formals(FUN), alist(INDEX=) ) + } + if( ! "NAMES" %in% fnames ){ + formals(FUN) <- append( formals(FUN), alist(NAMES=) ) + } + mapply(FUN,X,INDEX=index, NAMES=namesX,MoreArgs=list(...)) + } > ancestors<-do.call(rbind,yapply(ancestorsList,function(.ele){cbind(rep(NAMES,length(.ele)),.ele)})) > colnames(ancestors)<-c("entrez_id","go_term") > library("ChIPpeakAnno") #使用用其中addGeneIDs可以很方便地在种类ID之间进行转换 > anno<-addGeneIDs(as.character(genelist.entrez.id),"org.Hs.eg.db",c("symbol","genename"),"entrez_id") > go.anno<-merge(anno,ancestors,by="entrez_id") > c.cell.death<-lapply(cell.death.GO.terms,function(.ele,x){as.character(x[x$go_term==.ele,"symbol"])},go.anno)
今天又得到高人指点，有了新的更简单的脚本
> library(org.Hs.eg.db) > help("org.Hs.egGO2ALLEGS") > e<-mget(cell.death.GO.terms,org.Hs.egGO2ALLEGS) > t<-lapply(e,function(.ele,x){intersect(as.character(x),.ele)},genelist.entrez.id)
Tags: bioconductor, 生物信息学, 编程
文氏图是一种非常常用的图示手段，主要用于显示组与组之间重叠的程度。
R当中可以画文氏图的包有好几个，使用起来各有特点。最原始的工具，来自于2004年的《Venn Diagrams in R》–Duncan J.Murdoch, Journal of Statistical Software. 但是，现在已经无法找到venn包了。
之后使用得较为广泛的有两个工具，一个是LIMMA内嵌的vennDiagram, 另一个是gplots当中的venn。前者最多可以对3组数据画文氏图，后者可以最多对5组数据画文氏图。首先来介绍使用LIMMA的vennDiagram.
> source("http://bioconductor.org/biocLite.R") > biocLite("limma") > library(limma) > hsb2< -read.table("http://www.ats.ucla.edu/stat/R/notes/hsb2.csv", sep=',', header=T) > attach(hsb2) > hw< -(write>=60) > hm< -(math >=60) > hr< -(read >=60) > c3< -cbind(hw, hm, hr) > a < - vennCounts(c3); a hw hm hr Counts [1,] 0 0 0 113 [2,] 0 0 1 18 [3,] 0 1 0 8 [4,] 0 1 1 8 [5,] 1 0 0 12 [6,] 1 0 1 8 [7,] 1 1 0 11 [8,] 1 1 1 22 attr(,"class") [1] "VennCounts" > vennDiagram(a, include = "both", names = c("High Writing", "High Math", "High Reading"), cex = 1, counts.col = "red")
LIMMA绘制文氏图
从上面的脚本我们可以知道，LIMMA在绘制文氏图的时候，先是对数据转换成bool类型的矩阵，而后进行分组计数，分组就包括所有可能的组合，得出数字后绘图。
[......]
Tags: bioconductor, 生物信息学, 编程
首先升级到最新的firmware,这一步是为了可以使用ssh登录。如果你是第一次登录的话，密码是空的。
qiuworld$ ssh root@192.168.2.2
root@192.168.2.2's password:
Last login: Thu Jan 1 10:08:45 1970 from 192.168.2.3
Linux (none) 2.6.10_mvl401_AG_NAS_3.1.1.exported #1 Thu Sep 6 11:57:38 JST 2007 armv5tejl GNU/Linux
Welcome to MontaVista(R) Linux(R) Professional Edition 4.0.1 (0502020).
root@(none):~#系统时间非常不正确。
root@(none):~# date --set "Wed Jan 18 23:21:21 JST 2012" Wed Jan 18 23:21:21 JST 2012
安装funplug，funplug的作用是安装编译好的软件。因为NAS一般都不带有编译环境，所以使用编译好[......]
1， 出现Error in readRDS(file) : error reading from the connection错误。
出现的原因，安装过程中死机，导致安装不完整。
解决办法，
在Ｒ当中
> .libPaths() [1] "/usr/local/lib64/R/library"
进入/usr/local/lib64/R/library下删除最新安装的所有library。将BiocInsta[......]
Tags: bioconductor
本文心得自：The Split-Apply-Combine Strategy for Data Analysis, Hadley Wickham, Journal of Statistical Software, April 2011, V.40.
引子：
我们常常会遇到这样的问题，数据量很大，并不需要依顺序来依次处理。合理分块处理，并最终整合起来是一个不错的选择。这也就是所谓的Split-Apply-Combine Strategy策略。这在速度上会有比做一个loop有优势，因为它可以并行处理数据。
什么时候我们需要使用到化整为零的策略呢？有以下三种情况：
- 数据需要分组处理
- 数据需要按照每行或者每列来处理
- 数据需要分级处理，和分组很类似，但是分级时需要考虑分级之间的关系。
化整为零策略有点类似于由Google推广的map-reduce策略。当然map-reduce策略的基础是网格，而这里的Split-Apply-Combine的基础完全可以是单机，甚至不支持并行处理的单机都可以。
然而，化整为零并不是一个很直观的编程过程。最直观的过程是使用Loop循环。这里使用一个例子来讲解一下如何实现化整为零策略。在plyr包中有数据ozone,它是一个三维矩阵(24X24X72)，其中最后一维72是指的6年12个月每个月的结果。也就是ozone是一个包括了连续72个月24Ｘ24的三维矩阵数据。三维分别是lat,long,time。我们需要由对时间robust linear model之后的残基residuals。
> library(plyr) # need for dataset ozone > library(MASS) # need for function rlm > month < - ordered(rep(1:12, length=72)) #set time sequence > #try one set > one < - ozone[1,1,] > model < - rlm(one ~ month - 1); model Call: rlm(formula = one ~ month - 1) Converged in 9 iterations Coefficients: month1 month2 month3 month4 month5 month6 month7 month8 month9 month10 month11 month12 264.3964 259.2036 255.0000 252.0052 258.5089 265.3387 274.0000 276.6724 277.0000 285.0000 283.6036 273.1964 Degrees of freedom: 72 total; 60 residual Scale estimate: 4.45 > deseas < - resid(model)
UTF8_E[......]
Tags: bioconductor, 教程, 编程
近期评论