Tag Archives: bioconductor

十 29

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

Posted on by

经常会有人问这样的问题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"))

八 27

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

Posted on by 
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对象
七 19

Bioconductor系列之GenomicAlignments

Posted on by

Bioconductor系列包的安装方法都一样

source("http://bioconductor.org/biocLite.R")biocLite(“GenomicAlignments”)

这个包设计就是用来处理bam格式的比对结果的，具体功能非常多，其实还不如自己写脚本来做这些工作，更方便一点，实在没必要学别人的语法。大家有兴趣可以看GenomicAlignments的主页介绍，三个pdf文档来介绍。

http://www.bioconductor.org/packages/release/bioc/html/GenomicAlignments.html

我们首先读取示例pasillaBamSubset包带有的bam文件

library(pasillaBamSubset)

un1 <- untreated1_chr4()  # single-end reads

library(GenomicAlignments)

reads1 <- readGAlignments(un1)

查看reads1这个结果，可以看到把这个bam文件都读成了一个数据对象GAlignments object，

Bioconductor系列之GenomicAlignments632

readGAlignments这个函数就是读取我们bam文件的，生成的对象可以用多个函数来继续查看信息，比如它的列名都是函数

Seqnames(),strand(),cigar(),qwidth(),start(),end(),width(),njunc() 这些函数对这个GAlignments对象处理后会得到比对的各个情况的向量。

我们也可以读取这个GenomicAlignments包自带的bam文件，而不是pasillaBamSubset包带有的bam文件

> fls <- list.files(system.file("extdata", package="GenomicAlignments"),+                   recursive=TRUE, pattern="*bam$", full=TRUE)

> fls

[1] "C:/Program Files/R/R-3.1.2/library/GenomicAlignments/extdata/sm_treated1.bam"  [2] "C:/Program Files/R/R-3.1.2/library/GenomicAlignments/extdata/sm_untreated1.bam"

reads1 <- readGAlignments(fls[1])

比如我随意截图 cigar(reads1)的结果，就可以看出，大部分都是75M这样的比对情况，可以对这个向量进行统计完全匹配的情况，总共有6077总比对情况。

unique(cigar(reads1))

Bioconductor系列之GenomicAlignments1315

也可以用table格式来查看cigar，可以看到大部分都是M

head(cigarOpTable(cigar(reads1)))

Bioconductor系列之GenomicAlignments1385

接下来我们再读取PE测序数据的比对bam结果，看看分析方法有哪些。

Bioconductor系列之GenomicAlignments1422

U3.galp <- readGAlignmentPairs(untreated3_chr4(), use.names=TRUE, param=param0)head(U3.galp)

也可以分开查看左右两端测序数据的比对情况，跟上面的head是对应关系，上面的

[169,  205] -- [ 326,  362]在下面就被分成左右两个比对情况了。

> head(first(U3.galp))

GAlignments object with 6 alignments and 0 metadata columns:      seqnames strand       cigar    qwidth     start       end     width     njunc         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>  [1]     chr4      +         37M        37       169       205        37         0  [2]     chr4      +         37M        37       943       979        37         0  [3]     chr4      +         37M        37       944       980        37         0  [4]     chr4      +         37M        37       946       982        37         0  [5]     chr4      +         37M        37       966      1002        37         0  [6]     chr4      +         37M        37       966      1002        37         0  -------  seqinfo: 8 sequences from an unspecified genome

> head(last(U3.galp))

GAlignments object with 6 alignments and 0 metadata columns:      seqnames strand       cigar    qwidth     start       end     width     njunc         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>  [1]     chr4      -         37M        37       326       362        37         0  [2]     chr4      -         37M        37      1086      1122        37         0  [3]     chr4      -         37M        37      1119      1155        37         0  [4]     chr4      -         37M        37       986      1022        37         0  [5]     chr4      -         37M        37      1108      1144        37         0  [6]     chr4      -         37M        37      1114      1150        37         0  -------  seqinfo: 8 sequences from an unspecified genome

用isProperPair可以查看pe测序数据比对情况的，哪些没有配对比对成功

> table(isProperPair(U3.galp))FALSE  TRUE 29518 45828

还可以用Rsamtools包对我们的bam文件进行统计，看下面的代码中fls[1]和un1都是上面得到的bam文件全路径。

> library(Rsamtools)

> quickBamFlagSummary(fls[1], main.groups.only=TRUE)

group |    nb of |    nb of | mean / max                                   of |  records |   unique | records per                              records | in group |   QNAMEs | unique QNAMEAll records........................ A |     1800 |     1800 |    1 / 1  o template has single segment.... S |     1800 |     1800 |    1 / 1  o template has multiple segments. M |        0 |        0 |   NA / NA      - first segment.............. F |        0 |        0 |   NA / NA      - last segment............... L |        0 |        0 |   NA / NA      - other segment.............. O |        0 |        0 |   NA / NA Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.Indentation reflects this.

> library(Rsamtools)

> quickBamFlagSummary(un1, main.groups.only=TRUE)

group |    nb of |    nb of | mean / max                                   of |  records |   unique | records per                              records | in group |   QNAMEs | unique QNAMEAll records........................ A |   204355 |   190770 | 1.07 / 10  o template has single segment.... S |   204355 |   190770 | 1.07 / 10  o template has multiple segments. M |        0 |        0 |   NA / NA      - first segment.............. F |        0 |        0 |   NA / NA      - last segment............... L |        0 |        0 |   NA / NA      - other segment.............. O |        0 |        0 |   NA / NA Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.Indentation reflects this.

 

下面讲最后一个综合应用

biocLite("RNAseqData.HNRNPC.bam.chr14")

#这是一个665.5MB的RNA-seq测序数据的比对结果

library(RNAseqData.HNRNPC.bam.chr14)bamfiles <- RNAseqData.HNRNPC.bam.chr14_BAMFILESnames(bamfiles)  # the names of the runslibrary(Rsamtools)quickBamFlagSummary(bamfiles[1], main.groups.only=TRUE)flag1 <- scanBamFlag(isFirstMateRead=TRUE, isSecondMateRead=FALSE,                     isDuplicate=FALSE, isNotPassingQualityControls=FALSE)param1 <- ScanBamParam(flag=flag1, what="seq") #这里是设置readGAlignments的读取参数library(GenomicAlignments)gal1 <- readGAlignments(bamfiles[1], use.names=TRUE, param=param1)

mcols(gal1)$seqoqseq1 <- mcols(gal1)$seqis_on_minus <- as.logical(strand(gal1) == "-")oqseq1[is_on_minus] <- reverseComplement(oqseq1[is_on_minus])

这样就还原得到了原始测序数据啦！

七 19

Bioconductor系列之pasillaBamSubset

Posted on by

这个包主要有bam文件测试数据
> biocLite("pasillaBamSubset")
BioC_mirror: http://bioconductor.orgUsing Bioconductor version 3.0 (BiocInstaller 1.16.5), R version 3.1.2.
Installing package(s) 'pasillaBamSubset'
trying URL 'http://bioconductor.org/packages/3.0/data/experiment/bin/windows/contrib/3.1/pasillaBamSubset_0.3.1.zip'
Content type 'application/zip' length 31514402 bytes (30.1 Mb)
打开pasillaBamSubset包的安装地址就可以看到里面有几个bam文件
Several functions are available for reading BAM files into R:
而且加载包的同时也引入了几个读取bam文件的函数
readGAlignments()
readGAlignmentPairs()
readGAlignmentsList()
scanBam()
Bioconductor系列之pasillaBamSubset448
加载包就可以看到用两个函数得到包自带的数据文件的地址，主要是有很多人不一定把包安装在C盘，所以用函数来定位文件更加安全一点
> library(pasillaBamSubset)
> untreated1_chr4()
[1] "C:/Program Files/R/R-3.1.2/library/pasillaBamSubset/extdata/untreated1_chr4.bam"
> untreated3_chr4()
[1] "C:/Program Files/R/R-3.1.2/library/pasillaBamSubset/extdata/untreated3_chr4.bam"

接下来我们就看看如何读取这些bam文件的
library(pasillaBamSubset)
un1 <- untreated1_chr4() # single-end reads library(GenomicAlignments) reads1 <- readGAlignments(un1) cvg1 <- coverage(reads1) 查看reads1这个结果，可以看到把这个bam文件都读成了一个数据对象GAlignments object， Bioconductor系列之pasillaBamSubset1142
针对着个数据对象有很多操作，其中一个coverage操作是来自于GenomicFeatures
或者GenomicAlignments函数的，可以算出测序覆盖情况。
可以看到这个bam文件里面的比对情况大多几种在4号染色体里面
> cvg1$chr4
integer-Rle of length 1351857 with 122061 runs
Lengths: 891 27 5 12 13 45 5 12 13 ... 5 1 1 3 10
Values : 0 1 2 3 4 5 4 3 2 ... 12 11 10 6
> mean(cvg1$chr4)
[1] 11.33746
> max(cvg1$chr4)[1] 5627
可以看到平均测序深度是11.3X，最大测序深度是5627X

七 18

Bioconductor系列之biomaRt

Posted on by

biomaRt是一个超级网络资源库，里面的信息非常之多，就是网页版的biomaRt的R语言接口。谷歌搜索 the biomart user’s guide filetype:pdf 这个关键词，就看到关于这个包的详细介绍以及例子，我这里简单总结一下它的用法。

这个包一直在更新，函数总是变化，而且需要联网才能使用，基本上等于废物一个，希望大家不要使用它~

包的安装

Continue reading →

七 18

Bioconductor系列之hgu95av2.db

Posted on by

这个包主要都是芯片数据处理方面的应用，本来我是懒得看的，但是我无意中浏览了一下readme，发现它挺适合被作为数据库的典范来学习。
source("http://bioconductor.org/biocLite.R")
biocLite("org.Hs.eg.db")
biocLite("hgu95av2.db")
还是老办法安装了hgu95av2.db之后用ls可以查看这个包里面有36个映射数据，都是把芯片探针ID号转为其它36种主流ID号的映射而已。
library(hgu95av2.db)
> ls("package:hgu95av2.db")
[1] "hgu95av2" "hgu95av2ACCNUM" "hgu95av2ALIAS2PROBE"
[4] "hgu95av2CHR" "hgu95av2CHRLENGTHS" "hgu95av2CHRLOC"
[7] "hgu95av2CHRLOCEND" "hgu95av2.db" "hgu95av2_dbconn"
[10] "hgu95av2_dbfile" "hgu95av2_dbInfo" "hgu95av2_dbschema"
[13] "hgu95av2ENSEMBL" "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"
[16] "hgu95av2ENZYME" "hgu95av2ENZYME2PROBE" "hgu95av2GENENAME"
[19] "hgu95av2GO" "hgu95av2GO2ALLPROBES" "hgu95av2GO2PROBE"
[22] "hgu95av2MAP" "hgu95av2MAPCOUNTS" "hgu95av2OMIM"
[25] "hgu95av2ORGANISM" "hgu95av2ORGPKG" "hgu95av2PATH"
[28] "hgu95av2PATH2PROBE" "hgu95av2PFAM" "hgu95av2PMID"
[31] "hgu95av2PMID2PROBE" "hgu95av2PROSITE" "hgu95av2REFSEQ"
[34] "hgu95av2SYMBOL" "hgu95av2UNIGENE" "hgu95av2UNIPROT"
但是用这个包自带的函数capture.output(hgu95av2())可以看到这些映射并没有囊括我们标准的hg19版本的2.3万个基因，也就说这个芯片设计的探针只有1.1万个左右。但它根据已有的org.Hs.eg.db来构建了它自己芯片独有的数据包，这样就显得更加正式规范化了，这启发我们研发人员在做一些市场产品的同时也可以把自己的研究成果通主流数据库数据形式结合起来，更易让他人接受。
> capture.output(hgu95av2())
[1] "Quality control information for hgu95av2:"
[2] ""
[3] ""
[4] "This package has the following mappings:"
[5] ""
[6] "hgu95av2ACCNUM has 12625 mapped keys (of 12625 keys)"
[7] "hgu95av2ALIAS2PROBE has 33755 mapped keys (of 103735 keys)"
[8] "hgu95av2CHR has 11540 mapped keys (of 12625 keys)"
[9] "hgu95av2CHRLENGTHS has 93 mapped keys (of 93 keys)"
[10] "hgu95av2CHRLOC has 11474 mapped keys (of 12625 keys)"
[11] "hgu95av2CHRLOCEND has 11474 mapped keys (of 12625 keys)"
[12] "hgu95av2ENSEMBL has 11460 mapped keys (of 12625 keys)"
[13] "hgu95av2ENSEMBL2PROBE has 9814 mapped keys (of 28553 keys)"
[14] "hgu95av2ENTREZID has 11543 mapped keys (of 12625 keys)"
[15] "hgu95av2ENZYME has 2125 mapped keys (of 12625 keys)"
[16] "hgu95av2ENZYME2PROBE has 786 mapped keys (of 975 keys)"
[17] "hgu95av2GENENAME has 11543 mapped keys (of 12625 keys)"
[18] "hgu95av2GO has 11245 mapped keys (of 12625 keys)"
[19] "hgu95av2GO2ALLPROBES has 17214 mapped keys (of 18826 keys)"
[20] "hgu95av2GO2PROBE has 12920 mapped keys (of 14714 keys)"
[21] "hgu95av2MAP has 11519 mapped keys (of 12625 keys)"
[22] "hgu95av2OMIM has 10541 mapped keys (of 12625 keys)"
[23] "hgu95av2PATH has 5374 mapped keys (of 12625 keys)"
[24] "hgu95av2PATH2PROBE has 228 mapped keys (of 229 keys)"
[25] "hgu95av2PMID has 11529 mapped keys (of 12625 keys)"
[26] "hgu95av2PMID2PROBE has 372382 mapped keys (of 432400 keys)"
[27] "hgu95av2REFSEQ has 11506 mapped keys (of 12625 keys)"
[28] "hgu95av2SYMBOL has 11543 mapped keys (of 12625 keys)"
[29] "hgu95av2UNIGENE has 11533 mapped keys (of 12625 keys)"
[30] "hgu95av2UNIPROT has 11315 mapped keys (of 12625 keys)"
[31] ""
[32] ""
[33] "Additional Information about this package:"
[34] ""
[35] "DB schema: HUMANCHIP_DB"
[36] "DB schema version: 2.1"
[37] "Organism: Homo sapiens"
[38] "Date for NCBI data: 2014-Sep19"
[39] "Date for GO data: 20140913"
[40] "Date for KEGG data: 2011-Mar15"
[41] "Date for Golden Path data: 2010-Mar22"
[42] "Date for Ensembl data: 2014-Aug6"
这个hgu95av2.db所加载的36个包都是一种特殊的对象，但是可以把它当做list来操作，是一种类似于hash的对应表格，其中keys是独特的，但是value可以有多个。
既然是类似于list，那我就简单讲几个小技巧来操作这些数据对象。所有的操作都要用as.list()函数来把数据展现出来
> as.list(hgu95av2ENZYME[1])
$`1000_at`
[1] "2.7.11.24"
可以看到这样就提取出来了hgu95av2ENZYME的第一个元素，key是`1000_at`，它所映射的value是 "2.7.11.24"这个酶。
同理对list取元素的三个操作在这里都可以用
> as.list(hgu95av2ENZYME['1000_at'])
$`1000_at`
[1] "2.7.11.24"
> as.list(hgu95av2ENZYME$'1000_at')
[[1]]
[1] "2.7.11.24"
然而这只是list的操作，我们还有一个函数专门对这个对象来取元素，那就是get函数，取多个元素用mget，类似于as.list(hgu95av2ENZYME[1:10])
用函数mget(probes_id,hgu95av2ENZYME)就可以根据自己的probes_id这个向量来选取数据了，当然也要用as.list()来展示出来。
值得一提的是这里有个非常重要的函数是revmap，可以把我们的key和value调换位置。
因为我们的数据对象里面的对应关系不是一对一，而是一对多，例如一个基因往往有多个go通路信息，例如这个就有十几个
as.list(hgu95av2GO$'1000_at')
$`GO:0000165`
$`GO:0000165`$GOID
[1] "GO:0000165"

$`GO:0000165`$Evidence
[1] "TAS"

$`GO:0000165`$Ontology
[1] "BP"

$`GO:0000186`
$`GO:0000186`$GOID
[1] "GO:0000186"

$`GO:0000186`$Evidence
[1] "TAS"

$`GO:0000186`$Ontology
[1] "BP"
一层层的数据结构非常清晰，但是数据太多，所以它给了一个toTable函数来格式化，就是把一对多的list压缩成表格。
> toTable(hgu95av2GO[1])
probe_id go_id Evidence Ontology
1 1000_at GO:0000165 TAS BP
2 1000_at GO:0000186 TAS BP
3 1000_at GO:0000187 TAS BP
4 1000_at GO:0002224 TAS BP
5 1000_at GO:0002755 TAS BP
6 1000_at GO:0002756 TAS BP
7 1000_at GO:0006360 TAS BP
------------------------------------------------------------
81 1000_at GO:0004707 NAS MF
82 1000_at GO:0001784 IEA MF
83 1000_at GO:0005524 IEA MF
这样就非常方便的查看啦。
当然对这个数据对象还有很多实用的操作。length(),Llength(),Rlength(),keys(),Lkeys(),Rkeys()等等，还有count.mappedkeys(),mappedLkeys()，还有个比较特殊的isNA来判断是否有些keys探针没有对应任何信息。
而且以上这些函数不仅可以用来获取数据对象的信息，还可以修改这个数据对象。
Lkeys(hgu95av2CHR) 可以查看这个对象有12625个探针。
而> Rkeys(hgu95av2CHR)
[1] "19" "12" "8" "14" "3" "2" "17" "16" "9" "X" "6" "1" "7" "10" "11"
[16] "22" "5" "18" "15" "Y" "20" "21" "4" "13" "MT" "Un"
可以看到这些探针分布在不同的染色体上面，如果我这时候给
x=hgu95av2CHR
Rkeys(x) = c("1","2")
toTable(x)可以看到这时候x只剩下1946个探针了，就是1,2号染色体上面的基因了。
然后可以一个个看这些函数的用法，其实就是SQL的增删改查的基本操作而已。
> length(x)
[1] 12625
> Llength(x)
[1] 12625
> Rlength(x)
[1] 2
> head(keys(x))
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at" "1005_at"
> head(Lkeys(x))
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at" "1005_at"
> head(Rkeys(x))
[1] "1" "2"
> count
count.fields count.mappedLkeys countOverlaps counts<- count.links count.mappedRkeys countQueryHits countSubjectHits count.mappedkeys countMatches counts > count.mappedkeys(x)
[1] 1946
> count.mappedLkeys(x)
[1] 1946
> count.mappedRkeys(x)
[1] 2

五 28

org.Xx.eg.db系列包概述

Posted on by

在bioconductor的官网里面可以查到共有111个系列包，基本上跨越了我们常见的物种啦！

org.Xx.eg.db系列包介绍65

斑马鱼：Bioconductor - org.Dr.eg.db - /packages/release/data/annotation/html/org.Dr.eg.db.html

Details biocViews AnnotationData , Danio_rerio , OrgDb Version 3

拟南芥：Bioconductor - org.At.tair.db - /packages/release/data/annotation/html/org.At.tair.db.html

Details biocViews AnnotationData , Arabidopsis_thaliana , OrgDb Version 3

小鼠：Bioconductor - org.Mm.eg.db - /packages/release/data/annotation/html/org.Mm.eg.db.html

Details biocViews AnnotationData , Mus_musculus , OrgDb , mouseLLMappings Version 3

人类：Bioconductor - org.Hs.eg.db - /packages/release/data/annotation/html/org.Hs.eg.db.html

Details biocViews AnnotationData , Homo_sapiens , OrgDb , humanLLMappings Version 3

对这些系列包的函数都一样，包括以下几个：

columns(x)  keytypes(x)  keys(x, keytype, ...)  select(x, keys, columns, keytype, ...)  saveDb(x, file)  loadDb(file, dbType, dbPackage, ...)

 

 

这些包就是bioconductor已经做好的数据库，我们可以根据定义好的ID号来进行任意的基因转换，现在支持的信息有一下几种！

keytypes(org.Hs.eg.db)

[1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"      "ACCNUM"       "ALIAS"        "CHR"

[8] "CHRLOC"       "CHRLOCEND"    "ENZYME"       "MAP"          "PATH"         "PMID"         "REFSEQ"

[15] "SYMBOL"       "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "GENENAME"     "UNIPROT"

[22] "GO"           "EVIDENCE"     "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"  "OMIM"

[29] "UCSCKG"

这些包的确非常有用，大家可以看我博客里面关于它们的介绍！！！

 

五 28

菜鸟团第二次作业的部分答案

Posted on by

> library(org.Hs.eg.db)

载入需要的程辑包：AnnotationDbi载入需要的程辑包：stats4载入需要的程辑包：GenomeInfoDb载入需要的程辑包：S4Vectors载入需要的程辑包：IRanges载入程辑包：‘AnnotationDbi’The following object is masked from ‘package:GenomeInfoDb’:     species载入需要的程辑包：DBI

 

1、人共有多少个entrez id的基因呢？

x <- org.Hs.egENSEMBLTRANS

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(x)

[1] 47721

可知共有47721个基因都是有entrez ID号的

2、能对应转录本ID的基因有多少个呢？

length(xx)

[1] 20592

可以看到共有20592个基因都是有转录本的！

2、能对应ensembl的gene ID的基因有多少个呢？

x <- org.Hs.egENSEMBL

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 26019

可以看到只有26019是有ensembl的gene ID的

3、那么基因对应的转录本分布情况如何呢？

table(unlist(lapply(xx,length)))

菜鸟团第二次作业的部分答案863

可以看出绝大部分的基因都是20个转录本一下的，但也有极个别基因居然有高达两百个转录本，很可怕！

4、那么基因在染色体的分布情况如何呢？

x <- org.Hs.egCHR

# Get the entrez gene identifiers that are mapped to a chromosome

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 47232

可以看到有接近五百个基因居然是没有染色体定位信息的！！！

table(unlist(xx))

用barplot函数可视化一下，如图

 

菜鸟团第二次作业的部分答案1209

6、那么有多多少基因是有GO注释的呢？

x <- org.Hs.egGO

# Get the entrez gene identifiers that are mapped to a GO ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(xx)

[1] 18229

> length(x)

[1] 47721

可以看到只有18229个基因是有go注释信息的。

那么基因被注释的go的分布如何呢？

菜鸟团第二次作业的部分答案1477

可以看到大部分的基因都是只有30个go的，但是某些基因特别活跃，高达197个go注释。

还有kegg和omin数据库的我就不写了！

五 28

实战R语言bioconductor的seqinr包探究人的所有转录本的性质

Posted on by

首先安装这个包

source("http://bioconductor.org/biocLite.R")

biocLite("seqinr")

然后加载包，并读取我们的CDS.fa文件

library("seqinr")

human_cds=read.fasta("CDS.fa")

#这一个步骤非常耗时间，可能是因为我们的转录本文件有十万多个转录本的原因吧

str(human_cds) #查看可知读入了一个list，其中每个转录本都是list的一个元素

List of 100778

$ ENST00000415118:Class 'SeqFastadna'  atomic [1:8] g a a a ...

.. ..- attr(*, "name")= chr "ENST00000415118"

.. ..- attr(*, "Annot")= chr ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene tran"| __truncated__

$ ENST00000448914:Class 'SeqFastadna'  atomic [1:13] a c t g ...

.. ..- attr(*, "name")= chr "ENST00000448914"

.. ..- attr(*, "Annot")= chr ">ENST00000448914 havana_ig_gene:known chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene tran"| __truncated__

对list的每个元素都有几种函数可以处理得到信息：

Length,table,GC,count

其中count函数很有趣，数一数序列里面的这些组合出现的次数

count(dengueseq, 1)

count(dengueseq, 2)接下来我们随机取human_cds这个list的一个元素用这几个函数对它处理一下

> tmp=human_cds[[1]]

> tmp

[1] "g" "a" "a" "a" "t" "a" "g" "t"

attr(,"name")

[1] "ENST00000415118"

attr(,"Annot")

[1] ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene"

attr(,"class")

[1] "SeqFastadna"

再看看函数的结果

> length(tmp)

[1] 8

> table(tmp)

tmp

a g t

4 2 2

> GC(tmp)

[1] 0.25

> count(tmp,1)

 

a c g t

4 0 2 2

> count(tmp,2)

 

aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt

2  0  1  1  0  0  0  0  1  0  0  1  1  0  0  0

>

还是挺好用的，接下来我们应用R的知识来对着十万多个转录本进行一些简单的总结

human_cds_length=unlist(lapply(human_cds,length))

human_cds_gc=unlist(lapply(human_cds,GC))

这样就得到了所有转录本的长度和GC含量信息

然后我们简单统计一下，并画几个图表吧！

> summary(human_cds_length)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

3     366     699    1132    1425  108000

> summary(human_cds_gc)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

0.1467  0.4577  0.5285  0.5264  0.5932  0.8917

可以看到还是有很多很极端的转录本的存在！

最长的转录本也不过10k，而我记得最长的基因高达8M，看了内含子远大于外显子呀。

但是GC含量有很多高于80%，这些基因在二代测序的研究中是一个盲区。

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2075

这些极端基因可以通过biomaRt等包进行注释，得到gene名和功能信息。

 

hist(human_cds_gc)

hist(log10(human_cds_length))

GC含量分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2177

长度分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2186

附表：

http://www.bioinformatics.org/sms/iupac.html 所有字符的碱基氨基酸意义表格

 

五 21

Biostrings包简介

Posted on by

首先讲讲它的对象

有下面几个字符串对象BString, DNAString, RNAString and AAString可以通过以下代码构造它们：

b <- BString("I am a BString object")

d <- DNAString("TTGAAAA-CTC-N")

这两个对象的区别是DNAstring对象对字符串的要求严格一些，只有IUPAC字符和+-字符可以。

对构造好的对象可以通过下标来取子字符串对象，也可以通过subseq来取，但是子字符串仍然是数据对象，只有通过toString函数才能把它们转化成字符串。

用length(dd2)和nchar(toString(dd2))都可以找到我们Biostrings对象的长度。但是后者速度会很慢。

Views(RNAString("AU"), start=0, end=2)这个函数可以把string对象任意截取成list

start, end and width可以作用于我们截取的list，判断list里面的元素在原来的string对象上面的起始终止及长度信息。

 

接下来讲这个包带有的一个比对函数！

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] succ--eed subject: [1] supersede score: -33.99738

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",type = "local")

Local PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su subject: [1] su score: 5.578203

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",gapOpening = 0, gapExtension = 1)

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su-cce--ed subject: [1] sup--

可以看出这个比对函数可以调整的参数实在是太多了，而且改变参数之后比对情况大不一样，还有很多参数就不一一细讲了。

这个比对结果可以赋值给一个变量，保存比对的对象。

psa1 <- pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

class(psa1)

summary(psa1)

class(pattern(psa1))

class(summary(psa1))

score(psa2)

还可以自己构建打分矩阵来进行比对。

submat <-

+ matrix(-1, nrow = 26, ncol = 26, dimnames = list(letters, letters))

diag(submat) <- 0

Biostrings包简介1454

psa2 <-pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",substitutionMatrix = submat, gapOpening = 0, gapExtension = 1)

我们的包还自带了两个非常流行的氨基酸比对矩阵PAM和BLOSUM

ls("package:Biostrings")可以查看这个包所有的对象。

data(package="Biostrings")可以查看这个包所有的数据对象

还有很多其它函数

还可以去除adaptor，挺好玩的

既然有配对比对函数，那么就有多重比对函数！

我们可以读取clustaW, Phylip and Stolkholm这几种不同的比对结果文件来构造多重比对对象。

library(Biostrings)这个包里面自带了两个文件，我们可以示范一下构建对象。

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata", "msx2_mRNA.aln", package="Biostrings"), format="clustal")

phylipMAlign <- readAAMultipleAlignment(filepath = system.file("extdata","Phylip.txt", package="Biostrings"),format="phylip")

 

对构造好的多重比对对象就可以构建进化树啦，代码如下！

sdist <- stringDist(as(origMAlign,"DNAStringSet"), method="hamming")

> clust <- hclust(sdist, method = "single")

> pdf(file="badTree.pdf")

> plot(clust)

> dev.off()

Biostrings包简介2345

五 05

Bioconductor的数据包library(biomaRt)简介

Posted on by

 

这是发布在bioconductor平台上面的一个数据库文件，可以通过R里面下载安装并使用，非常方便。其实在ensembl数据库里面也有一个biomart，我之前也讲过这个平台，非常好用，可以把任意的数据库之间的ID号进行转换。

为了更好的理解和掌握biomaRt，我们可以先通过在线资源来了解一下它的原型biomart （http://www.biomart.org）。 biomart是为生物科研提供数据服务的免费软件，它为数据下载提供打包方案。它有许多成功的应用实例，比如欧洲生物信息学中心（The European Bioinformatics Institute ，EBI）维护的Ensembl数据库（http://www.ensembl.org/）就使用biomart提供数据批量下载服务， 还有COSMIC, Uniprot, HGNC, Gramene, Wormbase以及dbSNP等。

这个就是一个R平台的biomart而已，但是非常好用！

> library(biomaRt)

> head(listMarts(), 3)

biomart                           version

1    ensembl      ENSEMBL GENES 79 (SANGER UK)

2        snp  ENSEMBL VARIATION 79 (SANGER UK)

3 regulation ENSEMBL REGULATION 79 (SANGER UK)

这是这个biomart最具有代表性的三个数据库，用listMarts()可以查看得知，它总共有58个数据库。

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这是创建了人的ensembl数据库对象

> head(listFilters(ensembl), 3)

name     description

1 chromosome_name Chromosome name

2           start Gene Start (bp)

3             end   Gene End (bp)

可以看到对人的数据库ensembl来说，有多种字段可以来挑选自己感兴趣的东西，最常用的的当然是染色体号及起始终止坐标啦，用listFilters(ensembl),以查看得知，它总共有284中挑选感兴趣数据的方式。

既然 chromosome_name是其中一个挑选字段，那么我们就可以看看，是如何进行挑选的

用filterOptions(myFilter, ensembl)可以看到它挑选参数非常之多，远不止我们所认为的染色体号码。

染色体号一般是1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y,MT

还有一堆稀奇古怪的标志，LRG_101,LRG_102,LRG_103,LRG_104，因为我们组装好的人的标准基因组还有很多小的片段不被计入染色体中。

然后还可以看到人的ensembl数据库对象，有很多的属性，最常见的当然是基因ID和转录本ID和蛋白的ID号啦！

> head(listAttributes(ensembl), 3)

name           description

1       ensembl_gene_id       Ensembl Gene ID

2 ensembl_transcript_id Ensembl Transcript ID

3    ensembl_peptide_id    Ensembl Protein ID

用listAttributes(ensembl),,以查看得知，它总共有1166个ID号，太恐怖了，我实在是没有想到！

 

那么接下来我简单讲讲这个包的几个应用吧

首先是根据entrez ID号来找

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这样就得到了人的信息，然后我探究以下两个基因的其它信息。

entrzID=c("672","1")

getBM(attributes=c("entrezgene","hgnc_symbol","ensembl_gene_id"), filters = "entrezgene", values =entrzID, mart=ensembl )

entrezgene hgnc_symbol ensembl_gene_id

1          1        A1BG ENSG00000121410

2        672       BRCA1 ENSG00000012048

3        672       BRCA1         LRG_292

 

其实这个函数很简单，就是根据自定义的entrzID这个变量来找到一些数据，数据的属性是我自己定义的entrezgene","hgnc_symbol","ensembl_gene_id"，所以它就显示这个信息给我，在我之前弄好的人的数据库里面寻找！listAttributes(ensembl),,以查看得知，它总共有1166个ID号，就是说，你可以挑选你想要的基因的1166种信息，包罗万象！！！

其它功能也是很简单的啦，自己多看帮助文档！

 

从上面的操作来看，使用biomaRt只需要两步，1，指定mart数据库，2，使用getBM获得注释。但是首先，我们如何知道有哪些服务器，以及这些服务器上哪些数据库呢？其次，我们如何获阳getBM中attributes，filters的正确设置呢？

关于第一个问题，我们可以使用biomaRt中的listMarts以及listDatasets两个函数来解决。

> marts <- listMarts(); head(marts) #查看当前可用的数据源 ，总共有58个数据源。

> ensembl <- useMart("ensembl") #使用ensembl数据源

> datasets <- listDatasets(ensembl); datasets[1:10,] #查看ensembl中可用数据库，共有69个物种的数据库！

对于第二个问题，我们使用biomaRt中的listFilters以及listAttributes两个函数来解决。

> mart <- useMart("ensembl", "hsapiens_gene_ensembl")  #首先使用人的数据库

>listAttributes(ensembl) #,以查看得知，它总共有1166个ID号，就是说，人的数据库可供挑选的信息多达1166种。

> filters <- listFilters(mart); filters[grepl("entrez", filters[,1]),] #总共有284中挑选感兴趣数据的方式。

最后的问题是，biomaRt会被如何使用呢？我们做注释的时候，怎么就想到要使用biomaRt呢？因为在注释上，各种ID，symbol, name之间的转换都可以考虑使用biomaRt来做。更重要的是，biomaRt还会有很多SNP, alternative splicing, exon, intron, 5’utr, 3’utr等等信息。当然，只要能做也数据库并使用SQL访问的数据都可以使用biomaRt来获取。所以我们的思路可以更加发散一些。

 

 

 

五 05

Bioconductor的数据包library(org.Hs.eg.db)简介

Posted on by

 

这是发布在bioconductor平台上面的一个数据库文件，可以通过R里面下载安装并使用，非常方便。而且用的是数据库存储方式，所以搜索起来也是非常快速。

这个包里面有28个主流数据资料文件，这样我们可以用select函数根据我们自己的ID在这28个数据库里面随意转换自己想要的信息！！！

当然我本人是比较喜欢直接下载原文件，然后写脚本自己进行各种数据直接的转换。

首先我们加载这个数据包，可以看到这个数据包依赖于很多其它的包，如果是第一次安装。会耗时很长！

Bioconductor的数据包org.Hs.eg.db269

用这个函数，可以看到这个org.Hs.eg.db数据对象里面包含着各大主流数据库的数据，一般人都比较熟悉的entrez ID 和ensembl 数据库的ID。

keytypes(org.Hs.eg.db)

##  [1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"

##  [5] "ACCNUM"       "ALIAS"        "ENZYME"       "MAP"

##  [9] "PATH"         "PMID"         "REFSEQ"       "SYMBOL"

##  [13] "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"

##  [17] "GENENAME"     "UNIPROT"      "GO"           "EVIDENCE"

##  [21] "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"

##  [25] "OMIM"         "UCSCKG"

然后，我们用select函数，就可以把任意公共数据库的数据进行一一对应了。

ensids <- c("ENSG00000130720", "ENSG00000103257", "ENSG00000156414",

"ENSG00000144644", "ENSG00000159307", "ENSG00000144485")

cols <- c("SYMBOL", "GENENAME")

select(org.Hs.eg.db, keys=ensids, columns=cols, keytype="ENSEMBL")

比如说，我们有几个ensembl的基因ID号。然后我们想找它所对应的gene名和缩略词简称，就通过select函数来搞定即可！

Bioconductor的数据包org.Hs.eg.db1158

select(org.Hs.eg.db, keys="BRCA1", columns=c("ENSEMBL","UNIGENE","ENTREZID","CHR","GO","GENENAME"), keytype="SYMBOL")

这样得到了这个BRCA1基因的大部分信息，只是它的GO条目太多了，看得有点乱。

Bioconductor的数据包org.Hs.eg.db1318

 

 

 

五 05

Bioconductor简介

Posted on by

主页：http://www.bioconductor.org/

文字介绍我懒得写了，具体大家参考

http://www.bioconductor.org/about/

http://blog.csdn.net/shmilyringpull/article/details/8542607

这是一个R语言进行生信分析的流程发布平台，每个包都解决生信的一个流程问题。到目前为止2015年5月5日10:57:29已经有了1024个包，所以大家可以看到生信分析是一个海量的任务了。每一个包都有着详尽的说明文档及脚本代码，还附带着数据，非常容易弄懂，接下来我会花一个月的时间好好学习这些包！

这1024个虽然还是R语言的包，但是它的安装方式与常规的R语言包已经有所区别了。

需要用一下代码来安装

source("http://bioconductor.org/biocLite.R")biocLite()

biocLite(c("GenomicFeatures", "AnnotationDbi"))

也是非常easy的。

现在这个平台上面有1024个包，241个实验数据，917个数据库文件！！！

We are pleased to announce Bioconductor 3.1,

consisting of 1024 software packages,

241 experiment data packages,

and 917 up-to-date annotation packages.

在MOOC上面有很多关于这个的公开课

http://bioconductor.org/help/course-materials/

 

这里面有很多生信方向的分析流程，包括了我之前提到了snp-calling，RNA-seq，CHIP-seq等等，当然最主要的还是芯片数据的处理。

Workflows »

Common Bioconductor workflows include:

这些流程基本上涉及到了现在生物信息学的主流方向，所以基本上掌握了这些包，就是一个合格的生物信息学人才啦！

更重要的是它有着917个数据库文件，里面的信息分门别类，几乎可以算作是生物信息学的百科全书啦！

主要的数据库包括以下。

 

Package Description
AnnotationHub Ensembl, Encode, dbSNP, UCSC data objects
biomaRt Ensembl and other annotations
PSICQUIC Protein interactions
uniprot.ws Protein annotations
KEGGREST KEGG pathways
SRAdb Sequencing experiments.
rtracklayer genome tracks.
GEOquery Array and other data
ArrayExpress Array and other data

 

 

 

 

 

 

 

 

 

 

 

 

 

五 04

topGO简单使用

Posted on by

首先载入这个包

source("http://bioconductor.org/biocLite.R")

biocLite("topGO")

biocLite("ALL")

library(topGO)

library(ALL)

data(ALL)

data(geneList)

data(GOdatat)

这样就载入了很多变量， ls()查看如下

[1] "affyLib"      "ALL"          "geneList"     "topDiffGenes"

其中ALL这个数据我在另一篇日志里面重点介绍了一下。

然后我简单提一下"geneList"

head(geneList)

1095_s_at   1130_at   1196_at 1329_s_at 1340_s_at 1342_g_at

1.0000000 1.0000000 0.6223795 0.5412240 1.0000000 1.0000000

str(geneList) 是一个向量，有323个数字。

Named num [1:323] 1 1 0.622 0.541 1 ...

- attr(*, "names")= chr [1:323] "1095_s_at" "1130_at" "1196_at" "1329_s_at" ...

然后简单查询该包的安装地址和一些文件

system.file(package = "topGO")

[1] "C:/Program Files/R/R-3.1.1/library/topGO"

在这个目录下面可以找到文件"examples/geneid2go.map"

里面的内容格式如下，第一列是gene的ID号，一般是entrez ID ，第二列是该基因所对应的GO所有的条目，用逗号分隔。

068724 GO:0005488, GO:0003774, GO:0001539, GO:0006935, GO:0009288

119608 GO:0005634, GO:0030528, GO:0006355, GO:0045449, GO:0003677, GO:0007275

此处省略一万行。

readMappings(file = system.file("examples/geneid2go.map", package = "topGO"))

这个函数可以读取我们的文件，返回一个list。是gene到go的映射，每个基因都有一个或者多个go条目。

这个list可以用inverseList这个函数反转，变成每个go条目到基因的映射。

构建topGO这个大全，需要的数据包括：

  1. 基因identifier，可以附上某种分数以便后面施用某种统计处理，分数可以是t检验的p值或者与某个表型的correlation等；
  2. identifier和GO term间的map，如果是芯片数据的话BioC里有多种注释包，声明包的名称即可。至于我等蛋白界苦人，也能自己构建map，见下；
  3. GO的层级结构，由GO.db提供，目前这个包只支持GO.db提供的结构：GOslim就再说了。

topGOdata对象构建函数的参数包括：

  1. ontology，可指定要分析的GO term的类型，即BP、CC之类；
  2. description：topGOdata的描述，可选；
  3. allGenes：基因identifier的原始列表，和后面的geneSelectionFun联合作用，得出参与分析的基因，可以是numeric或factor。
  4. geneSelectionFun：基因选择函数，如果前面allGenes是numeric的话就必须得指明此参数；
  5. nodeSize：被认为富集的GO term辖下基因的最小数目（>=），默认为1。
  6. annotationFun：基因identifier map到GO term的函数。

代码如下

BPterms <- ls(GOBPTerm)

geneID2GO=readMappings(file = system.file("examples/geneid2go.map", package = "topGO"))

直接使用系统自带的data(GOdata)数据，自己构建太麻烦了！

其实就是这就对ALL这个数据集来进行分析！！！

GOdata包含实例topGOdata对象。它可以用来直接运行富集分析。

topGOdata对象构建好后，即可利用这个包里的各种方法和函数做分析。

numGenes(GOdata) 查看对象包含的基因的数目

[1] 318

> description(GOdata)

[1] "Simple topGOdata object"

genes(GOdata)可以得到该对象里面所有的318个基因

geneScore() 可以得到想318个基因的分数

函数sigGenes()返回一个character vector，为各显著变化基因identifier。函数numSigGenes()则用于查看显著变化基因的数目。

函数updateGenes()可以修改topGOdata对象里所包含的基因。

想要看全部基因都对应上了哪些GO term，可用函数usedGO()得到一个character

 

基因集富集分析（gene set enrichment analysis）

首先看看GSEA的三种方式：

  1. 基于count，即仅要求输入一组基因，此种方式最为流行，可用Fisher's exact test、Hypegeometric  test和binomial test进行检验；
  2. 基于基因的score或rank，可用Kolmogorov-Smirnov like tests（即05年那篇PNAS的GSEA文章所用方法），Gentleman's Category、t-test等方法；
  3. 基于基因的表达，可从expression matrix直接分析，如Goeman's globaltest，以及GlobalAncova。

topGO提供两种分析方法，一种自由度更高但上手不易，本菜鸟还是跟着第二种走吧，较为用户友好但集成度较高。简单来说，就是用runTest()这个函数，要求三个主要的argument，一个是之前构建好的topGOdata类实例，第二个参数algorithm用于指定生成GO graph的方法，而参数statistic用于指定所用的检验方法，比如：

> resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

> resultWeight <- runTest(GOdata, algorithm = "weight", statistic = "fisher")

> resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")

> resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks.elim")

runTest这一锤子买卖敲定后就能开始解读和展示结果了，得到的结果是topGOresult类的一个实例，其组成很简单，为对象的基本信息，以及各基因的分数（p值或者其他统计参数

 

 

我这里随便挑一个富集结果来看看

resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

 

-- Classic Algorithm --

 

the algorithm is scoring 590 nontrivial nodes

parameters:

test statistic:  fisher

 

resultWeight <- runTest(GOdata, algorithm = "weight", statistic = "fisher")

 

-- Weight Algorithm --

 

The algorithm is scoring 590 nontrivial nodes

parameters:

test statistic:  fisher : ratio

然后我们对这两种富集方式来画图

pvalFis=score(resultFis) 得到矫正的P值

pvalWeight <- score(resultWeight , whichGO = names(pvalFis))

返回resultFis和resultWeight共有的基因在resultWeight中的分数。有了这两组分数，可以做一些比较，比如关联分析：

cor(pvalFis, pvalWeight)

[1] 0.370151

library(lattice)

xyplot(pvalWeight ~ pvalFis) 画了一个散点图

 

五 04

R语言里面的一个数据集ALL(Acute Lymphoblastic Leukemia)简介

Posted on by

这个数据内容太多了，我感觉自己也理解的不是很清楚！

非常多的R的bioconductor包都是拿这个数据集来举例子的，所以我简单的介绍一下这个数据集。

这个数据集是对ALL这个病的研究数据，共涉及到了128个ALL病人，其中95个是B细胞的ALL，剩余33个是T细胞的ALL。

是一个芯片数据，同时还包含有其它的病人信息。

大家首先要在R里面安装这个数据集

source("http://bioconductor.org/biocLite.R")

biocLite("ALL")

library(ALL)

data(ALL)

data(geneList)

在R里面输入str(ALL)可以看到这个数据的具体信息，但是非常多!

ALL

ExpressionSet (storageMode: lockedEnvironment)

assayData: 12625 features, 128 samples 

element names: exprs

protocolData: none

phenoData

sampleNames: 01005 01010 ... LAL4 (128 total)

varLabels: cod diagnosis ... date last seen (21 total)

varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

 pubMedIds: 14684422 16243790 

Annotation: hgu95av2

我们首先它的BT变量记录的是什么

R语言里面的一个数据集ALL750

可以看出它记录的是数据病人的分组信息。

bcell = grep("^B", as.character(ALL$BT))通过这句话可以挑选出B细胞病人

然后我们看看它的ALL$mol.biol变量记录是是什么

R语言里面的一个数据集ALL857

可以看到它记录的是这些病人的几种突变情况（molecular biology testing）

types = c("NEG", "BCR/ABL")

moltyp = which(as.character(ALL$mol.biol) %in% types)

用这个命令就能挑选出我们想研究的两组突变的病人。

然后我们还可以把刚才的两个标准用来从ALL数据集里面取想要的子集

ALL_bcrneg = ALL[, intersect(bcell, moltyp)]

 

 

同理我们可以查看这个数据集的非常多的变量信息。

包括sex，age，cod，diagnosis，等等，这个'data.frame':共有128 obs. of  21 variables:

R语言里面的一个数据集ALL1190

 

我们除了可以查看这个ALL数据集自带的变量，还可以通过一些方法来访问它的其它信息。

Exprs这个方法可以查看它的表达数据，可以看到有128个变量，12625行的探针数据。

str(exprs(ALL))

num [1:12625, 1:128] 7.6 5.05 3.9 5.9 5.93 ...

- attr(*, "dimnames")=List of 2

..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...

..$ : chr [1:128] "01005" "01010" "03002" "04006" ...

 

还有很多很多函数都可以操作这个数据集，这样可以得到非常多的信息！我就不一一列举了

对这个数据的一系列操作可以画热图，见下面的教程！！！

http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/heatmap/