(13)Hg19基因组的一些分析-生信菜鸟团博客2周年精选文章集
搞生信研究的，大部分数据都是针对于人类的，那么人类的参考基因组就不得不知了！
与hg19的突变相关的一些数据解释。
Hg19基因组的分析
R的bioconductor包TxDb.Hsapiens.UCSC.hg19.knownGene详解
下载地址我就不贴了，随便谷歌一下即可！
Genome Reference Consortium Human —》 GRCh3
Feb. 2009 (hg19, GRCh37)这个是重点
Mar 2006 assembly = hg18 = NCBI36.
May 2004 assembly = hg17 = NCBI35.
July 2003 assembly = hg16 = NCBI34
以前的老版本就不用看啦，现在其实都已经有hg38出来啦，GRCh38 (NCBI) and hg38(UCSC)
参考：http://age.wang.blog.163.com/blog/static/119252448201092284725460/
http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/
人的hg19基因组是3G的大小，因为一个英文字符是一个字节，所以也是30亿bp的碱基。
包括22条常染色体和X,Y性染色体及M线粒体染色体。
查看该文件可以看到，里面有很多的N，这是基因组里面未知的序列，用N占位，但是觉得部分都是A.T.C.G这样的字符，大小写都有，分别代表不同的意思。
然后我用linux的命令统计了一下里面这个文件的行数，
perl -lne 'END { print $. }’ hg19.fa
awk 'END { print NR }’ hg19.fa
wc -l hg19.fa
然后我写了一个脚本统计每条染色体的长度，42秒钟完成任务！
看来这个服务器的性能还是蛮强大的，读取文件非常快！
[perl]
while(<>){
chomp;
if (/>/){
if (exists $hash_chr{$key} ){
$len = length $hash_chr{$key};
print "$key => $len\n";
}
undef %hash_chr;
$key=$_;
}
else {
$hash_chr{$key}.=$_;
}
}
[/perl]
然后我用seed统计了一下hg19的词频（我不知道生物信息学里面的专业描述词语是什么）
我的程序耗费了42分钟才跑完，感觉我写的程序应该是没有问题的，让我吃惊的是总共竟然只有105万条独特的10bp短序列。然后我算了一下4的10次方，(⊙o⊙)…悲剧，原来只有1048576，之所以出现这种情况，是因为里面有N这个字符串，不仅仅是A.T.C.G四个字符。我用grep -v N seed10.txt |wc -l命令再次统计了一下，发现居然就是1048576，也就是说，任意A.T.C.G四个字符组成的10bp字符串短序列在人的基因组里面都可以找到！！！
然后我测试了一下，还是真是这样的，真是一个蛮有意思的现象。虽然我无法解释为什么，但是根据这个结果我们可以得知连续的A或者T在人类基因组里面高频出现，而连续的G或者C却很少！
如果我们储存这个10bp字符串的同时，也储存着它们在基因组的位置，那么就可以根据这个seed来进行比对，这就是blast的原理之一！
#这里是下载人类的已知基因的信息 (35.4 MB)
source(“http://bioconductor.org/biocLite.R”)
biocLite(“TxDb.Hsapiens.UCSC.hg19.knownGene”)
然后查看我们下载的这个包里面所包含的信息
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
> txdb
TxDb object:
| Db type: TxDb
| Supporting package: GenomicFeatures
| Data source: UCSC
| Genome: hg19
| Organism: Homo sapiens
| UCSC Table: knownGene
| Resource URL: http://genome.ucsc.edu/
| Type of Gene ID: Entrez Gene ID
| Full dataset: yes
| miRBase build ID: GRCh37
| transcript_nrow: 82960
| exon_nrow: 289969
| cds_nrow: 237533
| Db created by: GenomicFeatures package from Bioconductor
| Creation time: 2014-09-26 11:16:12 -0700 (Fri, 26 Sep 2014)
| GenomicFeatures version at creation time: 1.17.17
| RSQLite version at creation time: 0.11.4
| DBSCHEMAVERSION: 1.0
可以看到这个UCSC数据库发布的经典的hg19版本基因组所有的基因信息，共有237533个CDS，共有289969个外显子
然后我们可以用几个简单的函数提取信息
>genes(txdb)
可以看到genes函数可以提取23056个基因信息，还是一个Granges对象
>exons(txdb)
而用exons函数可以提取这个txdb对象的exons信息，共289969个exon
同理还有 transcripts函数可以提取转录本信息，共82960个转录本
还有cds函数，提取到237533个cds信息
#我们可以提取外显子的GRanges对象来具体看看，也可以用genes,transcripts,cds等函数
exon_txdb=exons(txdb)
seqnames(exon_txdb)返回一个class 'Rle’ [package “S4Vectors”] with 4 slots，有93个染色体信息，以及每条染色体上面有多少个外显子信息
ranges(exon_txdb)返回外显子的起始终止位点，长度，以及其它信息，也是一个对象class 'IRanges’ [package “IRanges”] with 6 slots
还有很多函数
strand(exon_txdb)返回外显子的正负链信息，要么在正链要么在负链
mcols(exon_txdb)返回exon的id编号，1到27750个
seqlengths(exon_txdb)返回每条染色体的长度信息
names,length
GRanges对象还有很多其它类型的操作，非常好玩的，split,shift,resize,flank,reduce,gaps,disjoin,coverage
其它求交集并集和都可以用,union,intersect,setdiff,pintersect,psetdiff
然后我们再回头看看我们的重点，txdb对象
transcriptsBy(txdb,by=”gene”)
这个是把我们对象按照gene来对转录本分组，可以看到，分成了23459个元素的list，其中第一个基因有两个转录本，也有一些基因只有一个转录本，甚至有些基因会有非常多的转录本，也可以用exonsBy,cdsBy来对它进行处理
每一个元素都是GRangesList对象，就是前面的GRanges对象，
seqnames(x)
ranges(x)
strand(x)
mcols(x, use.names=FALSE)
elementMetadata(x)
values(x)
seqinfo(x)
seqlevels(x)
seqlengths(x)
isCircular(x)
genome(x)
seqnameStyle(x)
seqnames(x)
score(x)
还有很多关于它的介绍
http://web.mit.edu/r_v3.0.1/lib/R/library/GenomicRanges/html/GRangesList-class.html
ttp://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的核苷酸序列