十二 15

CpG Islands记录文件下载的4种方式

这个也是读者来信最多的,关于基因组某些区域的起始终止坐标的下载问题,genomic feature的问题,一般是gtf文件或者bed文件,比如人类hg19上面的所有外显子的坐标记录文件,所有基因的坐标记录文件,所有lncRNA,rRNA等等,我这里拿CpG Islands记录文件下载的4种方式举例子给大家说明一下: Continue reading

十二 11

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

这是系列文章,请先看:

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

ncbi现有的GPL已经过万了,但是bioconductor的芯片注释包不到一千,虽然bioconductor可以解决我们大部分的需要,比如affymetrix的95,133系列,深圳1.0st系列,HTA2.0系列,但是如果碰到比较生僻的芯片,bioconductor也不会刻意为之制作一个bioconductor的包,这时候就需要自行下载NCBI的GPL信息了,也可以通过R来解决:

##本质上是下载一个文件,读进R里面,然后解析行列式,得到芯片探针与基因的对应关系,看下面的代码,你就能理解了。 Continue reading

十二 01

吐血推荐snpedia数据库,非常丰富的snp信息记录

正好,我拿到了自己的全基因组测序数据,而前些天看到朋友圈推送的文章提到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,我就想顺便看看自己在这两个位点的变异情况。一般的流程是先找完变异位点,然后用vep/snpEFF对变异位点进行注释,然后看看有没有这两个位点。但我仅仅是想查看这两个位点,所以我会根据它的rsID来找到它的基因组坐标,再直接call这个位置的变异情况。以前我都是用dnSNP来查看rsID的基因组坐标的,
mkdir -p ~/annotation/variation/human/dbSNP
cd ~/annotation/variation/human/dbSNP
## https://www.ncbi.nlm.nih.gov/projects/SNP/
## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh38p2/
## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/
nohup wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz &
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz.tbi

Continue reading

十一 23

用R的bioconductor里面的stringDB包来做PPI分析

PPI本质上是根据一系列感兴趣的蛋白质或者基因(可以是几百个甚至上千个)来去PPI数据库里面找到跟这系列蛋白质或者基因的相互作用关系!

本次的主角是stringDB,顾名思义用得是大名鼎鼎的string数据库,
本来还以为需要自己上传自己的基因给这个数据库去做分析,没想到他们也开发了R包,主页见: http://www.bioconductor.org/packages/release/bioc/html/STRINGdb.html 而我比较喜欢用编程来解决问题,所以就学了一下这个包,非常好用!
它只需要一个3列的data.frame,分别是logFC,p.value,gene ID,就是标准的差异分析的结果。
然后用string_db$map函数给它加上一列是 string 数据库的蛋白ID,然后用string_db$add_diff_exp_color函数给它加上一列是color。
用string_db$plot_network函数画网络图,只需要 string 数据库的蛋白ID,如果需要给蛋白标记不同的颜色,需要用string_db$post_payload来把color对应到每个蛋白,然后再画网络图。
也可以直接用get_interactions函数得到所有的PPI数据,然后写入到本地,再导入到cytoscape进行画图

Continue reading

07

阅读文献下载原始reads之pacbio全基因组数据

目录

一:文献解读

二:NCBI搜索

三:构建脚本下载所有reads压缩包

四:用sra-toolkit解压

正文

一:文献解读

NCBI的原始reads数据里面最重要的就是SRP023104这样的字符串,能定位下载数据,也能定位原始文献。

文献地址如下

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4036802/

 

该研究小组还做了很多关于这个足部疾病的工作,其它文献地址如下

https://www.hgsc.bcm.edu/human/charcot-marie-tooth-project

二:NCBI搜索

http://www.ncbi.nlm.nih.gov/sra/SRX852869[accn]

SRX852869: PacBio whole genome sequencing of Homo sapien personal sample HS-1011
88 PACBIO_SMRT (PacBio RS II) runs: 14.4M spots, 89.6G bases, 319.9Gb downloads

这些数据都是可以下载的!!!

wxid_81t1aioe6cx122_1475841674027_6

如果有多余服务器资源,可以考虑下载玩玩这些数据!

这是很久以前的一个草稿帖子了,居然忘记发了,正好我开始有新的服务器了,一定要抽空玩一下三代测序啦!

04

生信人必学ftp站点之 dbsnp

这个数据库我也不想多解释了,也是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数据库

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

Continue reading

16

假基因资源中心

假基因是原来的能翻译成蛋白的基因经过各种突变导致丧失功能的基因。
比如
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包对芯片数据做差异分析

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

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

 

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语言,完全可以自己来!

05

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

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

随便一个位点,注释了如此多的信息!~~~

1       861389  .       C       T       5621.53 PASS    AC=4;AC_AFR=0;AC_AMR=0;AC_Adj=4;AC_EAS=0;AC_FIN=0;AC_Het=4;AC_Hom=0;AC_NFE=3;AC_OTH=1;AC_SAS=0;AF=3.300e-05;AN=121216;AN_AFR=10212;AN_AMR=11516;AN_Adj=119730;AN_EAS=8606;AN_FIN=6594;AN_NFE=65414;AN_OTH=890;AN_SAS=16498;BaseQRankSum=2.78;ClippingRankSum=-2.380e-01;DP=1488042;FS=7.913;GQ_MEAN=62.49;GQ_STDDEV=14.73;Het_AFR=0;Het_AMR=0;Het_EAS=0;Het_FIN=0;Het_NFE=3;Het_OTH=1;Het_SAS=0;Hom_AFR=0;Hom_AMR=0;Hom_EAS=0;Hom_FIN=0;Hom_NFE=0;Hom_OTH=0;Hom_SAS=0;InbreedingCoeff=-0.0004;MQ=59.70;MQ0=0;MQRankSum=0.198;NCC=409;QD=15.11;ReadPosRankSum=0.561;VQSLOD=0.392;culprit=FS;DP_HIST=373|361|219|102|34981|16744|5493|1367|498|210|121|54|32|18|13|9|3|3|3|4,0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|1|1|1;GQ_HIST=26|352|26|24|472|62|71|34|23|29|34|16|44468|8058|2176|2147|1116|370|365|739,0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4;CSQ=T|ENSG00000187634|ENST00000420190|Transcript|missense_variant|157|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|||ENSP00000411579||Q5SV95_HUMAN&I7FV93_HUMAN&A6PWC8_HUMAN|UPI000155D47C|deleterious(0)|probably_damaging(0.999)|2/7|||ENST00000420190.1:c.68C>T|ENSP00000411579.1:p.Pro23Leu||||||||||||||||||,T|ENSG00000268179|ENST00000598827|Transcript|missense_variant|211|211|71|G/R|Ggg/Agg||1||-1|AL645608.1|Clone_based_ensembl_gene||protein_coding|YES||ENSP00000471152||M0R0C9_HUMAN|UPI0000D61E05||probably_damaging(0.98)|6/6|||ENST00000598827.1:c.211G>A|ENSP00000471152.1:p.Gly71Arg||||||||||||||||||,T|ENSG00000187634|ENST00000437963|Transcript|missense_variant|128|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|||ENSP00000393181||Q5SV95_HUMAN&I7FV93_HUMAN|UPI000155D47B|deleterious(0)|probably_damaging(0.999)|2/5|||ENST00000437963.1:c.68C>T|ENSP00000393181.1:p.Pro23Leu||||||||||||||||||,T|ENSG00000187634|ENST00000342066|Transcript|missense_variant|151|68|23|P/L|cCg/cTg||1||1|SAMD11|HGNC|28706|protein_coding|YES|CCDS2.2|ENSP00000342313|SAM11_HUMAN|Q5SV95_HUMAN&I7FV93_HUMAN&A6PWC8_HUMAN|UPI0000D61E04|deleterious(0)|probably_damaging(0.999)|2/14|||ENST00000342066.3:c.68C>T|ENSP00000342313.3:p.Pro23Leu||||||||||||||||||,T||ENSR00000528850|RegulatoryFeature|regulatory_region_variant|||||||1||||||regulatory_region|||||||||||||||||||||||||||||||

头文件里面有对每一列的详细介绍,包括突变的标准格式

HGVS.c   --》ENST00000420190.1:c.68C>T

HGVS.p –》ENSP00000411579.1:p.Pro23Leu

还有该突变对蛋白功能的影响,包括sift和polyphen的打分~~~

不多说了,直接介绍该软件如何使用吧!

 

软件安装:

最新版是84:http://useast.ensembl.org/info/docs/tools/vep/script/vep_download.html

然后进入目录用perl的形式来安装这个软件:perl INSTALL.pl 即可

安装时其实有很多参数可以选择的,请仔细阅读介绍;http://useast.ensembl.org/info/docs/tools/vep/script/vep_download.html

 

前提是你已经安装好了两个模块!

perl -e 'use DBD::mysql'

perl -e 'use Archive::Extract'

如果不报错,就证明你已经安装过这些模块,如果报错,去搜索我以前关于perl模块的博客吧,不是很简单的事情。

By default the script will install the cache files in the ".vep" subdirectory of the user's home area. Using this option users can configure where cache files are installed.

我不想把cache文件放在默认的$HOME/.vep/下面,所以我安装的时候稍微做了更改

 

下载完了软件,接下来就要下载注释用的数据库啦!

它支持非常多的物种的注释,我这里拿人类做例子咯:ftp://ftp.ensembl.org/pub/release-82/variation/VEP/

我下载的是ftp里面的82 版本: wget ftp://ftp.ensembl.org/pub/release-82/variation/VEP/homo_sapiens_refseq_vep_82_GRCh37.tar.gz

有6.1G,所以会有点耗时~

下载完毕后直接用tar –zxvf解压即可使用啦!

我安装软件的时候指定了cache目录,而不是默认的$HOME/.vep/

 

Download the archive file for your species 

Extract the archive in your cache directory. By default the VEP uses $HOME/.vep/ as the cache directory, where $HOME is your UNIX home directory. 

mv homo_sapiens_vep_84.tar.gz ~/.vep/ cd ~/.vep/tar xfz homo_sapiens_vep_84.tar.gz

Run the VEP with the--cache option

所以要把下载的6.1G数据库放在我自己的cashe目录

如果你安装VEP的时候用的默认安装参数,就需要把自己下载的6.1G文件放在  ~/.vep/ 目录下面

参考:http://davetang.org/wiki2/index.php?title=VEP

输入数据:

它支持好几种输入格式数据:

  • BED: a simple tab-delimited format containing 3-12 columns of data. The first 3 columns contain the coordinates of the feature. If available, the VEP will use the 4th column of the file as the identifier of the feature.
  • GFF: a format for describing genes and other features. If available, the VEP will use the "ID" field as the identifier of this feature.
  • GTF: treated in an identical manner to GFF.
  • VCF: a format used to describe genomic variants. The VEP will use the 3rd column of the file as the identifier.
  • bigWig: a format for storage of dense continuous data. The VEP uses the value for the given position as the "identifier". Note that bigWig files contain their own indices, and do not need to be indexed by tabix.

Any other files can be easily converted to be compatible with the VEP; the easiest format to produce is a BED-like file containing coordinates and an (optional) identifier:

其实重点就是给出你的突变的坐标即可,在哪条染色体,什么位置!

我们可以拿snpEFF里面的example文件夹里面的数据来做测试。

 

 

运行命令:

可以直接进入安装目录(VEP_ensembl/ensembl-tools-release-84/scripts/variant_effect_predictor)运行那个主程序

variant_effect_predictor.pl -i example_GRCh37.vcf \

--cache --assembly GRCh37 \

--offline --force_overwrite

或者用全路径的形式去调用这个程序

参数非常复杂,详细介绍见:http://useast.ensembl.org/info/docs/tools/vep/script/vep_options.html

一般用标准参数就好啦,而且还有一些插件,其中我比较喜欢dbNSFP and LOFTEE plugins,这也是EXAC里面用过的。

3

 

结果解读:

这个非常复杂,对结果理解了多少,就是我们对软件理解了多少。

具体大家看readme吧,注释信息太多了,按需索取:

直接看EXAC(ExAC.r0.3.sites.vep.vcf.gz)文件里面近一亿条突变记录也能慢慢理解!

 

参考:http://gemini.readthedocs.io/en/latest/content/functional_annotation.html

05

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

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

官网是:http://snpeff.sourceforge.net/

1       889455  .       G       A       .       .        ## 假设我们的vcf文件里面记录的突变是这个,那么我们可以用snpEFF进行注释,注释得到的信息非常完全!

信息用|符号分割,所有很容易用脚本提取需要的信息

ANN=A|stop_gained|HIGH|NOC2L|ENSG00000188976|transcript|ENST00000327044|protein_coding|7/19|c.706C>T|p.Gln236*|756/2790|706/2250|236/749||,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000487214|processed_transcript||n.*865C>T|||||351|,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000469563|retained_intron||n.*878C>T|||||4171|,A|non_coding_exon_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000477976|retained_intron|5/17|n.2153C>T||||||;LOF=(NOC2L|ENSG00000188976|6|0.17);NMD=(NOC2L|ENSG00000188976|6|0.17)

包括突变类型是:non_coding_exon_variant

突变在各种转录本上面,在每个转录本的第几个碱基呀,哪个氨基酸的改变呀,氨基酸第几位呀!

标准突变表示形式是:

突变发生在NOC2L这个基因上面,它的ensembl 数据库ID是ENSG00000188976

 

其余的看头文件自己慢慢理解:

"Functional annotations: 'Allele | Annotation | Annotation_Impact | Gene_Name | Gene_ID | Feature_Type | Feature_ID | Transcript_BioType | Rank | HGVS.c | HGVS.p | cDNA.pos / cDNA.length | CDS.pos / CDS.length | AA.pos / AA.length | Distance | ERRORS / WARNINGS / INFO'

 

软件安装:

选择最新版软件下载:https://sourceforge.net/projects/snpeff/files/

wget https://sourceforge.net/projects/snpeff/files/snpEff_latest_core.zip

因为是java软件,unzip 解压之后就可以直接使用,当然前提是你有java平台。

1

输入数据:

首先下载用来做注释的数据库:java -jar snpEff.jar download GRCh37.75(自己选择需要的版本)

1

软件下载很快,但是数据库下载就需要一定时间啦,去喝杯咖啡吧。

然后软件本身会提供example文件,里面就是一堆各种各样的vcf数据,而且还提供了运行命令,非常简单(examples.sh) ,这些就是我们的输入数据啦!

运行命令:

运行也很简单:java -Xmx4G -jar snpEff.jar -i vcf -o vcf GRCh37.75 example.vcf > example_snpeff.vcf

指定输入输出格式都是vcf,然后指定刚才下载的必备数据库,然后输入输出文件即可!

也可以调用全路径,如果你写在脚本里面的话!

java -Xmx4G -jar path/to/snpEff/snpEff.jar \

-c path/to/snpEff/snpEff.config \

GRCh37.69 \

path/to/example.vcf > example_snpeff.vcf

 

结果解读:

这个非常复杂,对结果理解了多少,就是我们对软件理解了多少。

具体大家看readme吧,注释信息太多了,按需索取:

  1. chromosome_number_variation
  2. exon_loss_variant
  3. frameshift_variant
  4. stop_gained
  5. stop_lost
  6. start_lost
  7. splice_acceptor_variant
  8. splice_donor_variant
  9. rare_amino_acid_variant
  10. missense_variant
  11. inframe_insertion
  12. disruptive_inframe_insertion
  13. inframe_deletion
  14. disruptive_inframe_deletion
  15. 5_prime_UTR_truncation+exon_loss_variant
  16. 3_prime_UTR_truncation+exon_loss
  17. splice_branch_variant
  18. splice_region_variant
  19. splice_branch_variant
  20. stop_retained_variant
  21. initiator_codon_variant
  22. synonymous_variant
  23. initiator_codon_variant+non_canonical_start_codon
  24. stop_retained_variant
  25. coding_sequence_variant
  26. 5_prime_UTR_variant
  27. 3_prime_UTR_variant
  28. 5_prime_UTR_premature_start_codon_gain_variant
  29. upstream_gene_variant
  30. downstream_gene_variant
  31. TF_binding_site_variant
  32. regulatory_region_variant
  33. miRNA
  34. custom
  35. sequence_feature
  36. conserved_intron_variant
  37. intron_variant
  38. intragenic_variant
  39. conserved_intergenic_variant
  40. intergenic_region
  41. coding_sequence_variant
  42. non_coding_exon_variant
  43. nc_transcript_variant
  44. gene_variant
  45. chromosome

http://snpeff.sourceforge.net/VCFannotationformat_v1.0.pdf

 

28

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

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如何构建

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

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

基因组各种版本对应关系

我是受到了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

我前面写到过如何把数据库写到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数据库

在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

画基因结构图!

一个基因有一个甚至多个转录本!根据转录及翻译现象可以把一个基因人为的定义成多种结构!
首先自己查资料搞明白转录本,外显子和内含子,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

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

但是重要而且常用的芯片并不多!
一般分析芯片数据都需要把探针的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

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

华盛顿大学把所有的变异数据都用自己的方法注释了一遍,然后提供下载:
文献是: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的压缩包数据,我想,一般的公司或者单位都不会去用这个数据了!