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

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-06

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信技能树

生物信息学技能面试题(第5题)-根据GTF画基因的多个转录本结构

可以下载各种gtf,从NCBI,ENSEMBL,UCSC,GENCODE都可以!(记住,你下载什么样的gtf就需要修改成什么样的代码!!!)本文来源于我的个人博...

3268
来自专栏hotqin888的专栏

request+goquery+mahonia实现自动抓取网页数据

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/hotqin888/article/det...

544
来自专栏前端小栈

CSS3 黑科技 - 内凹圆角 - 径向渐变实现

可以拿个白色圆盒子盖住方形盒子的大半边实现,但是这样,是不透明的,背景发生改变时,就要改遮盖盒子的颜色,或者背景是渐变,改起来更麻烦,或背景是图片等等,就直接不...

491
来自专栏菩提树下的杨过

“AS3.0高级动画编程”学习:第三章等角投影(上)

什么是等角投影(isometric)? 刚接触这个概念时,我也很茫然,百度+google了N天后,找到了一些文章: [转载]等角(斜45度)游戏与数学 [转载]...

1975
来自专栏Golang语言社区

Golang中image/jpeg包和image/png包用法

jpeg包实现了jpeg图片的编码和解码 func Decode(r io.Reader) (image.Image, error) //Decode读取一...

3944
来自专栏琦小虾的Binary

二维码生成原理及解析代码

二维码生成原理及解析代码 自从大街小巷的小商小贩都开始布满了腾讯爸爸和阿里爸爸的二维码之后,我才感觉到我大天朝共享支付的优越性。最近毕业论文写的差不多了,在入职...

5176
来自专栏吴小龙同學

Android之美化控件Shape

除了使用drawable这样的图片外,随着图片资源增多,程序也越来越大,今天我们自定义图形shape的方法。 1 2 3 4 5 6 7 8 9 10 11 1...

3606
来自专栏人工智能

手把手带你用机器学习写unity AI

2017unity机器学习社区挑战赛参加地址:https://connect.unity.com/challenges/ml-agents-1 打开下载的uni...

2087
来自专栏CDA数据分析师

翻译 | 简单而有效的EXCEL数据分析小技巧

介绍 我一直很欣赏EXCEL蕴藏的巨大能量。这款软件不仅具备基本的数据运算,还能使用它对数据进行分析。EXCEL被广泛运用到很多领域,例如:金融建模和商业预测。...

16810
来自专栏CDA数据分析师

【技能get】简单而有效的 EXCEL 数据分析小技巧

作者 CDA 数据分析师 我一直很欣赏 EXCEL 蕴藏的巨大能量。这款软件不仅具备基本的数据运算,还能使用它对数据进行分析。EXCEL 被广泛运用到很多领域...

2349

扫描关注云+社区