转录组分析<一>之Hisat2不完全使用指南一

一 本系列的开启原因

转录组分析已经成为当今生物信息分析领域最最最常见的一种组学分析了。并且近些年因为其价格的大幅度下降以及分析手段的日渐成熟，已经成为打开每个生物学问题的所必须经历的第一道门槛了。虽然这是一种常见的分析，但是很多人还是对这些问题不够明晰。正好，本公众号最近需要重新分析一批转录组数据，这里就一边分析一边总结。希望这一个系列能够给其他读者展示我已经踩过的一些坑，从而帮助读者开始完成这一分析，进入更深一步的功能验证环境。

那么，我们要做转录组分析，不可避免的就是要使用到hisat2、Stringtie和edgeR这三个软件，并且这三个软件也是目前转录组分析的主流流程。那么本系列前几篇将先介绍一下这三个软件的使用。考虑到原文档太过繁琐，我会尽量简化一下hisat的使用文档。但是，如果读者想要知道细节，或者觉得我介绍的不够准确，可以去官网查看。虽然尽可能的确保准确，但本公众号仍对介绍的准确性不做任何保证。不过，如果有错误，后续也会进行相应的更正。那么，接下来就让我们先了解一下在使用hisat2软件前需要做哪些准备吧。

一 建立索引

现如今，随着基因组的大量出现，转录组分析也更多的偏向于有参转录组，因此，hisat2首先提供了一个hisat2-build选项，用于对所使用的基因组进行处理，即建立索引，方便分析。

1 主参数

<reference_in>: 以逗号分隔的FASTA文件列表，如每个染色体的fa文件。如果只有一个文件，那么不需要任何逗号。

<ht2_base>: 要写入的索引文件的前缀名（会产生许多以此前缀名为开头的文件）。

2 options选项

-f: 用于比对的输入文件是FASTA文件。

-c: 在命令行中给出参考序列，即<reference_in>是一个用逗号分隔的序列列表，而不是FASTA文件的列表。

--large-index: 强制hisat2-build建立一个大的索引，即使参考序列的长度小于~40亿bp核苷酸长度。

-a/--noauto: 禁用默认行为，即hisat2-build根据可用内存自动选择。

--bmax <int>: 一个区块中允许的最大后缀数。

--bmaxdivn <int>: 一个区块中允许的最大后缀数，以参考长度的分数表示。

--dcv <int>: 使用<int>作为差异覆盖样本的周期。

--nodc: 禁用差值覆盖样本的使用。

-r/--noref: 不建立索引的NAME.3.ht2和NAME.4.ht2部分，而这些部分包含参考序列的位包版本，可以用于成对端比对。

-3/--justref: 只建立索引中NAME.3.ht2和NAME.4.ht2的部分，其中包含参考序列的位包版本，用于配对端比对。

-o/--offrate <int>: 为了将比对结果映射到参考序列上的位置，有必要用基因组上的相应位置来注释（"标记"）部分或全部Burrows-Wheeler行。

-t/--ftabchars<int ftab>: 是用于计算与查询的第一个<int>字符有关的初始Burrows-Wheeler范围的查找表。

--localoffrate<int>: 这个选项控制了在本地索引中标记多少行：索引器将每2^<int>行标记一次。

--localftabchars<int>: 本地ftab是本地索引中的查找表。

-p <int>: 构建并行线程（默认：1）。

--snp <path>: 提供SNPs的列表格式如下。

--haplotype <path>: 提供一个单倍型列表的文件，格式如下。

--ss <path>: 注意这个选项应该和下面的--exon选项一起使用，提供一个剪接点的列表。

--exon <path>: 注意这个选项应该和上面的--ss选项一起使用。提供一个外显子列表，可以使用hisat2_extract_exons.py（HISAT2软件包），从GTF文件中提取外显子。

--seed <int>: 使用<int>作为伪随机数发生器的种子。

--cutoff <int>: 只索引参考序列的前<int>个碱基（跨序列累积），而忽略其余部分。

-q/--quiet: hisat2-build在默认情况下是静默的。使用这个选项，hisat2-build将只打印错误信息。

-h/--help: 打印使用信息并退出。

--version: 打印版本信息并退出。

3 示例

上面的命令表明的含义有以下几点：线程为2，基因组序列为genome.fa，生成索引文件的前缀为genome。

二 生成可变剪切的文件信息

注意：最好使用使用注释的转录本建立的索引（如genome_tran或genome_snp_tran），这比使用这个选项效果更好。提供已经包含在索引中的剪接位点没有影响。

hisat2_extract_splice_sites.py: 包含在HISAT2软件包中，

genes.gtf: 基因注释文件

spliceites.txt: 提供给HISAT2的剪接点列表。

三 比对gap

本部分主要介绍一下比对过程中该如何根据测序长度和读长来定义双端测序中两个read中所应该产生的gap，方便在下一章中介绍如何设定参数进行比对。

gap计算方法：如果指定了-I 60，并且一个配对端排列由两个20bp的排列组成，在适当的方向上，它们之间有20bp的间隙，那么这个排列被认为是有效的（只要-X也被满足）。在这种情况下，19bp的间隙是无效的。如果同时使用了修剪选项-3或-5，-I约束将被应用于未修剪的配体。

gap大小与比对性能：-I和-X之间的差异越大，HISAT2的运行速度就越慢。这是因为-I和-X之间的差异越大，HISAT2就需要扫描一个更大的窗口来确定是否存在一致的排列。对于典型的片段长度范围（200到400个核苷酸），HISAT2是非常有效的。

六 惯例小结

通过构建索引、生成可变剪切文件，我们就可以获得除测序文件外，所有应该被准备好的数据了。那么，接下来，就轮到正式比对啦。

本公众号开发的相关软件，Multi-omics Hammer软件和Multi-omics Visual软件欢迎大家使用。文末是本公众号在其他平台的账户，也欢迎大家关注并多提意见。

简书：WJ的生信小院

公众号：生信小院

博客园：生信小院

最后，也欢迎各位大佬能够在本平台上：1传播和讲解自己发表的论文；2：发表对某一科研领域的看法；3：想要达成的合作或者相应的招聘信息；4：展示自己以寻找博后工作或者博士就读的机会;5:博导提供博后工作或者博士攻读机会，都可以后台给笔者留言。希望本平台在进行生信知识分享的同时，能够成为生信分析者的交流平台，能够实现相应的利益互补和双赢（不一定能实现，但是梦想总得是有的吧）。

另外，怎么说呢，投币也可，不强求，但奢求。