(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 条评论
登录 后参与评论

相关文章

来自专栏岑玉海

hbase源码系列(十三)缓存机制MemStore与Block Cache

这一章讲hbase的缓存机制,这里面涉及的内容也是比较多,呵呵,我理解中的缓存是保存在内存中的特定的便于检索的数据结构就是缓存。 之前在讲put的时候,put是...

3287
来自专栏张高兴的博客

张高兴的 Windows 10 IoT 开发笔记:BH1750FVI 光照度传感器

2916
来自专栏Golang语言社区

Go 语言实现 MapReduce 框架

MapReduce 是 Google 提出的一个软件架构,用于大规模数据集(大于1TB)的并行运算。简而言之,就是将任务切分成很小的任务然后一个一个区的执行最后...

3566
来自专栏Java技术分享

Java基础常见英语词汇

Java基础常见英语词汇(共70个) OO:object-oriented ,面向对象

2047
来自专栏流媒体

命令模式

这里我们用一个场景来描述:去川湘阁饭店点一份剁椒鱼头和宫保鸡丁。 这里我们把整个关键流程写出来

874
来自专栏Crossin的编程教室

【每周一坑】特殊回文数

Thank God It’s Friday! again! 又到了惬意的周末了,祝大家周末愉快,闲暇之余记得练习编程技能啊!另外,上海的朋友可以安心在家写代码了...

3287
来自专栏前端笔记

【 测试文章 第一篇章 】Hello World 初始程序

Hello World ,中文意思:你好,世界。世界上的第一个程序就是Hello World,由Brian Kernighan创作。 中文名:你好,世界。 ...

30410
来自专栏小鹏的专栏

堆排序

堆排序排序是优秀的算法,但是在实际应用中,快速排序的性能一般会优于堆排序, 尽管如此,堆排序仍然有很多应用,例如:作为高效的优先队列,最大优先队列应用于共享计算...

1896
来自专栏安恒网络空间安全讲武堂

Writeup丨国赛线上初赛解题第二波~

1154
来自专栏数据结构与算法

1001 舒适的路线 2006年

1001 舒适的路线 2006年 时间限制: 2 s 空间限制: 128000 KB 题目等级 : 钻石 Diamond 题目描述 Descrip...

2364

扫码关注云+社区