1.genome survey
-
数据过滤
去除测序原始数据中可能包含低质量、接头污染以及含 N 过高的 reads
NT比对
通过BLAST对下机数据过滤后的有效数据进行 NT 比对评估，如果有较高比例的序列同时比对到非近源物种的基因组上，可能是样品存在污染引起的Kmer分析
通过GenomeScope进行Kmer分析，Kmer分析可以初步判断样本的基因组大小 、杂合情况和重复序列信息。SOAP denovo 组装结果
用SOAP denovo对高通量测序数据进行初步组装,获取拼接结果和基因组大小-
GC 分布
GC百分比是一个物种基因组的重要特征之一 ,同时 GC 含量分布的集中性有时也可以反映杂合、重复以及是否有污染等特征
Genome survey
2. genome assembly
- PacBio测序组装
1.1 基于 PacBio 平台的基因组测序
PacBio测序平台基于独特的单分子实时测序技术（ Single Molecule Real Time SMRT ），其应用了边合成边测序的原理，以 SMRT 芯片为测序载体，利用芯片
上纳米级别的零模波导孔（ zero mode waveguides, ZMWs ）和荧光标记的核苷酸焦磷酸链（ Phospholinked nucleotides ），每个零模波导孔中都能够包含一个 DNA聚合酶及一条 DNA 样品链进行单分子测序，并实时检测插入碱基的荧光信号，连续不断地读取该 DNA 样品链的碱基信息，因此可以极大的提高所获得的测序片段长度。
1.2 基于PacBio的基因组组装
1.2.1 基因组组装
Canu 在加载 reads 后将对 k mer 进行计数，用于计算序列间的 overlap 。 Canu 分为纠错、修整和组装三个步骤，在 reads 纠错时从 overlap中挑选一致性序列替换原始的噪声 reads ，修整时使用 overlap 区分并 确定 reads中 的 高质量区域 和 需要修整 的低质量区域 ，保留单个最高质量的序列块。最终组
装时， Canu 将根据一致的 overlap 对序列进行排列 layout ，得到 contig 的 组装结果 。
1.2.2 基因组纠错
Pilon 以 FASTA和 BAM 文件作为输入，根据比对结果对输入的参考基因组进行纠错。整个过程包括了比对、标记重复、过滤高质量比对的 read 、 polish 。 - 组装评估
2.1 GC分布评估
使用BWA 软件将过滤后的测序数据比对到基因组的组装结果上，使用 soap.coverage 软件将所有碱基的 soap 比对结果进行统计，得到基因组的单碱基深度信息；以 10kb 为滑窗长度在基因组上无重复前进，对每个滑窗长度内的碱基平均测序深度与 GC 含量进行统计。
2.2 深度/覆盖度评估
为了评估组装的准确性，选取小片段文库reads 采用 BWA 软件比对到组装的基因组上，统计 reads 的比对率、覆盖基因组的程度及深度的分布情况，评估组装的完整性和测序的均匀性。
2.3 SNP评估
单核苷酸多态性指在基因组上单个核苷酸变异形成的遗传标记，其数量很多，多态性丰富。我们利用 Samtools 等工具对 BWA 比对结果经过染色体坐标排序、去掉重复的 reads 等处理，进行 SNPCalling ，并对原始结果进行过滤。
2.4 BUSCO评估
BUSCO（Benchmarking sets of Universal Single Copy Orthologs ）评估是利用单拷贝直系同源基因，抽样了数百个基因组，从中选择单拷贝直系同源大于 90%的基因作为直系同源基因集，并对基因组组装结果中对同源基因的情况进行比对，以此评估基因组组装的完整性。
3.genome annotation
基因组注释主要包括四个方面内容 ：重复序列注释、基因结构预测 、基因功能注释、 ncRNA注释。
1 .重复序列注释
基因组重复序列
TRF（Tandem Repeat Finder）；RepeatMasker；RepeatProteinMask；De novo
2 .基因结构预测
de novo预测 使用软件 Augustus；homolog注释(近缘物种）；Transcript预测（RNA-seq数据）
使用Glean软件对上述三种证据集进行整合， 然后过滤掉部分基因，得到近缘物种基因结构统计结果。使用BUSCO软件对基因集得完整性进行评估。
3 .基因功能注释
用基因结构预测得到的蛋白质序列与Interpro、 KEGG、 Swissprot、 Tremble等已知蛋白库进行比对。
4 .非编码RNA注释
通过与已知ncRNA 库进行比对
- 重复序列注释
重复序列可分为串联重复序列(Tandem repeat)和散在重复序列 (Interspersed repeat)两大类。其中串联重复序列包括有微卫星序列，小卫星序列等；散在重复序列又称转座子元件，包括以 DNA-DNA方式转座的 DNA转座子和反转录转座子 (retrotransposon)。常见的反转录转座子类别有 LTR LINE和 SINE等。
使用两种方法进行重复序列注释：
基于RepBase (http://www.girinst.org/repbase) 的同源预测方法 (软件：RepeatMasker);
基于自身序列比对(软件 : RepeatModeler、 Piler、 RepeatScount)及重复序列特征 (软 : Trf和 LTR-FINDER)的 De novo预测方法。 - 基因结构预测
使用 GLEAN、 EVM、 Maker等软件对不同的证据集 进行整合，去除冗余，得到完整的基因集。其中，证据集的来源主要有三种：
Homolog预测 , 挑选 3~10个 近缘 物种，使用 Genewise软件进行 同源 预测；
De novo 预测，根据基因自身的结构特征，使用 Augustus, Genscan, Genemark, Glimmer, GeneID, SNAP等软件对基因组序列进行从头预测；
Transcript预测，有两种策略，一种是 先 使用 Trinity软件对 RNA-seq数据进行组装，然后使用 Blat软件将组装结果 与基因组序列进行比对，接着再用Transdecoder软件将比对结果去冗余；另一种是 利用 Hisat2软件直接将 RNA-seq数据与基因组序列进行比对，然后用 Stringtie软件 将比对结果转化为基因格式 - 基因功能注释
借助于外源蛋白数据库(SwissProt、 TrEMBL、 KEGG、 InterPro、 COG、NT、 NR和 GO)对基因集中的蛋白进行功能注释。 - 非编码RNA注释注释
非编码RNA（（Non-coding RNA）是指不编码蛋白质的）是指不编码蛋白质的RNA,包括包括rRNA，，tRNA，，snRNA和和miRNA等等。。这些这些RNA的共同特点是都能从基因组上转录而的共同特点是都能从基因组上转录而来，但是不翻译成蛋白，在来，但是不翻译成蛋白，在RNA 水平上就能行使生物学功能。水平上就能行使生物学功能。miRNA可降解靶可降解靶基因或抑制靶基因翻译成蛋白质，具有沉默基因的功能基因或抑制靶基因翻译成蛋白质，具有沉默基因的功能; tRNA、、rRNA直接参与直接参与蛋白质的合成蛋白质的合成; snRNA主要参与主要参与RNA前体的加工，是前体的加工，是RNA剪切体的主要成分。剪切体的主要成分。根据tRNA的结构特征，利用的结构特征，利用tRNAscan-SE软件来寻找基因组中的软件来寻找基因组中的tRNA序序列；由于列；由于rRNA具有高度的保守性，因此可以选择近缘物种的具有高度的保守性，因此可以选择近缘物种的rRNA序列作为参序列作为参考序列，通过考序列，通过BLASTN比对来寻找基因组中的比对来寻找基因组中的rRNA；另外，利用；另外，利用Rfam家族的家族的协方差模型，采用协方差模型，采用Rfam自带的自带的INFERNAL软件可预测基因组上的软件可预测基因组上的miRNA和和snRNA序列信息。序列信息。
4.Comparative genomics
与其近缘物种进行比较基因组学分析，主要包括基因家族聚类、系统发育树、分歧时间、基因家族的扩张和收缩、 基因家族功能富集等分析。
1.基因家族聚类
使用OrthoMCL软件 来进行基因家族聚类，使用BLASTP软件比对所有物种的蛋白序列， e值阈值设为 1e-5然后用OrthoMCL 软件对所有基因进行聚类。
2.系统发育
使用单拷贝基因家族构建系统发育树。首先使用 MUSCLE 比对单拷贝基因家族的蛋白序列 然后 基于比对结果，将蛋白序列反转录为 CDS 序列 提取每个比对的 4 倍简并位点串联成 super gene 然后 使用PhyML 和 Mrbayes 分别 进行构树 ，获得树形图文件最终使用 Figtree 将树形图文件图像化 。
3.分歧时间
使用PAML中的 MCMCTREE 来估计物种 分歧 时间 使用“ Correlated molecular clock ”分子钟模型和HKY85 ”核 酸 替换模型,校正点分歧时间来自TimeTree http://www.timetree. 。
4.基因家族扩张和收缩
根据基因家族聚类结果和物种间的系统发育关系，使用CAFE 进行基因家族扩张和收缩分析,对这些显著扩张的基因进行 KEGG 和 GO 富集分析.
使用PAML 中的 CodeML 进行正选择分析，选用“branch site” 模型得到受正选择基因(p<0.05).
5.LTR插入时间
逆转录转座子在插入宿主基因组时 两个 LTR 区域 通常是相同的。随着时间的推移，核苷酸的替换会导致两个 LTR 序列 出现 差异。 在核苷酸 替换率已知的情况下 ，可以根据两个 LTR 之间的差异 数估算插入时间,我们使用 LTR_FINDER 寻找基因组中的 LTR 区域，然后使用MUSCLE 进行多序列比对并使用 DISTMAT 计算距离矩阵，最后根据公式 T = K2P/2r 计算 LTR 插入时间 .
6.基因组共线性
我们使用BLASTP 来检测物种间的直系同源基因，选择最优的比对结果，使用 MCscan 识别同源基因区块，然后选择更长的同源基因区块进行下一步作图。
7.全基因组复制分析
由于同义突变在物种进化过程中不受自然选择， 所以在一定时间尺度下，同义突变的速率可以衡量物种进化的时间 也可以用来衡量物种全基因组复制时间发生的时间和次数 。四倍简并位点颠换率 (4dTv )分布 和同义替换率 (Ks) 分布 常被用来进行全基因组复制分析。 首先 使用 BLASTP 来检测物种内的旁系同源基因和物种间的直系同源基因 然后使用 MCscan 软件识别同源基因区块，计算同源基因区块的 4dTv 值,使用KaKs_caculator计算旁系同源基因的Ks值。
