查了好久的bug，终于搞清楚问题所在了！因为要对基因进行reads计数，所以要拿到基因在基因组上面的染色体起始终止坐标，结果发现了个十分诡异的现象，很多基因有多个坐标，比如下面这个PTPRS 在hg38这个基因组版本，居然有两个定位，因为我是写程序格式化得到的坐标，所以我check了我的程序，http://www.biotrainee.com/thread-472-1-1.html 感兴趣的同学可以点开看看我的代码！
十一
10
Tag Archives: ENSEMBL
三
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，那就有点麻烦了：
需要选择一系列参数：
1. Navigate to http://genome.ucsc.edu/cgi-bin/hgTables2. 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 browser3. 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
四
10
对snp进行注释并格式化输出
前面我已经讲了如何用annovar来把vcf格式的snp进行注释，注释之后大概是这样的，每个snp位点的坐标，已经在哪个基因上面，都标的很清楚啦，。而且该突变是在哪个基因的哪个转录本的哪个外显子都一清二楚，更强大的是，还能显示是第几个碱基突变成第几个，同样氨基酸的突变情况也很清楚。
但是这样不是很方便浏览具体突变情况，所以我写了一个脚本格式化该突变情况。
理论上是应该要做出上面这个样子，突变氨基酸前后各12个氨基酸都显示出来，突变的那个还要标红色突出显示！但是颜色控制很麻烦，我就没有做。效果如下
实现这样的格式化输出有三个重点，首先是NM开头的refseq的ID号要转换为ensembl数据库的转录本ID号，还有找到该转录本的CDS序列，这个都需要在biomart里面转换，或者自己写脚本，然后就用脚本爬取即可！
代码如下
[perl]
open FH1,"NM2ensembl.txt";
while(<FH1>){
chomp;
@F=split;
$hash_nm_enst{$F[4]}=$F[1] if $F[4];
}
open FH2,"ENST.CDS.fa";
while($line=<FH2>){
chomp $line;
if ($line=~/>/) {$key = (split /\|/,$line)[1];}
else {$hash_nucl{$key}.=$line;}
}
open FH3,"ENST.protein";
while($line=<FH3>){
chomp $line;
if ($line=~/>/) {$key = (split /\|/,$line)[1];}
else {$hash_prot{$key}.=$line;}
}
open FH4,"raw.mutiple.txt";
$i=1;
while(<FH4>){
chomp;
@F=split;
@tmp=split/:/,$F[1];
/:exon(\d+):/;$exon=$1;
/(NM_\d+)/; $nm=$1;
$enst=$hash_nm_enst{$nm};
print "$i. $tmp[0] $F[0] the $exon -th exon(s) of $enst \n";
$i++;
$tmp[3]=~/(\d+)/;$num_nucl=$1;
$tmp[3]=~/>([ATCG])/;$mutation_nucl=$1;
$tmp[4]=~/(\d+)/;$num_prot=$1;
$sequence=$hash_nucl{$enst};
$num_up=3*$num_prot-39;
$out_nucl=substr($sequence,$num_up,75);
print "WT:$out_nucl\n ";
for(my $j=0; $j < (length($out_nucl) - 2) ; $j += 3)
{print ' ';print $codon{substr($out_nucl,$j,3)} ;print ' ';}
print "\n";
$mutation_pos=$num_nucl-$num_up-1;
substr($out_nucl,$mutation_pos,1,$mutation_nucl) if ((length $out_nucl) == 75 );
print "MU:$out_nucl\n ";
for(my $j=0; $j < (length($out_nucl) - 2) ; $j += 3)
{print ' ';print $codon{substr($out_nucl,$j,3)} ;print ' ';}
print "\n";
print "\n";
print "\n";
}
[/perl]
三
31
Ensembl数据库在线网页工具biomart简单教程
这个工具主要是针对不会bioperl不会API调取数据的生信纯菜鸟准备的，主要是方便大家批量研究某些感兴趣的基因，需要准备的数据就是基因名或者基因的ID号，能从该网站获取的资料非常多，可以是关于你的输入的基因名的各种数据库有的信息。
http://www.ensembl.org/biomart/
第一步：选取数据库，我一般选取人的
第二步，选择上传数据的格式
这个下拉框里面可以选取很多种格式，你随便张贴进去哪一种格式的基因ID都可以，也可以把做好的ID文件上传进去，批量获取基因信息。
我这里输入的是几个免疫基因。
第三步，选择下载数据的格式
首先可以选择你上传的gene的可以转换的各种ID
然后可以选择你上传的gene的各种序列
可以选择的信息非常多，基本上可以想到的转换在这里都能做！！！
但是，始终没有脚本方便，只适合不太懂编程的菜鸟使用！
然后点击result即可，看到结果还可以导出成txt文档，点击右上角的GO即可
三
26
一个基因的生信之旅
感觉大家对很多生物信息学的术语都不甚了解，我这里简单的从一个基因开始，扩展开来讲一讲生信数据库，及它相关的一些术语！
我要讲的基因是BRCA1，这是一个与乳腺癌以及卵巢癌都息息相关的基因。而BRCA1是它的英文缩写简称，也是通常学者们进行交流十它的名字。它的全称是breast cancer 1,每个基因都会有一个简称，比如下面这些，在human里面这些简称多大47732个，正常人都不会认识它们所有，只需要碰到了去数据库搜索即可，但是搞医疗健康的，必须熟悉癌症50基因。
这样的缩写简称其实弊端很多，单词毕竟是有限的，而且缩写也没有语义。所以NCBI给每个基因都定义了一个entrez ID号，是整数的排序，具体大家可以去看NCBI发的一篇文献，专门讲解了entrez ID号的好处。
1 A1BG
2 A2M
3 A2MP1
9 NAT1
10 NAT2
11 NATP
12 SERPINA3
13 AADAC
14 AAMP
这里我们来找一下我们的BRCA1这个基因在生物信息数据库里面的其它信息，在NCBI的ftp里面有一个文件是Homo_sapiens.gene_info里面包含着人类所以基因的全部信息
9606 首先这个基因在human上面的，而human被NCBI定义的taxid是9606
672 然后这个基因的被NCBI定义的entrez ID号是672
BRCA1 这个当然就是这个基因的英文缩写名称啦
- 这个表明这个基因在负链什么
BRCAI|BRCC1|BROVCA1|FANCS|IRIS|PNCA4|PPP1R53|PSCP|RNF53
这个可能是基因以前的名称，或者是在其它研究领域的一些名称。MIM:113705|HGNC:HGNC:1100|Ensembl:ENSG00000012048|HPRD:00218|Vega:OTTHUMG00000157426
这里面包含在它在其它数据库的信息，我们的NCBI用entrez ID号672来标识它，相应的ensembl数据用ensembl ID号ENSG00000012048来标识它，还有什么MIM数据库，HGNC数据库，Vega数据库我就不详细讲啦
17 17q21 这个说明它在human的17号染色体的位置信息
下面一堆都是这个基因的描述，它的功能等等。
breast cancer 1, early onset protein-coding BRCA1 breast cancer 1, early onset
O BRCA1/BRCA2-containing complex, subunit 1|Fanconi anemia, complementation group S|RING finger protein 53|breast and ovarian cancer susceptibility protein 1|breast and ovarian cancer sususceptibility protein 1|breast cancer type 1 susceptibility protein|protein phosphatase 1, regulatory subunit 53
20150201
这样我们就把好几个数据库给串起来了，也大致了解了一个基因的各种信息，但是，这样肯定是不够的。
接下来我们就不用BRCA1来称呼这个基因了，我们统一用NCBI定义entrez ID号672来称呼这个基因，当然用ensembl ID号ENSG00000012048也可以，它们都是比较通用的。
ENSG00000012048 672 这个基因在GO数据库里面可以找到67个功能信息，分别是以下
GO:0000151 GO:0000724 GO:0000724 GO:0000794 GO:0003677 GO:0003684 GO:0003713 GO:0003723 GO:0004842 GO:0005515 GO:0005634 GO:0005654 GO:0005694 GO:0005737 GO:0005886 GO:0006260 GO:0006281 GO:0006301 GO:0006302 GO:0006302 GO:0006349 GO:0006357 GO:0006359 GO:0006633 GO:0006915 GO:0006974 GO:0006978 GO:0007059 GO:0007098 GO:0008270 GO:0008274 GO:0008630 GO:0009048 GO:0010212 GO:0010575 GO:0010628 GO:0015631 GO:0016567 GO:0016874 GO:0019899 GO:0030521 GO:0030529 GO:0031398 GO:0031436 GO:0031572 GO:0031625 GO:0035066 GO:0035067 GO:0042127 GO:0042981 GO:0043009 GO:0043234 GO:0043627 GO:0044030 GO:0044212 GO:0045717 GO:0045739 GO:0045766 GO:0045892 GO:0045893 GO:0045893 GO:0045944 GO:0045944 GO:0046600 GO:0050681 GO:0051571 GO:0051572 GO:0051573 GO:0051574 GO:0051865 GO:0070512 GO:0070531 GO:0071158 GO:0071356 GO:0071681 GO:0085020 GO:1902042 GO:2000378 GO:2000617 GO:2000620
由于GO太多了，我简单讲几个
ubiquitin ligase complex
double-strand break repair via homologous recombination
double-strand break repair via homologous recombination
condensed nuclear chromosome
DNA binding
damaged DNA binding
transcription coactivator activity
RNA binding
ubiquitin-protein transferase activity
protein binding
都是描述这个基因的功能的。
到这里我们大致了解了这个基因的功能，但是还不够。
然后可以查到它有一下6个转录本，都有二十多个外显子。
NR_027676
NM_007300
NM_007299
NM_007298
NM_007297
NM_007294
在hg19这个参考基因组的起始终止坐标，还有各个外显子的起始终止坐标都能找到。
41196311,41199659,41201137,41203079,41209068,41215349,41215890,41219624,41222944,41226347,41228504,41234420,41242960,41243451,41247862,41249260,41251791,41256138,41256884,41258494,41267742,41276033,41277198
41197819,41199720,41201211,41203134,41209152,41215390,41215968,41219712,41223255,41226538,41228631,41234592,41243049,41246877,41247939,41249306,41251894,41256278,41256973,41258550,41267796,41276132,41277340
http://asia.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000012048;r=17:43044295-43125483
在ensembl里面关于这个基因的描述如下。
breast cancer 1, early onset [Source:HGNC Symbol;Acc:HGNC:1100]
BRCC1, FANCS, PPP1R53, RNF53
Chromosome 17: 43,044,295-43,125,483 reverse strand.
chromosome:GRCh38:CM000679.2:43044295:43125483:1
This gene has 29 transcripts (splice variants), 63 orthologues, is a member of 4 Ensembl protein families and is associated with 11 phenotypes.
RefSeq Gene ID 672
Uniprot identifiers: P38398
而且ensembl里面可以可视化这个基因的所有信息。
然后简单检索一下关于这个BRCA1基因的文献发表状况，居然多达2111篇文献，看来这个基因很火呀！！！
awk '{if ($1==9606 && $2==672) print }' gene2pubmed |wc
9606 672 1676470
9606 672 2001833
9606 672 2270482
9606 672 4506230
9606 672 7481765
9606 672 7545954
9606 672 7550349
9606 672 7795652
9606 672 7894491
9606 672 7894492
第三列1676470等编号是pubmed数据库的文献编号，可以直接找到关于这个基因的文献发表情况。
而直接在NCBI的pubmed数据库里面可以搜到多达11339篇文献。
esearch -db pubmed -query 'BRCA1'
Esearch这个程序是NCBI提供的，挺好用的，希望大家可以熟悉一下。
esearch -db pubmed -query 'BRCA1' | efetch -format docsum | xtract -pattern DocumentSummary -present Author -and Title -element Id -first "Author/Name" -element Title >BRCA1.pubmed
用这个代码，可以找到所有关于这个BRCA1基因的文献的作者及标题，这样可以统计在这个基因领域的研究者最出名的是谁。
至于这个基因的序列，及其转录本翻译的蛋白我就不列了，太长了，而且占位子
三
23
与hg19的突变相关的一些数据解释。
http://statgenpro.psychiatry.hku.hk/limx/kggseq/download/resources/
这个网站收集了大部分资料，我们就用它的，如果它倒闭了，大家再想办法去搜索吧。
其实这些文件都是基于NCBI以及UCSC和ensembl数据库的文件用一些脚本转换而来的，都是非常简单的脚本。
首先我们看看humandb/hg19_refGene.txt 这个文件，总共2.5万多个基因的共5万多个转录本。
19 可能是entrez ID，但是又不像。
NM_001291929 参考基因名
chr11 染色体
-
89057521
89223909
89059923
89223852
17 89057521,89069012,89070614,89073230,89075241,89088129,89106599,89133184,89133382,89135493,89155069,89165951,89173855,89177302,89182607,89184952,89223774, 89060044,89069113,89070683,89073339,89075361,89088211,89106660,89133247,89133547,89135710,89155150,89166024,89173883,89177400,89182692,89185063,89223909,
0
NOX4 基因的英文简称，通俗名
cmpl
cmpl
2,0,0,2,2,1,0,0,0,2,2,1,0,1,0,0,0,
然后我们看看hg19_snp141.txt这个文件
1 10229 A - .
1 10229 AACCCCTAACCCTAACCCTAAACCCTA - .
1 10231 C A .
1 10231 C - .
1 10234 C T .
1 10248 A T .
1 10250 A C .
1 10250 AC - .
1 10255 A - .
1 10257 A C .
1 10259 C A .
1 10291 C T .
1 10327 T C .
1 10329 ACCCCTAACCCTAACCCTAACCCT - .
1 10330 C - .
1 10390 C - .
1 10440 C A .
1 10440 C - .
1 10469 C G .
1 10492 C T .
1 10493 C A .
1 10519 G C .
1 10583 G A 0.144169
1 10603 G A .
1 10611 C G 0.0188246
1 10617 CGCCGTTGCAAAGGCGCGCCG -
里面记录了以hg19为参考的所有的snp位点。
585
ENST00000518655 基因的ensembl ID号
chr1 + 11873 14409 14409 14409
4 基因有四个外显子
11873,12594,13402,13660, 12227,12721,13655,14409, 在基因的四个外显子的坐标
0
DDX11L1 基因的通俗英文名
none none -1,-1,-1,-1,
CTTGCCGTCAGCCTTTTCTTT·····gene的核苷酸序列
三
10
