Tag Archives: R

十二 29

生信分析人员如何系统入门R？

Posted on by
R语言的重要性我也就不再赘述了，它不仅在生物信息数据处理中发挥着重要作用，也是其它主流数据处理人士的首选工具。现在非常多的半路出家自学生物信息学的小伙伴必须学而且有可能学的就是R，所以写一个R的系统性入门指导是非常有必要的。这本来应该是我下面的生信分析人员如何系统入门编程语言的姊妹篇的，但是因为时隔太久，我的感悟可能发生了变化，所以这个R跟前面的两个看起来总结指引模式不太一样的。

生信分析人员如何系统入门perl？

生信分析人员如何系统入门linux？

我作为老一辈的生信工程师，所以喜欢perl一点，排斥python，其实呢，我也稍微看过一些python的语法，个人认为R和python几乎是一模一样的。R的特点就是内置了大量的函数，基本上你认识的英文单词都可以是一个函数，即使不是，你也可以自定义为函数。搞清楚了函数和变量，就可以看懂大部分的R代码了。

下面是生信菜鸟团QQ群管理员赵云对这3种编程语言的心得体会！

python跟perl都是高级语言， 两个开发的目的不同， perl更面向过程一些，优势是严谨，快。 python主流面向对象编程， 这个跟R类似， 数据结构等方面有些不同，但可以互相调用。 实际上以上三者之间可以互相调用部分功能。python的语法并不是很严谨，个人感觉，越偏向自然语言的编程语言越通俗但不严谨，以上，是跟C比较的。

R本身起源于S语言，是主要针对统计的， 也是面向对象的， 本质上，，是把一个比excel功能强大的软件归零化成了命令行吧.excel高级应用也是要编程的， 所以R的初级应用可以当成是没有用户交互界面的excel，细心一点， 把示例代码都打对，当功能强大但不好使版的excel吧， 这样至少心理上不会畏难跟抵触.

内部集成的越多，用户需要做的越少， 你用C画个图累死你， 用python得写几行， R一行就行了！
PS:菜鸟发言，如有误导，概不负责！

Continue reading →

十一 28

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

Posted on by

经过热心的小伙伴的提醒，我才知道我以前写的R语言画网络图三部曲竟然漏掉了最基础的一个包，就是igraph，不了解这个，后面的两个也是无源之水。

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

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

其实包括了3个包：igraph/RBGL/Rgraphviz
用到了一个测试数据，是构建好的PPI网络对象：We will first analyse a curated data set of protein-protein interactions in the yeast Saccharomyces cerevisiae extracted from published papers. This data set comes from with an R package called “yeastExpData”, which calls the data set “litG”. This data was first described in a paper by Ge et al (2001) in Nature Genetics (http://www.nature.com/ng/journal/v29/n4/full/ng776.html).

Continue reading →

十 15

R一大利器之对象的操作函数查询

Posted on by

对于生物出身的部分生物信息学工程师来说，很多计算机概念让人很头疼，尤其是计算机语言里面的高级对象。我以前学编程的时候，给我一个变量，一个数据，一个hash，我就心满意足了，可以解决大部分我数据处理问题，可事情远比想象之中复杂。因为很多高手喜欢用封装，代码复用，喜欢用高级对象。在R的bioconductor里面尤其是如此，经常会遇到各种包装好的S3，S4对象，看过说明书，倒是知道一些对象里面有什么，可以去如何处理那些对象，提取我们想要的信息，比如我就写过一系列的帖子：

Continue reading →

六 04

R的shiny 服务器管理-入门

Posted on by

如果你已经安装好了shiny 服务器，(安装教程)要开始使用了，掌握一些基础知识是必须的。这里我简单学习了一些入门资料，分享给大家，慢慢的我会写一个进阶资料。安装成功之后，系统会增加4个目录，是一定要掌握的：

1、这个目录只存放关键配置文件：/etc/shiny-server/shiny-server.conf   初始状态只有一个文件，记录着非常多的默认信息，默认的网站目录是根目录下的srv的shiny-server目录，端口是3838
2、网站运行log日子存放：/var/log/shiny-server  初始状态下该目录为空
3、程序存放目录是：/srv/shiny-server 初始状态，有一个测试程序：
4、最后是/opt/shiny-server/ 目录，这里面也有一个配置文件：/opt/shiny-server/config/default.config

Continue reading →

六 04

安装自己的shiny服务器-实战指南

Posted on by

个人比较欣赏R shiny制作的网页，入门简单，上手极快，多看点例子，制作复杂逻辑的网页也不是问题。这篇实战指南有四个步骤：

至少需要root权限的linux系统  （我测试了阿里云）
安装R   （一般安装最新版，）
在R中安装shiny模块   （一般还可以多安装一些模块）
下载并且安装shiny server安装包    （根据系统选择）

Continue reading →

五 23

读书笔记（R语言）

Posted on by

R与ASReml-R统计分析教程（林元震）中国林业出版社

1-3章简单介绍了R的基本语法，然后第4章着重讲了各种统计方法，第5章讲R的绘图，最后一张讲ASReml-R这个包
语法重点:

1,install.packages(),library(),help(),example(),demo(),length(),attribute(),class(),mode(),dim(),names(),str(),head(),
tail()

2,rep,seq,paste,array,matrix,data.frame,list,c(),factor(),

3,缺失值处理（na.omit,na.rm=T），类型转换（as.numeric(),as.character(),as.factor(),as.logical()）

Continue reading →

四 12

R包精讲第四篇：4种R包安装方式

Posted on by

请先看：R包精讲第一篇：如何查看你已经安装了和可以安装哪些R包？

第一种方式，当然是R自带的函数直接安装包了，这个是最简单的，而且不需要考虑各种包之间的依赖关系。

对普通的R包，直接install.packages()即可，一般下载不了都是包的名字打错了，或者是R的版本不够，如果下载了安装不了，一般是依赖包没弄好，或者你的电脑缺少一些库文件，如果实在是找不到或者下载慢，一般就用repos=来切换一些镜像。

> install.packages("ape")  ##直接输入包名字即可
Installing package into ‘C:/Users/jmzeng/Documents/R/win-library/3.1’
(as ‘lib’ is unspecified)  ##一般不指定lib，除非你明确知道你的lib是在哪里
trying URL 'http://mirror.bjtu.edu.cn/cran/bin/windows/contrib/3.1/ape_3.4.zip'
Content type 'application/zip' length 1418322 bytes (1.4 Mb)
opened URL   ## 根据你选择的镜像，程序会自动拼接好下载链接url
downloaded 1.4 Mb

package ‘ape’ successfully unpacked and MD5 sums checked  ##表明你已经安装好包啦

The downloaded binary packages are in  ##程序自动下载的原始文件一般放在临时目录，会自动删除
	C:\Users\jmzeng\AppData\Local\Temp\Rtmpy0OivY\downloaded_packages
>

对于bioconductor的包，我们一般是

source("http://bioconductor.org/biocLite.R") ##安装BiocInstaller

#options(BioC_mirror=”http://mirrors.ustc.edu.cn/bioc/“) 如果需要切换镜像
biocLite("ggbio")

或者直接BiocInstaller::biocLite('ggbio') ## 前提是你已经安装好了BiocInstaller

某些时候你还需要卸载remove.packages("BiocInstaller") 然后安装新的

第二种方式，是直接找到包的下载地址，需要进入包的主页

packageurl <- "http://cran.r-project.org/src/contrib/Archive/ggplot2/ggplot2_0.9.1.tar.gz"
packageurl <- "http://cran.r-project.org/src/contrib/Archive/gridExtra/gridExtra_0.9.1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
#packageurl <- "http://www.bioconductor.org/packages/2.11/bioc/src/contrib/ggbio_1.6.6.tar.gz"
#packageurl <- "http://cran.r-project.org/src/contrib/Archive/ggplot2/ggplot2_1.0.1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")

这样安装的就不需要选择镜像了，也跨越了安装器的版本！

第三种是，先把包下载到本地，然后安装：

download.file("http://bioconductor.org/packages/release/bioc/src/contrib/BiocInstaller_1.20.1.tar.gz","BiocInstaller_1.20.1.tar.gz")
##也可以选择用浏览器下载这个包
install.packages("BiocInstaller_1.20.1.tar.gz", repos = NULL)
## 如果你用的RStudio这样的IDE，那么直接用鼠标就可以操作了
或者用choose.files()来手动交互的选择你把下载的源码BiocInstaller_1.20.1.tar.gz放到了哪里。

这种形式大部分安装都无法成功，因为R包之间的依赖性很强！

第四种是：命令行版本安装

如果是linux版本，命令行从网上自动下载包如下：
sudo su - -c \
"R -e \"install.packages('shiny', repos='https://cran.rstudio.com/')\""
如果是linux，命令行安装本地包，在shell的终端
sudo R CMD INSTALL package.tar.gz
window或者mac平台一般不推荐命令行格式，可视化那么舒心，何必自讨苦吃

 

四 12

R包精讲第三篇：如何切换镜像？

Posted on by

这个技巧很重要，一般来说，R语言自带的install.packages函数来安装一个包时，都是用的默认的镜像！

如果你是用的Rstudio这个IDE，你的默认镜像就是： https://cran.rstudio.com/

如果你直接用的R语言，那么就是："http://cran.us.r-project.org" 但是一般你安装的时候会提醒你选择。

Continue reading →

四 12

R包精讲第二篇：如何安装旧版本的包？

Posted on by 
既然你点进来看，肯定是有需求咯！
一般来说，R语言自带的install.packages函数来安装一个包时，都是默认安装最新版的。
但是有些R包的开发者他会引用其它的一些R包，但是它用的是人家旧版本的功能，但他自己来不及更新或者疏忽了。
而我们又不得不用他的包，这时候就不得不卸载最新版包，转而安装旧版本包。

Continue reading →

四 11

R包精讲第一篇：如何查看你已经安装了和可以安装哪些R包？

Posted on by

最近经常出现一个错误，类似于package ‘airway’ is not available (for R version 3.1.0)

就是某些包在R的仓库里面找不到，这个错误非常普遍，stackoverflow上面非常详细的解答：

http://stackoverflow.com/questions/25721884/how-should-i-deal-with-package-xxx-is-not-available-for-r-version-x-y-z-wa

在阅读这个答案的时候，我发现了一个非常有用的函数！available.packages()可以查看自己的机器可以安装哪些包！

Continue reading →

三 03

用R语言包从EBI的arrayexpress数据库里面下载芯片数据

Posted on by
这个包跟GEOquery区别不是很大，只不过一个是正对NCBI的GEO数据库，一个是针对EBI的arrayexpress数据库，只有对写自动化脚本的人来说才有需求，一般个人分析者都是自己去数据库主页里面查找，然后拿到下载链接，一个个下载。
从EBI的arrayexpress数据库里面下载芯片数据：
update to 2016-3-1 11:41:27
63890 experiments
1912744 assays
40.53 TB of archived data 数据量还是蛮大的
所有的data，都可以在ftp服务器里面下载：ftp://ftp.ebi.ac.uk/pub/databases/arrayexpress/data/experiment/BUGS/
根据ID号很整齐的储存着。
也可以用一个R语言包：ArrayExpress R package
2009年，那个时候R语言用的人很少，这个简单的包都可以发文章，现在看来简直不可思议！
其实大部分数据都是跟GEO数据库对应的：比如https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-55645/  对应于：GEO - GSE55645
比如对NASH表达数据查找：https://www.ebi.ac.uk/arrayexpress/search.html?query=NASH++expression  30条结果里面只有4条是arrayexpress数据库独有的！
biocLite("ArrayExpress")
library(ArrayExpress)
1
如果用R语言，搜索如下：
可以用sets = queryAE(keywords = "NASH+expression", species = "homo+sapiens")
2
效果是一样的！
下载数据用：
back = getAE("E-MEXP-3291")
下载其实也就是里面存储了链接，直接调用R语言的下载函数即可！
3
一般没必要下载原始测序文件，直接用下面这个函数就可以得到一个数据对象，可以直接得到表达矩阵和实验的metadata

rawset = ArrayExpress("E-MEXP-3291")

二 28

在R里面操作SQLite

Posted on by
我前面写到过如何把数据库写到mysql，但是发现其实msyql并不方便，需要连接数据库什么的，如果发布一个离线小网页，这时候sqlite的优点就显示出来了！
基础代码很简单：
library(RSQLite)
sqlite    <- dbDriver("SQLite")
con <- dbConnect(sqlite,"hg19_bioconductor.sqlite") # makes a new file
suppressMessages(library(org.Hs.eg.db))
kegg2ID=toTable(org.Hs.egPATH)
#[1] "gene_id" "path_id"
dbWriteTable(con,'keggID2geneID',kegg2ID,row.name=F,overwrite=T)
具体代码，可以看我的github主页：https://github.com/jmzeng1314/my-R/blob/master/3-get-hg19-gene-mapping/get-hg19-gene-mapping-bioconductor.R
做出来的数据，如下，就是几个table存储在文件里面！
 2
最后这些数据都保存在了当前工作目录下的hg19_bioconductor.sqlite文件里面！
在其它程序里面就可以直接调用这个文件，而不需要加载一大堆的包了！
1
library(KEGG.db)
library(GO.db)

library(org.Hs.eg.db )

 

二 15

用R获取芯片探针与基因的对应关系三部曲-bioconductor

Posted on by

现有的基因芯片种类不要太多了！

但是重要而且常用的芯片并不多！
一般分析芯片数据都需要把探针的ID切换成基因的ID，我一般喜欢用基因的entrez ID。
一般有三种方法可以得到芯片探针与gene的对应关系。
金标准当然是去基因芯片的厂商的官网直接去下载啦！！！
一种是直接用bioconductor的包

一种是从NCBI里面下载文件来解析好！
首先，我们说官网，肯定可以找到，不然这种芯片出来就没有意义了！
然后，我们看看NCBI下载的，会比较大
这两种方法都比较麻烦，需要一个个的来！
所以我接下来要讲的是用R的bioconductor包来批量得到芯片探针与gene的对应关系！
一般重要的芯片在R的bioconductor里面都是有包的，用一个R包可以批量获取有注释信息的芯片平台，我选取了常见的物种，如下：
        gpl           organism                  bioc_package
1     GPL32       Mus musculus                        mgu74a
2     GPL33       Mus musculus                        mgu74b
3     GPL34       Mus musculus                        mgu74c
6     GPL74       Homo sapiens                        hcg110
7     GPL75       Mus musculus                     mu11ksuba
8     GPL76       Mus musculus                     mu11ksubb
9     GPL77       Mus musculus                     mu19ksuba
10    GPL78       Mus musculus                     mu19ksubb
11    GPL79       Mus musculus                     mu19ksubc
12    GPL80       Homo sapiens                        hu6800
13    GPL81       Mus musculus                      mgu74av2
14    GPL82       Mus musculus                      mgu74bv2
15    GPL83       Mus musculus                      mgu74cv2
16    GPL85  Rattus norvegicus                        rgu34a
17    GPL86  Rattus norvegicus                        rgu34b
18    GPL87  Rattus norvegicus                        rgu34c
19    GPL88  Rattus norvegicus                         rnu34
20    GPL89  Rattus norvegicus                         rtu34
22    GPL91       Homo sapiens                      hgu95av2
23    GPL92       Homo sapiens                        hgu95b
24    GPL93       Homo sapiens                        hgu95c
25    GPL94       Homo sapiens                        hgu95d
26    GPL95       Homo sapiens                        hgu95e
27    GPL96       Homo sapiens                       hgu133a
28    GPL97       Homo sapiens                       hgu133b
29    GPL98       Homo sapiens                     hu35ksuba
30    GPL99       Homo sapiens                     hu35ksubb
31   GPL100       Homo sapiens                     hu35ksubc
32   GPL101       Homo sapiens                     hu35ksubd
36   GPL201       Homo sapiens                       hgfocus
37   GPL339       Mus musculus                       moe430a
38   GPL340       Mus musculus                     mouse4302
39   GPL341  Rattus norvegicus                       rae230a
40   GPL342  Rattus norvegicus                       rae230b
41   GPL570       Homo sapiens                   hgu133plus2
42   GPL571       Homo sapiens                      hgu133a2
43   GPL886       Homo sapiens                     hgug4111a
44   GPL887       Homo sapiens                     hgug4110b
45  GPL1261       Mus musculus                    mouse430a2
49  GPL1352       Homo sapiens                       u133x3p
50  GPL1355  Rattus norvegicus                       rat2302
51  GPL1708       Homo sapiens                     hgug4112a
54  GPL2891       Homo sapiens                       h20kcod
55  GPL2898  Rattus norvegicus                     adme16cod
60  GPL3921       Homo sapiens                     hthgu133a
63  GPL4191       Homo sapiens                       h10kcod
64  GPL5689       Homo sapiens                     hgug4100a
65  GPL6097       Homo sapiens               illuminaHumanv1
66  GPL6102       Homo sapiens               illuminaHumanv2
67  GPL6244       Homo sapiens   hugene10sttranscriptcluster
68  GPL6947       Homo sapiens               illuminaHumanv3
69  GPL8300       Homo sapiens                      hgu95av2
70  GPL8490       Homo sapiens   IlluminaHumanMethylation27k
71 GPL10558       Homo sapiens               illuminaHumanv4
72 GPL11532       Homo sapiens   hugene11sttranscriptcluster
73 GPL13497       Homo sapiens         HsAgilentDesign026652
74 GPL13534       Homo sapiens  IlluminaHumanMethylation450k
75 GPL13667       Homo sapiens                        hgu219
76 GPL15380       Homo sapiens      GGHumanMethCancerPanelv1
77 GPL15396       Homo sapiens                     hthgu133b
78 GPL17897       Homo sapiens                     hthgu133a
这些包首先需要都下载
gpl_info=read.csv("GPL_info.csv",stringsAsFactors = F)
### first download all of the annotation packages from bioconductor
for (i in 1:nrow(gpl_info)){
  print(i)
  platform=gpl_info[i,4]
  platform=gsub('^ ',"",platform) ##主要是因为我处理包的字符串前面有空格
  #platformDB='hgu95av2.db'
  platformDB=paste(platform,".db",sep="")
  if( platformDB  %in% rownames(installed.packages()) == FALSE) {
    BiocInstaller::biocLite(platformDB)
    #biocLite(platformDB )
  }
}
下载完了所有的包， 就可以进行批量导出芯片探针与gene的对应关系！
for (i in 1:nrow(gpl_info)){
  print(i)
  platform=gpl_info[i,4]
  platform=gsub('^ ',"",platform)
  #platformDB='hgu95av2.db'
  platformDB=paste(platform,".db",sep="")
  if( platformDB  %in% rownames(installed.packages()) != FALSE) {
    library(platformDB,character.only = T)
    #tmp=paste('head(mappedkeys(',platform,'ENTREZID))',sep='')
    #eval(parse(text = tmp))
###重点在这里，把字符串当做命令运行
    all_probe=eval(parse(text = paste('mappedkeys(',platform,'ENTREZID)',sep='')))
    EGID <- as.numeric(lookUp(all_probe, platformDB, "ENTREZID"))
##自己把内容写出来即可
  }
}

 

 

一 22

用TCGA数据做cox生存分析的风险因子（比例风险模型）

Posted on by
再次强调一下，R里面实现生存分析非常简单！
用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')构建生存曲线。
用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)来做某一个因子的KM生存曲线。用 survdiff(my.surv~type, data=dat)来看看这个因子的不同水平是否有显著差异，其中默认用是的logrank test 方法。
用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。

Continue reading →

一 08

别人写的代码运行真快！！！

Posted on by
最近需要做十几万个差异基因分析，每个分析都是对大约5万个探针，200个样本的数据量进行批量T检验计算P值
然后发现自己无论怎么用R来写，每个分析都要耗时半分钟左右，因为我必须循环所有的探针，即使不用for，而用R推荐的apply系列函数，也快不到哪里去，但是我搜索时候发现有一个package里面自带了矩阵T检验，直接对5万个探针进行T检验，而不需要循环处理它们
看下面代码！
dat=matrix(rnorm(10000000),nrow = 50000)
dim(dat) #50000   200
system.time(
  apply(dat,1,function(x){
    t.test(x[1:100],x[101:200])$p.value
  })
)
#用户  系统  流逝
#29.29  0.04 30.64
library(pi0)
system.time(matrix.t.test(dat,1,100,100))
#用户 系统 流逝
#0.48 0.03 0.53
差距真的是非常的明显呀！！！
然后，我解析了它的代码，发现里面调用了C写的代码，我想这就是问题所在咯，可是他们到底怎么写，才能把速度搞这么快？？？
  tmp = .C("tstatistic", dat = x, n1 = n1, n2 = n2, ntests = ntests, 
        MARGIN = MARGIN, pool = pool, tstat = rep(0, ntests), 
        df = rep(0, ntests), PACKAGE = "pi0")

源码在这个package的github里面可以找到，有兴趣的童鞋可以研究一下

 
 

 

一 06

R语言软件的各种旧版本下载

Posted on by
这种编程语言，我之前以为只有包，或者模块什么的才烦人，没想到最近还碰到了版本的问题。
一般来说， 去R的官网，可以下载到基于各种操作系统的最新版R软件，但是在某些特别的情况下，我们可能需要下载以前的旧版本。这时候，想用百度这个破东西搜索出来下载地址简直是痴心妄想。
比如，我想搜索2.15.2这个版本，在Google里面很容易找到了
如果进入bin目录，可以看到各种操作系统的软件https://cran.r-project.org/bin/ ，里面都是可执行的文件，有下面几种操作系统
[DIR] linux/ 23-Jan-2008 19:47 -
[DIR] macos/ 19-Apr-2005 09:45 -
[DIR] macosx/ 12-Dec-2015 09:04 -
[DIR] windows/ 24-Feb-2012 18:41 -
大多数情况下面我们会用window这个平台的，毕竟，可视化嘛！
很多时候，我们也会需要linux平台的https://cran.r-project.org/bin/linux/ 但是linux平台本身就比较复杂
[DIR] debian/ 15-Dec-2015 02:06 -
[DIR] redhat/ 27-Jul-2014 21:12 -
[DIR] suse/ 16-Feb-2012 15:09 -
[DIR] ubuntu/ 06-Jan-2016 04:05 -
我用的Ubuntu比较多，即使选择了ubuntu，里面还有好几个选项，因为Ubuntu也有各种版本！
[DIR] precise/ 06-Jan-2016 04:03 -
[DIR] trusty/ 06-Jan-2016 04:04 -
[DIR] vivid/ 06-Jan-2016 04:04 -
[DIR] wily/ 06-Jan-2016 04:05 -

所以如果，大家是想在linux系统里面安装旧版本的R，建议大家直接下载c源码，直接三部曲就可以安装啦！

[DIR] R-0/ 04-Oct-2004 10:20 -
[DIR] R-1/ 04-Oct-2004 19:02 -
[DIR] R-2/ 01-Mar-2013 09:11 -
[DIR] R-3/ 10-Dec-2015 09:13 -
windows系统里面一般是是exe的可执行文件，直接安装，不需要源码变异，只需要不停的下一步即可
R 3.2.2 (August, 2015)
R 3.2.1 (June, 2015)
R 3.2.0 (April, 2015)
R 3.1.3 (March, 2015)
R 3.1.2 (October, 2014)
R 3.1.1 (July, 2014)
R 3.1.0 (April, 2014)
R 3.0.3 (March, 2014)
R 3.0.2 (September, 2013)
R 3.0.1 (May, 2013)
R 3.0.0 (April, 2013)
R 2.15.3 (March, 2013)
R 2.15.2 (October, 2012)
R 2.15.1 (June, 2012)
R 2.15.0 (March, 2012)
R 2.14.2 (February, 2012)
R 2.14.1 (December, 2011)
R 2.14.0 (November, 2011)
R 2.13.2 (September, 2011)
R 2.13.1 (July, 2011)
R 2.13.0 (April, 2011)
R 2.12.2 (February, 2011)
R 2.12.1 (December, 2010)
R 2.12.0 (October, 2010)
R 2.11.1 (May, 2010)
R 2.11.0 (April, 2010)
R 2.10.1 (December, 2009)
R 2.10.0 (October, 2009)
R 2.9.2 (August, 2009)
R 2.9.1 (June, 2009)
R 2.9.0 (April, 2009)
R 2.8.1 (December, 2008)
R 2.8.0 (October, 2008)
R 2.7.2 (August, 2008)
R 2.7.1 (June, 2008)
R 2.7.0 (April, 2008)
R 2.6.2 (February, 2008)
R 2.6.1 (November, 2007)
R 2.6.0 (October, 2007)
R 2.5.1 (July, 2007)
R 2.5.0 (April, 2007)
R 2.4.1 (December, 2006)
R 2.4.0 (October, 2006)
R 2.3.1 (June, 2006)
R 2.3.0 (April, 2006)
R 2.2.1 (December, 2005)
R 2.2.0 (October, 2005)
R 2.1.1 (June, 2005)
R 2.1.0 (April, 2005)
R 2.0.1 (November, 2004)
R 2.0.0 (October, 2004)
R 1.9.1 (June, 2004)
R 1.8.1 (November, 2003)
R 1.7.1 (June, 2003)
R 1.6.2 (January, 2003)
Installer for R 1.5.1 (June, 2002)
Installer for R 1.4.1 (January, 2002)
Installer for R 1.3.1 (September, 2001)
Binary files for R 1.2.2 (March, 2001)
Binary files for R 1.0.0 (February, 2000)

 

十 30

R语言也可以用hash啦

Posted on by

本来想着是如何习惯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

 

八 31

文献笔记-2010-R-softeware-identify-cancer_driver_genes

Posted on by

我们用188 non-small cell lung tumors数据来测试了一个R语言程序，find driver genes in cancer ~
软件地址如下：http://linus.nci.nih.gov/Data/YounA/software.zip
这是一个R语言程序，里面有readme，用法很简单。
准备好两个文件，分别是silent_mutation_table.txt and nonsilent_mutation_table.txt ，它们都是普通文本格式数据，内容如下，就是把找到的snp格式化，根据注释结果分成silent和nonsilent即可。
#Ensembl_gene_id Chromosome Start_position Variant_Type Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2 Tumor_Sample_Barcode
#ENSG00000122477 1 100390656 SNP G G A TCGA-23-1022-01A-02W-0488-09
然后直接运行程序包里面的主程序，在R语言里面source("main_R_script.r")
We reanalyzed sequence data for 623 candidate genes in 188 non-small cell lung tumors using the new method.
to identify genes that are frequently mutated and thereby are expected to have primary roles in thedevelopment of tumor
To find these driver genes, each gene is tested for whether its mutation rate is significantly higher than the background (or passenger) mutation rate.

Some investigators (Sjoblom et al., 2006) further divide mutations into several types according to the nucleotide and the neighboring nucleotides of the mutations.

Ding et al. (2008)的方法的三个缺点：
1、different types of mutations can have different impact on proteins.（越影响蛋白功能的突变，越有可能是driver mutation）
2、different samples have different background mutation rates. （在高突变背景的样本中的突变，很可能是高突变背景的原因，而不是因为癌症）
3、a different number of non-silent mutations can occur at each base pair according to the genetic code.（比如Tryptophan仅仅只有一个密码子，而arginine高达6个密码子）

我们提出的方法的4个优点是：
1，我们对非同义突变根据它们对蛋白功能的影响进行了评级打分。
2，我们允许不同的样品有着不同的BMR
3，that whether the mutation is non-silent or silent depends on the genetic code
4，we take into account uncertainties in the background mutation rate by using empirical Bayes methods

还有5个需要改进的地方：
1，However, the functional impact is also dependent on the position in which a mutation occurs.（我们仅仅考虑了突变对氨基酸的改变）
2，the current scoring system which assigns mutation scores in the order: missense mutation<inframe indel<mutation in splice sites<frameshift indel=nonsense mutation may be biased toward identifying tumor suppressor genes over oncogenes.
3，we may refine our background mutation model in Table 1 so that all six types of mutations, A:T→G:C, A:T→C:G, A: T→T :A,G:C→A:T, G:C→T :A, G:C→C:G have separate mutation rates.
4，we did not take into account correlations among mutations in identifying driver genes.
5，one might combine both copy number variation and sequencing data to identify driver genes.

HGNC定义的gene Symbol转为ensemble数据库的ID，的R语言代码：
library(biomaRt)
ensembl=useMart("ensembl",dataset = "hsapiens_gene_ensembl")
all.gene.table = read.table("all_gene.symbol", header=F)
convert=getBM(attributes = c("chromosome_name","ensembl_gene_id","hgnc_symbol"),filters =c("hgnc_symbol"),values=all.gene.table[,1],mart=ensembl)
chromosome=c(1:22,"X","Y","M")
convert=convert[!is.na(match(convert[,1],chromosome)),2:3] #remove names whose matching chromosome is not 1-22, X, or Y.
convert=convert[rowSums(convert=="")==0,]
write.table(convert,"ensembl2symbol.list",quote = F,row.names =F,col.names =F)
write.table(convert,"all_gene_name.txt",quote = F,row.names =F,col.names =F)

一个gene Symbol可能对应着多个ensemble ID号，但是在每个染色体上面是一对一的关系。
有些gene Symbol可能找不到ensemble ID号，一般情况是因为这个gene Symbol并不是纯粹的HGNC定义的，或者是比较陈旧的ID。
比如下面的TIGAR ，就很可能被写作是C12orf5
Aliases for TIGAR Gene
TP53 Induced Glycolysis Regulatory Phosphatase 2 3
TP53-Induced Glycolysis And Apoptosis Regulator 2 3 4
C12orf5 3 4 6
Probable Fructose-2,6-Bisphosphatase TIGAR 3
Fructose-2,6-Bisphosphate 2-Phosphatase 3
Chromosome 12 Open Reading Frame 5 2
Fructose-2,6-Bisphosphatase TIGAR 3
Transactivated By NS3TP2 Protein 3
EC 3.1.3.46 4
FR2BP 3
External Ids for TIGAR Gene
HGNC: 1185 Entrez Gene: 57103 Ensembl: ENSG00000078237 OMIM: 610775 UniProtKB: Q9NQ88
Previous HGNC Symbols for TIGAR Gene
C12orf5
Export aliases for TIGAR gene to outside databases

七 21

R语言用hclust进行聚类分析

Posted on by

聚类的基础就是算出所有元素两两间的距离，我们首先做一些示例数据，如下：

x=runif(10)

y=runif(10)

S=cbind(x,y)                                 #得到2维的数组

rownames(S)=paste("Name",1:10,"")             #赋予名称，便于识别分类

out.dist=dist(S,method="euclidean")           #数值变距离

这个代码运行得到的S是一个矩阵，如下

> S

x         y

Name 1   0.41517985 0.4697017

Name 2   0.35653781 0.1132367

Name 3   0.52253349 0.3680286

Name 4   0.80558684 0.9834687

Name 5   0.04564145 0.8560690

Name 6   0.11044397 0.2988598

Name 7   0.34984447 0.8515141

Name 8   0.28097709 0.1260050

Name 9   0.81771888 0.5976135

Name 10 0.40700158 0.5236567

可以看出里面共有10个点，它们的X,Y坐标均已知，我们有6总方法可以求矩阵

image001

注释：在聚类中求两点的距离有：

1，绝对距离：manhattan

2，欧氏距离：euclidean 默认

3，闵科夫斯基距离：minkowski

4，切比雪夫距离：chebyshev

5，马氏距离：mahalanobis

6，蓝氏距离：canberra

用默认的算法求出距离如下

算出距离后就可以进行聚类啦！

out.hclust=hclust(out.dist,method="complete") #根据距离聚类

注释：聚类也有多种方法：

1，类平均法：average

2，重心法：centroid

3，中间距离法:median

4，最长距离法：complete 默认

5，最短距离法：single

6，离差平方和法：ward

7，密度估计法：density

接下来把聚类的结果图画出来

plclust(out.hclust)                           #对结果画图

image003

rect.hclust(out.hclust,k=3)                   #用矩形画出分为3类的区域

image005

out.id=cutree(out.hclust,k=3)                 #得到分为3类的数值

这里的out.id就是把每个点都分类了的分类数组，1,2,3.