1. Office2003使用
1.1 显示大纲工具栏
点击“插入”——“引用”——“索引和目录”，弹出一个框，点击“目录”，点击“显示大纲工具栏”按钮，则会显示出大纲工具栏。
大纲工具栏用于设置大纲，指定大纲的级别，以利于目录的自动更新。
Microsoft Word使用笔记
1
分类目录归档：未分类
Pfam
1. simple introduction
The Pfam database is a large collection of protein families, each represented by multiple sequence alignments and hidden Markov models (HMMs).
There are two components to Pfam: Pfam-A and Pfam-B. Pfam-A entries ara high quality, manually curated families. Pfam-B families are of lower quality, but can be useful for identifying functionally conserved regions when no Pfam-A entries are found.
Pfam also generates higher-level groupings of related families, known as clans. A clan is a collection of Pfam-A entries which are related by similarity of sequence, structure of profile-HMM.
2. citation
The latest version of Pfam is release 27.0, which contains 14,831 manually curated protein families.
3.
454测序与组装
软件安装
$ tar zxf DataAnalysis_2.9_All_20130530_1559.tgz $ cd DataAnalysis_2.9_All/ $ sudo yum install zlib-devel.i686 \ libXi-devel.i686 libXtst-devel.i686 \ libXaw-devel.i686 在CentOS6 x86_64系统中需要安装这些。 $ ./setup.sh 这一步产生图形化界面，安装很简单。
OrthoMCL介绍
1. OrthoMCL的用途
基于序列的相似性，OrthoMCL能将一组proteins（比如全基因组的proteins）归类到ortholog groups、in-paralogs groups和co-orthologs。
2. OrthoMCL-DB
OrthoMCL-DB包含了很多proteins，这些proteins来自一些已经完全测序的真核或原核生物的基因组。OrthoMCL-DB将这些proteins进行了聚类，分成很多的ortholog groups。
2010.5.31，发布了OrthoMCL-DB第4版，包含 116,536个ortholog groups、1,270,853个proteins、88个真核生物基因组、16个古菌基因组、34个细菌基因组。
2011.5.31，发布了OrthoMCL-DB第5版，包含 124,740个ortholog groups、1,398,546 个proteins、150个基因组。
2013年末即将发布OrthoMCL-DB第6版。
3. OrthoMCL的两种使用方法
1. OrthoMCL-DB的官网已经将数据中的proteins进行了ortholog的聚类，其网站提供了一个工具，用于接收上传的基因组proteins，再将这些proteins group到相应的ortholog groups中。官网提供的工具Assign your proteins to OrthoMCL Groups用于进行分析。
2. 如果要对多个基因组的proteomes进行聚类，则可以使用OrthoMCL单机版的软件来进行运算。其用法详见：OrthoMCL的使用。
4. OrthoMCL算法
1. 将多个proteomes转换成orthomcl兼容的FASTA文件。
2. 移除低质量的序列。
3. All-versus-All BLASTP with 1e-5 cutoff。即使用这些proteomes的protein sequences构建blast数据库，再将所有的这些序列和数据库进行BLASTP比对，取evalue小于1e-5的比对结果。
4. Filter by percent match length。计算比对结果的percent match length ( 所有hsp中比对上序列的长度之和 / 两条序列中短的那条序列的长度 )。取50%的cutoff值。
5. 寻找不同物种间potential ortholog pairs(两两物种的protein序列相互是best hits)；寻找同一物种内in-paralog pairs（相互之间是better hits，即对于2个序列之中的任意一条序列，和其in-paralog序列之间的evalue值 <= 这条序列和其它物种比对的evalue值）.
6. 根据上一步结果寻找co-ortholog pairs(pairs connected by orhthology and in-paralog,并且pairs之间的evalue值低于1e-5).
7. 对所有的pairs进行E-values的Normalization，以利于下一步MCL的计算。见下一部分内容，或参考OrthoMCL Algorithm Document。
8. 将所有的ortholog,in-paralog和co-ortholog pairs，以及它们的标准化后的weight值输入到MCL程序中，来进行聚类分群。MCL documentation
5. pairs的evalue计算和标准化
pairs的evalue计算：pairs的两条序列相互blast后有两个evalue值，这两个值常常不相等。但是为了计算需要，于是pairs的之间的两个evalue值要进行一个计算，得到pairs weight，weight= ( -log10(evalue1) + -log10(evalue2) ) / 2 。
pairs的evalue的标准化：1. 对于in-paralog pairs，在某一个基因组中，取两条序列中任意一条序列有ortholog的in-paralog pairs为有效in-paralog pairs。若在这个基因组没有这样的pairs，则该基因组所有的in-paralog pairs都为有效的in-paalog pairs。最后得到所有基因组所有有效的in-paralog pairs。然后取这些有效in-paralog pairs的weight的平均值。最后，每个in-paralog pair的evalue标准化后的值为其weight除以average weight。 2. 对于ortholog或co-ortholog pairs则简单很多，求所有weight的平均值，然后使用各个pair的weight除以average weight，则将其标准化了。
6. 网络版的OrthoMCL的使用
OrthoMCL-DB已经对150个proteomes进行了OrthoMCL的分析，对orthologs进行了聚类。这个过程由于数据量大，因此，在好几百的CPU资源下也需要好几个星期才能做完。
在OrthoMCL-DB上上载 a set of proteins,服务器则会将所有上传的proteins比对到OrthoMCL-DB中所有的proteins上；选取evalue 1e-5和50% match的cutoff；然后将protein归类到其top hit所对应的protein的类上；如果top hit所对应的proteins没有group，则该protein归类到NO_GROUP。
然后，再对上一步cutoff掉的proteins来使用OrthoMCL-DB的in-paralog算法来创建in-paralogs pairs,然后再进行MCL的聚类。
使用该方法，最后将a set of proteins进行了同源基因的聚类，但是缺点如下:
1. 这种方法是单向最佳，根据protein比对的最佳结果去归类到已有的group中去。但是反过来，最佳比对结果对应的protein不一定和query protein是最佳的。这和OrthoMCL的算法是有出入的，所以该方法省了时间，但是结果和真正的结果是有一定差别的。
2. 只使用cutoff后剩下的proteins进行in-paralog分析，而没有进行所有query proteins之间的in-paralog分析。
3. 没有ortholog pairs和co-ortholog pairs的信息，没法进行单拷贝同源基因的提取与分析。
7. 本地OrthoMCL的使用
对指定的a set of proteomes进行同源基因分析，则使用本地的OrthoMCL进行分析。而官网不提供这种服务，因为消耗的计算机资源过大。
8. 注意事项
1. 序列都是使用protein序列，而不是nucleotide序列，是因为protein序列更精确。
2. proteomes中的序列要去除可变剪切，只留取alternative proteins中长度最长的。否则在有alternative proteins存在的情况下，则会造成pseudo-in-paralogs（即alternative proteins称为in-paralogs），给后续的分析造成麻烦。
3. paralog分为in-paralog和out-paralog。in-paralog是指同一个物种的paralogs的分化发生在物种分化之后，这样的话，代表in-paralogs之间的序列相似性比其orthologs的相似性要高；通过OrthoMCL的原理可以看出，很好得到分析。而out-paralog是指paralogs的分化发生在物种分化以前，这代表out-paralogs的序列之间的相似性比某一个物种的orthologs的相似性要低；这样是很不好分析的，因为不好定阈值，或者得到的结果不易得不到大众的认可；OrthoMCL也没进行out-paralog的分析；当然，也可以只将一个query proteome输入到orthomcl来进行分析，得到的是所有paralog分析结果，包含了out-paralog。
Variant calling 软件
1. SNP和INDEL caller
1.1 GATK
1.2 SAMtools
1.3 DISCOVAR
DISCOVAR的官方博客：http://www.broadinstitute.org/software/discovar/blog/。该软件和ALLPATHS-LG是同一个组织开发的，是在ALLPATHS-LG的基础上开发的，用于基因组的组装和Variants的发掘。2013.3 公布出了该软件。
1.4 PyroHMMvar
PyroHMMvar: a sensitive and accurate method to call short INDELs and SNPs for Ion Torrent and 454 data. 2013.8.28发表在bioinformatics上。
samtools常用命令详解
samtools的说明文档：http://samtools.sourceforge.net/samtools.shtml
samtools是一个用于操作sam和bam文件的工具合集。包含有许多命令。以下是常用命令的介绍
1. view
view命令的主要功能是：将sam文件转换成bam文件；然后对bam文件进行各种操作，比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。
bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。
view命令中，对sam文件头部的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]
默认情况下不加 region，则是输出所有的 region.
Options: -b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ ,
are taken into account.
-1 fast compression (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-s FLOAT fraction of templates to subsample; integer part as seed [-1]
-? longer help
例子：
将sam文件转换成bam文件 $ samtools view -bS abc.sam > abc.bam $ samtools view -b -S abc.sam -o abc.bam 提取比对到参考序列上的比对结果 $ samtools view -bF 4 abc.bam > abc.F.bam 提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可 $ samtools view -bF 12 abc.bam > abc.F12.bam 提取没有比对到参考序列上的比对结果 $ samtools view -bf 4 abc.bam > abc.f.bam 提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式 $ samtools view abc.bam scaffold1 > scaffold1.sam 提取scaffold1上能比对到30k到100k区域的比对结果 $ samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam 根据fasta文件，将 header 加入到 sam 或 bam 文件中 $ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
2. sort
sort对bam文件进行排序。
Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix> -m 参数默认下是 500,000,000 即500M（不支持K，M，G等缩写）。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间。 -n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。
例子：
$ samtools sort abc.bam abc.sort
$ samtools view abc.sort.bam | less -S
3.merge
将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件不需要则是已经sort过了的。
Usage: samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam>[...]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to <out.bam> [in1.bam]
Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users
must provide the correct header with -h, or uses Picard which properly maintains
the header dictionary in merging.
4.index
必须对bam文件进行默认情况下的排序后，才能进行index。否则会报错。
建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。比如samtool的tview命令就需要；gbrowse2显示reads的比对图形的时候也需要。
Usage: samtools index <in.bam> [out.index]
例子：
以下两种命令结果一样
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
5. faidx
对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列
Usage: samtools faidx <in.bam> [ [...]] 对基因组文件建立索引 $ samtools faidx genome.fasta 生成了索引文件genome.fasta.fai,是一个文本文件，分成了5列。第一列是子序列的名称； 第二列是子序列的长度；个人认为“第三列是序列所在的位置”，因为该数字从上往下逐渐变大， 最后的数字是genome.fasta文件的大小；第4和5列不知是啥意思。于是通过此文件，可以定 位子序列在fasta文件在磁盘上的存放位置，直接快速调出子序列。 由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列 $ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta
6. tview
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
Usage: samtools tview <aln.bam> [ref.fasta] 当给出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。 按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第 10号scaffold的第1000个碱基位点处。 使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。 使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。 Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离 可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示； 20～30黄色；10～20绿色；0～10蓝色。 使用点号'.'切换显示碱基和点号；使用r切换显示read name等 还有很多其它的使用说明，具体按 ？ 键来查看。
7. flagstat
给出BAM文件的比对结果
Usage: samtools flagstat <in.bam> $ samtools flagstat example.bam 11945742 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数 0 + 0 duplicates 7536364 + 0 mapped (63.09%:-nan%) #总体上reads的匹配率 11945742 + 0 paired in sequencing #有多少reads是属于paired reads 5972871 + 0 read1 #reads1中的reads数 5972871 + 0 read2 #reads2中的reads数 6412042 + 0 properly paired (53.68%:-nan%) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值 6899708 + 0 with itself and mate mapped #paired reads中两条都比对到参考序列上的reads数 636656 + 0 singletons (5.33%:-nan%) #单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。 469868 + 0 with mate mapped to a different chr #paired reads中两条分别比对到两条不同的参考序列的reads数 243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一个，只是其中比对质量>=5的reads的数量
7. depth
得到每个碱基位点的测序深度,并输出到标准输出。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.bam> [...]
8. 其它有用的命令
reheader 替换bam文件的头
$ samtools reheader <in.header.sam> <in.bam>
cat 连接多个bam文件，适用于非sorted的bam文件
$ samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
idxstats 统计一个表格，4列，分别为”序列名，序列长度，比对上的reads数，unmapped reads number”。第4列应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数。
$ samtools idxstats <aln.bam>
9. 将bam文件转换为fastq文件
有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序：http://www.hudsonalpha.org/gsl/information/software/bam2fastq
$ wget http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz $ tar zxf bam2fastq-1.1.0.tgz $ cd bam2fastq-1.1.0 $ make $ ./bam2fastq <in.bam>
10. mpileup
samtools还有个非常重要的命令mpileup，以前为pileup。该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件，在samtools的安装文件夹中可以找到。
最常用的参数有2： -f 来输入有索引文件的fasta参考序列； -g 输出到bcf格式。用法和最简单的例子如下
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
$ samtools mpileup -f genome.fasta abc.bam > abc.txt
$ samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf
$ samtools mpileup -guSDf genome.fasta abc.bam | \
bcftools view -cvNg - > abc.vcf
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：
scaffold_1 2841 A 11 ,,,...,.... BHIGDGIJ?FF scaffold_1 2842 C 12 ,$,,...,....^I. CFGEGEGGCFF+ scaffold_1 2843 G 11 ,,...,..... FDDDDCD?DD+ scaffold_1 2844 G 11 ,,...,..... FA?AAAA<AA+ scaffold_1 2845 G 11 ,,...,..... F656666166* scaffold_1 2846 A 11 ,,...,..... (1.1111)11* scaffold_1 2847 A 11 ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG %.+....-..) scaffold_1 2848 N 11 agGGGgGGGGG !!$!!!!!!!! scaffold_1 2849 A 11 c$,...,..... !0000000000 scaffold_1 2850 A 10 ,...,..... 353333333
mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。其中第5列比较复杂,解释如下：
1 ‘.’代表与参考序列正链匹配。
2 ‘,’代表与参考序列负链匹配。
3 ‘ATCGN’代表在正链上的不匹配。
4 ‘atcgn’代表在负链上的不匹配。
5 ‘*’代表模糊碱基
6 ‘^’代表匹配的碱基是一个read的开始；’^’后面紧跟的ascii码减去33代表比对质量；这两个符号修饰的是后面的碱基，其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。
7 ‘$’代表一个read的结束，该符号修饰的是其前面的碱基。
8 正则式’\+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基；比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。
9 正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基；
pileup具体的参数如下：
输入参数 -6 Assume the quality is in the Illumina 1.3+ encoding. -A Do not skip anomalous read pairs in variant calling. -B Disable probabilistic realignment for the computation of base alignment quality (BAQ). BAQ is the Phred-scaled probability of a read base being misaligned. Applying this option greatly helps to reduce false SNPs caused by misalignments. -b FILE List of input BAM files, one file per line [null] -C INT Coefficient for downgrading mapping quality for reads containing excessive mismatches. Given a read with a phred-scaled probability q of being generated from the mapped position, the new mapping quality is about sqrt((INT-q)/INT)*INT. A zero value disables this functionality; if enabled, the recommended value for BWA is 50. [0] -d INT At a position, read maximally INT reads per input BAM. [250] -E Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt specificity a little bit. -f FILE The faidx-indexed reference file in the FASTA format. The file can be optionally compressed by razip. [null] -l FILE BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null] -M INT cap mapping quality at INT [60] -q INT Minimum mapping quality for an alignment to be used [0] -Q INT Minimum base quality for a base to be considered [13] -r STR Only generate pileup in region STR [all sites] 输出参数 -D Output per-sample read depth (require -g/-u) -g Compute genotype likelihoods and output them in the binary call format (BCF). -S Output per-sample Phred-scaled strand bias P-value (require -g/-u) -u Similar to -g except that the output is uncompressed BCF, which is preferred for piping. Options for Genotype Likelihood Computation (for -g or -u): -e INT Phred-scaled gap extension sequencing error probability. Reducing INT leads to longer indels. [20] -h INT Coefficient for modeling homopolymer errors. Given an l-long homopolymer run, the sequencing error of an indel of size s is modeled as INT*s/l. [100] -I Do not perform INDEL calling -L INT Skip INDEL calling if the average per-sample depth is above INT. [250] -o INT Phred-scaled gap open sequencing error probability. Reducing INT leads to more indel calls. [40] -P STR Comma dilimited list of platforms (determined by @RG-PL) from which indel candidates are obtained. It is recommended to collect indel candidates from sequencing technologies that have low indel error rate such as ILLUMINA. [all]
11. 使用bcftools
bcftools和samtools类似，用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件，后者为其二进制文件。
bcftools使用简单，最主要的命令是view命令，其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为：
$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]
[-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]
[-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]
[-T trioType] in.bcf [region]
$ bcftools view -cvNg abc.bcf > snp_indel.vcf
生成的结果文件为vcf格式，有10列，分别是：1 参考序列名；2 varianti所在的left-most位置；3 variant的ID（默认未设置，用’.’表示)；4 参考序列的allele；5 variant的allele(有多个alleles，则用’,’分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分号隔开；9 FORMAT of the genotype fields, separated by colon (optional)； 10 SAMPLE genotypes and per-sample information (optional)。
例如：
scaffold_1 2847 . A AACGGTGAAG 194 . INDEL;DP=11;VDB=0.0401;AF1=1;AC1=2;DP4=0,0,8,3;MQ=35;FQ=-67.5 GT:PL:GQ 1/1:235,33,0:63 scaffold_1 3908 . G A 111 . DP=13;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,5,7;MQ=42;FQ=-63 GT:PL:GQ 1/1:144,36,0:69 scaffold_1 4500 . A G 31.5 . DP=8;VDB=0.0034;AF1=1;AC1=2;DP4=0,0,1,3;MQ=42;FQ=-39 GT:PL:GQ 1/1:64,12,0:21 scaffold_1 4581 . TGGNGG TGG 145 . INDEL;DP=8;VDB=0.0308;AF1=1;AC1=2;DP4=0,0,0,8;MQ=42;FQ=-58.5 GT:PL:GQ 1/1:186,24,0:45 scaffold_1 4644 . G A 195 . DP=21;VDB=0.0198;AF1=1;AC1=2;DP4=0,0,10,10;MQ=42;FQ=-87 GT:PL:GQ 1/1:228,60,0:99 scaffold_1 4827 . NACAAAGA NA 4.42 . INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=40;FQ=-37.5 GT:PL:GQ 0/1:40,3,0:3 scaffold_1 4854 . A G 48 . DP=6;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,2,1;MQ=41;FQ=-36 GT:PL:GQ 1/1:80,9,0:16 scaffold_1 5120 . A G 85 . DP=8;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,5,3;MQ=42;FQ=-51 GT:PL:GQ 1/1:118,24,0:45
第8列中显示了对variants的信息描述，比较重要，其中的 Tag 的描述如下：
Tag Format Description AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele DP int Raw read depth (without quality filtering) DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise MQ int Root-Mean-Square mapping quality of covering reads PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2 PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias QCHI2 int Phred-scaled PCHI2 RP int # permutations yielding a smaller PCHI2 CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint UGT string Most probable genotype configuration without the trio constraint CGT string Most probable configuration with the trio constraint
bcftools view 的具体参数如下：
Input/Output Options: -A Retain all possible alternate alleles at variant sites. By default, the view command discards unlikely alleles. -b Output in the BCF format. The default is VCF. -D FILE Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null] -F Indicate PL is generated by r921 or before (ordering is different). -G Suppress all individual genotype information. -l FILE List of sites at which information are outputted [all sites] -N Skip sites where the REF field is not A/C/G/T -Q Output the QCALL likelihood format -s FILE List of samples to use. The first column in the input gives the sample names and the second gives the ploidy, which can only be 1 or 2. When the 2nd column is absent, the sample ploidy is assumed to be 2. In the output, the ordering of samples will be identical to the one in FILE. [null] -S The input is VCF instead of BCF. -u Uncompressed BCF output (force -b). Consensus/Variant Calling Options: -c Call variants using Bayesian inference. This option automatically invokes option -e. -d FLOAT When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0] 当有多个sample用于variants calling时，比如多个转录组数据或多个重测序 数据需要比对到参考基因组上，设置该值，表明至少有该<float 0~1>比例的 samples在该位点都有覆盖才计算入variant.所以对于只有一个sample的情况 下，该值设置在0～1之间没有意义，大于1则得不到任何结果。 -e Perform max-likelihood inference only, including estimating the site allele frequency, testing Hardy-Weinberg equlibrium and testing associations with LRT. -g Call per-sample genotypes at variant sites (force -c) -i FLOAT Ratio of INDEL-to-SNP mutation rate [0.15] -p FLOAT A site is considered to be a variant if P(ref|D) -t FLOAT Scaled muttion rate for variant calling [0.001] -T STR Enable pair/trio calling. For trio calling, option -s is usually needed to be applied to configure the trio members and their ordering. In the file supplied to the option -s, the first sample must be the child, the second the father and the third the mother. The valid values of STR are ‘pair’, ‘trioauto’, ‘trioxd’ and ‘trioxs’, where ‘pair’ calls differences between two input samples, and ‘trioxd’ (‘trioxs’) specifies that the input is from the X chromosome non-PAR regions and the child is a female (male). [null] -v Output variant sites only (force -c) Contrast Calling and Association Test Options: -1 INT Number of group-1 samples. This option is used for dividing the samples into two groups for contrast SNP calling or association test. When this option is in use, the following VCF INFO will be outputted: PC2, PCHI2 and QCHI2. [0] -U INT Number of permutations for association test (effective only with -1) [0] -X FLOAT Only perform permutations for P(chi^2)
使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行尤为重要，提供了4个数据：1 比对结果和正链一致的reads数、2 比对结果和负链一致的reads数、3 比对结果在正链的variant上的reads数、4 比对结果在负链的variant上的reads数。可以设定 （value3 + value4）大于某一阈值，才算是variant。比如：
$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ && ($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf > snp_indel.final.vcf
12. samtools rmdup
NGS上机测序前需要进行PCR一步，使一个模板扩增出一簇，从而在上机测序的时候表现出为1个点，即一个reads。若一个模板扩增出了多簇，结果得到了多个reads，这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。使用rmdup命令即可完成.
Usage: samtools rmdup [-sS] -s 对single-end reads。默认情况下，只对paired-end reads -S 将Paired-end reads作为single-end reads处理。 $ samtools input.sorted.bam output.bam
GBrowse的安装和使用
1. GBrowse的安装
GBrowse安装说明文档：http://gmod.org/wiki/GBrowse_2.0_Install_HOWTO
GBrowse的安装很少有能顺利安装成功的。需要不断的摸索，看文档，并搜索相关错误，google看别人是怎么解决的。有管一些我安装过程遇到的困难如下：
1.1 安装 io-lib和Bio-SCF, io-lib是安装Bio-SCF所必须的。在这一步容易出问题，导致无法安装。
$ wget http://downloads.sourceforge.net/project/staden/io_lib/1.13.0/io_lib-1.13.0.tar.gz $ make -j 8 $ sudo make install $ wget http://search.cpan.org/CPAN/authors/id/L/LD/LDS/Bio-SCF-1.03.tar.gz $ tar zxf Bio-SCF-1.03.tar.gz $ cd Bio-SCF-1.03 $ perl Makefile.PL $ make $ sudo make install
1.2 安装Bio::Graphics
当使用CPAN安装的时候test错误，使用手工安装，不用进行test。
2. Gbrowse的配置
2.1 Gbrowse文件
安装好Gbrowse后，须知道几个主要的文件存放位置，默认如下：
GBrowse-2.54/ 解压的GBrowse安装目录，其中bin目录下有很多重要的程序，比如创建Gbrowse2的用户的程序等。
/etc/gbrowse2/ 存放Gbrowse的配置文件，有全局的配置文件 GBrowse.conf 和 自己建立的一个物种的配置文件 species.conf;
/var/www/html/gbrowse2/ 存放gbrowse2的一些网页文件，比如gbrowse2的使用教程等；
/var/www/cgi-bin/gb2/ gbrowse2的web程序可执行文件，
/var/lib/gbrowse2/ gbrowse2的数据库目录和用户目录等。需要修改数据库目录的用户拥有者，以便于导入数据。
$ sudo chown chenlianfu /var/lib/gbrowse2/databases/
2.2 /etc/gbrowse2/GBrowse.conf 几个可能需要需要定制的配置：
1. max_render_processes = 12 设置渲染需要的最大CPU线程数 2. image widths = 450 640 800 1024 1280 1440 设置基因组浏览器显示的宽度，可选的值，在perference项中进行设置时可选的值 default width = 1024 默认的值 3. show sources = 1 默认下开启下拉菜单 数据源 ，以利于转移到其它物种的基因组浏览中。 4. #include "themes/warm_colors" 基因组浏览器的3个主题，此为默认的主题。 # #include "themes/transparent_colors" # #include "themes/solid_gray_colors" 5. default source = yeast 设置基因组浏览器中默认的物种，即http://host/cgi-bin/gb2/gbrowse/默认所访问的物种。 6. [lentinula_edodes] 设置所要建立的物种的基因组浏览器的配置文件路径 description = Lentinula edodes Genome path = lentinula_edodes.conf 7. 在GBrowse.conf同目录下有个languages的文件夹，其中gbrowse2默认使用zh.pm模块，于是GBrowse的显示文字默认是繁体中文，可以使用其目录下的其它语言模块代替。
2.3 Data Source Sections的配置
对所需要浏览基因组的物种，则需要建立一个专门的配置文件，比如本文中的 /etc/gbrowse2/lentinula_edodes.conf 文件。
[GENERAL] 全局设置
restrict = require user chenlianfu hzaumycology 设置该基因组浏览器的访问权限，只有chenlianfu和hzaumycology能访问。使用命令gbrowse_create_account.pl来创建gbrowse用户
description = Lentinula edodes Genome Database 对数据库的描述
database = gene_Prediction_EVM 数据库的名字，这个一定要有，不然会提示错误。当然如果将只有一个数据库，并放入到GENERAL中，则不需要该项
initial landmark = scaffold_1:20000..40000 初始访问数据库时候显示的区域
default tracks = Genes_EVM 初始访问数据库时候显示的tracks
Genes_Augustus
metadata = 对数据库的一些描述
-description Lentinula edodes Genome (strain: W1-26)
-maintainer Lianfu Chen
-created 2013-05-28
-modified 2013-05-28
-authority hzaumycology
-coordinates http://www.hzaumycology.com/
-coordinates_version 1
-source Scaffold
-testrange scaffold_1:103361..105454
-species Lentinula edodes W1-26
example = scaffold_1 给出的几个例子用于选择
scaffold_1000:2164..4463
#################################
# database definitions 数据库设置。
#################################
[gene_Prediction_EVM:database] 第一个数据库。该数据库稍微特殊些，在全局中使用该数据库，需要将fasta文件和scaffold信息导入到该数据库。
db_adaptor = Bio::DB::SeqFeature::Store 数据库的读取方法
db_args = -adaptor DBI::mysql 使用mysql数据库
-dsn lentinula_edodes_EVM mysql数据库中的数据库名
-user chenlianfu mysql数据库的可写用户
-password 1234567 用户的密码
search options = default 该数据库中的搜索选项
[gene_Prediction_Augustus:database] 另外的一个数据库。将不同的基因注释放入不同的数据库，然后放入不同的tracks，有利于阅读和使用基因组。
db_adaptor = Bio::DB::SeqFeature::Store
db_args = -adaptor DBI::mysql
-dsn lentinula_edodes_Augustus
-user chenlianfu
-password 1234567
search options = default
########################
# Default glyph settings 默认的glyph设置
########################
[TRACK DEFAULTS]
glyph = generic glyph的默认类型
height = 10 glyph的高度
r = black
font2color = blue
label density = 25 当labels比该数目多的时候，则关闭labels的显示以节约空间
bump density = 100 当features的数目多于该值的时候， 则不在垂直方向上显示features，它们都被限制在了一个水平线上。
link = AUTO 点击feature的时候，链接到feature的信息文档中。
################## TRACK CONFIGURATION ####################
# the remainder of the sections configure individual tracks 设置track
###########################################################
#[Genes_EVM:overview] 将该track置于overview中，而不是detail中，此时，分类自动成为overview。好处是在整体上看到该track的特征，但是细节上就没法放大了。
#[Genes_EVM:region] 同上。
[Genes_EVM] track名
database = gene_Prediction_EVM track所用到的数据库
feature = gene track所用的feature
glyph = gene track的glyph
starnd_arrow = 1 glyph具有方向性；有些glyph内在就具有或不具方向性，设置该值则不影响。
bgcolor = peachpuff 颜色
decorate_introns = 1 intron显示方法
label_transcripts = 1
draw_translation = 1
category = Genes track所属的分类，对应着gbrowse2的“Select Tracks"的分类
label_transcripts = 1
visible = 1 初始访问数据库时候显示该tracks
key = Genes Predictions Intergrated by EVM track在浏览器中的名称
citation = EVM was used to integrate the genes prediction results of Augustus, SNAP and GeneMarkES ;PASA was used to add UTR annotations and Alternatively spliced isoforms. 该track的介绍。
[Genes_Augustus]
database = gene_Prediction_Augustus
#feature = gene:AUGUSTUS feature为gff3文件的type：source。这样做的话，就可以不必建多个mysql数据库，只需要把source设置好即可。
feature = gene
glyph = gene
bgcolor = peachpuff
decorate_introns = 1
label_transcripts = 1
draw_translation = 1
category = Genes
label_transcripts = 1
key = Genes Predicted by Augustus
3. Gbrowse的数据的导入
3.1 EVM将多个基因组预测结果进行融合后，使用PASA加上5’和3’端非翻译区后，得到最终基因组注释结果。将该文件导入到mysql数据库 gene_Prediction_EVM 中。由于该最终的注释结果文件中mRNA的Nama的值包含的字符过长，需要进行缩短，以利于阅读；该文件中的genes也需要进行排序。因此编写程序提取出适合于gbrowse导入的gff3文件。可用于导入到gbrowse2的gff3文件有一些特点：
首先，每个scaffold、chromosome或contig之前要有一行指定其feature和name。比如：
scaffold_1 . scaffold 1 322871 . . . Name=scaffold_1
这样gbrowse才能识别scaffod_1是属于scaffold类型，并有个Name是scaffold_1。才会将其在基因组浏览器中显现出来。而正常的gff3文件是没有这样一行的。
gff3文件中mRNA中的Name有些太长，在基因组浏览器中的图片中占空间太大，需要重命名得简洁些；gff3文件中intron的可以去掉；gff3文件中将feature为transcript的改变为mRNA;去掉注释行等。以上这些都会影响gborows的显示结果。
使用如下命令来将基因预测信息导入到数据库中：
$ perl parse_evm_pasa_gff3.pl LEdodesGGTrinity.gene_structures_post_PASA_updates.26576.gff3 genome.fasta 该perl程序生成适合于导入gbrowse2的gff3文件gbrowse.gff3 和 protein.fasta文件.后者为预测的蛋白组文件。 $ mysql -h localhost -u root -p mysql > CREATE DATABASE gene_Prediction_EVM; 创建一个名为 gene_Prediction_EVM 的 mysql数据库 EOF mysql > Bye $ /usr/local/bin/bp_seqfeature_load.pl -c -a DBI::mysql -d gene_Prediction_EVM -u root -p password genome.fasta gbrowse.gff3 该程序能将gff3文件或fasta文件导入到数据库。其参数： -c 清空数据库 -a 导入的数据库类型 -d mysql数据库对应的数据库名称 -u mysql数据库用户名 -p mysql数据库密码 该程序导入的时间有点长，依据feature的数目，时间长短不一。对于1.2万个基因的注释，需要约10分钟导入完成。
3.2 将Augustus的基因预测信息导入
perl prepare_Augustus_gff3_for_gbrowse2.pl Agustus.gff3 > gbrowse2.gff3 $ mysql -h localhost -u root -p mysql > CREATE DATABASE gene_Prediction_Augustus; 创建一个名为 gene_Prediction_Augustus 的 mysql数据库 EOF mysql > Bye $ /usr/local/bin/bp_seqfeature_load.pl -c -a DBI::mysql -d gene_Prediction_Augustus-u root -p password genome.fasta gbrowse2.gff3
3.3 其它SNAP和PASA的基因预测信息导入和上面2中一致。不赘述。
4. NGS数据的导入
4.1 安装Bio::DB::Sam。需要有samtools安装，并且该samtools的安装和正常安装不一样。
$ wget http://garr.dl.sourceforge.net/project/samtools/samtools/0.1.19/samtools-0.1.19.tar.bz2 $ tar jxf samtools-0.1.19.tar.bz2 $ cd samtools-0.1.19 $ perl -p -i -e 's/CFLAGS.*/CFLAGS= -g -Wall -O2 -fPIC #-m64 #-arch ppc/' Makefile $ make -j 8 $ make clean 如果之前已经安装过samtools的话，需要修改makefile，再重新安装 $ make install $ wget http://search.cpan.org/CPAN/authors/id/L/LD/LDS/Bio-SamTools-1.38.tar.gz $ tar zxf Bio-SamTools-1.38.tar.gz $ cd Bio-SamTools-1.38 $ perl Build.pl $ ./Build $ sudo ./Build install
4.2 修改data source配置文件
[NGS_Genome:database]
db_adaptor = Bio::DB::Sam # 数据库的读取方法
db_args = -fasta genome.fasta # 基因组的fasta文件
-bam bowtie2.bam # NGS reads的比对结果
search options = default
[GenomeReadCoverageXyplot]
feature = coverage # 基因组测序的reads的coverage
glyph = wiggle_xyplot # 使用峰图来显示reads的覆盖度
database = NGS_Genome
height = 50
fgcolor = black
bicolor_pivot = 20 # 设定一个颜色变换的coverage值
pos_color = blue # 当coverage > 以上设置的数值，s使用蓝色
neg_color = red # 当coverage < 以上设置的该数值，s使用红色
category = GenomeReads
label = 0 # labels on wiggle tracks are redundant
key = Coverage (xyplot) of Genome NGS data
[GenomeReadCoverageDensity]
feature = coverage
glyph = wiggle_density # 使用密度来显示reads的覆盖度，覆盖度越高则线条颜色越深
database = NGS_Genome
height = 30
bgcolor = blue
bicolor_pivot = 5 # 小于该值，则线条无颜色
pos_color = blue
neg_color = red
category = GenomeReads
label = 0
key = Coverage (density plot) of Genome NGS data
[GenomeReads] # reads比对到基因组的图形显示
feature = match
glyph = segments
draw_target = 1
show_mismatch = 1
mismatch_color = red
database = NGS_Genome
bgcolor = blue
fgcolor = white
height = 5
label density = 50
bump = fast
category = GenomeReads
key = Reads of Genome NGS data
[GenomeReadsPair] # reads pair比对到基因组的图形显示
feature = read_pair
glyph = segments
database = NGS_Genome
draw_target = 1
show_mismatch = 1
bgcolor = sub {
my $f = shift;
return $f->attributes('M_UNMAPPED') ? 'red' : 'green';
}
fgcolor = green
height = 3
label = sub {shift->display_name}
label density = 50
bump = fast
connector = dashed
balloon hover = sub {
my $f = shift;
return '' unless $f->type eq 'match';
return 'Read: '.$f->display_name.' : '.$f->flag_str;
}
category = GenomeReads
key = Read Pairs of Genome NGS data
[GenomeReadsMappingQuality] # reas比对到基因组的Mapping质量图，高质量使用深蓝色表示，低质量使用浅蓝色显示
feature = match
glyph = segments
draw_target = 1
show_mismatch = 1
mismatch_color = red
database = NGS_Genome
bgcolor = sub {
my $feature = shift;
my $blueness = 255 - int($feature->qual * 2.40);
my $colour = chr(35) . sprintf("%X", $blueness) .
sprintf("%X", $blueness) . "FF";
return $colour;
}
fgcolor = black
height = 5
label density = 50
bump = fast
category = GenomeReads
key = Reads' Mapping Quality of Genome NGS data
4.3 将Bam文件放置到配置文件中对应的位置。
若在track中提示该错误：
Track rending error: No index file for bam file; try opeing file with -autoindex at /usr/local/lib64/perl5/Bio/DB/Sam.pm line 2064
则表示缺少bam文件对应的index文件。该文件以bai为后缀，使用samtools生成
$ samtools index geonme.bam
GFF3格式
GFF3的官方介绍：Generic Feature Format Version 3 （GFF3）
1. GFF3文件格式描述
GFF3格式文件为文本文件，分为9列，以TAB分开。控制符使用 RFC 3986 Percent-Encoding 编码。比如：%20 代表着ASCII的空格。
9列文件依次是：
1. seqid：参考序列的id。该id的取名不能以’>’开头，不能包含空格。
2. source ：注释的来源。如果未知，则用点（.）代替。一般指明产生此gff3文件的软件或方法。
3. type ：属性的类型。建议使用符合SO惯例的名称（sequence ontology，参看[[Sequence Ontology Project]]) ,如gene，repeat_region，exon，CDS等。
4. start position ：属性对应片段的起点。从1开始计数。
5. end position ：属性对应片段的终点。一般比起点的数值要大。
6. score ：得分，对于一些可以量化的属性，可以在此设置一个数值以表示程度的不同。如果为空，用点（.）代替。
7. strand ：“＋”表示正链，“－”表示负链，“.”表示不需要指定正负链。
8. phase ：步进。对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。可以是0，1或2，表示到达下一个密码子需要跳过的碱基个数。
对于其它属性，则用点（.）代替。
9. attributes ：属性
一个包含众多属性的列表。格式为“标签＝值”（tag=value）。不同属性之间以分号相隔。可以存在空格，不过若有“,=;”则用URL转义（URL escaping rule），同时TAB也需要转换为“%09”表示。所有以大写字幕开头的标签被保留，用于大众认可的用途，而以小写字母开头的标签则根据自己安排随意应用。
常用的标签有：
ID
Feature的标识。该ID具有唯一性。
Name
Feature的展示名称。Name的值在可视化的时候得到展示。因此，Name可以根据自己展示的需要随意取值。
Alias
Feature的第2个Name。
Parent
指明feature所从属的上一级ID。用于将exons聚集成transcript，将transripts聚集成gene。
Target
指明比对的目标区域，一般用于表明序列的比对结果。格式为”target_id start end [strand]”,其中strand是可选的(“+”或”-“), target_id中如果包含空格，则要转换成’%20’。
Gap
比对结果的gap信息，和Target一起，用于表明序列的比对结果。
Note
文本描述
Is_circular
表明featrue是否为环化的。用于环状基因组序列。
同一个tag如果有多个值，则多个值之间使用逗号隔开，比如：
Parent=AF2312,AB2812,abc-3
Alias=M19211,gna-12,GAMMA-GLOBULIN
能够使用多个值的tag有：Parent, Alias, Note, Dbxref and Ontology_term。
2. GFF3文件检测
检验GFF3格式文件： GFF3 Validator
Augustus training
1. Convert GFF file to Genbank format file
$ $AugustusHome/scripts/gff2gbSmallDNA.pl PASA.gff genome.fa 1000 genes.raw.gb
将gff文件转换成genebank格式，左右侧翼各加1000bp序列。gff文件可以由PASA将RNA-Seq的转录子比对到genome得到。而PASA得到的gff文件是有5’端非翻译区注释的，这样的信息会被trainig忽略。it is sufficient to have only the coding parts of the gene structure (CDS).
当然，genebank文件也可以使用NCBI的nucleotide数据库进行检索得到。
2. remove these problematic locis from genes.raw.gb
$ $AugustusHome/bin/etainig --species=SPECIES --stopCodonExcludedFromCDS=false genes.raw.gb 2> train.err $ cat train.err | perl -pe 's/.*in sequence (\S+): .*/$1/' > badgenes.lst $ $AugustusHome/scripts/filterGenes.pl badgenes.lst genes.raw.gb > genes.gb
第一条命令用于输出trainig过程中的错误信息，根据错误信息找到 badgenes，然后在去掉这些badgenes，剩下的genes用于training。
值得注意的是：
1. 至少有200个gene structures用于training，才能得到不错的结果。越多的gene，则training的效果越好；当然，达到1000个genes的时候，提升的效果就很小了。
2. 当有多于1000个基因的时候，则需要注重基因的质量，而不是数量了。要保证multi-exon genes的数目要多，这样用于tain the introns。并且gene structures越精确越好。
3. gene set should be non-redundant.如果2个不同的基因序列绝大部分的amino acid sequence是一致的，则去掉其中一个。推荐的条件是：gene set里面任意两个gene在amino acid level上的identity要不高于80%。可以使用blast来解决，由于80%的阈值算是比较高的，一般也就需要去除掉20多个基因。
3. Split gene structure set into training and test set
$ $AugustusHome/scripts/randomSplit.pl genes.gb 100
将genes.gb分隔成了genes.gb.test和genes.gb.train两个文件。其中前者为genes.gb中随机取出的100个genes，后者为剩下的genes。后者将用于不停地traning。
4. CREATE A META PARAMETERS FILE FOR YOUR SPECIES
$ $AugustusHome/scripts/new_species.pl --species=lentinula_edodes
假如我们要建立香菇物种的traning参数，则上命令建立了其参数文件和文件夹，不过文件内容是初始的。
注意的是，用于training的gene的最后一个CDS的最后3个碱基若不是终止密码子，则需要手动修改Lentinula_edodes_parameters.cfg文件，将其中的stopCodonExcludedFromCDS由默认的false改为true。
5. MAKE AN INITIAL TRAINING
$ $AugustusHome/bin/etrainig --species=lentinula_edodes genes.gb.train $ $AugustusHome/bin/augustus --species=lentinula_edodes genes.gb.test | tee firsttest.out $ grep -A 22 Evalustion firsttest.out
使用genes.gb.train做一次trainig，然后使用genes.gb.test来检测training的精确性。分别在nucleotide，exon和gene level上检测其sensitivity和specificity。
sensitivity表示被被检测出来的百分率；specificity表示检测出来的nucleotide,exon或gene和test set中的完全一致的百分率。
6. RUN THE SCRIPT optimize_augustus.pl
$ $AugustusHome/scripts/optimize_augustus.pl --species=lentinula_edodes --cpus=8 genes.gb.train $ $AugustusHome/bin/etrainig --species=lentinula_edodes genes.gb.train $ $AugustusHome/bin/augustus --species=lentinula_edodes genes.gb.test
1. optimize_augustus.pl所做的事情：
默认情况下，optimize_augustus.pl将genes.gb.train中的genes随机分成8等份，然后使用其中的7个等份的genes做training，另外的1个做精确性评估。这样相互下来，共有8个方案，每个方案取1个等份用于精确性评估，另外7个用于training。
进行一次随机分配后再运行10次training和精确性评估，即为一次预测，得到一个target value。该值是 base，exon和gene level上sensitivities和specificities的权重值。
每次预测，如果得到更高的target value，则修正参数文件中的值：lentinula_edodes_parameters.cfg。
默认下参数文件中有28项参数需要按一定顺序进行优化；一般情况下每个参数最多设置5个值各进行一次预测(即对一项参数而言，这设置的5个值其中可能有1个值是用于之前的预测，故每个参数优化需要运行最多5次预测)，取最大的target value对应的值为参数的值；对所有的参数进行优化一次是一轮，则5轮参数优化完毕后程序会停止运行(以1800左右个genes来进行training，则每次augustus对200多个gene进行预测需要1min，那么每个参数优化需要28*4*8*1min=896min=15h，5轮参数的优化总共需要75h，即3天)，或如果在一轮参数优化中没有improvements则提前停止运行。当然，如果等不及，也能手动停止程序运行。由于optimize_augustus.pl运行时间太长，最好使用screen来运行。
如果了解了上述运行原理，则可以视情形终止其运行，或保存配置文件后接着运行。
2. 在optimize_augustus.pl完成或中断之后，需要(re)train AUGUSTUS with genes.gb.train。然后在使用genes.test.gb进行预测的精确性检测，如果gene level sensitivity低于20%，则表明training set不够大，或者质量不够好，或者物种somehow special。
7. Training AUGUSTUS UTR parameters
这部分的Training则需要5’和3’端的UTR都存在的gene structure。