28

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

在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

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

> 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包探究人的所有转录本的性质

首先安装这个包

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 所有字符的碱基氨基酸意义表格

 

25

用R画GO注释二级分类统计图

群里有朋友问这个图怎么画,我想了想,这肯定是ggplot完成的,非常简单,但是菜鸟们缺乏实践,可能会困惑,所以我模拟数据画了一个!

图片1

首先构造数据

dat=data.frame(name=LETTERS[1:21],
 number=abs(rnorm(21)*10),
 type=c(rep("BP",7),rep("CC",7),rep("MF",7))
)
# 请务必自己查看dat是一个什么数据,print出来即可
# 然后对这个数据画图,一行代码即可!!!
library(ggplot2)
ggplot(dat,aes(x=name,y=number,fill=type))+geom_bar(stat="identity")+coord_flip()

看起来是不是很像回事啦!细节我就懒得调控啦!

图片2

其实自己搜索即可!坐标轴和主题都是可以控制的

http://rstudio-pubs-static.s3.amazonaws.com/3364_d1a578f521174152b46b19d0c83cbe7e.html

http://docs.ggplot2.org/0.9.3.1/coord_flip.html

21

Biostrings包简介

首先讲讲它的对象

有下面几个字符串对象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

21

Bioconductor的DO.db包介绍

Bioconductor的包都是同样的安装方法:

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

还有GO.bd包是完全一模一样的规则!!!

加载这个包可以发现它依赖于好几个其它的包,这也是我比较喜欢R的原因,它会自动把它需要的包全部安装加载进来,不需要自己一个个调试!

> library(DO.db)

载入需要的程辑包:AnnotationDbi

载入需要的程辑包:stats4

载入需要的程辑包:GenomeInfoDb

载入需要的程辑包:S4Vectors

载入需要的程辑包:IRanges

载入需要的程辑包:DBI

> help(DO.db)

> ls("package:DO.db")

[1] "DO"          "DO_dbconn"   "DO_dbfile"   "DO_dbInfo"   "DO_dbschema" "DOANCESTOR"  "DOCHILDREN"  "DOID"        "DOMAPCOUNTS"

[10] "DOOBSOLETE"  "DOOFFSPRING" "DOPARENTS"   "DOSYNONYM"   "DOTERM"      "DOTerms"     "Secondary"   "show"        "Synonym"

[19] "Term"

这个包里面有19个数据对象!都是比较高级的S4对象。

比如我们可以拿DOTERM[1:10]这个小的数据对象来做例子!example=DOTERM[1:10]

因为example是一个高级对象,所以无法直接查看,需要用as.list方法来查看

> as.list(example)

$`DOID:0001816`DOID: DOID:0001816Term: angiosarcomaSynonym: DOID:267Synonym: DOID:4508Synonym: "hemangiosarcoma" EXACT []Secondary: DOID:267Secondary: DOID:4508

~~~~~~~~~~~~共十个DO条目

对每一个DO条目来说都有DOID,Term,Synony这些函数可以取对应的值。

下面是对DO的有向无环图的数据解读

xx <- as.list(DOANCESTOR)可以查看每个DO与它所对应的上级条目DO,每个DO都会有不止一个的上级DO。

xx <- as.list(DOPARENTS)可以查看每个DO与它所对应的父条目DO,每个DO都有且只有一个父DO。

xx <- as.list(DOOFFSPRING)可以查看每个DO与它所对应的下级DO的关系列表,大多数DO都不止一个子条目DO,所有的下级DO都会列出。

xx <- as.list(DOCHILDREN)以查看每个DO与它所对应的子条目DO的关系列表,大多数DO都不止一个子条目DO。

还有Lkeys(DOTERM)可以查看数据库里面的所有的DO条目的ID号

> head(keys(DOTERM))

[1] "DOID:0000000" "DOID:0001816" "DOID:0002116" "DOID:0014667" "DOID:0050004" "DOID:0050012"

dbmeta(GO_dbconn(), "GOSOURCEDATE")

可以查看这个DO库的制备时间

> dbmeta(DO_dbconn(), "DOSOURCEDATE")

[1] "20140417"

05

Bioconductor的数据包library(biomaRt)简介

 

这是发布在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)简介

 

这是发布在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简介

主页: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

 

 

 

 

 

 

 

 

 

 

 

 

 

05

国外最出名的R语言大会-useR

这是2014年的会议报告以及ppt,但是好像很多ppt都是需要翻墙才能下载

http://user2014.stat.ucla.edu/#tutorials

Morning Tutorials Monday, 9:15

Room Presenter Title
Palisades Salon A+B Max Kuhn Applied Predictive Modeling in R
Palisades Salon C+F Winston Chang Interactive graphics with ggvis
Palisades Salon D+E Yihui Xie Dynamic Documents with R and knitr [Slides] [Examples]
Hermosa Romain Francois C++ and Rcpp11 for beginners [slides]
Venice Bob Muenchen Managing Data with R
Sproul-Landing building, 3rd floor Matt Dowle Introduction to data.table [Tutorial] [Talk]
Sproul-Landing building, 4th floor Virgilio Gomez Rubio Applied Spatial Data Analysis with R
Sproul-Landing building, 5th floor Martin Morgan Bioconductor

Afternoon Tutorials Monday, 14:00

Room Presenter Title
Palisades Salon A+B Hadley Wickham Data manipulation with dplyr
Palisades Salon C+F Garrett Grolemund Interactive data display with Shiny and R
Palisades Salon D+E Drew Schmidt Programming with Big Data in R
Hermosa S繪ren H繪jsgaard Graphical Models and Bayesian Networks with R
Venice John Nash Nonlinear parameter optimization and modeling in R [slides]
Sproul-Landing building, 3rd floor Dirk Eddelbuettel An Example-Driven Hands-on Introduction to Rcpp [slides]
Sproul-Landing building, 4th floor Ramnath Vaidyanathan Interactive Documents with R
Sproul-Landing building, 5th floor Thomas Petzoldt Simulating differential equation models in R

 

然后2015年的也要开始了,有兴趣的朋友可以June 30 - July 3, 2015
Aalborg, Denmark看看,有很多干货分享!

http://user2015.math.aau.dk/#BN

2015的内容如下

 

04

topGO简单使用

首先载入这个包

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)简介

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

非常多的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/

 

01

Perl及R及python模块碎碎念

老实说,模块其实是一个很讨厌的东西,但是它也实实在在的节省了我们很多时间,也符合我的理念:避免重复造轮子!此教程可能过期了,请直接看最新版(perl模块安装大全)

1,perl的那些模块

如果有root权限,用root权限

进入cpan然后install ExtUtils::Installed模块

这样就可以执行instmodsh这个脚本了,可以查看当前环境下所有的模块 Continue reading

01

R的包(package)

关于R语言包的一些操作,挺重要的!!!

R的包(package)通常有两种:
1 binary package:这种包属于即得即用型(ready-to-use),但是依赖与平台,即Win和Linux平台下不同。
2 Source package: 此类包可以跨平台使用,但用之前需要处理或者编译(compiled)。

以下一些常用的包相关的函数:
.libPaths():查看包的安装目录

ls('package:ggplot2')可以查看该包里面所有的函数
library():查看已经安装的包目录
library(mypackage):载入mypackage包

getOption("defaultPackages"):查看启动R时自动载入的包。
help(package = 'mypackage'):查看‘mypackage’的帮助
args(function):查看函数的参数
example(function):自动运行该函数帮助文档中的例子,很赞!
demo("package"):展示一些包中demostration,需要再看下??
vignette('mypackage'):有的包,特别是bioconductor的包有vignette,用函数查看
openVignette('mypackage'):这个函数也可以查看vignette,更好用一些
RSiteSearch("helpinfor"):搜索R网站上的“helpinfor”相关信息
help.start():查看已经安装包的详细HTML文档,这个命令非常爽
更新:
search():查看当前载入的包

sessionInfo():查看R中载入的包
methods():查看某个S3泛型函数中所有的方法或者一个类中所有的方法(S3:S version 3)

showMethods(class = "myClass"):查看S4类的方法

findMethods("myMethods"):查看method的代码

class(myObject):查看某个对象的类
getClass(“class/package”):查看某个class或者包的具体内容

getSlots("class"):查看某个class的slot

slotNames(MyObject):查看某个对象的slot。

可以使用Myobject@slotNames访问对象的slot值,这个@设计实在是太爽了,可以连续用。
查询包内信息:1. ?function/method:查看某个“函数”或者“方法”的详细内容
2. class?graph::graph:查看“组”的详细内容的一个例子。这个例子的来源是查询graph包时候,查看其中class的信息,输入??graph后出现一个graph::graph-class
ls("package:mypackage"):查看"mypackage"中的所有对象。

安装source package方法

1 在终端输入 # R CMD INSTALL /.../mypackage.tar.gz
使用此方法,需要解决包依赖问题,即安装此包所依赖的包,安装过程有提示

2 也可以使用R的install.packages()函数安装
回答:可以使用install.packages()函数安装,而且比较简便,即联网即可装,装了就可用。
# R
> install.packages('mypackage')

回答2:可以使用install.packages()安装本地下载的包,尤其适用于在服务器上安装包

$ R

> install.packages( c("XML_0.99-5.tar.gz", "http://www.cnblogs.com/Interfaces/Perl/RSPerl_0.8-0.tar.gz"), repos = NULL, configure.args = c(XML = '--with-xml-config=xml-config', RSPerl = "--with-modules='IO Fcntl'"))
3 Bioconductor的安装方法
> source("http://bioconductor.org/biocLite.R")
> biocLite("mypackage")

 

4 卸载package

remove.packages("mypackage")

 

5 查看R及其package的version

R version: version 或者 R.version

R package version:

 

6 更新包

update.packages( )  可以定期执行以下

 

7 使用别人安装的包

修改.bashrc文件,添加环境变量R的lib路径

export R_LIBS=/home/.../R/lib64/R/library

R中用.libPaths()函数查看lib路径,如果有多个lib,install.packages()默认是安装在第一个目录下

 

17

转录组cummeRbund操作笔记

转录组cummeRbund操作笔记

这是跟tophat和cufflinks套装紧密搭配使用的一个R包,能出大部分文章要求的标准化图片。

一:安装并加装该R包

安装就用source("http://bioconductor.org/biocLite.R") ;biocLite("cummeRbund")即可,如果安装失败,就需要自己下载源码包,然后安装R模块。

转录组cummeRbund操作笔记220

然后把cuffdiff输出的文件目录拷贝到R的工作目录,或者自己设置工作目录

 

二:读取FN目录下面的所有文件。

转录组cummeRbund操作笔记239

可以看到把cuffdiff下面的文件夹所有的文件都读取到了,里面有如下文件,包括genes,isoforms,cds,tss这四种差异情况都读取了。

转录组cummeRbund操作笔记316

 

三:表达水平分布图

转录组cummeRbund操作笔记328

转录组cummeRbund操作笔记330
四、表达水平箱线图

csBoxplot(genes(cuff_data))

转录组cummeRbund操作笔记371
五、画基因表达差异热图

转录组cummeRbund操作笔记386

画出热图如下

转录组cummeRbund操作笔记396

 

六、得到差异的genes,isoforms,TSS,CDS等等

 

  • 得到上调下调基因列表

diffData <- diffData(myGenes )

转录组cummeRbund操作笔记430

转录组cummeRbund操作笔记474

 

只有一百个有表达差异的基因

转录组cummeRbund操作笔记490

 

 

最后贴出一个综合性的代码,算了,太浪费空间了,把整个空间搞得不好看,就不贴了。

这个代码可以自动运行出图;

转录组cummeRbund操作笔记3781

16

转录组edgeR分析差异基因

转录组edgeR分析差异基因

edgeR是一个研究重复计数数据差异表达的Bioconductor软件包。一个过度离散的泊松模型被用于说明生物学可变性和技术可变性。经验贝叶斯方法被用于减轻跨转录本的过度离散程度,改进了推断的可靠性。该方法甚至能够用最小重复水平使用,只要至少一个表型或实验条件是重复的。该软件可能具有测序数据之外的其他应用,例如蛋白质组多肽计数数据。可用性:程序包在遵循LGPL许可证下可以从Bioconductor网站。

一:下载安装该软件

下载安装edgeR这个R包,因为这是一次讲R包的下载,我就啰嗦一点,这种生物信息学的包不同于普通的R包,是需要用biocLite来安装的,命令如下

转录组edgeR分析差异基因304

 

Continue reading

15

仿写fastqc软件的一些功能-R代码

仿写fastqc软件的一些功能(下)

文件来自于上面perl代码的输出文件,好像算法有点问题,26G的文件居然处理近一个小时才出数据!

仿写fastqc软件的一些功能-下-R代码263

R语言本身自带的画图工具都很丑,懒得说了,可以用ggplot2来重新画一个,不是项目要求没有报酬我就懒得画了,大家面前看看画图原理即可。

Continue reading