Category Archives: 基础数据库

主要包含了NCBI，ensembl，UCSC，uniprot，IMGT，KEGG，OMIN，TIGR，GO等常见数据库

八 04

生信人必学ftp站点之 dbsnp

Posted on by
这个数据库我也不想多解释了，也是host在NCBI上，不仅有常见的模式生物已经被研究过的所有variation位点信息，还有很多其它物种的数据，主站点是：ftp://ftp-trace.ncbi.nih.gov/snp/organisms/
人类是物种ID是9606，可以看到variation位点信息有基于hg19和hg38的两种下载方式，如果还有其它需求，可以自己用基因组坐标转换工具。在NCBI的snp页面也有对各种物种的variation位点信息记录文件的统计：http://www.ncbi.nlm.nih.gov/snp/   http://www.ncbi.nlm.nih.gov/SNP/同时也是NCBI做好的一个网页版查询工具，因为下载一个 variation位点信息记录文件 动辄就是十几个G，一般人也不会处理那个文件，不知道从里面应该如何提取需要的信息，这时候学习它的网页版查询工具也挺好的。

Continue reading →

七 08

基因组标准注释文件-Gencode数据库

Posted on by
Gencode数据库是ENCODE计划的衍生品，也是由大名鼎鼎的sanger研究所负责整理和维护，主要记录了基因组的功能注释，比如基因组每条染色体上面有哪些编码蛋白的基因，哪些假基因，哪些lncRNA的基因，它们坐标是什么，基因上面的外显子内含子坐标是什么，UTR区域坐标是什么？我以前通常是在EBI的ENSEMBL的FTP服务器下载，后来才发现了这个Gencode数据库，现在以这个为金标准啦！

Continue reading →

五 16

假基因资源中心

Posted on by
假基因是原来的能翻译成蛋白的基因经过各种突变导致丧失功能的基因。
比如
PTEN-->PTENP1
KRAS-->KRASP1
NANOG-->NANOGP1
很好理解，一般来说看到结尾是P1,等字眼的都是假基因，现在共有一万多假基因，我一般以http://www.genenames.org/cgi-bin/statistics （人类基因命名委员会）为标准参考。
研究的时候可能需要更全面一点，所以我又谷歌了一下，发现了一个还算比较全面的收集。
就是 http://pseudogene.org/Human/  （中心网站）
现在主要是 ENCODE计划的GENCODE 21. 和 耶鲁大学的Ensembl genome release 79.
Human Pseudogene Annotation

GENCODE Annotation

- Data: The current human pseudogene annotation is in GENCODE 21. .

- Description: The GENCODE annotation of pseudogenes contains models that have been created by the Human and Vertebrate Analysis and Annotation (HAVANA) team, an expert manual annotation team at the Wellcome Trust Sanger Institute. This is informed by, and checked against, computational pseudogene predictions by thePseudoPipe and RetroFinder pipelines.

PseudoPipe Output

- Data: The current PseudoPipe results are on Ensembl genome release 79. .

- Description: Genome-wide human pseudogene annotation predicted by PseudoPipe. PseudoPipe is a homology-based computational pipeline that searches a mammalian genome and identifies pseudogene sequences.

- Reference:

Other Human Pseudogene Sets

- Data: .

- Description: Archived pseudogene annotation on previous human genome releases from PseudoPipe. Genome-wide annotation or specific subset.
五 05

用samr包对芯片数据做差异分析

Posted on by

本来搞差异分析的工具和包就一大堆了，而且limma那个包已经非常完善了，我是不准备再讲这个的，正好有个同学问了一下这个包，我就随手测试了一下，顺便看看它跟limma有什么差异没有！手痒了就记录了测试流程！

学习一个包其实非常简单，就是找到包的官网看看说明书即可！说明书链接

 

Continue reading →

五 05

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

Posted on by

第一次听说这个软件，是一个香港朋友推荐的：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语言，完全可以自己来！

五 05

用VEP对vcf格式的突变数据进行注释

Posted on by

VEP是国际三大数据库之一的ENSEMBL提供的，也是非常主流和方便，但它是基于perl语言的，所以在模块方面可能会有点烦人。跟snpEFF一样，也是对遗传变异信息提供更具体的注释，而不仅仅是基于位点区域和基因。如果你熟悉外显子联盟这个数据库EXAC(ExAC.r0.3.sites.vep.vcf.gz)，你可以下载它所有的突变记录数据，看看它对每个变异位点到底注释了些什么，它就是典型的用VEP来注释的。 Continue reading →

五 05

用snpEFF对vcf格式的突变数据进行注释

Posted on by

这个软件比较重要，尤其是对做遗传变异相关研究的，很多人做完了snp-calling后喜欢用ANNOVAR来进行注释，但是那个注释还是相对比较简单，只能得到该突变位点在基因的哪个区域，那个基因这样的信息，如果想了解更具体一点，就需要更加功能化的软件了，snpEFF就是其中的佼佼者，而且是java平台软件，非常容易使用！而且它的手册写的非常详细：http://snpeff.sourceforge.net/SnpEff_manual.html Continue reading →

四 28

下载最新的蛋白相互作用数据库-STRING

Posted on by

string数据库是PPI领域里面最完备已经最受欢迎的数据库了。如果直接在谷歌里面搜索PPI，映入眼帘就是string的官网，它们的主页现在是html5啦，比较精美： http://string-db.org/

1

写的很霸气，近两亿的记录，不过一般大家只会关心一个物种，比如人，其实还不到一千万！

我们直接进入下载界面，找到人类的数据，人类的物种ID是9606.

2

需要一定许可才能下载完整版本，我这里测试最上面那个公开版本数据！

数据很简单，就是protein+protein+score，共八百多万行记录，记录着string数据库搜集的所有可能以及可信的蛋白相互作用！但是它的蛋白ID是ENSEMBL的ID，所以需要转换成基因的ID，才能被大多数人使用，因为大家的研究单位一般是基因，所以蛋白相互作用略等于基因相互作用。

基因ID转换，我推荐用org.Hs.eg.db这个R的包，很容易就可以实现的！

> tmp=toTable(org.Hs.egENSEMBLPROT)
> dim(tmp)
[1] 110916      2
> head(tmp)
  gene_id         prot_id
1       1 ENSP00000263100
2       1 ENSP00000470909
3       2 ENSP00000443302
4       2 ENSP00000323929
5       2 ENSP00000438599
6       2 ENSP00000445717
>

有约500多个蛋白ID是无法转换成对应的基因的，这个很正常，毕竟这种ID本来就不稳定，很多用着用着就失效了！

转换好之后就可以上传到数据库啦，然后可以供其它可视化或者分析程序使用！

 

四 09

生信数据库的schema如何构建

Posted on by

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

这时候可能会用到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)
);

 

 

 

三 15

基因组各种版本对应关系

Posted on by 
我是受到了SOAPfuse的启发才想到整理各种基因组版本的对应关系，完整版！！！
以后再也不用担心各种基因组版本混乱了，我还特意把所有的下载链接都找到了，可以下载任意版本基因组的基因fasta文件，gtf注释文件等等！！！
首先是NCBI对应UCSC，对应ENSEMBL数据库：
GRCh36 (hg18): ENSEMBL release_52.
GRCh37 (hg19): ENSEMBL release_59/61/64/68/69/75.
GRCh38 (hg38): ENSEMBL  release_76/77/78/80/81/82.
可以看到ENSEMBL的版本特别复杂！！！很容易搞混！
但是UCSC的版本就简单了，就hg18,19,38, 常用的是hg19，但是我推荐大家都转为hg38
看起来NCBI也是很简单，就GRCh36,37,38，但是里面水也很深！
Feb 13 2014 00:00    Directory April_14_2003
Apr 06 2006 00:00    Directory BUILD.33
Apr 06 2006 00:00    Directory BUILD.34.1
Apr 06 2006 00:00    Directory BUILD.34.2
Apr 06 2006 00:00    Directory BUILD.34.3
Apr 06 2006 00:00    Directory BUILD.35.1
Aug 03 2009 00:00    Directory BUILD.36.1
Aug 03 2009 00:00    Directory BUILD.36.2
Sep 04 2012 00:00    Directory BUILD.36.3
Jun 30 2011 00:00    Directory BUILD.37.1
Sep 07 2011 00:00    Directory BUILD.37.2
Dec 12 2012 00:00    Directory BUILD.37.3
可以看到，有37.1,   37.2，  37.3 等等，不过这种版本一般指的是注释在更新，基因组序列一般不会更新！！！
反正你记住hg19基因组大小是3G，压缩后八九百兆即可！！！
如果要下载GTF注释文件，基因组版本尤为重要！！！
对于ensembl：
变幻中间的release就可以拿到所有版本信息：ftp://ftp.ensembl.org/pub/
对于UCSC，那就有点麻烦了：
需要选择一系列参数：
2. Select the following options:
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Genes and Gene Predictions
track: UCSC Genes
table: knownGene
region: Select "genome" for the entire genome.
output format: GTF - gene transfer format
output file: enter a file name to save your results to a file, or leave blank to display results in the browser
3. Click 'get output'.
 现在重点来了，搞清楚版本关系了，就要下载呀！
UCSC里面下载非常方便，只需要根据基因组简称来拼接url即可：
或者用shell脚本指定下载的染色体号：
for i in $(seq 1 22) X Y M;
do echo $i;
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
## 这里也可以用NCBI的：ftp://ftp.ncbi.nih.gov/genomes/M_musculus/ARCHIVE/MGSCv3_Release3/Assembled_Chromosomes/chr前缀
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fasta;
done
rm -fr chr*.fasta
二 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 )

 

二 28

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

Posted on by
在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)
用这个代码，可以在数据库里面多创建两个表，会有用的！
二 21

画基因结构图！

Posted on by
一个基因有一个甚至多个转录本！根据转录及翻译现象可以把一个基因人为的定义成多种结构！
首先自己查资料搞明白转录本，外显子和内含子，5端UTR区域和CDS区域，还有3端UTR区域的概念。
真核生物的基因含有外显子和内含子，是区别原核生物的特征之一，能转录的就是外显子，不能转录的就是内含子！
内含子（英语：Intron）是一个基因中非编码DNA片段，它分开相邻的外显子。更精确的定义是：内含子是阻断基因线性表达的序列。DNA上的内含子会被转录到前体RNA中，但RNA上的内含子会在RNA离开细胞核进行转译前被剪除。在成熟mRNA被保留下来的基因部分被称为外显子。
通常我们看到基因结构示意图会像下面这样！
gene-structure
但是对于这个图，我之前有非常多的疑问，直到我自己写程序去仔细统计，画图探究后才真正搞明白！
从图中可以看到5端UTR区域后面接着的是第一外显子，但是事实上不是！
外显子和内含子是转录本上面的概念，是基于转录这个行为的定义。
而5端UTR区域和CDS区域，还有3端UTR区域,是基于翻译这种行为的定义！
把它们画成这样，是严重的干扰信息。
实际上，我们需要更正很多概念，我检验了一下，得到下面几个正确的：
1，如果基因有多个转录本，基因的起始坐标，就是该基因所有转录本的第一个外显子的起始坐标的最小值，同理基因的终止坐标，就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。
2，通过这个概念，可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。（这个是我乱说的，并没有这个定义）
3，如果基因只有一个转录本，那么基因的起始终止坐标，就是转录本的起始终止坐标！
4，一个基因的一个转录本的5'utr区域可以包括多个外显子区域，前者是翻译行为，后者是转录行为
5，起始密码子和终止密码子是CDS的起止处，是基于翻译的概念
6，一个基因的多个转录本的外显子坐标不一定会排列整齐，每个转录本的剪切位点并不一定要比其它转录本一致！
比如对这个ANXA1基因来说，非常多的转录本，但是基因的起始终止坐标，是所有转录本起始终止坐标的极大值和极小值！同时，它是一个闭合基因，因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的，虽然多个转录本会共享某些外显子，但是也存在某些外显子比同区域其它外显子长的现象！
 anxa1-gene-structure
再比如下面这个例子：对DDX11L11这个基因来说，前两个外显子都不会翻译，直到第三个外显子才开始翻译，构成CDS。
有些转录本是没有utr的，所以该转录本的起始坐标，就是CDS的起始坐标
5UTR包括多个exon行为
首先把gtf文件格式化导入mysql数据库
我用的是broadinstitute提供的GENCODE的gtf文件，是hg19版本的
Modified GENCODE GTF file for human with contig names of form ("1","2", etc)
##description: evidence-based annotation of the human genome (GRCh37), version 7 (Ensembl 62)
##provider: GENCODE
##format: gtf
##date: 2011-03-23
很简单，下载数据，parse好之后，用R导入mysql即可！
下面是mysql代码：
--mysql code
use test;
drop table  if exists hg19_gtf;
create table hg19_gtf (
    gene_name VARCHAR(30),
    transcript_name VARCHAR(30) ,
    record  VARCHAR(15) NOT NULL ,
    chr VARCHAR(2) NOT NULL ,
    start INT NOT NULL ,
    end INT NOT NULL ,
    source VARCHAR(10) NOT NULL ,
    strand VARCHAR(1) NOT NULL ,
    gene_id VARCHAR(30) NOT NULL ,
    transcript_id VARCHAR(30) NOT NULL ,
    gene_status VARCHAR(30) ,
    gene_type VARCHAR(30)  ,
    transcript_type VARCHAR(30) ,
    transcript_status VARCHAR(30)
);
--我的网页好像不支持mysql代码高亮，大家凑合着看吧，反正就是一个简单的建表语句！
select * from hg19_gtf limit 100;
select * from hg19_gtf where gene_name='DMD';
select count(*) from hg19_gtf where gene_name='DMD' and record='start_codon';  --18 start condon
select count(distinct(transcript_name)) from hg19_gtf where gene_name='DMD' ;  --34 transcript
select count(distinct(transcript_name)) c ,gene_name from hg19_gtf where record='transcript' group by gene_name  order by c desc;
我是随意设计的一个表，主要是为了画图方便！
接下来是R code来具体的对基因进行画图！

[perl]
suppressMessages(library(ggplot2))
suppressMessages(library(RMySQL))
con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")
dbSendQuery(con, "USE test")
gene='SOX10'
#gene='DDX11L11'
if (T){
query=paste("select * from hg19_gtf where gene_type='protein_coding' and gene_name=",shQuote(gene),sep="")
structure=dbGetQuery(con,query)
tmp_min=min(c(structure$start,structure$end))
structure$new_start=structure$start-tmp_min
structure$new_end=structure$end-tmp_min
tmp_max=max(c(structure$new_start,structure$new_end))
num_transcripts=nrow(structure[structure$record=='transcript',])
tmp_color=rainbow(num_transcripts)
x=1:tmp_max;y=rep(num_transcripts,length(x))
#x=10000:17000;y=rep(num_transcripts,length(x))
plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))
title(main = gene,sub = paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
j=0;
tmp_legend=c()
for (i in 1:nrow(structure)){
tmp=structure[i,]
if(tmp$record == 'transcript'){
j=j+1
tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
}
if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}
# legend('topleft',legend=tmp_legend,lty=1,lwd = 4,col = tmp_color);

}
[/perl]

通过这个代码，，就能简单的得到我上面显示的基因结构图啦！

 

二 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"))
##自己把内容写出来即可
  }
}

 

 

一 14

华盛顿大学把所有的变异数据都用自己的方法注释了一遍，然后提供下载

Posted on by
华盛顿大学把所有的变异数据都用自己的方法注释了一遍，然后提供下载：
文献是：Kircher M, Witten DM, Jain P, O'Roak BJ, Cooper GM, Shendure J. 

A general framework for estimating the relative pathogenicity of human genetic variants.
Nat Genet. 2014 Feb 2. doi: 10.1038/ng.2892.
PubMed PMID: 24487276.

文中的观点是：现在大多的变异数据注释方法都非常单一，通常是看看该位点是否保守，对蛋白功能的改变，在什么domain上面等等。
但这样是远远不够的，所以他们提出了一个新的注释方法，用他们自己的CADD方法把现存的一些公共数据库的变异位点（约86亿的位点）都注释了一下，并对每个位点进行了打分。
C scores correlate with allelic diversity, annotations of functionality, pathogenicity, disease severity, experimentally measured regulatory effects and complex trait associations, and they highly rank known pathogenic variants within individual genomes.
总之，他们的方法是无与伦比的！
所有他们已经注释好的数据下载地址是：http://cadd.gs.washington.edu/download
这些数据在很多时候非常有用，尤其是想跟自己得到的突变数据做交叉验证，或者做一下统计分析的时候！
clipboard
人的基因组才300亿个位点，他们就注释了86亿！！！
所以有三百多G的压缩包数据，我想，一般的公司或者单位都不会去用这个数据了！
一 14

蛋白质相互作用（PPI）数据库大全

Posted on by
最近遇到一个项目需要探究一个gene list里面的基因直接的联系，所以就想到了基因的产物蛋白的相互作用关系数据库，发现这些数据库好多好多！
一个比较综合的链接是：A compendium of PPI databases can be found in http://www.pathguide.org/.

里面的数据库非常多，仅仅是对于人类就有

Your search returned 207 results in 9 categories with the following search parameters:

人类的六个主要PPI是：Analysis of human interactome PPI data showing the coverage of six major primary databases (BIND, BioGRID, DIP, HPRD, IntAct, and MINT), according to the integration provided by the meta-database APID.
BIND the biomolecular interaction network database died link
DIP the database of interacting proteins http://dip.doe-mbi.ucla.edu/ 
MINT the molecular interaction database http://mint.bio.uniroma2.it/mint/ 
STRING Search Tool for the Retrieval of Interacting Genes/Proteins http://string-db.org/  
HPRO Human protein reference database http://www.hprd.org/ 
BioGRID The Biological General Repository for Interaction Datasets http://thebiogrid.org/ 
这些数据库大部分都还有维护者，还在持续更新，每次更新都可以发一篇paper，而数据库收集的paper引用一般都上千，如果你做了一个数据库，才十几个人引用，那就说明你是自己在跟自己玩。
其中比较好用的是宾夕法尼亚州匹兹堡的大学的一个：http://severus.dbmi.pitt.edu/wiki-pi/
(a) PPI definition; a definition of a protein-to-protein interaction compared to other biomolecular relationships or associations.
(b)PPI determination by two alternative approaches: binary and co-complex; a description of the PPIs determined by the two main types of experimental technologies.
(c) The main databases and repositories that include PPIs; a description and comparison of the main databases and repositories that include PPIs, indicating the type of data that they collect with a special distinction between experimental and predicted data.
(d) Analysis of coverage and ways to improve PPI reliability; a comparative study of the current coverage on PPIs and presentation of some strategies to improve the reliability of PPI data.
(e) Networks derived from PPIs compared to canonical pathways; a practical example that compares the characteristics and information provided by a canonical pathway and the PPI network built for the same proteins. Last, a short summary and guidance for learning more is provided.
现在的蛋白质相互作用数据库的数据都很有限，但是在持续增长，一般有下面四种原因导致数据被收录到数据库
There are four common approaches for PPI data expansions:
1) manual curation from the biomedical literature by experts;
2) automated PPI data extraction from biomedical literature with text mining methods;
3) computational inference based on interacting protein domains or co-regulation relationships, often derived from data in model organisms; and
4) data integration from various experimental or computational sources.
Partly due to the difficulty of evaluating qualities for PPI data, a majority of widely-used PPI databases, including DIP, BIND, MINT, HPRD, and IntAct, take a "conservative approach" to PPI data expansion by adding only manually curated interactions. Therefore, the coverage of the protein interactome developed using this approach is poor.
In the second literature mining approach, computer software replaces database curators to extract protein interaction (or, association) data from large volumes of biomedical literature . Due to the complexity of natural language processing techniques involved, however, this approach often generates large amount of false positive protein "associations" that are not truly biologically significant "interactions".
The challenge for the integrative approach is how to balance quality with coverage.
In particular, different databases may contain many redundant PPI information derived from the same sources, while the overlaps between independently derived PPI data sets are quite low .
参考：
2009年发表的HIPPI数据库：http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-10-S1-S16#CR6_2544 （是对HPRD [11], BIND [20], MINT [21], STRING [26], and OPHID数据库的整合）
一 14

居然可以下载千人基因组计划的所有数据bam，vcf数据

Posted on by
它有两个ftp站点存储所有的数据！
直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了！
也可以直接看最新版的vcf文件，记录了这两千多人的所有变异位点信息！
可以直接看到所有的位点，具体到每个人在该位点是否变异！
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610
它有两个ftp站点存储所有的数据！
直接看最新版的数据，共有NA编号开头的1182个人，HG开头的1768个人！
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了！
也可以直接看最新版的vcf文件，记录了这两千多人的所有变异位点信息！
可以直接看到所有的位点，具体到每个人在该位点是否变异！
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的，具体原理见：http://www.ncbi.nlm.nih.gov/pubmed/23093610
clipboard

 

十二 29

所有TCGA收集的mRNA表达数据集数据集-GSE62944

Posted on by
无意中看一篇文献，有提到这个数据集，我简单看了一下，居然还真的这么全面！！！
它处理了目前（大概是2015年6月）TCGA收集的所有癌症样本的mRNA表达数据，并且统一处理成了count和RPKM两种表达量形式。
GEO地址：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Title Alternatively processed and compiled RNA-Sequencing and clinical data for thousands of samples from The Cancer Genome Atlas
Organism Homo sapiens
Experiment type Expression profiling by high throughput sequencing
Summary We reprocessed RNA-Seq data for 9264 tumor samples and 741 normal samples across 24 cancer types from The Cancer Genome Atlas with "Rsubread". Rsubread is an open source R package that has shown high concordance with other existing methods of alignment and summarization, but is simple to use and takes significantly less time to process data. Additionally, we provide clinical variables publicly available as of May 20, 2015 for the tumor samples where the TCGA ids are matched.
这样就非常方便大家使用了。
可以直接拿来进行聚类，看看不同癌症种类如何聚集在一起，也可以直接拿来跟某个临床指标进行关联分析.
比如，如果我们想用DEseq来做差异分析，那么需要的数据就是基因的count，这里就有近一万个癌症样本的基因表达的count值，随便都可以分类看看差异情况，再富集分析分析。非常适合初学者练手。
癌症种类见官网：http://gdac.broadinstitute.org/
clipboard

 

十二 24

使用R包cgdsr来下载TCGA的数据

Posted on by

前面我讲到TCGA的数据可以在5个组织机构可以获取，他们都提供了类似的接口来供用户下载数据

每个接口都有使用教程，比如http://firebrowse.org/tutorial/FireBrowse-Tutorial.pdf

非常详细！！！

有的还专门写了软件接口：https://confluence.broadinstitute.org/display/GDAC/Download

或者写了R的接口：http://www.cbioportal.org/cgds_r.jsp

接下来我们要讲的就是cbioportal网站提供的一个R接口，非常好用，只需记住4个函数即可！！！ Continue reading →