hisat2构建索引流程(hisat2建立索引太耗内存)

扫码添加渲大师小管家，免费领取渲染插件、素材、模型、教程合集大礼包！

大家好，今天来介绍hisat2构建索引流程的问题，以下是渲大师小编对此问题的归纳和整理，感兴趣的来一起看看吧！

hisat2索引构建命令

命令行：

hisat2-build[options]*

注意：

       如果使用--snp，--ss和/或--exon选项，hisat2-build将需要大约200GB的运行内存来满足人类基因组规模大小的基因组的索引构建，因为建立索引涉及到graph construction。否则，就可以用8GB的运行内存在个人电脑构建索引了。

主要参数

       以逗号分隔的FASTA文件列表，其中包含要比对的参考序列，例如chr1.fa，chr2.fa，chrX.fa，chrY.fa；如果指定了-c，则可以具体的GGTCATCCT，ACGGGTCGT，CCGTTCTATGCGGCTTA序乎洞列。

       要写入的索引文件的basename；默认情况下，hisat2-build写入文件名为NAME.1.ht2, NAME.2.ht2, NAME.3.ht2, NAME.4.ht2, NAME.5.ht2, NAME.6.ht2, NAME.7.ht2, NAME.8.ht2的文件。就是文件前缀NAME。

选项

1.-f

参考基因组输入文件(指定为_)是FASTA文件(通常具有.fa、.mfa、.fna或类似的扩展名)。

2.-c

在命令行上给出参考序列。也就是说是一个逗号分隔的序列列表，而不是一个FASTA文件列表。

3.--large-index

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

4.-a/--noauto

禁用hisat2-build根据可用内存自动选择--bmax，--dcv和[--packed]参数的这一默认行为。相反，用户可以为这些参数指定值。如果内存在构建索岁圆枯引期间耗尽，将输出错误信息；由用户决定是否尝试新的参数。

5.--bmax

block 中允许的最大后缀数。允许每个block使用更多的后缀可以加快索引构建速度，但会增加内存的峰值使用。设置此选项将覆盖以前对--bmax或--bmaxdivn的任何设置。--bmaxdivn默认值是4。这是默认自动配置的；使用-a/- noauto则可以手动配置。

6.--bmaxdivn

block中允许的最大后缀数，表示为参考序列长度的一部分。设置此选项将覆盖以前对--bmax或--bmaxdivn的任何设置。默认值: --bmaxdivn 4。这是默认自动配置的；使用-a/- noauto则可以手动配置。

7.--dcv

使用作为 difference-cover 样本的period。较大的period产生较少的内存开销，但可能使后缀排序变慢，特别是如果存在重复。必须是2的整数幂且必须不大于4096。默认值:1024。这是默认自动配置的；使用-a/- noauto手动配置。

8.--nodc

禁用difference-cover样本的使用。在最坏的情况下(最坏的情况是极度重复的参考序列)，后缀排序变成了二次排序。默认值:关闭。

9.-r/--noref

不要构建名称为NAME.3.ht2 、NAME.4.ht2的索引部分，它包含参考序列的bitpacked版本， 用于双端测序的比对。

10.-3/--justref

只构建名称为NAME.3.ht2 、NAME.4.ht2的索引部分，它包含引用序列的bitpacked版本，用于双端测序的比对。

11.-o/--offrate

为了将比对结果映射回参考序列上，有必要用基因组上相应位置标注(标记)部分或全部的Burrows-Wheeler

rows。-o/- offrate统计有多少行被标记:索引器将标记每2^行。标记更多的行可以使序列-位置查找更快，但是需要更多的内存来在运行时保存注释。默认值为4(每16行标记一次；对于人类基因组来说，注释大约有680兆字节)。

12.-t/--ftabchars

ftab是用于计算的第一个字符的初始Burrows-Wheeler范围的查找表。较大的将生成较大的查找表，但查腔乎询时间更快。ftab的大小为4^(+1)字节。默认设置为10 (ftab为4MB)。

13.--localoffrate

这个选项统计在本地索引中标记多少行:索引器将标记每2^行。标记更多的行可以使引用位置查找更快，但是需要更多的内存来在运行时保存注释。默认值为3(每标记第8行，每个本地索引大约占用16KB)。

14.--localftabchars

本地ftab是本地索引中的查找表。默认设置为6 (ftab为每个本地索引8KB)。

15.-p

并行运算线程数(默认:1)。

16.--snp

提供一个snp列表(HISAT2自己的格式)，如下(五列)。

SNPID snp type (single, deletion, or insertion) chromosomename zero-offset based genomic position of a SNP alternative base (single), the length of SNP (deletion), or insertion sequence(insertion)

例如：rs58784443，single，13，18447947，T

hisat2_extract_snps_haplotypes_UCSC.py(在HISAT2包中)从dbSNP文件(例如http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp144Common.txt.gz)中提取SNPs和haplotypes。或者hisat2_extract_snps_haplotypes_VCF.py从VCF文件中提取SNPs和haplotypes(例如 ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 /supporting/GRCh38_

positions/ALL.chr22.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz)。

17.--haplotype

提供一个单倍型列表(使用HISAT2自己的格式)，如下所示(五列)。

HaplotypeID chromosome name zero-offset based left coordinate ofhaplotype zero-offset based right coordinate of haplotype a comma separated list of SNP ids in the haplotype

例如：ht35，13，18446877，18446945，rs12381094,rs12381056,rs192016659,rs538569910

关于如何提取单倍型，请参阅上面的-snp选项。这个选项不是必需的，但是单倍型信息可以防止索引构造激增，并大幅减少索引大小。

18.--ss

注意， 这个选项应该与下面的--exon选项一起使用。提供一个 剪切位点 列表(HISAT2自己的格式)，如下(四列)：

chromosomename zero-offset based genomic position of the flanking base on theleft side of an intron zero-offset based genomic position of theflanking base on the right strand

可以用hisat2_extract_splice_sites.py 脚本从GTF注释文件中提取剪切位点文件。

19.--exon

注意 ，这个选项应该与上面的--ss选项一起使用。提供一个外显子列表(HISAT2自己的格式)，如下(三列)：

chromosomename zero-offset based left genomic position of an exon zero-offset based right genomic position of an exon

可以用hisat2_extract_exons.py脚本从GTF注释文件中提取外显子文件。

20. --seed

伪随机数生成器的种子

21.--cutoff

22.-q/--quiet

静默输出，只会输出错误信息。

23.-h/--help

帮助文档

24.--version

软件版本信息

hisat2构建基因组index

选用：gencode的GRCm38.fasta和gtf

第渣信岩一步：转录本的index：

第如御二步：基因坦茄组的index：

插件 Hisat2+StringTie 本地界面化（Win/Mac）点点点完成转录组数据分析

早前，我已经通过插件的方式，让所有 TBtools 用户，都能完成 RNAseq 数据分析，从测序原始数据到基因表达量，使用的是一个曲线救国的策略，即直接使用 kallisto，跳过读段回帖，直接进行读段计数。
目前，更为常用的 RNAseq 上游数据分析流程，应该还是读段回帖之后进行读段计数。一般情况下，使用的软件是：star / hisat2。前者对内存要求高，而后者做了专门的层级索引设计，可以在个人电脑甚至是笔记本（比如我的笔记本 8G 内存）上完成绝大多数物种的转录组读段回帖。
于是，前几天对应的插件都开发出来了，即 hisat2-build 和 hisat2-align。走到这里，我们还能更进一步，做更有意义的事情。
早前的Kallisto本身是依赖于基因组基因结构注释的，其准确程度颇受已有注释的影响，而hisat2等基于回帖的，我们可以进一步做注释“自动校正”以及新转录本或基因挖掘。更为全面一些。这些，则往往常用的软件凳猛凯是 Stringtie。
Stringtie目前为止，并没有人编译windows版本（有点像 MCScanX 当初的情况），于是我做了尝试，调整了源码，并编译了（注：苹果用户 Mac 直接有可用程序，不存在这个问题）。折腾折腾，现在我们可以直接在 TBtools 里面进行转录组的有参考组装以及基于读段回帖的表达量估计。
于是，有必要整理一个教程，理清四个插件的使用，步骤如下：

插件直接从 TBtools 插件商店获取。主要到推荐从高速商店获取，参考前述推文《Plugin 高知桐速版插件商店！我又有一个绝妙的 idea》。

设置基因组序列文件，用于建立索引

点击Start，并等待即可

如此，即完成了索引构建。

总的来说，基本没什么特别要注意的，除非数据是链特异的，那么最好设置一下。另外是，是否很关注多匹配的reads，如repeat区域，那么可以考虑提高max hits。
恩，Threads 参数控制的是并行任务数目，而不是stringtie运行时的线程数。简单来说，假设输入的是 6 个样品，Threads设置为 2 ，那么同时会有最多两个样品在进行组装（即并行）。
输出结果会放置在输出目录下，

大体如下，

可能唯一需要注意的就是....并行任务数，可参考前述推文，其实常常也无需修改，一般按照电脑有多少个线程，保留2个，剩下的都可以用上试试。

示例数据只有一个样品，所以只组装出一个XXXX.assembly.gtf。无论有多少个输入样品，最终每个样品都会被独立组装，最后合并成一个 merged.stringtie.gtf。这个文件，可用于后续任何分析（亦即，枣唤完成了转录本组装）。

Stringtie 除了进行组装，还可以估算转录本以及基因的表达量。

按照要求设置文件即可，可能需要调整的就是read length，如果你想要得到 read counts，用于下一步差异表达分析的话。
运行后，可以看到在输出目录增加了 6 个文件。

插件均已上传至高速商店，

感兴趣地同样参考前述推文《Plugin 高速版插件商店！我又有一个绝妙的 idea》

今天是大年初一 ~~~
新年新气象，
祝所有 TBtools 用户朋友，
牛年大吉！

STEP4: 得到表达矩阵的流程

这是RNA-Seq 上游分析的大致流程，比对+定量。当然实验目的若只需要定量已知基因，也可以选择free-alignment 的流程工具如kallisto/Salmon/Sailfish，其优点是可用于RNA-seq的基因表达的快速定量，但是对于小RNA和表达量低的基因分析效果并不好（2018年刚发表的一篇文章对free-alignment 的工具尺肆进行了质量评估，doi: https://doi.org/10.1101/246967 ）。基于比对的流程，比对工具也有很多选择，如Hisat，STAR,Tophat(hista可以替代tophat),bowtie等, 还有据说速度超快的Subread。根据比对工具的选择，定量软件也可以配套选择，如
STAR/bowtie—RSEM；
Hisat —featureCounts/HTseq；
Subread—featureCounts。
由于后续要鉴定新的转录本，因此需要对转录本组装，组装可以选择cufflinks,或者stringtie。文章是用的Tophat+cufflinks。我的后续分析打算用Hisat2比对，stringtie组装，featureCounts定量，同时打算尝试下Subread+featureCounts的流程。
Hisat2+stringtie+featureCounts；
Subread+featureCounts 。

sra格式的数据需要先用fastq-dump转换, --split-3 表示双端测序，--gzip将生成的fastq文件压缩。

质控检测采用fastqc和multiqc联合使用， multiqc有利于多个样本同时展示质量检测结果。FastQC的检测报告左侧会给出测序质量的一个summary,通常并不是所有的检测项都会是绿色通过，会有腔蔽一些警告和fail的项目。主要可以从Per base sequence quality 看一下测序碱基质量，Per sequence GC content 看一下GC含量，如果实际的GC含量（红线）出现双峰，且导致后期的序列比对很低时，需要引起注意，有可能是存在样品污染；再者就是看一下Adapter是否去除及去除的是否干净。这里的Adapter虽然显示通过，但是去除的并不干净，所以后续还会进一步去除Adapter。

Trimmomatic 支持多线程，处理数据速度快，主要用来去除 Illumina 平台的 fastq 序列中的接头，并根据碱基质量值对 fastq 进行修剪。软件有两种过滤模式，分别对应 SE 和 PE 测序数据，同时支持 gzip 和 bzip2 压缩文件。另外也支持 phred-33 和 phred-64 格式互相转化。

Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关，通常的过滤步骤如下：
ILLUMINACLIP : 过滤 reads 中的 Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2, 该引物序列陵圆轿可以在Trimmomatic软件的安装目录下找到，双端通常选择TruSeq3-PE-2。
SLIDINGWINDOW: 从 reads 的 5’ 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗。
MAXINFO : 一个自动调整的过滤选项，在保证 reads 长度的情况下尽量降低测序错误率，最大化 reads 的使用价值。
LEADING : 从 reads 的开头切除质量值低于阈值的碱基。
TRAILING : 从 reads 的末尾开始切除质量值低于阈值的碱基。
CROP : 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
HEADCROP : 从 reads 的开头切掉指定数量的碱基。
MINLEN : 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
AVGQUAL : 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
TOPHRED33 : 将 reads 的碱基质量值体系转为 phred-33。
TOPHRED64 : 将 reads 的碱基质量值体系转为 phred-64。
参考 ： NGS 数据过滤之 Trimmomatic 详细说明

Hisat2 比对，首先用hisat2-build 构建index，然后比对，输出sam格式，再用samtools将sam转为bam格式，并排序构建索引。
featureCounts 目前已经整合在subread，subread是一个综合的软件，还包括比对的工具和对应的R包， featrueCounts的定量效果和HTSeq差不多，但是featureCounts的优点非常快！

二进制版本的软件解压即可使用

（PS：存储不够了，后续选择两组数据做流程分析。）

SRR4042230和SRR4042231的比对率分别为88.02%和88.81%，总比对时间将近5小时。

hisat_counts.txt.summary 包含一些基本的统计信息。

可以看出subjunc比对不到2小时，比hisat2快3个多小时。

根据第1列是Geneid，第7,8列是counts数，用awk提取出geneID和counts。

一个植物转录组项目的实战
史上最快的转录组流程-subread
转录组定量-FEATURECOUNTS

HISAT2 建立索引警告和比对时报错解决方案

tags： HISAT2 RNA-seq

HISAT2 发表的文章中强调了它的速度很快，我就测试了一下这个工具。

HISAT2 建立索引：

然而没多久就看到这样的警告：

只是警告，并没有报错。

HISAT2 建参考索引很慢，等 HISAT2 建完索引（建索引花了 16 个小时），然后用 HISAT2 比对 RNA-seq 数据测试。

几分钟之后报错了：

github 上有人在 HISAT2 项目中报告过这个错误，虽然没有最终讨论出解决办法，但是都觉得跟建索引不完整有关，或许与建索引时候的警告有关。我查了一些资料，综合 biostar 和 SEQanswer 中的讨论，建立索引时遇到的警告是由于参考序列中存在大段的 n 碱基导致的，例如其中一条 fasta 中 n （我遇到的是小写的 n）太多。解决办法也很简单，过滤掉参考序列中长度小于 50bp 的 contig 和序列中连续 n 碱基超过 40bp 的contig 。然后重新建索引，就没有任何警告了，但是会损失部分参考序列。这样处理之后建的索引再睁唤庆用 HISAT2 比对 RNA-seq 数据，就没有问题了。

HISAT2 比对速度确实很快，一个样本转录组数据比对 2G 的核糖体 RNA 参考基因组约 25 分钟，bowtie2 需要 170 分钟（线程数都是 4）。 bowtie2 的 local 模式比对出 21% 的 rRNA 污染链宴，而 HISAT2 比对 2% 的 rRNA 污染，差异也挺大的……

但是，HISAT2 的线程数不能设置太高，用户手册中建议对人类参考基因组进行比对时，线程数设置在 1-8 之间。我用文献悉握 Transcript-level expression analysis of RNA-seq
experiments with HISAT, StringTie and Ballgown 中的数据测试也发现当线程数超过 12 时整个比对步骤消耗的时间反而会增加。