生物信息学大实验实验指导.doc

生物信息学大实验生物信息学大实验实验指导适用专业：生物信息学生物与制药大类编写：解增言生物信息学院2014年4月目录实验1 基因组序列组装（软件CAP3的使用）3实验2 基因组序列组装（软件velvet的使用）7实验3 原核生物基因识别（软件Glimmer的使用）9实验4 真核生物基因识别（软件GlimmerM的使用）17实验5 HAP3基因家族的多序列比对分析（软件ClustalX与ClustalW的使用）25实验6 HAP3基因家族的分子进化分析（软件PHYLIP和MEGA的使用）28实验7 HAP3基因家族的分子进化分析（软件MrBayes的使用）33实验8 HAP3蛋白的结构分析（软件RasMol的使用）42实验1 基因组序列组装（软件CAP3的使用）一、实验目的1 了解基因组测序原理和主要策略；2 掌握CAP3序列组装软件的使用方法。二、实验原理基因组测序常用的两种策略是克隆法（clone-based strategy）和全基因组鸟枪法（whole genome shotgun method）。克隆法先将基因组DNA打成大的片段，连到载体上，构建DNA文库；再对每一个大片段（克隆）打碎测序。序列组装时先组装成克隆，再组装成染色体。克隆测序法的好处在于序列组装时可以利用已经定位的大片段克隆, 所以序列组装起来较容易, 但是需要前期建立基因组物理图谱, 耗资大, 测序周期长。全基因组鸟枪法测序无需构建各类复杂的物理图谱和遗传图谱，采用最经济有效的实验设计方案，直接将整个基因组打成不同大小的DNA片段构建Shotgun文库，再用传统Sanger测序法或Solexa等新一代测序技术对文库进行随机测序。最后运用生物信息学方法将测序片段拼接成全基因组序列。该方法具有高通量、低成本优势。序列组装时，先把把单条序列（read）组装成叠连群（contig）、再把叠连群组装成“支架”（scaffold），最后组装成染色体。本实验将练习在Linux环境下用CAP3软件组装流感病毒基因组。1CAP3序列组装程序简介Huang Xiaoqiu. 和 Madan，A. 开发的一套用于序列拼接的软件，此软件适用于小的数据集或 EST 拼接，它有如下特征：1. 应用正反向信息更正拼接错误、连接contigs。2. 在序列拼接中应用 reads 的质量信息。3. 自动截去 reads5端、3端的低质量区。4. 产生 Consed 程序可读的ace 格式拼接结果文件。5. CAP3 能用于Staden软件包的中的GAP4 软件。2下载此软件可以免费下载，下载地址：http：//download.html。填写基本信息表格，即可下载。CAP3 详细参考文档可见：http：//sas.html。3安装（1）上传cap3 的压缩包到本地linux/unix 运算服务器；（2）解压缩： bash-2.05b$ tar xvf cap3.tar CAP3/ CAP3/README CAP3/cap3 CAP3/doc CAP3/aceform CAP3/formcon（3）查看解压缩后的文件： bash-2.05b$ ls l total 240 -rwxr-xr-x 1 soft bgi 25844 Sep 2 2002 formcon* -rwxr-xr-x 1 soft bgi 169836 Sep 2 2002 cap3* -rw-r- 1 soft bgi 513 Aug 22 2002 README -rw- 1 soft bgi 18448 Aug 22 2002 aceform -rw-r- 1 soft bgi 18922 Jun 21 2002 doc4. 使用程序运行命令行：cap3 options cap3.out5输入：输入序列是普通的FASTA格式，如果序列文件名为“xyz”，则质量文件应命名为“xyz.qual”，约束文件应命名为“xyz.con”。在命令行中只需输入序列文件，程序会自动在相应的目录中寻找相应的质量文件和约束文件。 “xyz”格式如下：Sequence1ACGTGCGCGATCGCCTGCTAGGCGTACGTCGCAGGCGATCGATGTGCTAGATCAGATGACASequence2GGGCTAGATTAGCACCACATACATCGCTCA“xyz.qual”格式如下：R16 8 8 8 15 17 17 17 12 12 20 20 29 31 34 34 38 38 40 40 49 49 37 33 3333 33 30 31 24 24 34 45 45 45 45 38 38 38 45 40 40 40 40 40 40 40 40 40 4033 33 33 33 33 33 40 37 40 40 45 45 45 40 40 40 45 45 45 45 49 49 49 49 4540 43 43 43 40 40 40 37 40 49 49 40 40 37 37 37 42 45 4049 45 45 45 45 4036 36 36 36 33 33 27 27 21 19 19 27 33 33 34 36 36 36 36 38 36 36 40 33 35R298 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 98 9837 37 37 37 37 37 37 37 37 37 37 37 37 34 34 34 34 37 37 37 37 34 34 37 3834 37 34 37 37 37 37 37 45 37 37 37 37 37 37 37 40 37 37 32 45 41 45 45 41约束文件“xyz.con”中每一行都以如下格式指定了正反向的约束： ReadA ReadB MinDistance MaxDistance其中“ReadA”和“ReadB”是两个 reads 的名称；“MinDistance”、“MaxDistance”是最小、最大距离(bp)。约束文件*.con 可由此软件包中的 formcon 程序生成，用法： formcon 序列文件 最小长度 最大长度 此处最小、最大长度指克隆的长度限制，单位为 bp。克隆长度限制要与插入片段长度相差1000bp 到 1500bp 左右，如：插入片段为 2kb 到 3kb，建议 500 为克隆最小长度，4000 为克隆最大长度。 输入的序列文件中一对正反向的reads 名称在第一个句点前要保持相同。6输出输出文件格式：1. xyz.cap.ace：ace格式文件。注意：reads 的 5、3的低质量区没有被显示在 ace 格式中。2. xyz.cap.contigs：生成的contigs 序列文件。3. xyz.cap.contigs.qual：生成的contigs 质量文件。4. xyz.cap.singlets：没有用于拼接的reads 文件。5. ：关于拼接的额外信息文件。6. cap3.out：拼接的结果文件。7参数参数选项(默认值)： -a N specify band expansion size N 10 (20) -b N specify base quality cutoff for differences N 15 (20) 去除低质量时的质量值 N 5 (12) -c N -d N specify max qscore sum at differences N 20 (200) -e N specify clearance between no. of diff N 10 (30) 重叠部分最大 gap 长度 N 1 (20) -f N gap 罚分 N 0 (6) -g N -h N specify max overhang percent length N 2 (20) 比对分值 N 0 (2) -m N 不匹配的分值 N 20 (40) -p N specify overlap percent identity cutoff N 65 (80) -r N specify reverse orientation value N = 0 (1) -s N specify overlap similarity score cutoff N 400 (900) 匹配得最大长度 N 30 (300) -t N 用于修正得最小约束数目 N 0 (3) -u N 用于连接得最小约束数目 N 0 (2) -v N 序列去除信息的文件名 (none) -w N 输出文件名称的前缀 (cap) -x N 去除碱基范围 N 5 (250) -y N -z N min no. of good reads at clip pos N 0 (3)三、实验内容（步骤）1登陆，在个人目录下新建文件夹：seq_assembl_cap3：$ cd$ mkdir seq_assembl_cap32将服务器目录/home/pub/genome/virus/influenza_a_virus/raw_data下的所有文件复制到seq_assembl_cap3目录下：$ cd seq_assembl_cap3$ cp /home/pub/genome/virus/influenza_a_virus/raw_data/* . #注意最后的点“.”及其前面的空格3. 将序列文件解压：$zcat fasta.influenza_a_virus_a_new_york_ur07_0093_2008_h3n2_.001.gz fasta4. 运行CAP3程序：$ cap3 fasta5. 检查结果四、实验报告1运算环境（包括操作系统和软件），实验步骤，结果文件记录；2实验中遇到的问题，如何解决的。五、参考文献Huang, X. and Madan, A. (1999) CAP3: A DNA Sequence Assembly Program. Genome Research, 9: 868-877.实验2 基因组序列组装（软件velvet的使用）一、实验目的1 了解新一代测序技术所测序列的特点2 掌握软件velvet的使用方法二、实验原理新一代测序技术又称第二代测序技术，主要包括SOLiD，Solexa和454测序技术，其特点是通量高、测序成本低、测序速度快、周期短，缺点是所测序列较短，给后续的序列组装带来了许多问题，适用于第一代测序技术的组装软件如CAP3，phrap等不能组装新一代测序技术尤其是SOLiD和Solexa所测的序列。因此，人们又开发了适用于新一代测序技术的软件，如NextGENe（ABI开发，商业软件），ELAND（Solexa），EDENA，SSAKE及velvet等。本实验将学习使用velvet软件组装Solexa基因组测序结果。Velvet软件主要有两个程序组成：velveth和velvetg。(1)velveth的输入默认是fasta格式的序列文件，也能识别fastq、fasta.gz、fastq.gz、sam、bam、eland和gerald文件。序列类型默认是short，也可以是shortPaired、short2、shortPaired2、long或longPaired。命令格式为：$ ./velveth output_directory hash_length -file_format-read_type filename例如：$ ./velveth output_directory/ 21 -fasta -short solexa1.fa solexa2.fa solexa3.fa long capillary.fa 也可以写成：$ ./velveth output_directory/ 21 solexa*.fa -long capillary.fa（21是hash表的字长度）velveth运行的结果生成一个hash表，并输出3个文件，其中Roadmaps和Sequences文件是下一步velvetg程序运行必须的：Log：日志文件Roadmaps：路线图文件Sequences：序列文件，包含所有输入的序列(2)velvetg是velvet的核心程序，其命令格式为：$ ./velvetg output_directory/ -cov_cutoff-max_coverage运行的结果输出以下文件：contigs.fa：fasta格式的组装好的片段，长度大于2k（k为velveth运行时用的字长）PreGraph：中间组装图LastGraph：最后组装图Graph：最后组装图stats.txt：统计信息三、实验内容（步骤）1 Velvet软件编译：$ cd$ mkdir seq_assemble_velvet$ cp /home/pub/software/velvet_0.7.62.tar.gz seq_assemble_velvet/$ cd seq_assemble_velvet$ tar xzf velvet_0.7.62.tar.gz #解压$ cd velvet_0.7.62$ make#编译2. 运行velveth：$ cd data$ mkdir output$ ./velveth output 21 test_reads.fa3. 运行velvetg：$ ./velvetg output4. 检查结果。5. 改变选项和参数，看结果有无不同。四、实验报告1运行环境（包括操作系统和软件），实验步骤，结果文件记录；2实验中遇到的问题，如何解决的。五、参考文献Daniel Zerbino. （2008）“Velvet Manual - version 0.7”实验3 原核生物基因识别（软件Glimmer的使用）一、实验目的1 了解原核生物基因识别的特点2 掌握软件Glimmer的使用方法二、实验原理 Glimmer 预测系统先用 build-icm 程序对该物种已知的基因序列生成一个马尔可夫模型参数集合，glimmer再应用这个参数集对DNA 序列进行基因预测。此软件适合对原核生物进行预测。 本软件使用分两步： 第一步为编码序列的可能模型的建立,通过 build-icm 来实现,此训练数据集来源有以下三种:1. 某基因组的已知基因2. 通过 long-orfs 产生的长的无重叠的orfs3. 高度相似的物种的基因命令行：build-icm options train.model input_filebuild-icm 程序参数 -d Set depth of model to ,默认为 7 -F Ignore input strings with in-frame stop codons -h Print this message -p Set period of model to ,默认为 3 -r Use the reverse of input strings to build the model -t Output model as text (for debugging only) -v Set verbose level; higher is more diagnostic printouts -w Set length of model window to ,默认为 12 输入文件input_file为 fasta 格式的序列文件, 输出文件 train.model 做为下一步运行 glimmer 的输入,为马尔可夫模型参数集合文件(train.model缺省为二进制格式文件)。 第二步：运行 glimmer3命令行:glimmer3 options my.seq train.model resultglimmer3 参数：-A -start_codons Use comma-separated list of codons as start codons Sample format: -A atg, gtg Use -P option to specify relative proportions of use. If -P not used, then proportions will be equal -b -rbs_pwm Read a position weight matrix (PWM) from to identify the ribosome binding site to help choose start sites -C -gc_percent Use as GC percentage of independent model Note: should be a percentage, e.g., -C 45.2 -E -entropy Read entropy profiles from . Format is one header line, then 20 lines of 3 columns each. Columns are amino acid, positive entropy, negative entropy. Rows must be in order by amino acid code letter -f -first_codon Use first codon in orf as start codon -g -gene_len Set minimum gene length to -h -help Print this message -i -ignore specifies regions of bases that are off limits, so that no bases within that area will be examined -l -linear Assume linear rather than circular genome, i.e., no wraparound -L -orf_coords Use to specify a list of orfs that should be scored separately, with no overlap rules -M -separate_genes is a multifasta file of separate genes to be scored separately, with no overlap rules -o -max_olap Set maximum overlap length to . Overlaps this short or shorter are ignored. -P -start_probs Specify probability of different start codons (same number & order as in -A option). If no -A option, then 3 values for atg, gtg and ttg in that order. Sample format: -P 0.6,0.35,0.05 If -A is specified without -P, then starts are equally likely. -q -ignore_score_len Do not use the initial score filter on any gene or more base long -r -no_indep Dont use independent probability score column -t -threshold Set threshold score for calling as gene to n. If the in-frame score = , then the region is given a number and considered a potential gene. -X -extend Allow orfs extending off ends of sequence to be scored -z -trans_table Use Genbank translation table number for stop codons -Z -stop_codons Use comma-separated list of codons as stop codons Sample format: -Z tag,tga,taa输入输入文件 my.seq 就是普通的 FASTA 格式文件,如下:Sequence1ACGTGCGCGATCGCCTGCTAGGCGTACGTCGCAGGCGATCGATGTGCTAGATCAGATGACASequence2GGGCTAGATTAGCACCACATACATCGCTCA输出glimmer3 输出两个文件,一个是包含运行命令、参数及预测信息等的*.detail 文件,另一个为预测结果文件*.predict:1. *.detail 文件,此文件内容包含三部分,第一行为程序命令行,第二部分为详细参数信息,第三部分为每个输入序列的预测信息。% more out.detailCommand: ././bin/glimmer3 tpall.fna train.model outSequence file = tpall.fnaNumber of sequences = 1ICM model file = train.modelExcluded regions file = noneList of orfs file = noneInput is NOT separate orfsIndependent (non-coding) scores are usedCircular genome = trueTruncated orfs = falseMinimum gene length = 100 bpMaximum overlap bases = 30Threshold score = 30Use first start codon = falseStart codons = atg,gtg,ttgStart probs = 0.600,0.300,0.100Stop codons = taa,tag,tgaGC percentage = 52.8%Ignore score on orfs longer than 799gi|15638995|ref|NC_000919.1| Treponema pallidum subsp. pallidum str. Nichols, complete genome Sequence length = 1138011. .对于第三部分预测信息各列解释如下:ID:预测基因的编号。Frame:Orf的阅读框,“+”代表正链,“-”代表反链,由终止密码子最左侧碱基所在位置决定。Start of Orf:Orf的第一个碱基位置。Start of Gene:Gene起始密码子的第一个碱基位置。Stop:终止密码子最后一个碱基位置。Length of Orf:Orf的碱基数。Length of Gene:基因的碱基数(不包含终止密码子)。Scores:各个预测分析所得分值,详细如下:Raw Score:This is 100 times the per-base log-odds ratio of the in-frame coding ICM score to the independent (i.e., non-coding) model score. It gives a rough quantification to how well an orf scores that can be compared between any two orfs.InFrm Score:The normalized (to the range 0 . . . 99) score of the gene in its reading frame. This is just the appropriate-frame value among the next six scores.Frame Scores:The normalized (to the range 0 . . . 99) score of the gene in each reading frame. A “-” indicates the presence of a stop codon in that reading frame. The normalization compares only scores without stop codons and the independent (non-coding) NC score. If the orf is sufficiently long, i.e., longer than the value stated in “Ignore score on orfs longer than. . . ”, the score is not used.NC Score:The normalized independent (i.e., non-coding or intergenic) model score. This model is adjusted for the fact that the orf, by definition, has no inframe stop codons.EDR Score:An additional column of scores is produced if the -E option is specified. This is the entropy-distance ratio, i.e., the ratio of the distance of the amino-acid distribution from a positive model to the distance from a negative model. Scores below 1.0 are more likely to be genes; scores above 1.0 less likely to be genes. It is not currently used in the scoring process.2. *.predict 文件,此文件内容为最后的基因预测结果,每条序列的结果以 fasta 头开始,然后接着每行代表一个基因。% more out.predictgi|15638995|ref|NC_000919.1| Treponema pallidum subsp. pallidum str. Nichols,complete genomeorf00001 4 1398 +1 2.96orf00003 1641 2756 +3 2.89orf00004 2776 3834 +1 2.89orf00009 4391 6832 +2 2.98orf00010 9446 8340 -3 2.95orf00011 9449 9553 +2 2.04orf00012 10396 12378 +1 2.97orf00013 13029 13133 +3 0.97orf00014 13966 15780 +1 2.97orf00016 15758 18316 +2 2.98基因的各列信息分别为:Column 1 预测基因编号,此编号和*.detail文件里编号一致。Column 2 基因的开始位置。Column 3 基因的结束位置。为终止密码子的最后一个碱基位置,也就是说包含终止密码子。Column 4 阅读框。Column 5 基因的“raw”分值。三、实验内容预测 E. coli 的基因,与已经注释好的做比对。1. 训练及测试数据：/home/pub/genome/bacteria/Escherichia_coli/data_for_glimmer3/。在自己家目录下新建目录，并将该目录中的文件复制到其中。2. Glimmer程序安装在：/usr/local/bin/glimmer下，将该路径添加的环境变量PATH，即可直接运行。源程序及使用手册在/home/pub/software/下。四、实验报告1运行环境（包括操作系统和软件），实验步骤，结果文件记录；2实验中遇到的问题，如何解决的。五、参考文献1. A.L. Delcher, D. Harmon, S. Kasif, O. White, and S.L. Salzberg. Improved microbial gene identification with GLIMMER, Nucleic Acids Research 27:23 (1999), 4636-4641.2. S. Salzberg, A. Delcher, S. Kasif, and O. White. Microbial gene identification using interpolated Markov models, Nucleic Acids Research 26:2 (1998), 544-548.3. A.L. Delcher, K.A. Bratke, E.C. Powers, and S.L. Salzberg. Identifying bacterial genes and endosymbiont DNA with Glimmer. Bioinformatics (Advance online version) (2007).4. 北京华大基因研究中心. 生物信息学培训教材实验4 真核生物基因识别（软件GlimmerM的使用）一、实验目的1 了解真核生物基因识别的特点2 掌握软件GlimmerM的使用方法二、实验原理简介GlimmerM是 TIGR 开发的用于真核生物基因预测的软件。该软件包可以从/pub/software/GlimmerM/GlimmerM-2.5.1.tar.gz免费下载。安装1、上传 GlimmerM 的压缩包到本地 linux/unix 运算服务器;2、解压缩:gzip d GlimmerM-2.5.1.tar.gztar xvf GlimmerM-2.5.1.tar3、查看解压缩后的文件:bash-2.05b$ ls -ltotal 32drwxr-x- 2 soft bgi 4096 Jul 24 2003 sources/drwxr-x- 3 soft bgi 4096 Jun 3 2003 train/drwxr-x- 2 soft bgi 67 May 15 2003 bin/-rw-r- 1 soft bgi 5428 May 15 2003 modifyStartStop.txtdrwxr-x- 4 soft bgi 42 May 15 2003 trained_dir/-rw-r- 1 soft bgi 1311 May 15 2003 README.firstdrwxr-x- 5 soft bgi 72 May 15 2003 Malaria/-rwxr-x- 1 soft bgi 114 May 15 2003 COPYRIGHT*-rwxr-x- 1 soft bgi 5207 May 15 2003 LICENSE*4、编译源程序:(1)进入 sources 目录(cd sources)在命令行键入 make编译完成后,可将执行文件拷贝到 bin 目录下面。(2)进入train目录(cd train) ,键入make完成编译。如果perl不在/usr/local/bin下,需将trainGlimmerM程序首行的perl的位置,改为相应的位置。使用此软件分两步运行:第一步为训练数据集,程序命令行: trainGlimmerM optional_parameters第二步基于训练结果进行基因预测,程序命令行: glimmerm options输入1. 需要预测的基因组序列,为 Fasta 格式的核酸序列:AK058311_Chr03_11499998_11507414AATATCACCAGTGTCTTATAAGCAATTACACCAATTTTATTGGGGGTGTTCCTCAATAGTAACACGATTTTACAATGCCCCTTAACCAATTACATAAAATTTGAATATTCTATAACAAATTTTGCATTTTAGGGTTCCGTAGGTCTTCTAGATACACTAAAGCATCTGCGTAAGATGAGCTCGACGATCATGATGTCATGTCACCGAGGAATGAAAACGGAGCAGATAATTCTCTAATAGTTCAGATACGGATATGATTATGGATATTTGCTCTCGGATACGAATACAGGTATGATGTCATGGTTTCCACTGGATACGGATATCCGATGAGCAGTGCTGTTCGGATATCCGCTGCGAACATGAGACATACTCCAATTCTTTATAACTCTACCACCTTCAATATACCTTTTAGATTTTATAATGAGTAAATATCAGAAAACTGTAACTTTAGTAACAAAACTATCAGTTAGCTGCAACATT.AAGCCCAAGATACTTGAAAAGAACAGAAAAACAGAAG2. 训练数据时使用的基因信息文件,格式为:序列 1 exon1 起始位置 exon1 终止位置序列 1 exon2 起始位置 exon2 终止位置序列 1 exon1 起始位置 exon1 终止位置序列 1 exon2 起始位置 exon2 终止位置序列 2 exon1 起始位置 exon1 终止位置序列 2 exon2 起始位置 exon2 终止位置.一条序列有多个基因,基因间用空行分隔,不完整基因用表示,如: seq1 5 15 seq1 20 34 seq1 50 48 seq1 45 36 seq2 234 seq3 105输出glimmerM 运行结果屏幕标准输出,可重定向到文件,格式如下: GlimmerM (Version 3.0) Sequence name: CNS06C8G Sequence length: 209982 bp Predicted genes/exons 参数详细的参数列表可以通过直接输入命令查看:bash-2.05b$ trainGlimmerM Train GlimmerM module. Usage: trainGlimmerM optional_parameters is a multifasta file containing the sequences for training with the usual format: is a file with the exon coordinates relative to the sequences contained in the ; different genes are separated by a blank line; I am assuming a format like below: seq1 5 15 seq1 20 34 seq1 50 48 seq1 45 36 seq2 17 20 In this example seq1 has two genes: one on the direct strand and another one on the complementary strandThe partial genes can be specified as in the following example: seq2 234 seq3 105 optional_parameters -a filter_value specifies the window length for filtering locally maximal acceptor sites (default=60) -d filter_value specifies the window length for filtering locally maximal acceptor sites (default=60) 运行完 trainGlimmerM, 在运行目录会产生一