十二 25

TCGA表达数据的多项应用之1–下载数据并且导入mysql

这个TCGA表达数据的多项应用系列帖子是应群里朋友的要求来写的,你们也可以继续提需求,我会接着写下去,其实从TCGA数据库里面下载到了数据之后,后面的所有分析都跟TCGA没有半毛钱关系了,大家要有这个想法,别三两句就问TCGA数据怎么分析,http://www.bio-info-trainee.com/?s=TCGA&submit=Search 本系列最后会形成一个shiny版本的交互式表达数据查询,处理,绘图,统计的网页APP。
我这里偷懒一下了,直接下载GEO里面的TCGA的表达数据,而不是去TCGA的官网里面下载:
它处理了目前(大概是2015年6月)TCGA收集的所有癌症样本的mRNA表达数据,并且统一处理成了count和RPKM两种表达量形式。 GEO地址:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Continue reading

十一 07

mysql的table居然有最大列限制

想着把TCGA的RPKM值矩阵表格写入到mysql,然后做一个查询网页给生物学家,我下载的是所有TCGA收集的mRNA表达数据集数据集-GSE62944 ,共9264个癌症样本,和741个正常组织的表达数据。当我想写入癌症表达矩阵的时候,报错了:
Error in .local(conn, statement, ...) :
could not run statement: Too many columns
简单搜索了一下,发现是mysql有最大列的限制,但是我不是很懂计算机,所以没太看明白该如何调整参数使得mysql列限制扩充:http://dev.mysql.com/doc/refman/5.7/en/column-count-limit.html 所以就把癌症表达矩阵根据癌症拆分了,癌种数量如下:

Continue reading

05

用GEMINI来探索vcf格式的突变数据

第一次听说这个软件,是一个香港朋友推荐的:http://davetang.org/muse/2016/01/13/getting-started-with-gemini/ 他写的很棒,但是我当初以为是一个类似于SQLite的数据库浏览模式,所以没在意。实际上,我现在仍然觉得这个软件没什么用!

软件官网有详细的介绍:https://gemini.readthedocs.io/en/latest/

而且提供丰富的教程:

We recommend that you follow these tutorials in order, as they introduce concepts that build upon one another.

  • Introduction to GEMINI, basic variant querying and data exploration. html pdf
  • Identifying de novo mutations underlying Mendelian disease html pdf
  • Identifying autosomal recessive variants underlying Mendelian disease html pdf
  • Identifying autosomal dominant variants underlying Mendelian disease html pdf
  • Other GEMINI tools html pdf

软件本身并不提供注释,虽然它的功能的确包括注释,号称可以利用(ENCODE tracks, UCSC tracks, OMIM, dbSNP, KEGG, and HPRD.)对你的突变位点注释,比如你输入1       861389  .       C       T       ,它告诉你这个突变发生在哪个基因,对蛋白改变如何?是否会产生某些疾病?

虽然它本身没有注释功能,但是它会调用snpEFF或者VEP进行注释,你需要自己先学习它们。

1

软件安装:

GEMINI是用python写的,有一个小脚本可以自动完成安装过程:

7.3K May  4 14:44 gemini_install.py

下载这个脚本,然后安装即可

wget https://github.com/arq5x/gemini/raw/master/gemini/scripts/gemini_install.py

python gemini_install.py $tools $data

PATH=$tools/bin:$data/anaconda/bin:$PATH

where $tools and $data are paths writable on your system.

我把$tools用的就是当前文件夹,$data也是当前文件夹下面的gemini文件夹。

这样就会在当前文件夹下面生成两个文件夹,bin是存储程序,gemini是存储数据用的,而且注意要把bin目录的全路径添加到环境变量!

输入数据:

我们可以直接下载软件作者提供的测试数据

首先是22号染色体的所有突变位点经过WEP注释的文件

然后是一个三口直接的突变ped格式数据

数据存放在亚马逊云,所有的教程pdf也在

http://s3.amazonaws.com/gemini-tutorials

如果是你自己的vcf文件,需要自己用VEP注释一下

1

运行命令:

2

结果解读:

产生是chr22.db就是一个数据库格式的文件,但是需要用gemini 来进行查询,个人认为,并没有多大意思!

你只要熟悉mySQL等SQL语言,完全可以自己来!

09

生信数据库的schema如何构建

大家分析生物信息学数据的时候必不可少的步骤就是利用各种公共资源对自己的数据进行注释。

这时候可能会用到mysql,把一些公共数据库本地化,方便使用,但是数据的下载已经存储到mysql等数据库中间会有很多值得玩味的事情。

我这里给大家指出一个还算比较标准的参考,就是bioconductor官方制作的数据库设计代码。

bioconductor官方注释方面的包(主要是各种ID的转换,KEGG或者GO这样的功能注释,基因信息注释,转录本,外显子起始终止等等)

目前为止,bioconductor是3.3版本,共896个包

大部分包都是以sqlite的数据库标准发布,所以建表语句是一样的。

所有代码见:https://github.com/Bioconductor-mirror/AnnotationDbi/blob/release-3.2/inst/DBschemas

部分代码如下:

CREATE TABLE metadata (
name VARCHAR(80) PRIMARY KEY,
value VARCHAR(255)
);
CREATE TABLE go_ontology (
ontology VARCHAR(9) PRIMARY KEY,               -- GO ontology (short label)
term_type VARCHAR(18) NOT NULL UNIQUE          -- GO ontology (full label)
);
CREATE TABLE go_term (
_id INTEGER PRIMARY KEY,
go_id CHAR(10) NOT NULL UNIQUE,               -- GO ID
term VARCHAR(255) NOT NULL,                   -- textual label for the GO term
ontology VARCHAR(9) NOT NULL,                 -- REFERENCES go_ontology
definition TEXT NULL,                         -- textual definition for the GO term
FOREIGN KEY (ontology) REFERENCES go_ontology (ontology)
);
CREATE TABLE sqlite_stat1(tbl,idx,stat);
CREATE TABLE go_obsolete (
go_id CHAR(10) PRIMARY KEY,                   -- GO ID
term VARCHAR(255) NOT NULL,                   -- textual label for the GO term
ontology VARCHAR(9) NOT NULL,                 -- REFERENCES go_ontology
definition TEXT NULL,                         -- textual definition for the GO term
FOREIGN KEY (ontology) REFERENCES go_ontology (ontology)
);

 

 

 

28

把bioconductor的gene mapping信息上传到mysql数据库

在R语言里面的bioconductor系列包里面有一个物种注释信息包,其中人类是org.Hs.eg.db
里面有基于人的hg19基因组版本的大部分基因信息之间的转换数据!
包括基因的entrez ID,symbol,name,locus,refseq,kegg pathway,GO ontology对应关系!
但是它们散布在该包的各个数据结构里面!
> ls("package:org.Hs.eg.db")
 [1] "org.Hs.eg"                "org.Hs.eg.db"            
 [3] "org.Hs.eg_dbconn"         "org.Hs.eg_dbfile"        
 [5] "org.Hs.eg_dbInfo"         "org.Hs.eg_dbschema"      
 [7] "org.Hs.egACCNUM"          "org.Hs.egACCNUM2EG"      
 [9] "org.Hs.egALIAS2EG"        "org.Hs.egCHR"            
[11] "org.Hs.egCHRLENGTHS"      "org.Hs.egCHRLOC"         
[13] "org.Hs.egCHRLOCEND"       "org.Hs.egENSEMBL"        
[15] "org.Hs.egENSEMBL2EG"      "org.Hs.egENSEMBLPROT"    
[17] "org.Hs.egENSEMBLPROT2EG"  "org.Hs.egENSEMBLTRANS"   
[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"         
[21] "org.Hs.egENZYME2EG"       "org.Hs.egGENENAME"       
[23] "org.Hs.egGO"              "org.Hs.egGO2ALLEGS"      
[25] "org.Hs.egGO2EG"           "org.Hs.egMAP"            
[27] "org.Hs.egMAP2EG"          "org.Hs.egMAPCOUNTS"      
[29] "org.Hs.egOMIM"            "org.Hs.egOMIM2EG"        
[31] "org.Hs.egORGANISM"        "org.Hs.egPATH"           
[33] "org.Hs.egPATH2EG"         "org.Hs.egPFAM"           
[35] "org.Hs.egPMID"            "org.Hs.egPMID2EG"        
[37] "org.Hs.egPROSITE"         "org.Hs.egREFSEQ"         
[39] "org.Hs.egREFSEQ2EG"       "org.Hs.egSYMBOL"         
[41] "org.Hs.egSYMBOL2EG"       "org.Hs.egUCSCKG"         
[43] "org.Hs.egUNIGENE"         "org.Hs.egUNIGENE2EG"     
[45] "org.Hs.egUNIPROT"        
其实比较重要的信息,我们需要把它们统一关联起来,做成一个表格,这样可以上传到数据库里面,做成网页,可视化展现,查找!
suppressMessages(library(org.Hs.eg.db))
all_EG=mappedkeys(org.Hs.egSYMBOL)
###下面的表格最后是基于entrez ID而关联起来的!
tmp=unlist(as.list(org.Hs.egSYMBOL))
EG2Symbol=data.frame(EGID=names(tmp),symbol=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egENSEMBL))
EG2ENSEMBL=data.frame(EGID=names(tmp),ENSEMBL=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egGENENAME))
EG2name=data.frame(EGID=names(tmp),name=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egMAP))
EG2MAP=data.frame(EGID=names(tmp),MAP=as.character(tmp))
     
###EG2GO and    using mySQL
EG2path=as.list(org.Hs.egPATH)
EG2path=lapply(EG2path, function(x) paste(x,collapse = ":"))
tmp=unlist(EG2path)
EG2path=data.frame(EGID=names(tmp),path=as.character(tmp))
 
#as.list(head(org.Hs.egGO2ALLEGS))
GO2allEG=as.list(org.Hs.egGO2ALLEGS)
tmp=unlist(GO2allEG)
names(tmp)=substring(names(tmp),1,10) ## change GO:0000002.IMP to GO:0000002 
EG2allGO <- tapply(tmp,tmp,function(x){names(x)})
EG2allGO=lapply(EG2allGO, function(x) paste(x,collapse = ":"))
tmp=unlist(EG2allGO)
EG2GO=data.frame(EGID=names(tmp),GO=as.character(tmp))
 
tmp=merge(EG2Symbol,EG2MAP,by='EGID',all=TRUE)
tmp=merge(tmp,EG2ENSEMBL,by='EGID',all=TRUE)
tmp=merge(tmp,EG2path,by='EGID',all=TRUE)
tmp=merge(tmp,EG2name,by='EGID',all=TRUE)
my_gene_mapping=merge(tmp,EG2GO,by='EGID',all=TRUE)
##[1] "EGID"    "symbol"  "MAP"     "ENSEMBL" "path"    "name"    "GO"
1
因为我用的merge的all=TRUE,所以最后记录会越来越多
最后的结果就是下面这样,各种ID转换都储存到了my_gene_mapping这个数据对象里面!
可以看到有些基因是没有ensemble数据库对应的ID的,非常多的基因没有被注释到KEGG通路数据库!
我这里如果一个基因对应多个通路,我用冒号连接起来了,成了一个字符串!在后面会有用!

2

现在需要把它上传到数据库里面,这样我就可以用R的shiny可视化网页来查找数据库,做ID转换啦!!!
suppressMessages(library(RMySQL))
con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")
dbSendQuery(con, "USE test")
dbWriteTable(con,'my_gene_mapping',my_gene_mapping)
dbDisconnect(con)
然后就可以在自己的mysql里面看到它啦!!! 
其实 "org.Hs.egOMIM"  也很重要,不过我这里只是举个例子,就不深究了!
还有一点,就是那个pathway ID,因为我要以entrez ID为主键,所以把多个pathway给合并了!
suppressMessages(library(org.Hs.eg.db))
kegg2ID=toTable(org.Hs.egPATH)
#[1] "gene_id" "path_id"
dbWriteTable(con,'kegg2ID_bioconductor',kegg2ID,row.name=F,overwrite=T)
go2id=toTable(org.Hs.egGO2ALLEGS)
## gene_id      go_id Evidence Ontology
dbWriteTable(con,'go2id_bioconductor',go2id,row.name=F,overwrite=T)
dbDisconnect(con)
用这个代码,可以在数据库里面多创建两个表,会有用的!
16

生物信息学学者学习mysql之路

我一直都知道mysql其实很有用的, 哪怕是在bioinformatics领域。也断断续续的看过不少mysql教程,只是苦于没有机会应用。毕竟应用才是最好的学习方法,正好这些天需要用了,我就又梳理了一遍作为一个生物信息学学者,该如何学习mysql数据库。
然后再搜搜一堆技巧
差不多就可以开始啦。
我们不拿数据库来做网页,所以需要的仅仅是查询公共数据库的数据,当然,一般人都会选择直接去网页可视化的查询,或者去ftp批量下载后自己写脚本来查询,我以前也是这样想的,所以感觉mysql没什么用,因为它能做的, 我写一个脚本都能做到。但是任何事物能发展到如此流行的程度毕竟还是有它的优点的。
而在我看来,mysql的优点就是,不需要存储大量的文件信息,随查随用,如果我们想把数据库备份到本地,就要建立一大堆的文件夹,存放各种refgene信息呀,entrez gene信息呀,转录本,外显子等等各个文件夹,每个文件夹下面还有一堆文件,而且还要分物种存储,总之就是很麻烦,但是在数据库就不一样啦。
比如我们可以连接UCSC的数据库 (前提是你的机器里面可以允许mysql这个命令,而且你可以联网)
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A
就这么简单, 你就用mysql远程登录了UCSC的数据库,可以show databases;或者use database hg19 ; 等等
里面有两百多个数据库,主要是多物种多版本,然后如果我们看hg19这个数据库,里面还有一万多个数据表,包含着hg19的全面信息。
还有很多其它的公共数据库可以练习
来自于:https://www.biostars.org/p/474/#9095

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29
ENSEMBL http://uswest.ensembl.org/info/data/mysql.html
GO http://www.geneontology.org/GO.database.shtml#mirrors

1000 Genomes: since June 16, 2011: http://www.1000genomes.org/public-ensembl-mysql-instance

mysql -h mysql-db.1000genomes.org -u anonymous -P 4272

Flybase has direct access to its postgres chado database.
http://flybase.org/forums/viewtopic.php?f=14&t=114
hostname: flybase.org port: 5432 username: flybase password: no password database name: flybase
e.g. psql -h flybase.org -U flybase flybase

mysql -h database.nencki-genomics.org -u public
mysql -h useastdb.ensembl.org -u anonymous -P 5306

你都可以登录进去看看里面有什么,也可以练习练习mysql的语法,但是增删改查种的查是可以用的
然后我们可以用R或者perl或者Python来连接数据库,也是蛮好用的, 我现在比较倾向于R
所以我就简单看了一下这个包的说明书,然后成功连接了
#Connect to the MySQL server using the command:
#mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A
#The -A flag is optional but is recommended for speed
library(RMySQL)
my.port="";
my.user="genome";
my.password="";
my.db="hg19";
#there are 203 databases,such as hg18,hg38,mm9,mm10,ce10
con <- dbConnect(MySQL(), host=my.host, user=my.user,dbname=my.db)
dbListTables(con) # there are 11016 tables in this hg19 database;

是不是很简单呀,只有你认真的学习,其实这些应用的东西都还是蛮简单的。

下面这本书也比较好,就讲了R或者perl或者Python来连接数据库,很全面

当然,如果想看mysql在bioinformatics方面的应用,下面还有很多学习资料
进阶版还可以看看具体事例,GO数据库的设计:http://geneontology.org/page/lead-database-schema
从这个来看,python要比perl 好很多http://www.personal.psu.edu/iua1/courses/files/2010/week15.pdf
01

ubuntu服务器解决方案第九讲-mysql和apache的安装

发现自己搞服务器遇到的困难还是蛮多的,所以记录了一下,给菜鸟们指个路。
ubuntu对生信菜鸟来说是最好用的linux服务器,没有之一,因为它有apt-get。

很多时候大家的服务器可能并不是想联网,只是想玩一下,或者只是因为生信的某些软件要求数据库,所以大家可能会单独安装mysql,或者想学习perl的CGI模块,需要apache。

ubuntu上安装mysql

非常简单只需要几条命令就可以完成。

1. sudo apt-get install mysql-server

2. sudo apt-get install mysql-client

3.  sudo apt-get install libmysqlclient-dev

安装过程中会提示设置密码什么的,注意设置了不要忘了,安装完成之后可以使用如下命令来检查是否安装成功:

sudo netstat -tap | grep mysql

通过上述命令检查之后,如果看到有mysql 的socket处于 listen 状态则表示安装成功。

登陆mysql数据库可以通过如下命令:

mysql -u root -p

-u 表示选择登陆的用户名, -p 表示登陆的用户密码,上面命令输入之后会提示输入密码,此时输入密码就可以登录到mysql。

Ubuntu上安装Apache

Ubuntu上安装Apache,有两种方式:1 使用开发包的打包服务,例如使用apt-get命令;2 从源码构建Apache。本文章将详细描述这两种不同的安装方式。

方法一:使用开发包的打包服务——apt-get

安装apache,在命令行终端中输入一下命令:

$ sudo apt-get install apache2

安装完成后,重启apache服务,在命令行终端中输入一下命令:

$ sudo /etc/init.d/apache2 restart

可能会出现的问题1: NameVirtualHost *:80 has no VirtualHosts,

出现上述问题的原因:定义了多个NameVirtualHost,故将/etc/apache2/ports.conf中的NameVirtualHost *:80注释掉即可。

可能会出现的问题2: Could not reliably determine the server's fully qualified domain name, using 127.0.1.1 for ServerName

原因:

根据提示,无法可靠的确定服务器的有效域名,使用127.0.1.1作为服务器域名。应此,在下面的测试中,应该使用127.0.1.1,而不是127.0.0.1!

解决:

$ vim /etc/apache2/httpd.conf,在文件中添加:

ServerName localhost:80,再次重启apache2,就可以使用127.0.0.1来访问web服务器啦!

测试:

在浏览器里输入http://localhost或者是http://127.0.0.1,如果看到了It works!,那就说明Apache就成功的安装了,Apache的默认安装,会在/var下建立一个名为www的目录,这个就是Web目录了,所有要能 过浏览器访问的Web文件都要放到这个目录里。