14

R语言画网络图三部曲之sna

如果只是画网络图,那么只需要把所有的点,按照算好的坐标画出来,然后把所有的连线也画出即可!

其中算法就是,点的坐标该如何确定?Two of the most prominente algorithms (Fruchterman & Reingold’s force-directed placement algorithm and Kamada-Kawai’s)
有一个networkD3的包可以直接画图,但是跳过了确定点的坐标这个步骤,我重新找了一包,可以做到!
作者只是用sna包来得到数据,其实用的是ggplot来画网络图!
需要熟悉network包里面的network对象的具体东西,如何自己构造一个,然后数学sna包如何计算layout即可
解读这个包,也可以自己画网络图,代码如下:
plot(plotcord)
text(x=plotcord$X1+0.2,y=plotcord$X2,labels = LETTERS[1:10])
for (i in 1:10){
  for (j in 1:10){
    if(tmp[i,j]) lines(plotcord[c language="(i,j),1"][/c],plotcord[c language="(i,j),2"][/c])
  }
}
当然,我们还没有涉及到算法,就是如何生成plotcord这个坐标矩阵的!
大家看下面这个示意图就知道网络图是怎么样画出来的了,首先我们有一些点,它们之间有联系,都存储在networData这个数据里面,是10个点,共9个连接,然后我用reshape包把它转换成连接矩阵,理论上10个点的两两相互作用应该有100条线,但是我们的数据清楚的说明只有9条,所以只有9个1,其余的0代表点之间没有关系。接下来我们用sna这个包对这个连接矩阵生成了这10个点的坐标(这个是重点),最后很简单了,把点和线画出来即可!

1

另外一个例子:
net=network(150, directed=FALSE, density=0.03)
m <- as.matrix.network.adjacency(net) # get sociomatrix  
# get coordinates from Fruchterman and Reingold's force-directed placement algorithm.
plotcord <- data.frame(gplot.layout.fruchtermanreingold(m, NULL)) 
# or get it them from Kamada-Kawai's algorithm: 
# plotcord <- data.frame(gplot.layout.kamadakawai(m, NULL)) 
colnames(plotcord) = c("X1","X2")  ###所有点的坐标,共150个点
edglist <- as.matrix.network.edgelist(net) ##所有点之间的关系-edge ##共335条线
edges <- data.frame(plotcord[edglist[,1],], plotcord[edglist[,2],]) ##两点之间的连线的具体坐标,335条线的起始终止点点坐标
原始代码如下:
library(network)
library(ggplot2)
library(sna)
library(ergm)
clipboard
大家可以试用这个代码,因为它用的ggplot,肯定比我那个简单R作图要好看多了
参考算法文献:
 

 

14

R语言画网络图三部曲之networkD3

首先,我们需要了解一下基础知识:
网络图是为了展示数据与数据之间的联系,在生物信息学领域,一般是基因直接的相互作用关系,或者通路之间的联系!
通俗点说,就是我有一些点,它们之间并不是两两相互联系,而是有部分是有连接的,那么我应该如何把这些点画在图片上面呢?因为这些都并没有X,Y坐标,只有连接关系,所以我们要根据一个理论来给它们分配坐标,这样就可以把它们画出来了,然后也可以对这些点进行连线,连完线后,就是网络图啦!!!
而给它们分配坐标的理论有点复杂,大概类似于物理里面的万有引力和洛仑磁力相结合来给它们分配一个位置,使得总体的能量最小,就是最稳定状态!而通常这个状态是逼近,而不是精确,所以我们其实对同样的数据可以画出无数个网络图,只需使得网络图合理即可!
看这个ppt:http://statmath.wu.ac.at/research/friday/resources_WS0708_SS08/igraph.pdf

基本上就都理解了
画网络图真的好复杂呀!!!
大部分人都是用D3JS来画图:http://bl.ocks.org/mbostock/2706022
一步步讲受力分析图是如何画的:https://github.com/mbostock/d3/wiki/Force-Layout

接下来, 我们直接看看R里面是如何画网络图的,我们首推一个包:networkD3/

这个包非常好用,只需要做好data,然后用它提供的几个函数即可!
重要的是熟悉输入数据是什么,还可以结合shinny包来展示数据,非常好用!!!
具体怎么安装这个R包,我就不讲了,它的github主页上面其实也有说明书,很容易看懂的
最简单的用法,就是构造一个联结矩阵,把所有的连接关系都存储起来,然后直接用一个函数就可以画图啦!
# Plot
library(networkD3)
simpleNetwork(networkData,fontSize=25)
可以看出网络图,就是把所有的点,按照算好的坐标画出来,然后把所有的连线也画出即可!
其中算法就是,点的坐标该如何确定?Two of the most prominente algorithms (Fruchterman & Reingold’s force-directed placement algorithm and Kamada-Kawai’s)

clipboard

但是这个软件有一个弊端就是,生成图片之后,并没有给出这些点的坐标,当然,这种坐标基本上也很少有人需要,可视化就够了!
如果这些点还分组了,而且连接还有权重,那么网络图就复杂一点,比如下面这个:
就不仅仅是需要提供所有的link的信息,link还多了一列,是value,而且还需提供所有的node信息,node多了一列是分组。
clipboard
当然,还有好几个图,大家可以自己慢慢用,自己体会!
sankeyNetwork
sankeyNetwork
radialNetwork
diagonalNetwork
dendroNetwork
也可以把生成的网络图直接保存成网页,动态显示!
十二 25

使用进度条和并行计算来模拟计算赌博输赢概率

http://www.zhihu.com/question/19605599

以前看到过一个赌徒「必胜方法」
有一个广泛流传于赌徒中的“必胜方法”,是关于赌徒悖论极好的说明:一个赌场设有“猜大小”的赌局,玩家下注后猜“大”或者猜“小”,如果输了,则失去赌注,如果赢的话,则获得本金以及0.9倍的利润。 “必胜法”是这么玩的:
(1)押100猜“大”;
(2)如果赢的话,返回(1)继续;
(3)如果输的话则将赌注翻倍后继续猜“大”,因为不可能连续出现“大”,总会有获胜的时候,而且由于赌注一直是翻倍,只要赢一次,就会把所有输掉的钱赢回;
(4)只要赢了,就继续返回(1)。

一个朋友说写了程序来验证正确与否,我感觉蛮有趣的,想想自己也学了不少,也应该实践一下啦!
写出程序不难,就一个递归而已,不过是应该用并行计算,不然算的太慢了
另外一个重点是,如何随机???
[R]
judge=function(l=left,c=count,n=number,w=win){
c=c+1;
if(rnorm(1)>0){
l=l+2^n
n=1 # win
w=w+1
}else{
l=l-2^n
n=n+1 #lose
}
return(c(l,c,n,w))
}
my_fun=function(x){
tmp=judge(0,1,1,0)
for(i in 1:x){
tmp=judge(tmp[1],tmp[2],tmp[3],tmp[4])
print(tmp)
}
return(tmp[1])
}
my_fun(10)
[/R]
第一个函数是根据现在的状态,接受你的本金,第几次赌了,以及第几次连续下注,最后,你总共赢了多少次,这四个参数,然后随机给一个大小概率,这样你的本金会变化,赌的次数会加1
第二个函数是传入我们需要赌多少次,看结果:
> my_fun(10)
[1] -2  2  2  0
[1] 2 3 1 1
[1] 0 4 2 1
[1] -4  5  3  1
[1] 4 6 1 2
[1] 6 7 1 3
[1] 8 8 1 4
[1] 6 9 2 4
[1]  2 10  3  4
[1] 10 11  1  5
[1]  8 12  2  5
[1] 8
如果我们赌11次,可以看到,我最后会剩余8块钱,每次输赢的情况都反应在里面了,可以自己模拟多次看看!
因为我只赌了11次,所以很快,如果我赌1000次,而且还想检验一下10000次模拟结果,就会比较慢了!
我首先使用进度条模拟一下结果,代码如下:
2
还是比较慢的##Time difference of 1.861802 mins
我用了apply,好像时间是节省了一些,不过聊胜于无!
> library(pbapply)
> start=Sys.time()
> results=pbsapply(1:10000,function(i){my_fun(1000)})
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
> Sys.time()-start  ##9.909524 secs
Time difference of 1.301539 mins
那接下来我们分析一下模拟的结果吧
> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -64580     974     998     985    1020    1110
可以看到,我们平均下来是会赢钱的,而且赢面很大,赢的次数非常多!!!
但是,我们看下面的图就知道,一个很诡异的结果!
而且,我这里用的模拟胜负,并不是很好
3
这其实还只是小批量模拟,如果要模拟百亿次,首先我的笔记本肯定不行,cpu太破了,如果用服务器,就需要用并行计算啦!
##下面是用多核并行计算的代码,大家有兴趣可以自己去玩一下
library(parallel)
cl.cores <- detectCores() 
cl <- makeCluster(cl.cores)
clusterExport(cl, "judge") 
start=Sys.time()
results=parSapply( 
  cl=cl, 
  1:10000,
  my_fun(1000)
)
Sys.time()-start  ##4.260994 secs
stopCluster(cl)


十二 15

用超几何分布检验做富集分析

我们可以直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

image001

我们挑选logFC的绝对值大于0.5,并且P-value小雨0.05的基因作为差异基因,并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下:

image002

但是这样我们就忽略了其中的原理,我们不知道这些数据是如何算出来的,只是由别人写好的包得到了结果罢了。

事实上,这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理,我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单,球分成黑白两色,数量已知,那么你随机抽有限个球,应该抽多少白球的问题!
公式就是 exp_count=n*M/N
然后你实际上抽了多少白球,就可以计算一个概率值!
换算成通路的富集概念就是,总共有多少基因,你的通路有多少基因,你的通路被抽中了多少基因(在差异基因里面属于你的通路的基因),这样的数据就足够算出上面表格里面所有的数据啦!
tmp=toTable(org.Hs.egPATH)
GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#phyper(k-1,M, N-M, n, lower.tail=F)
#n*M/N
diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))
n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球
N=length(GeneID2Path) #5870  ##这里算出你总共有多少个球(这里是错的,有多少个球取决于背景基因!一般是两万个)
options(digits = 4)
for (i in names(Path2GeneID)){
 M=length(Path2GeneID[[i]])  ##这个算出你的所有的球里面,白球有多少个
 exp_count=n*M/N  ###这里算出你抽取的球里面应该多多少个是白色
 k=0         ##这个k是你实际上抽取了多少个白球
 for (j in diff_gene_has_path){
 if (i %in% GeneID2Path[[j]]) k=k+1
 }
 OddsRatio=k/exp_count
 p=phyper(k-1,M, N-M, n, lower.tail=F)  ##根据你实际上抽取的白球个数,就能算出富集概率啦!
 print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))
}
随便检查一下,就知道结果是一模一样的!

 

30

R语言也可以用hash啦

本来想着是如何习惯R里面的list对象,来实现hash的功能,无意中发现了居然还有这个包

https://cran.r-project.org/web/packages/hash/hash.pdf

我简单看了看文档,的确跟perl里面的hash功能类似,非常好用。

创建一个hash很简单

h <- hash( keys=letters, values=1:26 )

h <- hash( letters, 1:26 )

类似于下面的列表

L=setNames(as.list(LETTERS),1:26)

如果是列表要提取keys需要用names函数,如果需要提取values,需要用unlist函数,而用了hash之后,这两个函数就可以独自运行啦

还有很多其它函数,其实在list中也可以实现,主要是看你对哪种语法更熟悉,感觉自己现在编程能力差不多算是小有所成了,看什么都一样了,要是以前学R的时候看到了这个hash包,我肯定会很兴奋的!

clear signature(x = "hash"): Remove all key-value pairs from hash
del signature(x = "ANY", hash = "hash"): Remove specified key-value pairs from hash
has.key signature(key = "ANY", hash = "hash"): Test for existence of key
is.empty signature(x = "hash"): Test if no key-values are assigned
length signature(x = "hash"): Return number of key-value pairs from the hash
keys signature(hash = "hash"): Retrieve keys from hash
values signature(x = "hash"): Retrieve values from hash
copy signature(x = "hash"): Make a copy of a hash using a new environment.
format signature(x = "hash"): Internal function for displaying hash

 

30

R里面list对象和AnnDbBimap对象区别

AnnDbBimap对象是R里面的bioconductor系列包的基础对象,在探针数据里面会包装成ProbeAnnDbMap,跟go通路相关又是GOTermsAnnDbBimap对象。

但是它都是AnnDbBimap对象衍生过来的

主要存在于芯片系列的包和org系列的包,其实AnnDbBimap对象就是list对象的包装,比如下面这些例子:

ls("package:hgu133plus2.db")

 [1] "hgu133plus2"              "hgu133plus2.db"           "hgu133plus2_dbconn"
 [4] "hgu133plus2_dbfile"       "hgu133plus2_dbInfo"       "hgu133plus2_dbschema"
 [7] "hgu133plus2ACCNUM"        "hgu133plus2ALIAS2PROBE"   "hgu133plus2CHR"
[10] "hgu133plus2CHRLENGTHS"    "hgu133plus2CHRLOC"        "hgu133plus2CHRLOCEND"
[13] "hgu133plus2ENSEMBL"       "hgu133plus2ENSEMBL2PROBE" "hgu133plus2ENTREZID"
[16] "hgu133plus2ENZYME"        "hgu133plus2ENZYME2PROBE"  "hgu133plus2GENENAME"
[19] "hgu133plus2GO"            "hgu133plus2GO2ALLPROBES"  "hgu133plus2GO2PROBE"
[22] "hgu133plus2MAP"           "hgu133plus2MAPCOUNTS"     "hgu133plus2OMIM"
[25] "hgu133plus2ORGANISM"      "hgu133plus2ORGPKG"        "hgu133plus2PATH"
[28] "hgu133plus2PATH2PROBE"    "hgu133plus2PFAM"          "hgu133plus2PMID"
[31] "hgu133plus2PMID2PROBE"    "hgu133plus2PROSITE"       "hgu133plus2REFSEQ"
[34] "hgu133plus2SYMBOL"        "hgu133plus2UNIGENE"       "hgu133plus2UNIPROT"

那么我们可以随便挑选包中的一个数据分析一下

x <- hgu133plus2SYMBOL

xlist=as.list(x)

我们查看X对象,可知,它是object of class "ProbeAnnDbMap",而这个对象,就是常见的list对象,被包装了一下, 只有我们明白了它和list对象到底有什么区别,才算是真正搞懂了这一系列包

但是这个ProbeAnnDbMap对象,在GO等包里面会被包装成更复杂的对象-GOTermsAnnDbBimap,但是对他们的理解都大同小异。

我们先回顾一下在R语言里面的list的基础知识:

R 的 列表(list)是一个以对象的有序集合构成的对象。 列表中包含的对象又称为它的分量(components)。

分量可以是不同的模式或者类型, 如一个列表可以包括数值向量,逻辑向量,矩阵,复向量, 字符数组,函数等等

如果你会perl的话,可以把它理解成hash。

分量常常会被编号的(numbered),并且可以利用它来 访问分量

列表的分量可以被命名,这种情况下 可以通过名字访问。

特别要注意一下 Lst[[1]] 和 Lst[1] 的差别。 [[...]] 是用来选择单个 元素的操作符,而 [...] 是一个更为一般 的下标操作符。

因此前者得到的是列表 Lst 中的第一个对象, 并且含有分量名字的命名列表(named list)中分量的名字会被排除在外的。

后者得到的则是列表 Lst 中仅仅由第一个元素 构成的子列表。如果是命名列表, 分量名字会传给子列表的。

那么接下来我们就看看x和xlist的区别。

它们里面的数据都是affymetrix公司出品的人类的hgu133plus2芯片的探针ID与基因symbol的对应关系

如果我想拿到所有的探针ID,那么对于AnnDbBimap对象,需要用mappedkeys(x),对于普通的list对象,需要用names(xlist).

对于普通的list对象,如果我想看前几个元素,直接head就可以了,但是对于AnnDbBimap对象,数据被封装好了,需要先as.list,然后才能head

mapped_probes <- mappedkeys(x)

PID2=head(mapped_probes)

那么,如果我们想根据以下探针ID来查看它们在这些数据里面被对应着哪些基因symbol 呢?

> PID2 #这是一串探针ID,后面的操作都需要用的

[1] "1053_at"   "117_at"    "121_at"    "1255_g_at" "1316_at"   "1320_at"

如果是对于AnnDbBimap对象,我们可以用mget函数来操作,取多个探针的对应基因symbol

accnum <- mget(PID2, env=hgu133plus2ACCNUM);

gname <- mget(PID2, env=hgu133plus2GENENAME)

gsymbol <- mget(PID2, env=hgu133plus2SYMBOL)

mget函数返回的就是普通的list函数了,可以直接查看了。

如果是对于普通的list对象,我们想取多个探针的对应基因symbol也是非常简单的,xlist[PID2]即可。

那么我们不禁有问了,既然它们两个功能完全一样,何苦把它包装成一个对象了,我直接操作list对象不就好了,学那么多规矩干嘛?

所以,重点就来了

>  length(mapped_probes)

[1] 42125

> length(names(xlist))

[1] 54675

看懂了吗?

 

但是,事实上用处也不大,你觉得下面这两个有区别吗?

SYMBOL <- AnnotationDbi::as.list(hgu133plus2SYMBOL)

SYMBOL <- SYMBOL[!is.na(SYMBOL)];

 

x <- hgu133plus2SYMBOL

mapped_probes <- mappedkeys(x)

SYMBOL <- AnnotationDbi::as.list(x[mapped_probes])

 

PS,在R里面创建一个list对象是非常简单的

The setNames() built-in function makes it easy to create a hash from given key and value lists.(Thanks to Nick K for the better suggestion.)

Usage: hh <- setNames(as.list(values), keys)

Example:

players <- c("bob", "tom", "tim", "tony", "tiny", "hubert", "herbert")
rankings <- c(0.2027, 0.2187, 0.0378, 0.3334, 0.0161, 0.0555, 0.1357)
league <- setNames(as.list(rankings), players)

29

gene的各种ID转换终结者-bioconductor系列包

经常会有人问这样的问题I have list of 10,000 Entrez IDs and i want to convert the multiple Entrez IDs into the respective gene names. Could someone suggest me the way to do this?等等类似的基因转换,能做的基因转换的方法非常多,以前不懂编程的时候,都是用各种网站,而最常用的就是ensembl的biomart了,它支持的ID非常多,高达几百种,好多ID我到现在都不知道是什么意思。

现在学会编程了,我比较喜欢的是R的一些包,是bioconductor系列,一般来说,其中有biomart,org.Hs.eg.db,annotate,等等。关于biomart我就不再讲了,我前面的博客至少有七八篇都提到了它。本次我们讲讲简单的, 我就以把gene entrez ID转换为gene symbol 为例子把。

当然,首先要安装这些包,并且加载。

if("org.Hs.eg.db" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("org.Hs.eg.db")}
suppressMessages(library(org.Hs.eg.db))  我比较喜欢这样加载包

library(annotate) #一般都是这样加载包

如果是用org.Hs.eg.db包,首先你只需要读取你的待转换ID文件,构造成一个向量,tmp,然后只需要symbols <- org.Hs.egSYMBOL[as.character(tmp)]就可以得到结果了,返回的symbols是一个对象,需要用toTable这个函数变成数据框。但是这样转换容易出一些问题,比如如果你的输入数据tmp,里面含有一些无法转换的gene entrez ID,就会报错。

而且它支持的ID转换很有限,具体看看它的说明书即可:https://www.bioconductor.org/packages/release/data/annotation/manuals/org.Hs.eg.db/man/org.Hs.eg.db.pdf

org.Hs.eg.db
org.Hs.eg_dbconn
org.Hs.egACCNUM
org.Hs.egALIAS2EG
org.Hs.egCHR
org.Hs.egCHRLENGTHS
org.Hs.egCHRLOC
org.Hs.egENSEMBL
org.Hs.egENSEMBLPROT
org.Hs.egENSEMBLTRANS
org.Hs.egENZYME
org.Hs.egGENENAME
org.Hs.egGO
org.Hs.egMAP
org.Hs.egMAPCOUNTS
org.Hs.egOMIM
org.Hs.egORGANISM
org.Hs.egPATH
org.Hs.egPMID
org.Hs.egREFSEQ
org.Hs.egSYMBOL
org.Hs.egUCSCKG
org.Hs.egUNIGENE
org.Hs.egUNIPROT

如果是用annotate包,首先你还是需要读取你的待转换ID文件,构造成一个向量,tmp,然后用getSYMBOL(as.character(tmp), data='org.Hs.eg')这样直接就返回的还是以向量,只是在原来向量的基础上面加上了names属性。说明书:http://www.bioconductor.org/packages/3.3/bioc/manuals/annotate/man/annotate.pdf

然后你可以把转换好的向量写出去,如下:

1 A1BG
2 A2M
3 A2MP1
9 NAT1
10 NAT2
12 SERPINA3
13 AADAC
14 AAMP
15 AANAT
16 AARS

PS:如果是芯片数据,需要把探针的ID转换成gene,那么一般还需要加载特定芯片的数据包才行:

platformDB <- paste(eset.mas5@annotation, ".db", sep="") #这里需要确定你用的是什么芯片
cat("the annotation is ",platformDB,"\n")
if(platformDB %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");tmp=try(biocLite(platformDB))}
library(platformDB, character.only=TRUE)
probeset <- featureNames(eset.mas5)
rowMeans <- rowMeans(exprSet)

library(annotate) # lookUp函数是属于annotate这个包的
EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

29

使用GEOmetadb包来获取对应GEO数据的实验信息

理论上我前面提到的GEOquery包就可以根据一个GSE索引号来获取NCBI提供的所有关于这个GSE索引号的数据了,包括metadata,表达矩阵,soft文件,还有raw data

但是很多时候,那个metadata并不是很整齐,而且一个个下载太麻烦了,所以就需要用R的bioconductor的另一个神奇的包了GEOmetadb

它的示例:http://bioconductor.org/packages/devel/bioc/vignettes/GEOmetadb/inst/doc/GEOmetadb.R

里面还是很多数据库基础知识的
代码托管在github,它的示例代码是这样连接数据库的:
library(GEOmetadb)
if(!file.exists('GEOmetadb.sqlite')) getSQLiteFile()
file.info('GEOmetadb.sqlite')
con <- dbConnect(SQLite(),'GEOmetadb.sqlite')
dbDisconnect(con)
但是一般不会成功,因为这个包把它的GEOmetadb.sqlite文件放在了国外网盘共享,在国内很难访问,推荐大家想办法下载到本地
tmp2
用这个代码就会成功了,需要自己下载GEOmetadb.sqlite文件然后放在指定目录:/path/GEOmetadb.sqlite 需要自己修改
我们的diabetes.GEO.list文件内容如下:
GSE1009
GSE10785
GSE1133
GSE11975
GSE121
GSE12409
那么会产生的表格文件如下:共有32列数据信息,算是蛮全面的了
tmp
22

R语言中的排序,集合运算,reshape,以及merge总结

一直以来我都是随便看了点R的编程教程,因为我学了一点点C,所以还算有基础,现在基本上简单看看教程就能懂一门语言了,区别只是熟练度而已。R用得比较多,所以还算擅长,但是很多快捷应用的地方,我总是寄希望于到时候再查资料,所以没能用心的记住,这次花了点时间好好整理了一下R里面关于数据操作的重点,我想,以后再碰到类似的数据处理要求,应该很快能解决了把。

首先看看排序:

在R中,和排序相关的函数主要有三个:sort(),rank(),order()。

sort(x)是对向量x进行排序,返回值排序后的数值向量。

rank()是求秩的函数,它的返回值是这个向量中对应元素的“排名”。

order()的返回值是对应“排名”的元素所在向量中的位置。

其中sort(x)等同于x[order(x)]

下面以一小段R代码来举例说明:

> x<-c(97,93,85,74,32,100,99,67)

> sort(x)

[1]  32  67  74  85  93  97  99 100

> order(x)

[1] 5 8 4 3 2 1 7 6

> rank(x)

[1] 6 5 4 3 1 8 7 2

> x[order(x)]

[1]  32  67  74  85  93  97  99 100

其中比较有用的order,它可以用来给数据框进行排序

dat[order(dat[,1]),] 以该数据框的第一列进行排序

dat[order(dat[,1],dat[,2]),] 以该数据框的第一列为主要次序,第二列为次要序列进行排序

 

然后我们看看集合运算:

在R里面除了简单的对两个向量求交集并集补集之外,比较重要的就是match和 %in% 了,需要重点讲讲。

#首先对集合A,B,C赋值

> A<-1:10

> B<-seq(5,15,2)

> C<-1:5

> #求A和B的并集

> union(A,B)

[1]  1  2  3  4  5  6  7  8  9 10 11 13 15

> #求A和B的交集

> intersect(A,B)

[1] 5 7 9

> #求A-B

> setdiff(A,B)

[1]  1  2  3  4  6  8 10

> #求B-A

> setdiff(B,A)

[1] 11 13 15

> #检验集合A,B是否相同

> setequal(A,B)

[1] FALSE

> #检验元素12是否属于集合C

> is.element(12,C)

[1] FALSE

> #检验集合A是否包含C

> all(C%in%A)

[1] TRUE

> all(C%in%B)

从上面可以看到%in%这个操作符只返回逻辑向量TRUE 或者FALSE,而且返回值应该与%in%这个操作符前面的向量程度相等。也就是说它相当于遍历了C里面的一个个元素,判断它们是否在B中出现过,然后返回是或者否即可。

而match(C,B)的结果就很不一样了,它的返回结果同样与前面的向量等长,但是它并非返回逻辑向量,而是遍历了C里面的一个个元素,判断它们是否在B中出现过,如果出现就返回在B中的索引号,如果没有出现,就返回NA

>B

[1]  5  7  9 11 13 15

>C

[1] 1 2 3 4 5

>match(C,B)

[1] NA NA NA NA  1

>C%in%B

[1] FALSE FALSE FALSE FALSE  TRUE

 接下来我们看看reshape:

这是一个需要安装的包,起得就是R里面最经典的把长型数据变宽,和把宽数据拉长的作用。

其中melt函数是把很宽的数据拉长,它就是需要指定几列数据是保证不被融合的, 其余每一列数据都必须被融合到一列了,融合后的这一列数据每个元素旁边就用列名来标记该数据来自于哪一列。

# example of melt function
library(reshape)
mdata <- melt(mydata, id=c("id","time"))

融合后的数据如下:

id time variable value
1 1 x1 5
1 2 x1 3
2 1 x1 6
2 2 x1 2
1 1 x2 6
1 2 x2 5
2 1 x2 1
2 2 x2 4

可以看到variable列里面有两个不同的元素,说明是把旧数据中的两列给融合了,融合后的一个很长的列就是value

而cast函数的功能就是把刚才融合好的数据给还原。

# cast the melted data
# cast(data, formula, function)
subjmeans <- cast(mdata, id~variable, mean)
timemeans <- cast(mdata, time~variable, mean)

可以看到它们返回的结果是:

subjmeans

id x1 x2
1 4 5.5
2 4 2.5

timemeans

time x1 x2
1 5.5 3.5
2 2.5 4.5

可以看到cast函数比较复杂一点,formula公式右边的变量是需要拆开的variable,这一列变量有多少不重复元素,就新增多少列,左边的变量是行标记了,其余所有数据都会被计算一下再放在合适的位置。

在reshape2这个包里面还进化出来dcast和acast函数,功能大同小异。

 

最后我们来看看merge函数:

这个函数的功能非常强大,类似于SQL语句里面的join系列函数

测试数据如下,它们这两个表的连接是作者名

Capture1

 

如果要实现类似sql里面的inner join 功能,则用代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name")

如果要实现left join功能则用代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all.x=TRUE)

right join功能代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all.y=TRUE)

all join功能代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all=TRUE)

关于单变量匹配的总结就是这些,但对于多变量匹配呢,例如下面两个表,需要对k1,k2两个变量都相等的情况下匹配

x <- data.frame(k1 = c(NA,NA,3,4,5), k2 = c(1,NA,NA,4,5), data = 1:5)

y <- data.frame(k1 = c(NA,2,NA,4,5), k2 = c(NA,NA,3,4,5), data = 1:5)

匹配代码如下merge(x, y, by = c("k1","k2"))  #inner join

另外一个多行匹配的例子如下:

Capture2

我们的测试数据如上,这两个表的连接在于作者名。下面这两个语句等价

merge(authors, books, by=1:2)

merge(authors, books, by.x=c("FirstName", "LastName"),

by.y=c("AuthorFirstName", "AuthorLastName"),

all.x=TRUE)

都可以把两张表关联起来。

当然,在我搜索资料的时候,发现了另外一个解决问题的方法:

A[with(A, paste(C1, C2, sep = "\r")) %in% with(B, paste(C1, C2, sep="\r")), ]

 

 

参考:http://blog.sina.com.cn/s/blog_6caea8bf0100spe9.html

http://blog.sina.com.cn/s/blog_6caea8bf010159dt.html

http://blog.csdn.net/u014801157/article/details/24372441

http://www.statmethods.net/management/reshape.html

https://docs.tibco.com/pub/enterprise-runtime-for-R/1.5.0_may_2013/TERR_1.5.0_LanguageRef/base/merge.html

http://r.789695.n4.nabble.com/Matching-multiple-columns-in-a-data-frame-td890832.html

http://bbs.pinggu.org/thread-3234639-1-1.html

http://www.360doc.com/content/14/0821/13/18951041_403561544.shtml

http://developer.51cto.com/art/201312/423612.htm

22

用R和mysql把基因对应的最大表达量探针挑出来

懂点芯片数据分析的都应该知道,芯片设计的时候,对一个基因设计了多个探针,这样设计是为了更好的捕获某些难以发现的基因,或者重点研究某些基因。
但是对我们的差异分析不方便,所以我们只分析哪些有对应了entrez ID的探针,而且对每个entrez ID,我们只需要挑选它表达量最高的那个探针。
所以就演化为一个编程问题:分组求最大值,多公共列合并!
如果是在R语言里面,那么首先这个table的表示形式如下
> esetDataTable[1:10,c(7,8)]
EGID rowMeans
1000_at 5595 1840.04259751826
1001_at 7075 799.075414422572
1002_f_at 1557 50.4884096416177
1003_s_at 643 142.372008051308
1004_at 643 211.65300963049
1005_at 1843 4281.29318032004
1006_at 4319 38.5784289213085
1007_s_at NA 1489.98158531843
1008_f_at 5610 4013.576753977
1009_at 3094 3070.50648167305
我们首先看看一个R语言函数处理方式吧,这个是比较容易想到的算法,但是用到了循环,非常的不经济,计算量很大。

Capture1
因为里面涉及到了双重循环。我进行了人工计时,这个程序耗费了一分钟十二秒,程序里面的计时器有点问题。

然后我再讲一个精简版的算法
dat=esetDataTable[,c(7,8)]
dat=as.data.frame(apply(dat,2,as.numeric))
dat$prob=rownames(esetDataTable)
首先可以得到需要提取的数据所在行的两个值
match_row=aggregate(dat,by=list(dat[,1]),max)
类似于这句SQL语句:SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID
现在要根据match_row表去原表esetDataTable里面提取我们的探针ID数据。

Capture2

这样就完美解决了问题,我们的合并后的数据的prob列,就是前面那个函数计算了一分多钟的返回的非冗余探针ID向量,relevantProbesets,但是这次只花了不到一秒钟。
tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))
setdiff(as.character(relevantProbesets),as.character(tmp_prob$prob))
length(union(as.character(relevantProbesets),as.character(tmp_prob$prob)))
setdiff(as.character(tmp_prob$prob),as.character(relevantProbesets))
我顺便检查了一下上面那个复杂的R函数跟我这次精简版的结果,发现这次才是对的,上面那个错了。
你们能发现上面那个为什么错了吗?

如果是在mysql里面它的表现形式如下;
mysql> select row_names,EGID,rowMeans from esetDataTable order by EGID limit 10;
+------------+-------+------------------+
| row_names | EGID | rowMeans |
+------------+-------+------------------+
| 38912_at | 10 | 85.8820246154773 |
| 41654_at | 100 | 301.720067595558 |
| 907_at | 100 | 273.100008206028 |
| 2053_at | 1000 | 354.207060199715 |
| 2054_g_at | 1000 | 33.8472900312781 |
| 40781_at | 10000 | 425.007640082848 |
| 35430_at | 10001 | 152.885791914329 |
| 35431_g_at | 10001 | 181.915087187117 |
| 35432_at | 10001 | 184.869017764782 |
| 31366_at | 10002 | 44.9716205901791 |
+------------+-------+------------------+
10 rows in set (0.05 sec)

如果要用SQL来做同样的事情需要下面这个语句,这个就非常简单啦!
select b.*
from test.esetDataTable b
inner join (
SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID
) as c on b.EGID=c.EGID and b.rowMeans=c.val
结果是:8640 rows in set (4.56 sec) 看起来mysql还没有R语言快速,但是这个mysql语句很明显也不是很高效,但是我的mysql水平有限。
稍微解释一下这个mysql语句,其中SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID类似于R里面的match_row=aggregate(dat,by=list(dat[,1]),max),就是把每个entrez ID对应着的最大值提取出来成一个新的表
而inner join on b.EGID=c.EGID and b.rowMeans=c.val 就是我们的tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans")) 根据两列来合并两个数据框

14

用wget批量下载需要认证的网页或者ftp站点里面的pdf文档

这是一个终极命令,比写什么爬虫简单多了
wget -c -r -np -k -L -p  -A.pdf --http-user=CS374-2011 --http-passwd=AlgorithmsInBiology http://ai.stanford.edu/~serafim/CS374_2011/papers/
我这里的例子是斯坦福大学的生物信息学算法课程里面推荐阅读的的所有pdf格式的paper
课程的网址是:http://ai.stanford.edu/~serafim/CS374_2011/
可以看到,这个网站推荐的文献分成8大类,本身这个网站打开就需要登录用户名和密码
--http-user=CS374-2011 --http-passwd=AlgorithmsInBiology
每一篇文献的单独地址是http://ai.stanford.edu/~serafim/CS374_2011/papers/Miscellaneous_topics/Self-assembly_of_DNA/self_healing_and_proofreading.pdf 类似的格式。
我这里简单解释一下这些参数的意思:
-c -r -np -k -L -p  -A.pdf
-c 断点续传
-r 递归下载,下载指定网页某一目录下(包括子目录)的所有文件
-nd 递归下载时不创建一层一层的目录,把所有的文件下载到当前目录(特殊要求会选择这个参数)
-np 递归下载时不搜索上层目录,如wget -c -r www.xxx.org/pub/path/
没有加参数-np,就会同时下载path的上一级目录pub下的其它文件 (所以一定要加上这个参数,不然会下载太多东西的)
-k 将绝对链接转为相对链接,下载整个站点后脱机浏览网页,最好加上这个参数
-L 递归时不进入其它主机,如wget -c -r www.xxx.org/
-p 下载网页所需的所有文件,如图片等
-A 指定要下载的文件样式列表,多个样式用逗号分隔
至于最后的--http-user=CS374-2011 --http-passwd=AlgorithmsInBiology 就是登录该课程网站需要的用户名和密码
29

关于qsub和condor两种在集群上面提交任务的方式比对

毕竟不是计算机科班出身,很多计算机概念我其实并不清楚,很多时候对一个任务的解决只是凭着自己的悟性来模仿接触到的做法,大多数情况下并无可厚非,能完成任务即可,但是碰到多人协作的环境,往往就出了问题。

我一直都知道可以 cat /proc/cpuinfo 来查看cpu数量,用free -g 来查看内存的总量及消耗,用top来查看当前运行的任务情况。
以前都是三到五人的环境公用一个服务器,一般都有几十个cpu和几百G的内存,而且服务器还经常空闲着,所以不会面临资源的问题。
而最近接触的项目是多个小组公用服务器,所以面临计算资源的分配,接触了qsub和condor两种不同的任务提交模式。
qsub和condor其实都是用来运行一个脚本的,它与我们用bash或者sh来运行一个脚本的区别在于是否立即执行我们的脚本。
比如我写了一个脚本myscript.sh
我可以 bash myscript.sh 然后里面看到这个脚本被运行了
但是用 qsub myscript.sh, 就只是向我们的集群提交了一个任务而已,这个脚本什么时候运行,占用多少内存和cpu来运行,需要特定的设置好。
这里有个并行计算的ppt:http://www.marquette.edu/mugrid/bootcamp/2010Summer/Parallel.pdf 从计算机专业的角度来讲解了两者的区别。
其实只要看了这个https://wikis.nyu.edu/display/NYUHPC/Tutorial+-+Submitting+a+job+using+qsub 就能很容易明白qsub的用法,毕竟它只是一个命令而已,但是我们需要记住的是下面几个命令,因为经常用
PBS commands covered in this section
qsub
submit a job for execution
qstat
examine the status of a job (we have discussed what this status may be)
qhold
put a job on hold
qrls
release a job
qsig
send a signal to a job
qdel
delete a job

而condor其实很少见的。

An example Condor submission file is shown below:
Executable = ./myexecutable
Universe = vanilla
Log = condor.log
Output = out.log
Error = err.log
should_transfer_files = YES
when_to_transfer_output = ON_EXIT
transfer_input_files = mydata.txt
Save your script to a file with a name such as mysubmitfile.condor.
Once you have saved your submission script you can submit using the following command:
/dcs/condor/condor/bin/condor_submit mysubmitfile.condor

condor_q 可以用来查看任务提交情况

condor_rm 可以用来杀掉提交的任务。

任务提交模式还是挺有用的

25

Shell里面的各种括号的区别

[],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

[[ ]] double brackets

(())Double parentheses

{{}}double curly brackets

我们必须要记住的是下面

[] 相当于test,作逻辑判断

$( ) 与` ` (反引号) 都是用来做命令替换用

${ } 吧... 它其实就是用来作变量替换用的啦

(())就是用来计算的,相当于expr函数。

参考:http://sayle.net/book/

http://tldp.org/LDP/abs/html/index.html

 

我们首先看看一对的括号

首先[]是用来逻辑判断的,必须有空格

if [ -f binom.py ]

then

echo 'binom.py exists'

fi

或者

nub=$((i%4))

#echo $nub

if [ $nub == 0 ];then

echo "we need to sleep 4 hours"

sleep 14000

fi

这个[]操作符等价于test函数

if test $1 -gt 0
then
echo "$1 number is positive"
fi

但是都必须有空格!!!

参考:http://www.freeos.com/guides/lsst/ch03sec02.html

关于shell的test操作符还有很多http://tldp.org/LDP/abs/html/fto.html

( ) 将command group 置于 sub-shell 去执行,也称 nested sub-shell。

{ } 则是在同一个 shell 内完成,也称为non-named command group。

补充一个: {} 还可以做变量扩展 {5..9}  或者 {abcd}e, 自己运行一下就知道效果啦

这两个差异很小,而且一般用不着,就不讲了。

那么这一对的括号加上了$符号后又变成了上面鬼东西呢?

当然,只有$( ) ${ }才是合法的。

在 bash shell 中,$( ) 与` ` (反引号) 都是用来做命令替换用(command substitution)的。

在操作上,用$( ) 或` ` 都无所谓,用$( )的优点是:

1, ` ` 很容易与' ' ( 单引号)搞混乱,尤其对初学者来说

2, 在多层次的复合替换中,` ` 须要额外的跳脱( \` )处理,而$( ) 则比较直观

再让我们看${ } 吧... 它其实就是用来作变量替换用的啦。

一般情况下,$var 与${var} 并没有啥不一样。

但是用${ } 会比较精确的界定变量名称的范围,比方说:

[code][/code]

$ A=B

$ echo $AB

还可以用来截取变量,这个就很多花样啦

# 是去掉左边(在鉴盘上# 在$ 之左边)

% 是去掉右边(在鉴盘上% 在$ 之右边)

单一符号是最小匹配﹔两个符号是最大匹配

 

然后我们看看两对的括号:

nub=$((i%4)) 等价于$nub=`expr $i % 1` ;

((i++)) 等价于$i=`expr $i + 1` ;

所以(())就是用来计算的,而且里面的变量不需要$来标记啦

(在 $(( )) 中的变量名称,可于其前面加$ 符号来替换,也可以不用)

在(())前面加上$只是为了把计算结果给保存而已。

而两个中括号和两个大括号都是不合法的!

 

24

perl的模块组织方式

如何使用自己写的私人模块

模块通俗来讲,就是一堆函数的集合。

Personally I prefer to keep my modules (those that I write for myself or for systems I can control) in a certain directory, and also to place them in a subdirectory. As in:

/www/modules/MyMods/Foo.pm
/www/modules/MyMods/Bar.pm

And then where I use them:

use lib qw(/www/modules);useMyMods::Foo; 

useMyMods::Bar;

As reported by "perldoc -f use":

It is exactly equivalent to
BEGIN { require Module; import Module LIST; }
except that Module must be a bareword.

Putting that another way, "use" is equivalent to:

  • running at compile time,
  • converting the package name to a file name,
  • require-ing that file name, and
  • import-ing that package.

So, instead of calling use, you can call require and import inside a BEGIN block:

BEGIN{require'../EPMS.pm';
  EPMS->import();}

And of course, if your module don't actually do any symbol exporting or other initialization when you call import, you can leave that line out:

BEGIN{require'../EPMS.pm';}
比如我的一个模块如下,命名为my_stat.pm
package my_stat;
sub mean{
my $sum=0;
$sum+=$_ foreach @_;
$sum/($#_+1);
}
#print &mean(1..10),"\n";
sub stddev{
$avg=&mean(@_);
#print "$avg\n";
my $sum=0;
$sum+=($_-$avg)**2 foreach @_;
sqrt($sum/($#_));
#It will be different if you use $#_+1;
#sqrt($sum/($#_+1));
}
#print &stddev(1..10),"\n";
1;
里面有我定义好的两个函数 mean 和 stddev , 那么我就可以在我的其它perl程序里面直接引用这个模块,从而使用我的两个自定义函数。
use lib "./";  #取决于你把自定义模块my_stat.pm放在哪个目录
use my_stat;
print my_stat::stddev(1..10),"\n";
18

一个基因坐标定位到具体基因的程序的改进

这是为了回答以前的一个疑问:任意给定基因组的 chr:pos, 判断它在哪个基因上面?这个程序难吗?

基因的chr,start,end都是已知的

学术一点讲述这个问题:已知CNV数据在染色体上的position如chr1:2075000-2930999,怎样批量获取其对应的Gene Symbol呢(批量)

数据如下:

head gene_position.hg19 //21629

1 chr19 58858171 58874214 A1BG ENSG00000121410

2 chr12 9220303 9268558 A2M ENSG00000175899

3 chr12 9381128 9386803 A2MP1 ENSG00000256069

9 chr8 18027970 18081198 NAT1 ENSG00000171428

10 chr8 18248754 18258723 NAT2 ENSG00000156006

12 chr14 95058394 95090390  ENSG00000273259

13 chr3 151531860 151546276 AADAC ENSG00000114771

14 chr2 219128851 219134893 AAMP ENSG00000127837

15 chr17 74449432 74466199 AANAT ENSG00000129673

16 chr16 70286296 70323412 AARS ENSG00000090861

head pfam.df.hg19.bed   //340960

chr1  12190          12689          Helicase_C_2     0        +        12190          12689          255,255,0

chr1  69157          69220          7tm_4        0        +        69157          69220          255,255,0

chr1  69184          69817          7TM_GPCR_Srsx         0        +        69184          69817          255,255,0

chr1  69190          69931          7tm_1        0        +        69190          69931          255,255,0

chr1  69490          69910          7tm_4        0        +        69490          69910          255,255,0

现在需要对我们的pfam数据进行注释,根据每一行的chrpos来看看是属于哪一个基因

总共会有338879 pfam记录可以注释上基因。

注释之后应该是 head pfam.gene.df.hg19  这个样子

CDK11B       chr1  1571423     1573930     Pkinase        0        -         1571423     1573930     255,255,0

CDK11B       chr1  1572048     1573921     Pkinase_Tyr          0        -         1572048     1573921     255,255,0

CDK11B       chr1  1572120     1572823     Kinase-like  0        -         1572120     1572823     255,255,0

CDK11B       chr1  1572120     1572820     Kinase-like  0        -         1572120     1572820     255,255,0

CDK11B       chr1  1572120     1572817     Kinase-like  0        -         1572120     1572817     255,255,0

CDK11B       chr1  1573173     1573918     Kinase-like  0        -         1573173     1573918     255,255,0

CDK11B       chr1  1575747     1577317     Daxx  0        -         1575747     1577317     255,255,0

CDK11B       chr1  1576417     1577347     Nop14         0        -         1576417     1577347     255,255,0

CDK11B       chr1  1576423     1577332     Mitofilin      0        -         1576423     1577332     255,255,0

CDK11B       chr1  1576432     1577317     SAPS  0        -         1576432     1577317     255,255,0

我的第一个程序用的是全基因位点扫描到hash的方法。这样需要扫描13,1390,4974个位点,多于三分之一的基因组,这样是非常浪费内存的,尤其是keys需要多个字节。

我用了256G的服务器都没有运行完。

后来我取巧了把我的gene_position.hg19文件用split命令分成了25个,然后循环25次对pfam.df.hg19.bed  文件进行注释。

 

这样的确可以解决问了,而且只需要32G的内存的服务器即可,时间也很快,就十多分钟吧。

但这只是取巧的方法,应该要从算法上面优化,首先我仅仅做一个改动,就是不再扫描全基因的位点,对每个基因,我以1K的窗口来取位点进行扫描。这样我判断pfam的坐标时候,也以1K为最小单位进行判断。

这样只需要不到30s就可以出结果,总共注释了303474pfam记录,还不是最终的338879,因为我这次只注释了基因的1000整数倍基因区间,这样如果pfam记录落在一个基因的起始终止点不到1K位置时就不会被注释。这时候需要对代码进行继续优化。

 脚步懒得上传了,在我的有道云笔记里面。

http://note.youdao.com/share/?id=58e66d138e9434284ffa61c53b65abdc&type=note

28

R语言实现并行计算

前面我提到有一个大的运算任务需要很久才完成,所以用到了进度条来监控过程,但并不是改善了计算速度,所以需要用到并行计算,我又在网上找了找。

同样也是一个包,跟matlab的实现过程很像

library(parallel)

cl.cores <- detectCores() #检查当前电脑可用核数。

cl <- makeCluster(cl.cores) #使用刚才检测的核并行运算

#这里用clusterEvalQ或者par开头的apply函数族就可以进行并行计算啦

stopCluster(cl)

R-Doc里这样描述makeCluster函数:Creates a set of copies of R running in parallel and communicating over sockets. 即同时创建数个R进行并行运算。在该函数执行后就已经开始并行运算了,电脑可能会变卡一点。尤其在执行par开头的函数时。

在并行运算环境下,常用的一些计算方法如下:

1、clusterEvalQ(cl,expr)函数利用创建的cl执行expr。这里利用刚才创建的cl核并行运算expr。expr是执行命令的语句,不过如果命令太长的话,一般写到文件里比较好。比如把想执行的命令放在Rcode.r里:clusterEvalQ(cl,source(file="Rcode.r"))

2、par开头的apply函数族。这族函数和apply的用法基本一样,不过要多加一个参数cl。一般如果cl创建如上面cl <- makeCluster(cl.cores)的话,这个参数可以直接用作parApply(cl=cl,…)。当然Apply也可以是Sapply,Lapply等等。注意par后面的第一个字母是要大写的,而一般的apply函数族第一个字母不大写。

另外要注意,即使构建了并行运算的核,不使用parApply()函数,而使用apply()函数的话,则仍然没有实现并行运算。换句话说,makeCluster只是创建了待用的核,而不是并行运算的环境。

参考:http://www.r-bloggers.com/lang/chinese/1131

然后我模仿着用并行计算实现自己的需求

#it did work very fast
library(parallel)
cl.cores <- detectCores()
cl <- makeCluster(cl.cores)
clusterExport(cl, "all_dat_t")  #这里是重点,因为并行计算里面用到了自定义函数
clusterExport(cl, "all_prob_id") #但是这个函数需要用到这两个数据,所以需要把这两个数据加载到并行计算环境里面
prob_202723_s_at=parSapply( #我这里用的parSapply来实现并行计算
cl=cl,  #其中cl是我前面探测到的core数量,

deviation_prob, #deviation_prob是我待并行处理的向量

test_pro #这里其实应该是一个自定义函数,我这里就不写出来了,对上面的deviation_prob向量的每个探针都进行判断
)

28

R语言实现进度条

我也是临时在网上搜索到的教程,然后简单看了一下就实现了,其实就是就用到了一个名称为tcltk的包,直接查看函数tkProgressBar就可以知道怎么用啦!

下面是网上的一个小的示例代码(么有实际意义,仅为举例而已):

library(tcltk2)

u <- 1:2000

plot.new()

pb <- tkProgressBar("进度","已完成 %",  0, 100)

for(i in u) {

x<-rnorm(u)

points(x,x^2,col=i)

info <- sprintf("已完成 %d%%", round(i*100/length(u)))

setTkProgressBar(pb, i*100/length(u), sprintf("进度 (%s)", info), info)

}

close(pb)#关闭进度条

但是下面的代码是我模仿上面这个教程自己实现的。

[R]

# 以下是实现进度条
library(tcltk2)
plot.new()
pb <- tkProgressBar("进度","已完成 %", 0, 100)
prob_202723_s_at_value=rep(0,length(deviation_prob))
start_time=Sys.time() #这里可以计时,因为要实现进度条的一般都是需要很长运算时间
for (i in 1:length(deviation_prob)) {
tmp=test_pro(deviation_prob[i]) #test_pro是我自定义的一个函数,判断该探针是否符合要求。
if (length(tmp)!=0){prob_202723_s_at_value[i]=tmp}
info <- sprintf("已完成 %d%%", round(i*100/length(deviation_prob)))  #进度条就是根据循环里面的i来看看循环到哪一步了
setTkProgressBar(pb, i*100/length(deviation_prob), sprintf("进度 (%s)", info), info)
}
close(pb)#关闭进度条
end_time=Sys.time()
cat(end_time-start_time)

[/R]

28

R语言-比较数据框提取列的速度

结论:从数据框里面取某列数据,三种方法的时间消耗区别很大,直接用索引值,是最快的,而用$符号其次,用列名最慢。

我在R里面建立了一个表达量矩阵,列名是一个个样品,行是一个个探针,矩阵值是该探针在该样品测定的表达量。

那么,如果我要看看名为"202723_s_at"的探针的表达向量与其它所有探针的表达向量的相关系数,我可以用以下三种方法:

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t$"202723_s_at",x)))

user  system elapsed

22.93    0.03   23.03

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,"202723_s_at"],x)))

Timing stopped at: 92.02 5.32 97.66

太耗时间了,省去

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,grep(prob,names(all_dat_t))],x)))

Timing stopped at: 13.55 0.04 13.66

> prob_num=grep(prob,names(all_dat_t))

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,prob_num],x)))

user  system elapsed

8.14    0.01    8.17

可以看出,如果我首先根据探针名,grep出它在该表达量矩阵的列数,然后用列数来提取它的表达量是最快的,而且时间改善非常明显!

我们再探究一下cor函数的效率问题

探究的矩阵有54675个变量,每个变量均有189个观测值,如果取这个大矩阵的部分变量来求相关系数,结果如下!

> system.time(cor(all_dat_t[,1:10]))

user  system elapsed

0.001   0.000   0.001

> system.time(cor(all_dat_t[,1:100]))

user  system elapsed

0.003   0.000   0.003

> system.time(cor(all_dat_t[,1:1000]))

user  system elapsed

0.107   0.002   0.108

> system.time(cor(all_dat_t[,1:10000]))

user  system elapsed

11.102   0.849  11.983

> system.time(cor(all_dat_t)) 约六分钟也是可以搞定的

但是如果cor(all_dat_t),六分钟后得到的相关系数矩阵约32G,非常恐怖!

但是它很明显没有把这个32G相关系数矩阵存储到内存,因为我的机器本来就16G内存。我至今不能明白R具体实现机理

 

27

根据基因表达量对样品进行分类ConsensusClusterPlus

bioconductor系列的包都是一样的安装方式:
source("http://bioconductor.org/biocLite.R")
biocLite("ConsensusClusterPlus")

这个包是我见过最简单的包, 加载只有做好输入数据,只需要一句话即可运行,然后默认输出所有结果

读这个包的readme,很容易学会
就是做好一个需要来进行分类的样品的表达量矩阵。或者选择上一篇日志用GEOquery这个包下载的表达量矩阵也可以进行分析
因为这个包是用ALL数据来做测试的,所以可以直接加载这个数据结果,这样就能得到表达矩阵啦
library(ALL)
data(ALL)
d=exprs(ALL)
d[1:5,1:5]
可以看到数据集如下

> d[1:5,1:5]

             01005    01010    03002    04006    04007
1000_at   7.597323 7.479445 7.567593 7.384684 7.905312
1001_at   5.046194 4.932537 4.799294 4.922627 4.844565
1002_f_at 3.900466 4.208155 3.886169 4.206798 3.416923
1003_s_at 5.903856 6.169024 5.860459 6.116890 5.687997
1004_at   5.925260 5.912780 5.893209 6.170245 5.615210
> dim(d)
[1] 12625   128
共128个样品,12625个探针数据
也有文献用RNAs-seq的RPKM值矩阵来做
对上面这个芯片表达数据我们一般会简单的进行normalization ,然后取在各个样品差异很大的那些gene或者探针的数据来进行聚类分析
mads=apply(d,1,mad)
d=d[rev(order(mads))[1:5000],]

d = sweep(d,1, apply(d,1,median,na.rm=T))

#也可以对这个d矩阵用DESeq的normalization 进行归一化,取决于具体情况
library(ConsensusClusterPlus)
#title=tempdir() #这里一般改为自己的目录
title="./" #所有的图片以及数据都会输出到这里的
results = ConsensusClusterPlus(d,maxK=6,reps=50,pItem=0.8,pFeature=1,
 title=title,clusterAlg="hc",distance="pearson",seed=1262118388.71279,plot="png")
这样就OK了,你指定的目录下面会输出大于9个图片

clipboard

大家看看说明书就知道这个包的输出文件是什么了。
很多参数都是需要调整的,一般我们的maxK=6是根据实验原理来调整,如果你的样品应该是要分成6类以上,那么你就要把maxK=6调到一点。
查看结果results[[2]][["consensusClass"] 可以看到各个样品被分到了哪个类别里面去
results[[3]][["consensusClass"]
results[[4]][["consensusClass"] 等等
27

从GEO数据库下载矩阵数据-可以直接进行下游分析

bioconductor系列的包都是一样的安装方式:
source("http://bioconductor.org/biocLite.R")
    biocLite("GEOquery")

以前GEO数据库主要是microarray的芯片数据,后来有了RNA-seq,如果同时做多个样品的RNA-seq,表达量矩阵后来也可以上传到GEO数据库里面,只有看到文献里面有提到GEO数据库,都可以通过这个R包俩进行批量下载,其实就是网页版的一个API调用而已:

GEO数据库里面有四种数据

At the most basic level of organization of GEO, there are four basic entity types.
 The first three (Sample, Platform, and Series) are supplied by users;
the fourth, the dataset, is compiled and curated by GEO sta from the user-submitted data.
GEO accession number (GPLxxx). 
GEO accession number (GSMxxx)
GEO accession number (GSExxx). 
GEO DataSets (GDSxxx)
记住大小关系:一个GDS可以有多个GSM,一个GSM可以有多个GSE,至于GPL,一般不接触的
我们通常接触的都是GSE系列(一个GSE里面有多个GSM)的数据,而且这个包最重要的就是一个getGEO函数。
只要你通过文献确定了你的检索号,就可以通过这个函数来下载啦
检索号一般是A character string representing a GEO object for download and
          parsing.  (eg., 'GDS505','GSE2','GSM2','GPL96'

这个函数有很多参数,除非你需要下载的文件,那么就设置destdir到你喜欢的目录,如果只需要表达量数据就不用了。

 getGEO(GEO = NULL, filename = NULL, destdir = tempdir(), GSElimits=NULL,
     GSEMatrix=TRUE,AnnotGPL=FALSE)
例如:
gds <- getGEO("GDS10") 会返回一个对象,而且下载数据一般会在tmp目录下面,当然如果你需要保存这些文件,你可以自己制定下载目录及文件名。
gse2553 <- getGEO("GSE2553")
GDS2eSet函数可以把上面这个下载函数得到的对象(要确定是GDS而不是GSE)变成表达对象
pData和exprs函数都可以处理上面这个表达对象,从而分别得到样品描述矩阵和样品表达量矩阵
综合一起就是
g4100 <- GDS2eSet(getGEO("GDS4100"))
g4102 <- GDS2eSet(getGEO("GDS4102"))
e4102<-exprs(g4102)
e4100<-exprs(g4100)
这样的代码,这个e4100和e4102就都是一个数值矩阵啦,可以进行下游分析,但是如果是下载的GSM数据
就用下面这个代码,GSE26253_series_matrix.txt是通过GSEMatrix=TRUE这个参数特意下载到你的目录的
expr_dat=read.table("GSE26253_series_matrix.txt",comment.char="!",stringsAsFactors=F)
这样读取也是一个数值矩阵
具体大家可以看这个包的说明书
#Download GDS file, put it in the current directory, and load it:
gds858 <- getGEO('GDS858', destdir=".")
如果使用了GSEMatrix=TRUE这个参数,那么除了下载soft文件,还有表达量矩阵文件,可以直接用read.table读取那个文件。
#Or, open an existing GDS file (even if its compressed):
gds858 <- getGEO(filename='GDS858.soft.gz')
下面这个下载的是GSE对象,GDS对象还有大一点
GEOquery-下载结果gds对象