(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的核苷酸序列