|
搞生信研究的，大部分数据都是针对于人类的，那么人类的参考基因组就不得不知了！ 与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./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:///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./ | 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./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的核苷酸序列 |
|
|
