Category Archives: 生信基础

基础软件，基础数据库，基础数据格式

七 17

转载-VCF格式详解

Posted on by

CHROM(chromosome):染色体

POS - position:参考基因组variant碱基位置，如果是INDEL（插入缺失），位置是INDEL的第一个碱基位置

ID - identifier: variant的ID。比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’.'表示其为一个novel variant。

REF - reference base(s):参考碱基，染色体上面的碱基，必须是ATCGN中的一个，N表示不确定碱基

ALT - alternate base(s):与参考序列比较发生突变的碱基

QUAL - quality: Phred格式(Phred_scaled)的质量值，表 示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。

FILTER - _filter status: 使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为”PASS”;若variant不可靠，则该项不为”PASS”或”.”。

INFO - additional information:  这一行是variant的详细信息，具体如下：

DP-read depth：样本在这个位置的reads覆盖度。是一些reads被过滤掉后的覆盖度。          DP4:高质量测序碱基，位于REF或者ALT前后

MQ：表示覆盖序列质量的均方值RMS Mapping Quality

FQ：phred值关于所有样本相似的可能性

AF1： AF(Allele Frequency) 表示Allele的频率，AF1为第一个ALT allele 发生频率的可能性评估

AC1：AC表示Allele（等位基因）的数目,AC1为对第一个ALT allele count的最大可能性评估

AN：AN(Allele Number) 表示Allele的总数目

IS：插入缺失或部分插入缺失的reads允许的最大数量

AC：AC(Allele Count) 表示该Allele的数目

G3：ML 评估基因型出现的频率

HWE：chi^2基于HWE的测试p值和G3

CLR：在受到或者不受限制的情况下基因型出现可能性log值

UGT：最可能不受限制的三种基因型结构

CGT：最可能受限制三种基因型的结构

PV4：四种P值得误差，分别是（strand、baseQ、mapQ、tail distance bias）

INDEL：表示该位置的变异是插入缺失

PC2：非参考等位基因的phred（变异的可能性）值在两个分组中大小不同

PCHI2：后加权chi^2，根据p值来测试两组样本之间的联系

QCHI2：Phred scaled PCHI2.

PR：置换产生的一个较小的PCHI2

QBD：Quality by Depth，测序深度对质量的影响

RPB：序列的误差位置（Read Position Bias）

MDV：样本中高质量非参考序列的最大数目

VDB：Variant Distance Bias，RNA序列中过滤人工拼接序列的变异误差范围

GT：样品的基因型（genotype）。两个数字中间用’/'分 开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。因此： 0/0 表示sample中该位点为纯合的，和ref一致； 0/1 表示sample中该位点为杂合的，有ref和variant两个基因型； 1/1 表示sample中该位点为纯合的，和variant一致。

GQ：基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值，表示在该位点该基因型存在的可能性；该值越高，则Genotype的可能性越 大；计算方法：Phred值 = -10 * log (1-p) p为基因型存在的概率。

GL：三种基因型（RR RA AA）出现的可能性，R表示参考碱基，A表示变异碱基

DV：高质量的非参考碱基

SP：phred的p值误差线

PL：指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。和之前不一致，该值越大，表明为该种基因型的可能 性越小。 Phred值 = -10 * log (p) p为基因型存在的概率。

FORMAT 和 BC1-1-base.sorted.bam：这两行合起来提供了’ BC1-1-base′这个sample的基因型的信息。’ BC1-1-base′代表这该名称的样品，是由BAM文件中的@RG下的 SM 标签决定的。

七 17

转-基因突变种类大全

Posted on by

突变（Mutation, 即基因突变）:在生物学上的含义，是指细胞中的遗传基因（通常指存在于细胞核中的脱氧核糖核酸）发生的改变。它包括单个碱基改变所引起的点突变，或多个碱基的缺失、重复和插入。原因可以是细胞分裂时遗传基因的复制发生错误、或受化学物质、辐射或病毒的影响。

以功能分类:

失去功能的突变Loss-of-function mutations

：失去功能的突变是指发生的突变会造成基因完全地失去活性，原因可分成两类。一类是由于基因被删除或是调控基因表现的过程受到影响让基因不表现，另一种则是由于基因本身受到影响，使得基因的产物蛋白质失去功能。又称剔除突变（null mutations）或是敲除突变（knockout mutations）。

次形态突变Time form mutation：此种突变会使基因的表现或是基因产物的活性减弱，但不会消失。

超形态突变hypermorphic mutations：此种突变与次形态突变相反，会使基因的表现加强。

获得功能的突变gain-of-function mutation：获得功能的突变是指发生的突变让原本应该是不表现的基因产生活性，进而影响细胞功能，这样的突变多半需要染色体程度的突变较有可能产生，而最常发生获得功能的突变就是癌细胞。

以突变机理分类:

  1. 点突变point mutation:DNA序列中涉及单个核苷酸或碱基的变化称为点突变。 通常有两种情况：一是一种碱基或核苷酸被另一种碱基或核苷酸所替换；二是一个碱基的插入或缺失。

                   （1）沉默突变silent mutation

当点突变发生在基因及其调控序列之外，或使基因序列内一种密码子变成编码同一种氨基酸的另一种同义密码子时，不会改变生物个体的基因产物，因而不引起性状变异。不引起生物性状变异的突变称为沉默突变。

                   （2）错义突变missense mutation

指由于某个碱基对的改变，使编码一种氨基酸的密码子变成编码另外一种氨基酸的密码子，结果使构成蛋白质的数百上千个氨基酸中有一个氨基酸发生变化。（实例：镰刀形细胞贫血症）

                   （3）移码突变frameshift mutation

指在DNA链上，有时一个或几个非3的整数倍的碱基的插入或缺失，往往产生比碱基替换突变更严重的后果。 这种插入或缺失突变会造成阅读框的改变，翻译过程中其下游的三联密码子都被错读，产生完全错误的肽链或肽链合成提前终止。这种插入或缺失突变又称为移码突变。

                   （4）无义突变nonsense mutation

是指当点突变使一个编码氨基酸的密码子变成终止子时，则蛋白质合成进行到该突变位点时会提前终止，结果产生一个较短的多肽链或较小的蛋白质。

  1. 大突变

大突变是可能涉及整个基因以至多个基因的一长段DNA序列的改变，大突变常常导致染色体畸变。

（1）缺失:指DNA分子丢失一段碱基序列。（染色体缺失）(Deletion)

（2）插入:指DNA分子的正常序列中插入一段DNA序列。(Insertion②)

（3）重排:重排包括某段DNA序列的重复（duplication），倒位（inversion），易位（translocation）等。

 

七 13

使用Seq2HLA进行HLA分型

Posted on by

基于高通量测序数据进行HLA分型的软件挺多的，比较老的有三个，作者分别是Boegel et al.Kim et al.Major et al.，然后他们都被OptiType这个软件的作者被批评了，我这里先介绍Kim et al的seq2HLA使用方法，以下是它的一些链接。

功能概述

seq2HLA is a computational tool to determine Human Leukocyte Antigen (HLA) directly from existing and future short RNA-Seq reads. It takes standard RNA-Seq sequence reads in fastq format as input, uses a bowtie index comprising known HLA alleles and outputs the most likely HLA class I and class II types, a p-value for each call, and the expression of each class.

软件简介

Type of tool     Program

Nature of tool          Standalone

Operating system   Unix/Linux, Mac OS X

Language        Python, R

Article     (Boegel et al., 2013) HLA typing from RNA-Seq sequence reads. Genome medicine.

PubMed http://www.ncbi.nlm.nih.gov/pubmed/23259685

URL          https://bitbucket.org/sebastian_boegel/seq2hla

源代码，下载并安装

https://bitbucket.org/sebastian_boegel/seq2hla/src

http://tron-mainz.de/tron-facilities/computational-medicine/seq2hla/

第一版是这样的

image001

第二版是这样的

image002

只有第二版才支持gz压缩包格式的fastq，而且不需要指定length了

其中reference文件夹下面的是发布这个软件的团体已经制备好来的HLA库文件

image003

下载即可使用，前提是你的系统其它环境都OK

用法：

python seq2HLA.py -1 <readfile1> -2 <readfile2> -r "<runname>" [-p <int>]* [-3 <int>]**

image004

很简单，-1和-2指定我们的双端测序数据即可，可以是压缩包格式的（自动调用gzip），-r的输出目录，会输出7个文件，需要一个个解读，-p指定线程数给bowtie用的，-3是指定需要trim几个低质量碱基。

但是运行这个软件的要求非常多，需要安装好python和R，而且还有版本限制，需要安装好biopython而且还必须是双端测序，而且当前文件夹下面的reference文件夹下面必须有参考基因组的bowtie索引，而且系统必须安装好了bowtie，还需要在快捷方式里面！

我这里用的是第二版的

image006

所以，我用的也是第二版改进的命令。非常好用,我这里用的是一个外显子测序数据，是hiseq2500测的PE100

python seq2HLA.py -1 ../../6-exon/PC3-1.read1_Clean.fastq.gz -2 ../../6-exon/PC3-1.read2_Clean.fastq.gz -r PC3

貌似输出文件太多了一点

#Output:#The results are output to stdout and to textfiles. Most important are:

#i) <prefix>-ClassI.HLAgenotype2digits => 2 digit result of Class I

#ii) <prefix>-ClassII.HLAgenotype2digits => 2 digit result of Class II

#iii) <prefix>-ClassI.HLAgenotype4digits => 4 digit result of Class I

#iv) <prefix>-ClassII.HLAgenotype4digits => 4 digit result of Class II

#v) <prefix>.ambiguity => reports typing ambuigities (more than one solution for an allele possible)

#vi) <prefix>-ClassI.expression => expression of Class I alleles

#vii) <prefix>-ClassII.expression => expression of Class II alleles

根据文献，我简单看了一下，文件的确好复杂，不过我们只需要看输出日志即可

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

A       A*68   7.287148e-05   A*24   0.03680272

B       B*52   0.1717737       B*53   0.3952319

C       C*12   0.03009331     hoz("C*14")     0.6783964

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassI.bowtielog

A: 7.93 RPKM

C: 9.75 RPKM

B: 8.35 RPKM

The digital haplotype is written into BC1-1/BC1-1-ClassI.digitalhaplotype3

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!A     A*68:01 7.287148e-05   A*24:02 0.03680272

!B     B*52:01 0.1717737       B*53:01'       0.6542288

!C     C*12:02 0.03371717     C*12:02 0.6783964

上面的HLA的class I的数据结果

接下来是class II的数据结果，是不是很简单呀！

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

DQA     DQA1*01 0.1511134       DQA1*02 0

DQB     DQB1*02 0.02321615     DQB1*05 0.42202

DRB     DRB1*15 2.595144e-05   DRB1*07 0.321219

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassII.bowtielog

DQB1: 4.47 RPKM

DRB1: 5.59 RPKM

DQA1: 0.44 RPKM

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!DQA   DQA1*01:02'     0.1511134       DQA1*02:01     0.0

!DQB   DQB1*02:01'     0.02321615     DQB1*05:01     0.42202

!DRB   DRB1*15:02'     2.595144e-05   DRB1*07:01     0.321219

七 06

GATK使用注意事项

Posted on by

GATK这个软件在做snp-calling的时候使用率非常高，因为之前一直是简单粗略的看看snp情况而已，所以没有具体研究它。

这些天做一些外显子项目以找snp为重点，所以想了想还是用起它，报错非常多，调试了好久才成功。

所以记录一些注意事项!

GATK软件本身是受版权保护的，所以需要申请才能下载使用，大家自己去broad institute申请即可。

下载软件就可以直接使用，java软件不需要安装，但是需要你的机器上面有java，当然软件只是个开始，重点是你还得下载很多配套数据，https://software.broadinstitute.org/gatk/download/bundle（ps:这个链接可能会失效，下面的文件，请自己谷歌找到地址哈。），而且这个时候要明确你的参考基因组版本了！！！ b36/b37/hg18/hg19/hg38，记住b37和hg19并不是完全一样的，有些微区别哦！！！
Continue reading →

六 24

Genomemapper软件使用说明书

Posted on by

 我以前一直以为有了bwa跟bowtie，没什么必要用其它的alignment软件，直到我碰到了高插入删除的helicos三代测序数据，我才发现，这个古董软件genomemapper居然大有用武之地了。

一．下载并且安装该软件

这是最新版本了

Release 0.4.4 2012-10-30 source code including documentation

Wget http://1001genomes.org/data/software/genomemapper/genomemapper_0.4.4/genomemapper-0.4.4.tar.gz

这个软件安装很简单，解压进入目录，make一下即可

image001

看到make完了之后就会多了两个软件，其中一个是用来构建参考基因组索引，一个用来比对的！

二．准备数据

既然是比对软件，那么肯定是一个参考基因组，一个测序的fastq原始文件咯

当然这个软件比较奇葩，它还支持Multi-FASTA, FASTQ2 or SHORE flat file format,

三、比对命令

这里要分两步走，首先是构建参考基因组的索引，然后才是比对

/home/jmzeng/bio-soft/genomemapper-0.4.4/gmindex \

-i BRCA1.fa -x BRCA1.idx -t BRCA1.meta

首先构建索引，种子长度就用默认的12即可，然后构建完索引如下。

image002

然后进行比对即可

/home/jmzeng/bio-soft/genomemapper-0.4.4/genomemapper \

-i BRCA1.fa -q SRR258835.fastq -M 4 -G 2 -E 4 -o mapped_reads.fl -u unmapped_reads.fl

成功比对的都输出到了mapped_reads.fl -这个文件，未比对上的在unmapped_reads.fl

我有12344条序列，成功比对的只有5276条，但是如果我用精确比对的算法，只有一千五百条是可以比对的，所以用这个允许4个mismatch和2个gap的比对算法，大大提高了比对率。

然后我修改了比对参数可以达到5605,5654,5696的提升。但是没有质的飞跃，估计本身我的这种helicos测序数据错误率就太可怕了。

四，输出结果解读

image004

这个是很规则的tab键分割的文本字符，我就不解读了，大家看readme

五 28

菜鸟团第二次作业的部分答案

Posted on by

> library(org.Hs.eg.db)

载入需要的程辑包：AnnotationDbi载入需要的程辑包：stats4载入需要的程辑包：GenomeInfoDb载入需要的程辑包：S4Vectors载入需要的程辑包：IRanges载入程辑包：‘AnnotationDbi’The following object is masked from ‘package:GenomeInfoDb’:     species载入需要的程辑包：DBI

 

1、人共有多少个entrez id的基因呢？

x <- org.Hs.egENSEMBLTRANS

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(x)

[1] 47721

可知共有47721个基因都是有entrez ID号的

2、能对应转录本ID的基因有多少个呢？

length(xx)

[1] 20592

可以看到共有20592个基因都是有转录本的！

2、能对应ensembl的gene ID的基因有多少个呢？

x <- org.Hs.egENSEMBL

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 26019

可以看到只有26019是有ensembl的gene ID的

3、那么基因对应的转录本分布情况如何呢？

table(unlist(lapply(xx,length)))

菜鸟团第二次作业的部分答案863

可以看出绝大部分的基因都是20个转录本一下的，但也有极个别基因居然有高达两百个转录本，很可怕！

4、那么基因在染色体的分布情况如何呢？

x <- org.Hs.egCHR

# Get the entrez gene identifiers that are mapped to a chromosome

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 47232

可以看到有接近五百个基因居然是没有染色体定位信息的！！！

table(unlist(xx))

用barplot函数可视化一下，如图

 

菜鸟团第二次作业的部分答案1209

6、那么有多多少基因是有GO注释的呢？

x <- org.Hs.egGO

# Get the entrez gene identifiers that are mapped to a GO ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(xx)

[1] 18229

> length(x)

[1] 47721

可以看到只有18229个基因是有go注释信息的。

那么基因被注释的go的分布如何呢？

菜鸟团第二次作业的部分答案1477

可以看到大部分的基因都是只有30个go的，但是某些基因特别活跃，高达197个go注释。

还有kegg和omin数据库的我就不写了！

五 28

实战R语言bioconductor的seqinr包探究人的所有转录本的性质

Posted on by

首先安装这个包

source("http://bioconductor.org/biocLite.R")

biocLite("seqinr")

然后加载包，并读取我们的CDS.fa文件

library("seqinr")

human_cds=read.fasta("CDS.fa")

#这一个步骤非常耗时间，可能是因为我们的转录本文件有十万多个转录本的原因吧

str(human_cds) #查看可知读入了一个list，其中每个转录本都是list的一个元素

List of 100778

$ ENST00000415118:Class 'SeqFastadna'  atomic [1:8] g a a a ...

.. ..- attr(*, "name")= chr "ENST00000415118"

.. ..- attr(*, "Annot")= chr ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene tran"| __truncated__

$ ENST00000448914:Class 'SeqFastadna'  atomic [1:13] a c t g ...

.. ..- attr(*, "name")= chr "ENST00000448914"

.. ..- attr(*, "Annot")= chr ">ENST00000448914 havana_ig_gene:known chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene tran"| __truncated__

对list的每个元素都有几种函数可以处理得到信息：

Length,table,GC,count

其中count函数很有趣，数一数序列里面的这些组合出现的次数

count(dengueseq, 1)

count(dengueseq, 2)接下来我们随机取human_cds这个list的一个元素用这几个函数对它处理一下

> tmp=human_cds[[1]]

> tmp

[1] "g" "a" "a" "a" "t" "a" "g" "t"

attr(,"name")

[1] "ENST00000415118"

attr(,"Annot")

[1] ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene"

attr(,"class")

[1] "SeqFastadna"

再看看函数的结果

> length(tmp)

[1] 8

> table(tmp)

tmp

a g t

4 2 2

> GC(tmp)

[1] 0.25

> count(tmp,1)

 

a c g t

4 0 2 2

> count(tmp,2)

 

aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt

2  0  1  1  0  0  0  0  1  0  0  1  1  0  0  0

>

还是挺好用的，接下来我们应用R的知识来对着十万多个转录本进行一些简单的总结

human_cds_length=unlist(lapply(human_cds,length))

human_cds_gc=unlist(lapply(human_cds,GC))

这样就得到了所有转录本的长度和GC含量信息

然后我们简单统计一下，并画几个图表吧！

> summary(human_cds_length)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

3     366     699    1132    1425  108000

> summary(human_cds_gc)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

0.1467  0.4577  0.5285  0.5264  0.5932  0.8917

可以看到还是有很多很极端的转录本的存在！

最长的转录本也不过10k，而我记得最长的基因高达8M，看了内含子远大于外显子呀。

但是GC含量有很多高于80%，这些基因在二代测序的研究中是一个盲区。

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2075

这些极端基因可以通过biomaRt等包进行注释，得到gene名和功能信息。

 

hist(human_cds_gc)

hist(log10(human_cds_length))

GC含量分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2177

长度分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2186

附表：

http://www.bioinformatics.org/sms/iupac.html 所有字符的碱基氨基酸意义表格

 

五 21

Biostrings包简介

Posted on by

首先讲讲它的对象

有下面几个字符串对象BString, DNAString, RNAString and AAString可以通过以下代码构造它们：

b <- BString("I am a BString object")

d <- DNAString("TTGAAAA-CTC-N")

这两个对象的区别是DNAstring对象对字符串的要求严格一些，只有IUPAC字符和+-字符可以。

对构造好的对象可以通过下标来取子字符串对象，也可以通过subseq来取，但是子字符串仍然是数据对象，只有通过toString函数才能把它们转化成字符串。

用length(dd2)和nchar(toString(dd2))都可以找到我们Biostrings对象的长度。但是后者速度会很慢。

Views(RNAString("AU"), start=0, end=2)这个函数可以把string对象任意截取成list

start, end and width可以作用于我们截取的list，判断list里面的元素在原来的string对象上面的起始终止及长度信息。

 

接下来讲这个包带有的一个比对函数！

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] succ--eed subject: [1] supersede score: -33.99738

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",type = "local")

Local PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su subject: [1] su score: 5.578203

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",gapOpening = 0, gapExtension = 1)

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su-cce--ed subject: [1] sup--

可以看出这个比对函数可以调整的参数实在是太多了，而且改变参数之后比对情况大不一样，还有很多参数就不一一细讲了。

这个比对结果可以赋值给一个变量，保存比对的对象。

psa1 <- pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

class(psa1)

summary(psa1)

class(pattern(psa1))

class(summary(psa1))

score(psa2)

还可以自己构建打分矩阵来进行比对。

submat <-

+ matrix(-1, nrow = 26, ncol = 26, dimnames = list(letters, letters))

diag(submat) <- 0

Biostrings包简介1454

psa2 <-pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",substitutionMatrix = submat, gapOpening = 0, gapExtension = 1)

我们的包还自带了两个非常流行的氨基酸比对矩阵PAM和BLOSUM

ls("package:Biostrings")可以查看这个包所有的对象。

data(package="Biostrings")可以查看这个包所有的数据对象

还有很多其它函数

还可以去除adaptor，挺好玩的

既然有配对比对函数，那么就有多重比对函数！

我们可以读取clustaW, Phylip and Stolkholm这几种不同的比对结果文件来构造多重比对对象。

library(Biostrings)这个包里面自带了两个文件，我们可以示范一下构建对象。

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata", "msx2_mRNA.aln", package="Biostrings"), format="clustal")

phylipMAlign <- readAAMultipleAlignment(filepath = system.file("extdata","Phylip.txt", package="Biostrings"),format="phylip")

 

对构造好的多重比对对象就可以构建进化树啦，代码如下！

sdist <- stringDist(as(origMAlign,"DNAStringSet"), method="hamming")

> clust <- hclust(sdist, method = "single")

> pdf(file="badTree.pdf")

> plot(clust)

> dev.off()

Biostrings包简介2345

五 21

Bioconductor的DO.db包介绍

Posted on by

Bioconductor的包都是同样的安装方法：

source("http://bioconductor.org/biocLite.R");biocLite("DO.db")

还有GO.bd包是完全一模一样的规则！！！

加载这个包可以发现它依赖于好几个其它的包，这也是我比较喜欢R的原因，它会自动把它需要的包全部安装加载进来，不需要自己一个个调试！

> library(DO.db)

载入需要的程辑包：AnnotationDbi

载入需要的程辑包：stats4

载入需要的程辑包：GenomeInfoDb

载入需要的程辑包：S4Vectors

载入需要的程辑包：IRanges

载入需要的程辑包：DBI

> help(DO.db)

> ls("package:DO.db")

[1] "DO"          "DO_dbconn"   "DO_dbfile"   "DO_dbInfo"   "DO_dbschema" "DOANCESTOR"  "DOCHILDREN"  "DOID"        "DOMAPCOUNTS"

[10] "DOOBSOLETE"  "DOOFFSPRING" "DOPARENTS"   "DOSYNONYM"   "DOTERM"      "DOTerms"     "Secondary"   "show"        "Synonym"

[19] "Term"

这个包里面有19个数据对象！都是比较高级的S4对象。

比如我们可以拿DOTERM[1:10]这个小的数据对象来做例子！example=DOTERM[1:10]

因为example是一个高级对象，所以无法直接查看，需要用as.list方法来查看

> as.list(example)

$`DOID:0001816`DOID: DOID:0001816Term: angiosarcomaSynonym: DOID:267Synonym: DOID:4508Synonym: "hemangiosarcoma" EXACT []Secondary: DOID:267Secondary: DOID:4508

~~~~~~~~~~~~共十个DO条目

对每一个DO条目来说都有DOID,Term,Synony这些函数可以取对应的值。

下面是对DO的有向无环图的数据解读

xx <- as.list(DOANCESTOR)可以查看每个DO与它所对应的上级条目DO，每个DO都会有不止一个的上级DO。

xx <- as.list(DOPARENTS)可以查看每个DO与它所对应的父条目DO，每个DO都有且只有一个父DO。

xx <- as.list(DOOFFSPRING)可以查看每个DO与它所对应的下级DO的关系列表，大多数DO都不止一个子条目DO，所有的下级DO都会列出。

xx <- as.list(DOCHILDREN)以查看每个DO与它所对应的子条目DO的关系列表，大多数DO都不止一个子条目DO。

还有Lkeys(DOTERM)可以查看数据库里面的所有的DO条目的ID号

> head(keys(DOTERM))

[1] "DOID:0000000" "DOID:0001816" "DOID:0002116" "DOID:0014667" "DOID:0050004" "DOID:0050012"

dbmeta(GO_dbconn(), "GOSOURCEDATE")

可以查看这个DO库的制备时间

> dbmeta(DO_dbconn(), "DOSOURCEDATE")

[1] "20140417"

五 21

SAMStat软件使用说明书

Posted on by

这个软件是对我们的比对结果（通常是bwa,bowtie,tophat,hisat,star）bam或者sam来进行一个可视化的总结，类似于fastqc对我们的fastq测序结果做一个可视化总结，非常好用。

一．下载并安装该软件

软件主页是http://samstat.sourceforge.net/ 里面对该软件进行非常详细的说明

包括installation和usage，我这里简单的翻译一下。

Wget http://liquidtelecom.dl.sourceforge.net/project/samstat/samstat-1.5.tar.gz

解压开看里面的readme有介绍如何安装这个软件

Unpack the tarball:

bash-3.1$ tar -zxvf samstat-XXX.tar.gz

bash-3.1$ cd samstat

bash-3.1$ ./configure

bash-3.1$ make

bash-3.1$ make check

bash-3.1$ make install

如果用root命令就可以直接用samstat啦

如果没有root权限，安装的时候稍微有点不同

./configure  --prefix=/home/jmzeng/my-bin/

make

make install

很简单的

二，数据，就是我们的bam文件啦

三，运行命令

SAMStat软件使用说明书710

四，结果

简单看看samtools flagstat 740WT1.bam  的结果

19232378 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

18846845 + 0 mapped (98.00%:-nan%)

0 + 0 paired in sequencing

0 + 0 read1

0 + 0 read2

0 + 0 properly paired (-nan%:-nan%)

0 + 0 with itself and mate mapped

0 + 0 singletons (-nan%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

然后再看看我们的samstat的结果！

740WT1.bam.samstat.html

一个网页，非常丰富的内容

SAMStat软件使用说明书1168

内容太多了，我懒得解释了

见软件说明书http://davetang.org/wiki/tiki-index.php?page=SAMStat

五 08

GWAS研究现状及资源下载

Posted on by

 GWAS研究是非常火的，NHGIR还专门为它开辟了专栏来介绍，下面这个图片也是来自于NHGIR组织，是GWAS近年来发表文章的状况。

图片1

可以在该文章上面下载这个所有的数据

wget http://www.genome.gov/admin/gwascatalog.txt

截至目前为止。2015年5月8日21:08:34

这个文档有19603行的数据，但是只有2113篇pubmed文献，共涉及到七千多个基因

有293种杂志都发过GWAS的文章，总共有2113篇文献，发表关联分析突变位点最多的是这篇文献23251661 在PLoS One杂志上面，共 949个rs突变

杂志排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}' |cut -f 2 |perl -alne '{$hash{$_}++}END{print "$_\t$hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
Hum Genet 41
Am J Hum Genet 62
Mol Psychiatry 64
PLoS One 132
PLoS Genet 145
Hum Mol Genet 168
Nat Genet 397
文章的rs突变点排序
cut -f 2,5 gwascatalog.txt |perl -alne '{$hash{$_}++}END{print "$_ $hash{$_}" foreach sort {$hash{$a} <=> $hash{$b}} keys %hash}'
24324551 PLoS One 241
24097068 Nat Genet 245
24816252 Nat Genet 299
23382691 PLoS Genet 699
23251661 PLoS One 949

数据打开如下：

我取了表头和第一行数据，然后把它转置了，这样方便查看

Date Added to Catalog 10/22/2014
PUBMEDID 24528284
First Author Ji Y
Date 08/01/2014
Journal Br J Clin Pharmacol
Link http://www.ncbi.nlm.nih.gov/pubmed/24528284
Study Citalopram and escitalopram plasma drug and metabolite concentrations: genome-wide associations.
Disease/Trait Response to serotonin reuptake inhibitors in major depressive disorder (plasma drug and metabolite levels)
Initial Sample Description 300 European ancestry Escitalpram treated individuals, 130 European ancestry Citalopram treated individuals
Replication Sample Description NA
Region 17q25.3
Chr_id 17
Chr_pos 79831041
Reported Gene(s) CBX4
Mapped_gene CBX8 - CBX4
Upstream_gene_id 57332
Downstream_gene_id 8535
Snp_gene_ids
Upstream_gene_distance 33.93
Downstream_gene_distance 2.12
Strongest SNP-Risk Allele rs9747992-?
SNPs rs9747992
Merged 0
Snp_id_current 9747992
Context Intergenic
Intergenic 1
Risk Allele Frequency 0.086
p-Value 2.00E-07
Pvalue_mlog 6.698970004
p-Value (text) (S-DCT concentration)
OR or beta NR
95% CI (text) NR
Platform [SNPs passing QC] Illumina [7,537,437] (Imputed)
CNV N

 

 

上面这个文件是由tab键分割的，每一列的意义如下！

Note: The SNP data in the catalog has been mapped to dbSNP Build 142 and Genome Assembly,

GRCh38/hg37.p13.

DATE ADDED TO CATALOG: Date added to catalog

PUBMEDID: PubMed identification number

FIRST AUTHOR: Last name of first author

DATE: Publication date (online (epub) date if available)

JOURNAL: Abbreviated journal name

LINK: PubMed URL

STUDY: Title of paper (linked to PubMed abstract)

DISEASE/TRAIT: Disease or trait examined in study

INITIAL SAMPLE SIZE: Sample size for Stage 1 of GWAS

REPLICATION SAMPLE SIZE: Sample size for subsequent replication(s)

REGION: Cytogenetic region associated with rs number (NCBI)

CHR_ID: Chromosome number associated with rs number (NCBI)

CHR_POS: Chromosomal position associated with rs number (dbSNP Build 132,

NCBI)

REPORTED GENE (S): Gene(s) reported by author

MAPPED GENE(S): Gene(s) mapped to the strongest SNP (NCBI). If the SNP is

located within a gene, that gene is listed. If the SNP is intergenic, the upstream and

downstream genes are listed, separated by a hyphen. UPSTREAM_GENE_ID:

Entrez Gene ID for nearest upstream gene to rs number, if not within gene (NCBI)

DOWNSTREAM_GENE_ID: Entrez Gene ID for nearest downstream gene to rs

number, if not within gene (NCBI)

SNP_GENE_IDS: Entrez Gene ID, if rs number within gene; multiple genes

denotes overlapping transcripts (NCBI)

UPSTREAM_GENE_DISTANCE: distance in kb for nearest upstream gene to rs

number, if not within gene (NCBI)

DOWNSTREAM_GENE_DISTANCE: distance in kb for nearest downstream

gene to rs number, if not within gene (NCBI)

STRONGEST SNP-RISK ALLELE: SNP(s) most strongly associated with trait +

risk allele (? for unknown risk allele). May also refer to a haplotype.

SNPS: Strongest SNP; if a haplotype is reported above, may include more than one

rs number (multiple SNPs comprising the haplotype)

MERGED: denotes whether the SNP has been merged into a subsequent rs record

(0 = no; 1 = yes; NCBI)

SNP_ID_CURRENT: current rs number (will differ from strongest SNP when

merged = 1)

CONTEXT: SNP functional class (NCBI)

INTERGENIC: denotes whether SNP is in intergenic region (0 = no; 1 = yes;

NCBI)

RISK ALLELE FREQUENCY: Reported risk allele frequency associated with

strongest SNP

P-VALUE: Reported p-value for strongest SNP risk allele (linked to dbGaP

Association Browser)

PVALUE_MLOG: -log(p-value)

P-VALUE (TEXT): Information describing context of p-value (e.g. females,

smokers).

Note that p-values are rounded to 1 significant digit (for example, a published pvalue of 4.8 x 10-7 is rounded to 5 x 10-7).

OR or BETA: Reported odds ratio or beta-coefficient associated with strongest

SNP risk allele

95% CI (TEXT): Reported 95% confidence interval associated with strongest SNP

risk allele

PLATFORM (SNPS PASSING QC): Genotyping platform manufacturer used in

Stage 1; also includes notation of pooled DNA study design or imputation of

SNPs, where applicable

CNV: Study of copy number variation (yes/no)

Updated: January 13, 2015

 

 

 

 

五 05

Bioconductor的数据包library(biomaRt)简介

Posted on by

 

这是发布在bioconductor平台上面的一个数据库文件，可以通过R里面下载安装并使用，非常方便。其实在ensembl数据库里面也有一个biomart，我之前也讲过这个平台，非常好用，可以把任意的数据库之间的ID号进行转换。

为了更好的理解和掌握biomaRt，我们可以先通过在线资源来了解一下它的原型biomart （http://www.biomart.org）。 biomart是为生物科研提供数据服务的免费软件，它为数据下载提供打包方案。它有许多成功的应用实例，比如欧洲生物信息学中心（The European Bioinformatics Institute ，EBI）维护的Ensembl数据库（http://www.ensembl.org/）就使用biomart提供数据批量下载服务， 还有COSMIC, Uniprot, HGNC, Gramene, Wormbase以及dbSNP等。

这个就是一个R平台的biomart而已，但是非常好用！

> library(biomaRt)

> head(listMarts(), 3)

biomart                           version

1    ensembl      ENSEMBL GENES 79 (SANGER UK)

2        snp  ENSEMBL VARIATION 79 (SANGER UK)

3 regulation ENSEMBL REGULATION 79 (SANGER UK)

这是这个biomart最具有代表性的三个数据库，用listMarts()可以查看得知，它总共有58个数据库。

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这是创建了人的ensembl数据库对象

> head(listFilters(ensembl), 3)

name     description

1 chromosome_name Chromosome name

2           start Gene Start (bp)

3             end   Gene End (bp)

可以看到对人的数据库ensembl来说，有多种字段可以来挑选自己感兴趣的东西，最常用的的当然是染色体号及起始终止坐标啦，用listFilters(ensembl),以查看得知，它总共有284中挑选感兴趣数据的方式。

既然 chromosome_name是其中一个挑选字段，那么我们就可以看看，是如何进行挑选的

用filterOptions(myFilter, ensembl)可以看到它挑选参数非常之多，远不止我们所认为的染色体号码。

染色体号一般是1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y,MT

还有一堆稀奇古怪的标志，LRG_101,LRG_102,LRG_103,LRG_104，因为我们组装好的人的标准基因组还有很多小的片段不被计入染色体中。

然后还可以看到人的ensembl数据库对象，有很多的属性，最常见的当然是基因ID和转录本ID和蛋白的ID号啦！

> head(listAttributes(ensembl), 3)

name           description

1       ensembl_gene_id       Ensembl Gene ID

2 ensembl_transcript_id Ensembl Transcript ID

3    ensembl_peptide_id    Ensembl Protein ID

用listAttributes(ensembl),,以查看得知，它总共有1166个ID号，太恐怖了，我实在是没有想到！

 

那么接下来我简单讲讲这个包的几个应用吧

首先是根据entrez ID号来找

ensembl <-  useMart("ensembl", dataset = "hsapiens_gene_ensembl")

这样就得到了人的信息，然后我探究以下两个基因的其它信息。

entrzID=c("672","1")

getBM(attributes=c("entrezgene","hgnc_symbol","ensembl_gene_id"), filters = "entrezgene", values =entrzID, mart=ensembl )

entrezgene hgnc_symbol ensembl_gene_id

1          1        A1BG ENSG00000121410

2        672       BRCA1 ENSG00000012048

3        672       BRCA1         LRG_292

 

其实这个函数很简单，就是根据自定义的entrzID这个变量来找到一些数据，数据的属性是我自己定义的entrezgene","hgnc_symbol","ensembl_gene_id"，所以它就显示这个信息给我，在我之前弄好的人的数据库里面寻找！listAttributes(ensembl),,以查看得知，它总共有1166个ID号，就是说，你可以挑选你想要的基因的1166种信息，包罗万象！！！

其它功能也是很简单的啦，自己多看帮助文档！

 

从上面的操作来看，使用biomaRt只需要两步，1，指定mart数据库，2，使用getBM获得注释。但是首先，我们如何知道有哪些服务器，以及这些服务器上哪些数据库呢？其次，我们如何获阳getBM中attributes，filters的正确设置呢？

关于第一个问题，我们可以使用biomaRt中的listMarts以及listDatasets两个函数来解决。

> marts <- listMarts(); head(marts) #查看当前可用的数据源 ，总共有58个数据源。

> ensembl <- useMart("ensembl") #使用ensembl数据源

> datasets <- listDatasets(ensembl); datasets[1:10,] #查看ensembl中可用数据库，共有69个物种的数据库！

对于第二个问题，我们使用biomaRt中的listFilters以及listAttributes两个函数来解决。

> mart <- useMart("ensembl", "hsapiens_gene_ensembl")  #首先使用人的数据库

>listAttributes(ensembl) #,以查看得知，它总共有1166个ID号，就是说，人的数据库可供挑选的信息多达1166种。

> filters <- listFilters(mart); filters[grepl("entrez", filters[,1]),] #总共有284中挑选感兴趣数据的方式。

最后的问题是，biomaRt会被如何使用呢？我们做注释的时候，怎么就想到要使用biomaRt呢？因为在注释上，各种ID，symbol, name之间的转换都可以考虑使用biomaRt来做。更重要的是，biomaRt还会有很多SNP, alternative splicing, exon, intron, 5’utr, 3’utr等等信息。当然，只要能做也数据库并使用SQL访问的数据都可以使用biomaRt来获取。所以我们的思路可以更加发散一些。

 

 

 

五 05

Bioconductor的数据包library(org.Hs.eg.db)简介

Posted on by

 

这是发布在bioconductor平台上面的一个数据库文件，可以通过R里面下载安装并使用，非常方便。而且用的是数据库存储方式，所以搜索起来也是非常快速。

这个包里面有28个主流数据资料文件，这样我们可以用select函数根据我们自己的ID在这28个数据库里面随意转换自己想要的信息！！！

当然我本人是比较喜欢直接下载原文件，然后写脚本自己进行各种数据直接的转换。

首先我们加载这个数据包，可以看到这个数据包依赖于很多其它的包，如果是第一次安装。会耗时很长！

Bioconductor的数据包org.Hs.eg.db269

用这个函数，可以看到这个org.Hs.eg.db数据对象里面包含着各大主流数据库的数据，一般人都比较熟悉的entrez ID 和ensembl 数据库的ID。

keytypes(org.Hs.eg.db)

##  [1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"

##  [5] "ACCNUM"       "ALIAS"        "ENZYME"       "MAP"

##  [9] "PATH"         "PMID"         "REFSEQ"       "SYMBOL"

##  [13] "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"

##  [17] "GENENAME"     "UNIPROT"      "GO"           "EVIDENCE"

##  [21] "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"

##  [25] "OMIM"         "UCSCKG"

然后，我们用select函数，就可以把任意公共数据库的数据进行一一对应了。

ensids <- c("ENSG00000130720", "ENSG00000103257", "ENSG00000156414",

"ENSG00000144644", "ENSG00000159307", "ENSG00000144485")

cols <- c("SYMBOL", "GENENAME")

select(org.Hs.eg.db, keys=ensids, columns=cols, keytype="ENSEMBL")

比如说，我们有几个ensembl的基因ID号。然后我们想找它所对应的gene名和缩略词简称，就通过select函数来搞定即可！

Bioconductor的数据包org.Hs.eg.db1158

select(org.Hs.eg.db, keys="BRCA1", columns=c("ENSEMBL","UNIGENE","ENTREZID","CHR","GO","GENENAME"), keytype="SYMBOL")

这样得到了这个BRCA1基因的大部分信息，只是它的GO条目太多了，看得有点乱。

Bioconductor的数据包org.Hs.eg.db1318

 

 

 

五 05

Bioconductor简介

Posted on by

主页：http://www.bioconductor.org/

文字介绍我懒得写了，具体大家参考

http://www.bioconductor.org/about/

http://blog.csdn.net/shmilyringpull/article/details/8542607

这是一个R语言进行生信分析的流程发布平台，每个包都解决生信的一个流程问题。到目前为止2015年5月5日10:57:29已经有了1024个包，所以大家可以看到生信分析是一个海量的任务了。每一个包都有着详尽的说明文档及脚本代码，还附带着数据，非常容易弄懂，接下来我会花一个月的时间好好学习这些包！

这1024个虽然还是R语言的包，但是它的安装方式与常规的R语言包已经有所区别了。

需要用一下代码来安装

source("http://bioconductor.org/biocLite.R")biocLite()

biocLite(c("GenomicFeatures", "AnnotationDbi"))

也是非常easy的。

现在这个平台上面有1024个包，241个实验数据，917个数据库文件！！！

We are pleased to announce Bioconductor 3.1,

consisting of 1024 software packages,

241 experiment data packages,

and 917 up-to-date annotation packages.

在MOOC上面有很多关于这个的公开课

http://bioconductor.org/help/course-materials/

 

这里面有很多生信方向的分析流程，包括了我之前提到了snp-calling，RNA-seq，CHIP-seq等等，当然最主要的还是芯片数据的处理。

Workflows »

Common Bioconductor workflows include:

这些流程基本上涉及到了现在生物信息学的主流方向，所以基本上掌握了这些包，就是一个合格的生物信息学人才啦！

更重要的是它有着917个数据库文件，里面的信息分门别类，几乎可以算作是生物信息学的百科全书啦！

主要的数据库包括以下。

 

Package Description
AnnotationHub Ensembl, Encode, dbSNP, UCSC data objects
biomaRt Ensembl and other annotations
PSICQUIC Protein interactions
uniprot.ws Protein annotations
KEGGREST KEGG pathways
SRAdb Sequencing experiments.
rtracklayer genome tracks.
GEOquery Array and other data
ArrayExpress Array and other data

 

 

 

 

 

 

 

 

 

 

 

 

 

五 04

R语言里面的一个数据集ALL(Acute Lymphoblastic Leukemia)简介

Posted on by

这个数据内容太多了，我感觉自己也理解的不是很清楚！

非常多的R的bioconductor包都是拿这个数据集来举例子的，所以我简单的介绍一下这个数据集。

这个数据集是对ALL这个病的研究数据，共涉及到了128个ALL病人，其中95个是B细胞的ALL，剩余33个是T细胞的ALL。

是一个芯片数据，同时还包含有其它的病人信息。

大家首先要在R里面安装这个数据集

source("http://bioconductor.org/biocLite.R")

biocLite("ALL")

library(ALL)

data(ALL)

data(geneList)

在R里面输入str(ALL)可以看到这个数据的具体信息，但是非常多!

ALL

ExpressionSet (storageMode: lockedEnvironment)

assayData: 12625 features, 128 samples 

element names: exprs

protocolData: none

phenoData

sampleNames: 01005 01010 ... LAL4 (128 total)

varLabels: cod diagnosis ... date last seen (21 total)

varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

 pubMedIds: 14684422 16243790 

Annotation: hgu95av2

我们首先它的BT变量记录的是什么

R语言里面的一个数据集ALL750

可以看出它记录的是数据病人的分组信息。

bcell = grep("^B", as.character(ALL$BT))通过这句话可以挑选出B细胞病人

然后我们看看它的ALL$mol.biol变量记录是是什么

R语言里面的一个数据集ALL857

可以看到它记录的是这些病人的几种突变情况（molecular biology testing）

types = c("NEG", "BCR/ABL")

moltyp = which(as.character(ALL$mol.biol) %in% types)

用这个命令就能挑选出我们想研究的两组突变的病人。

然后我们还可以把刚才的两个标准用来从ALL数据集里面取想要的子集

ALL_bcrneg = ALL[, intersect(bcell, moltyp)]

 

 

同理我们可以查看这个数据集的非常多的变量信息。

包括sex，age，cod，diagnosis，等等，这个'data.frame':共有128 obs. of  21 variables:

R语言里面的一个数据集ALL1190

 

我们除了可以查看这个ALL数据集自带的变量，还可以通过一些方法来访问它的其它信息。

Exprs这个方法可以查看它的表达数据，可以看到有128个变量，12625行的探针数据。

str(exprs(ALL))

num [1:12625, 1:128] 7.6 5.05 3.9 5.9 5.93 ...

- attr(*, "dimnames")=List of 2

..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...

..$ : chr [1:128] "01005" "01010" "03002" "04006" ...

 

还有很多很多函数都可以操作这个数据集，这样可以得到非常多的信息！我就不一一列举了

对这个数据的一系列操作可以画热图，见下面的教程！！！

http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/heatmap/

 

五 03

生物信息入门视频推荐-新一代测序数据分析

Posted on by
  1. 新一代测序数据分析-在优酷里面可以搜索到，一下是配套视频的讲义及下载地址！
  2. Lecture Notes
  3. Lectures will appear below as they are presented. Homeworks are specified in each handout.
  4. Lecture 1 - slides, handouts. course information, homework and project information, introduction to computing, setting up you computer, basic unix command line usage, organizing your projects, homework 1.
  5. Lecture 2 - slides, handouts, The GFF format, sequence ontologies, basic Unix commands: wc, grep, cut, sort, redirecting input and output streams, piping commands, processing a tabular file with UNIX tools, homework 2
  6. Lecture 3 - slides, handouts. programming languages, download and install an proper editor, introduction to the AWK programming language, tabular file processing, filtering by feature types, Awk onliners explained, another collections of AWK oneliners, homework 3.
  7. Lecture 4 - slides, handouts, sequencing technologies, sequence representations, the FASTA format, processing FASTA files at the command line, homework 4.
  8. Lecture 5 - slides, handouts, string matching, edit distances, regular expressions, local and global alignments, homework 5.
  9. Lecture 6 - slides, handouts, introduction to using blast, legacy blast and blast+, preparing blast databases, performing a blastn query, formatting blast output, homework 6.
  10. Lecture 7 - slides, handouts, using blast, formatting databases, using the blastdbcmd, extract sequences, batch operations, formatting blast queries, homework 7.
  11. Lecture 8 - slides, handouts, blast score and E-values, search strategies, usage examples for blastn, blastp, blastx, tblastn, and tblastx, homework 8.
  12. Lecture 9 - slides, handouts, quality encodings, phred scales, the FASTQ format, homework 9.
  13. Lecture 10 - slides, handouts, file compression, gzip, zip, bz2, file archives, tarbombs, plotting fastq qualities homework 10.
  14. Lecture 11 - slides, handouts installing tools, quality control, adapter trimming, error corrections
  15. Lecture 12 - slides, handouts paired end sequencing, quality control for paired end sequencing, the bioawk language
  16. Lecture 13 - slides, handouts paired end sequencing, read stiching, automating tasks with shell scripts
  17. Lecture 14 - slides, handouts short read alignments, bwa, bowtie and other tools.
  18. Lecture 15 - slides, handouts the sequence alignment map SAM format
  19. Lecture 16 - slides, handouts the SAM/BAM format, sorting and indexing BAM files, using the samtools program
  20. Lecture 17 - slides, handouts aligning paired end reads, comparing and evaluating aligners, simulating sequencing reads with the wgsim tool
  21. Lecture 18 - slides, handouts read duplication, visualizing alignments with IGV and IGB
  22. Lecture 19, guest lecture by Nicholas Stoler - slides, the variant call format (VCF), calling variants with samtools mpileup
  23. Lecture 20,- slides, handouts origins of genome variations, more on SNP calling, successes and failures
  24. Lecture 21,- slides, handouts interval representation, BED and GFF formats, representing data
  25. Lecture 22,- slides, handouts interval operations: complement, extension, flanking, Using the BedTools package
  26. Lecture 23,- slides, handouts interval operations: intersect, window, selecting closest features
  27. Lecture 24,- slides, handouts an introduction to genome assembly, using the velvet assembler, evaluating genome assemblies with QUAST
  28. Lecture 25,- slides, handouts, meta.tar.gz (25MB) an introduction to metagenomics, software packages mothur, QIIME and MetaSim, online tools RDP, MG-RAST
  29. Lecture 26,- slides, handouts, lec26.tar.gz (25MB) an introduction to Chip-Seq technology, peak calling concepts, preprocessing and peak calling methods (part 1)
  30. Lecture 27,- slides, handouts, Chip-Seq peak calling sofware, preprocessing and peak calling methods (part 2)
  31. Lecture 28,- slides, handouts, lec28.tar.gz basic RNA-Seq data analysis concepts, split read mapping
  32. Lecture 29, slides, handouts, lec29.tar.gz RNA-Seq (part 2)
  33. Lecture 30, slides, handouts, bioinformatics beyond the command line: using R for data analysis
  34. Final Project 30, final-project, data for final project pony.tar.gz (17Mb) BMB 597D: Final project, 50% of the final grade, due 5pm Saturday Dec 14th, 2013
五 01

脚本作业-解读NCBI的ftp里面关于人的一些基因信息

Posted on by

为了感谢大家对我博客的关注，我在这里发布一个作业，适合菜鸟做的。里面有十几个类似的问题，大家可以下载数据自行处理，如果是问这些问题，我优先回答！

NCBI的ftp里面关于人的一些基因信息

我在NCBI的ftp服务器里面下载了这些数据，时间是2015年，大多是hg19系列的，文件名如下：

CDS.fa 这个是ensembl中人的CDS碱基序列文件，hg38

entrez2go.gene 这个是有go注释的基因情况，有一万八的基因都有go注释

entrez2name.gene 这个是NCBI的entrez ID号对应着基因名的文件

entrez2pubmed.gene 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号

entrez2refseq2ensembl.gene 这个是NCBI的entrez ID号对应着基因名的refseq的ID号和ensembl数据库的ID号

human_gene_info这个是基因的详细信息，包括基因的起始终止点坐标等等

Protein.fa 这个是ensembl中人的蛋白的氨基酸序列文件，有十万多个蛋白hg38

ref2ensembl.txt  这个是基因名的refseq的ID号和ensembl数据库的ID号

自行去NCBI的ftp服务器里面下载这些数据。

然后好好熟悉这些数据信息，回答一下几个问题：

人总的基因有多少个，它们分别分布在哪些染色体上面，基因的转录本分布情况如何，基因的长度分布如何，基因的外显子个数如何。

CD分子的基因有多少个，它们分别分布在哪些染色体上面，基因的转录本分布情况如何，基因的长度分布如何，基因的外显子个数如何。它们有没有氨基酸偏好性？？

MHC系列基因信息？CCL系列基因信息如何？CXCL系列信息如何？或者你感兴趣的基因家族信息？

现在研究最热门的基因是什么？发表文章最多的前十个基因是什么？

基因长度情况如何？最长的基因多长？最短的基因多少bp，可靠吗？

蛋白质长度情况如何？

每条染色体的基因分别情况？基因在染色体那个地方分别最多？

请用图形展示你的结论！！！

 

如果你能回答以上问题，证明你的脚本水平不错了。

如果找不到我，看旁边的公告，加入生信菜鸟群，我就在里面！！！

四 30

人的CD分子基因信息简介！

Posted on by

CD分子吧,它是Clusters of Differentiation的简写,是指一组分化抗原的家族,目前该家族已经有CD1——CD350甚至更多的成员.他们分布于T细胞等免疫细胞表面,参与免疫细胞各种表达,其中有整合素、受体、配体等蛋白分子,在免疫应答反应中参与识别、粘附和信号转导等功能.

我这里简单讲讲如何整理它们的基因信息，首先从NCBI里面下载的人的gene_info文件，然后通过脚本来查找CD分子信息。

perl -alne '{if (/\tCD\d+/ or /CD\d+\|/ ) {print}}' human_gene_info >CD.info

cut -f 2-5 CD.info >CD.table

再根据CD分子的排序把我们的信息重新排序

perl -alne '{/CD(\d+\w)/;$hash{$1}=$_}END{print $hash{$_} foreach sort {$a <=> $b}keys %hash}' CD.table >CD.table.sort

然后我发现了一个很有趣的问题，它们都是负义链上面的基因！

 

entrez ID gene symbol 正负链
911 CD1C - BDCA1|CD1|CD1A|R7
913 CD1E - CD1A|R2
909 CD1A - CD1|FCB6|HTA1|R4|T6
912 CD1D - CD1A|R3
910 CD1B - CD1|CD1A|R1
9266 CYTH2 - ARNO|CTS18|CTS18.1|PSCD2|PSCD2L|SEC7L|Sec7p-L|Sec7p-like
30011 SH3KBP1 - CD2BP3|CIN85|GIG10|HSB-1|HSB1|MIG18
23607 CD2AP - CMS
89886 SLAMF9 - CD2F-10|CD2F10|CD84-H1|CD84H1|SF2001
10849 CD3EAP - ASE-1|ASE1|CAST|PAF49
445347 TARP - CD3G|TCRG|TCRGC1|TCRGC2
915 CD3D - CD3-DELTA|IMD19|T3D
920 CD4 - CD4mut
922 CD5L - AIM|API6|PRO229|SP-ALPHA|Spalpha
925 CD8A - CD8|Leu2|MAL|p32
927 CD8BP - CD8B2
54675 CRLS1 - C20orf155|CLS|CLS1|GCD10|dJ967N21.6
3681 ITGAD - ADB2|CD11D
3683 ITGAL - CD11A|LFA-1|LFA1A
3684 ITGAM - CD11B|CR3A|MAC-1|MAC1A|MO1A|SLEB6
3687 ITGAX - CD11C|SLEB6
290 ANPEP - APN|CD13|GP150|LAP1|P150|PEPN
115708 TRMT61A - C14orf172|GCD14|Gcd14p|TRM61|hTRM61
2526 FUT4 - CD15|ELFT|FCT3A|FUC-TIV|FUTIV|LeX|SSEA-1
2215 FCGR3B - CD16|CD16b|FCG3|FCGR3|FCR-10|FCRIII|FCRIIIb
4055 LTBR - CD18|D12S370|LT-BETA-R|TNF-R-III|TNFCR|TNFR-RP|TNFR2-RP|TNFR3|TNFRSF3
930 CD19 - B4|CVID3

 

 

 

 

 

 

 

四 21

HGNC数据库简介

Posted on by

人类基因命名委员会(HUGO Gene Nomenclature Committee);人类基因组命名委员会!

其实有了NCBI的entrez ID，然后还有refseq里面的ID，还有ensembl的ID，还有基因本身的功能英文缩略简称，已经很麻烦了，又来了一个HGNC，唉，头疼！

The HGNC approves both a short-form abbreviation known as a gene symbol, and also a longer and more descriptive name.

可以下载整个数据，用脚本慢慢研究研究

wget ftp://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt

 

还是看看BRCA1这个基因，里面的信息挺多的，主要看HGNC:1100，就是这个数据库对它这个基因的编号

HGNC:1100

BRCA1  这个是基因名，需要得到该组织的认可！！！！

breast cancer 1, early onset protein-coding gene gene with protein product Approved

17q21.31 17q21.31

"RNF53|BRCC1|PPP1R53|FANCS" "BRCA1/BRCA2-containing complex, subunit 1|protein phosphatase 1, regulatory subunit 53|Fanconi anemia, complementation group S" "Ring finger proteins|Protein phosphatase 1 regulatory subunits" "58|694"

1991-02-20T00:00:00Z

2015-04-18T00:00:00Z

672     这里是entrez ID

ENSG00000012048  这里是ensembl的ID，

OTTHUMG00000157426 uc002ict.3 U14680 NM_007294

"CCDS11453|CCDS11454|CCDS11455|CCDS11456|CCDS11459|CCDS11455|CCDS11456|CCDS11459|CCDS11454" P38398 1676470 MGI:104537 RGD:2218

"Breast Cancer|http://research.nhgri.nih.gov/bic/|BRCA1 database at LOVD-China|http://genomed.org/LOVD/BC/home.php?select_db=BRCA1|LOVD - Leiden Open Variation Database|http://chromium.liacs.nl/LOVD2/cancer/home.php?select_db=BRCA1|LOVD - Leiden Open Variation Database|http://proteomics.bio21.unimelb.edu.au/lovd/genes/BRCA1|LRG_292|http://www.lrg-sequence.org/LRG/LRG_292"

BRCA1 113705 119068

 

数据结构大概就是这个样子的了！

这几个数据库的内容都是互相链接的！

 

然后我们看看HGNC数据库的一些统计信息

http://www.genenames.org/cgi-bin/statistics

总共有40392个基因信息

其中18990个是能编码蛋白产物的基因，它们大多有GO注释

其中5927个是non-coding RNA，是现在的研究热门。

还有12546个是假基因，挺复杂的

最后还有1188个免疫相关基因，位置基因，病毒基因等等

 

最后，送给大家一个彩蛋！还有十一个物种也是有一个命名委员会的！

类似于 Mouse Gene Nomenclature Committee (MGNC).  Please see the following links:

 

参考文献;

Gray KA, Yates B, Seal RL, Wright MW, Bruford EA. genenames.org: the HGNC resources in 2015. Nucleic Acids Res. 2015 Jan;43(Database issue):D1079-85. doi: 10.1093/nar/gku1071. PMID:25361968

四 18

Bowtie算法第五讲-index2tally

Posted on by

前面讲到了如何用笨方法进行字符串搜索，也讲了如何构建bwt索引，和把bwt索引还原成字符串！

原始字符串是ATGCGTANNGTC

排序过程是下面的

$ATGCGTANNGTC 12

ANNGTC$ATGCGT 6

ATGCGTANNGTC$ 0

C$ATGCGTANNGT 11

CGTANNGTC$ATG 3

GCGTANNGTC$AT 2

GTANNGTC$ATGC 4

GTC$ATGCGTANN 9

NGTC$ATGCGTAN 8

NNGTC$ATGCGTA 7

TANNGTC$ATGCG 5

TC$ATGCGTANNG 10

TGCGTANNGTC$A 1

现在讲讲如何根据bwt索引构建tally，并且用tally搜索方法来搜索字符串！

首先是bwt索引转换为tally

C 12

T 6

$ 0

T 11

G 3

T 2

C 4

N 9

N 8

A 7

G 5

G 10

A 1

这个其实非常简单的，tally就是增加四列计数的列即可

[perl]

$hash_count{'A'}=0;

$hash_count{'C'}=0;

$hash_count{'G'}=0;

$hash_count{'T'}=0;

open FH ,"<$ARGV[0]";

while(<FH>){

        chomp;

@F=split;

$last=$F[0]; # 读取上面的tally文件，分列，判断第一列，并计数

        $hash_count{$last}++;

   print  "$_\t$hash_count{'A'}\t$hash_count{'C'}\t$hash_count{'G'}\t$hash_count{'T'}\n";

}

[/perl]

输出的tally如下

C 12 0 1 0 0

T 6 0 1 0 1

$ 0 0 1 0 1

T 11 0 1 0 2

G 3 0 1 1 2

T 2 0 1 1 3

C 4 0 2 1 3

N 9 0 2 1 3

N 8 0 2 1 3

A 7 1 2 1 3

G 5 1 2 2 3

G 10 1 2 3 3

A 1 2 2 3 3

接下来就是针对这个tally的查询函数了