专栏首页生物信息云转录组分析 | 使用Hisat2进行序列比对

转录组分析 | 使用Hisat2进行序列比对

我前面已经对数据进行了质控:

转录组分析 | fastqc进行质控与结果解读

转录组分析 | 使用trim-galore去除低质量的reads和adaptor

接下来我们进行序列比对,利用的软件是Hisat2。

一.index文件下载

index文件直接去官网下载

(http://daehwankimlab.github.io/hisat2/download/),我测序的组织来自小鼠,所以我这里下载的是小鼠的。

下载完后上传到Linux服务器,解压后就可以直接用了。

我上传到了: /data/mouse_genome/ ,就是mm10_genome.tar.gz这个文件。

解压文件,解压过程中会在当前文件夹下创建mm10文件,解压后的文件就在mm10文件夹中。

 tar -zxvf /data/mouse_genome/mm10_genome_tar.gz

查看解压后的文件。

ll -h ./mm10

总共8个ht2格式文件,一个sh格式文件。

二.hisat2介绍

Hisat是一种高效的RNA-seq实验比对工具。它使用了基于BWT和Ferragina-manzini (Fm) index 两种算法的索引框架。使用了两类索引去比对,一类是全基因组范围的FM索引来锚定每一个比对,另一类是大量的局部索引对这些比对做快速的扩展。比对原理可阅读文献原文:HISAT: a fast spliced aligner with low memory requirements.

(一)背景

自2008年起,RNAseq已经成为研究基因表达、转录本结构、长链非编码RNA确定以及融合转录本的重要手段。随着测序深度的加深和read读长的延长,给比对工作带来很多困难。目前的工具如Tophat2和GSNAP等在对单个转录组实验的比对中耗时几天,而新型的STAR虽然很快,但是会吃掉大量的内存空间,如基于人类基因组的比对需要消耗28G左右,而hisat2是4.3G,小鼠3.8G左右。

(二)Hisat设计原则

1.优化了索引建立策略

hisat应用了基于bowtie2的方法去处理很多低水平的用于构建和查询FM索引的操作。但是与其它比对器不同的是,该软件应用了两类不同的索引类型:代表全基因组的全局FM索引和大量的局部小索引,每个索引代表64000bp。以人类基因组为例,创建了48000个局部索引,每一个覆盖1024bp,最终可以覆盖这个3 billion 的碱基的基因组。这种存在交叉(overlap)的边界可以轻松的比对那些跨区域的read(可变剪切体)。尽管有很多索引,但是hisat会把他们使用合适方法压缩,最终只占4gb左右的内存。

2.采用了新的比对策略

RNA-seq产生的reads可能跨长度比较大的内含子,哺乳动物中甚至最长能达到1MB,同时外显子比较短,read也比较短,会有很多read(模拟数据中大概34%)跨两个外显子的情况。为了更好的比对,将跨外显子的reads分成了三类:1)长锚定read,至少有16bp在两个外显子的每一个上 2)中间锚定read,有8-15bp在一个外显子上 3)短锚定read,只有1-7bp在一个外显子上。所以总的reads可以被划分为五类:1)不跨外显子的read 2)长锚定read 3)中间锚定read 4)短锚定read 5)跨两个外显子以上的read。在模拟的数据中,有25%左右的read是长锚定read,这种read在大多数情况下可以被唯一的定位到人的基因组上。5%为中间锚定read,对于这类,很多依赖于全局索引的算法就很难执行下去(需要比对很多次),而hisat,可以先将read中的长片段实现唯一比对,之后再使用局部索引对剩下的小片段进行比对(局部索引可以实现快速检索)。4.2%为短锚定read,因为这些序列特别短,因此只能通过在hisat比对其它read时发现的剪切位点或者用户自己提供的剪切位点来辅助比对。最后还有3%的是跨多个外显子的read,比对策略在hisat的online method中有介绍,文章中没有详解。比对过程中,中间锚定read、短锚定read、跨多个外显子read的比对占总比对时长的30%-60%,而且比对错误率很高!【引用:1】


官方手册:https://daehwankimlab.github.io/hisat2/manual/

(三).HiSat2进行比对的参数设置 【引用:2】

HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, www.ccb.jhu.edu/people/infphilo)
Usage: 
  hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

  <ht2-idx>  Index filename prefix (minus trailing .X.ht2).
             索引文件的前缀
  <m1>       Files with #1 mates, paired with files in <m2>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
             Paired-end测序的第一个文件,可以是压缩格式的(`.gz`或者`.bz2`)
  <m2>       Files with #2 mates, paired with files in <m1>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
             Paired-end测序的第二个文件
  <r>        Files with unpaired reads.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
             single-end测序的文件
  <sam>      File for SAM output (default: stdout)
             输出比对的结果(`.sam`)

  <m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
  specified many times.  E.g. '-U file1.fq,file2.fq -U file3.fq'.
  可以是多个序列文件用逗号分隔

Options (defaults in parentheses):

 Input:
  -q                 query input files are FASTQ .fq/.fastq (default)
                     输入检索序列为FASTQG格式
  --qseq             query input files are in Illumina's qseq format
                     输入检索序列是Illumina qseq格式
  -f                 query input files are (multi-)FASTA .fa/.mfa
                     输入文件是多个FASTA文件
  -r                 query input files are raw one-sequence-per-line
                     输入检索序列是每行一个原始序列
  -c                 <m1>, <m2>, <r> are sequences themselves, not files
                     直接输入序列,而不是文件
  -s/--skip <int>    skip the first <int> reads/pairs in the input (none)
                     忽略输入序列的前`<int>`个
  -u/--upto <int>    stop after first <int> reads/pairs (no limit)
                     只分析输入序列的前`<int>`个
  -5/--trim5 <int>   trim <int> bases from 5'/left end of reads (0)
                     剪去reads长度为`<int>`的5'端序列
  -3/--trim3 <int>   trim <int> bases from 3'/right end of reads (0)
                     剪去reads长度为`<int>`的3'端序列
  --phred33          qualities are Phred+33 (default)
                     序列质量打分为Phred+33
  --phred64          qualities are Phred+64
                     序列质量打分为Phred+64
  --int-quals        qualities encoded as space-delimited integers
                     序列质量打分为用空格分隔的整数

 Alignment:
  --n-ceil <func>    func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)
                     允许在比对结果中出现的最多非ACTG的计算函数
  --ignore-quals     treat all quality values as 30 on Phred scale (off)
                     将所有质量分值设为30
  --nofw             do not align forward (original) version of read (off)
                     不比对前向的read
  --norc             do not align reverse-complement version of read (off)
                     不比对后向的read

 Spliced Alignment:
  --pen-cansplice <int>              penalty for a canonical splice site (0)
                                     对已知剪接位点的罚分
  --pen-noncansplice <int>           penalty for a non-canonical splice site (12)
                                     对新的剪接位点的罚分
  --pen-canintronlen <func>          penalty for long introns (G,-8,1) with canonical splice sites
                                     对已知剪接位点中长内含子的罚分
  --pen-noncanintronlen <func>       penalty for long introns (G,-8,1) with noncanonical splice sites
                                     对未知剪接位点中长内含子的罚分
  --min-intronlen <int>              minimum intron length (20)
                                     最小允许的内含子长度
  --max-intronlen <int>              maximum intron length (500000)
                                     最大允许的内含子长度
  --known-splicesite-infile <path>   provide a list of known splice sites
                                     提供一个已知的剪接位点列表
  --novel-splicesite-outfile <path>  report a list of splice sites
                                     报告新的剪接位点
  --novel-splicesite-infile <path>   provide a list of novel splice sites
                                     提供新的剪接位点
  --no-temp-splicesite               disable the use of splice sites found
                                     不使用发现的新剪接位点
  --no-spliced-alignment             disable spliced alignment
                                     不允许剪接比对
  --rna-strandness <string>          specify strand-specific information (unstranded)
                                     指定链特异性的信息
  --tmo                              reports only those alignments within known transcriptome
                                     报告只包含已知转录组信息的比对结果
  --dta                              reports alignments tailored for transcript assemblers
                                     结果适合用于stringtie进行转录本的拼接
  --dta-cufflinks                    reports alignments tailored specifically for cufflinks
                                     结果适合于cufflinks的应用
  --avoid-pseudogene                 tries to avoid aligning reads to pseudogenes (experimental option)
                                     避免与假基因的比对
  --no-templatelen-adjustment        disables template length adjustment for RNA-seq reads
                                     不允许对reads的模板长度进行调整

 Scoring:
  --mp <int>,<int>   max and min penalties for mismatch; lower qual = lower penalty <6,2>
                     不匹配的最大与最小罚分值,低质量对应低罚分
  --sp <int>,<int>   max and min penalties for soft-clipping; lower qual = lower penalty <2,1>
                     soft-clipping的最大与最小罚分
  --no-softclip      no soft-clipping
                     不考虑Soft-clipping
  --np <int>         penalty for non-A/C/G/Ts in read/ref (1)
                     对read或参考序列中出现的非ATCG的罚分
  --rdg <int>,<int>  read gap open, extend penalties (5,3)
                     对read的gap-open, gap-extension的罚分
  --rfg <int>,<int>  reference gap open, extend penalties (5,3)
                     对参考序列的gap-open,gap-extension的罚分
  --score-min <func> min acceptable alignment score w/r/t read length
                     (L,0.0,-0.2)
                     相对read长度的最小接受比对分值

 Reporting:
  -k <int> (default: 5) report up to <int> alns per read
                     对每个read只报告`<int>`个比对结果

 Paired-end:
  -I/--minins <int>  minimum fragment length (0), only valid with --no-spliced-alignment
                     paired-end的最小片段长度
  -X/--maxins <int>  maximum fragment length (500), only valid with --no-spliced-alignment
                     最大片段长度
  --fr/--rf/--ff     -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)
                     `--fr==fw/rev`,`--rf==rev/fw`, `--ff==fw/fw`
  --no-mixed         suppress unpaired alignments for paired reads
                     不输出paired reads的不能匹配的比对结果
  --no-discordant    suppress discordant alignments for paired reads
                     不输出paired-reads的不统一的比对结果

 Output:
  -t/--time          print wall-clock time taken by search phases
                     输出搜索阶段所花费的wall-clock时间
  --un <path>           write unpaired reads that didn't align to <path>
                        输出没有与`<path>`比对的unpaired reads
  --al <path>           write unpaired reads that aligned at least once to <path>
                        输出至少与`<path>`有一次比对的unpaired reads
  --un-conc <path>      write pairs that didn't align concordantly to <path>
  --al-conc <path>      write pairs that aligned concordantly at least once to <path>
  (Note: for --un, --al, --un-conc, or --al-conc, add '-gz' to the option name, e.g.
  --un-gz <path>, to gzip compress output, or add '-bz2' to bzip2 compress output.)
                    将输出结果进行`gz`或者`bz2`压缩
  --summary-file     print alignment summary to this file.
                     将比对的总结输出到文件中
  --new-summary      print alignment summary in a new style, which is more machine-friendly.
                     将比对总结以新的形式输出
  --quiet            print nothing to stderr except serious errors
                     除了严重错误以外,不将任何结果输出到stderr
  --met-file <path>  send metrics to file at <path> (off)
                     将度量值i输出到文件
  --met-stderr       send metrics to stderr (off)
                     将度量值输出到stderr
  --met <int>        report internal counters & metrics every <int> secs (1)
                     每隔`<int>`秒输出当前的内部计数器和度量值
  --no-head          supppress header lines, i.e. lines starting with @
                     不输出header行,也就是所有以`@`开始的行
  --no-sq            supppress @SQ header lines
                     不输出`@SQ`行
  --rg-id <text>     set read group id, reflected in @RG line and RG:Z: opt field
                     设置`@RG`
  --rg <text>        add <text> ("lab:value") to @RG line of SAM header.
                     Note: @RG line only printed when --rg-id is set.
                     将`<text>`加到SAM文件的`@RG`行
  --omit-sec-seq     put '*' in SEQ and QUAL fields for secondary alignments.
                     将"*"置于SEQ和QUAL两列中,进行二次比对

 Performance:
  -o/--offrate <int> override offrate of index; must be >= index's offrate
                     重置索引的`offrate`,该值必须大于原有的`offrate`
  -p/--threads <int> number of alignment threads to launch (1)
                     比对所用的线程数
  --reorder          force SAM output order to match order of input reads
                     强制要求SAM输出顺序比对与输入reads一致
  --mm               use memory-mapped I/O for index; many 'hisat2's can share
                     使用内存映射的I/O作为索引

 Other:
  --qc-filter        filter out reads that are bad according to QSEQ filter
                     根据QSEQ filter过滤输入的reads
  --seed <int>       seed for random number generator (0)
                     设置RNG seeds的值
  --non-deterministic seed rand. gen. arbitrarily instead of using read attributes
                     不根据reads的属性产生随机数
  --remove-chrname   remove 'chr' from reference names in alignment
                     比对结果中的参考序列名称中去除`chr`
  --add-chrname      add 'chr' to reference names in alignment
                     在比对结果的参考序列名称中加入`chr` 
  --version          print version information and quit
                     输出版本信息
  -h/--help          print this usage message
                     输出帮助信息

hisat2可以自己建立index文件,如果没有现成的index的话,那就得自己去建立了,这个就比较麻烦而且耗内存和时间,下面是建立index的一些参数介绍。因为我们的index文件是直接下载的,所以我后面不需要建立index文件,关于这部分,后面有机会介绍。

HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, http://www.ccb.jhu.edu/people/infphilo)
Usage: hisat2-build [options]* <reference_in> <ht2_index_base>
    reference_in            comma-separated list of files with ref sequences
                            以逗号分隔的参考序列文件的列表
    hisat2_index_base       write ht2 data to files with this dir/basename
                            输出的索引文件的`basename`
Options:
    -c                      reference sequences given on cmd line (as
                            <reference_in>)
                            输入的参考序列
    --large-index           force generated index to be 'large', even if ref
                            has fewer than 4 billion nucleotides
                            强制要求产生的索引为`large`
    -a/--noauto             disable automatic -p/--bmax/--dcv memory-fitting
                            关闭-p/--bmax/--dcv等自动内存匹配
    -p                      number of threads
                            线程数
    --bmax <int>            max bucket sz for blockwise suffix-array builder
                            逐块后缀数组建立采用的最大bucket size
    --bmaxdivn <int>        max bucket sz as divisor of ref len (default: 4)
    --dcv <int>             diff-cover period for blockwise (default: 1024)
    --nodc                  disable diff-cover (algorithm becomes quadratic)
    -r/--noref              don't build .3/.4.ht2 (packed reference) portion
    -3/--justref            just build .3/.4.ht2 (packed reference) portion
    -o/--offrate <int>      SA is sampled every 2^offRate BWT chars (default: 5)
    -t/--ftabchars <int>    # of chars consumed in initial lookup (default: 10)
    --localoffrate <int>    SA (local) is sampled every 2^offRate BWT chars (default: 3)
    --localftabchars <int>  # of chars consumed in initial lookup in a local index (default: 6)
    --snp <path>            SNP file name
    --haplotype <path>      haplotype file name
    --ss <path>             Splice site file name
    --exon <path>           Exon file name
    --seed <int>            seed for random number generator
    -q/--quiet              verbose output (for debugging)
    -h/--help               print detailed description of tool and its options
    --usage                 print this usage message
    --version               print version information and quit

三. 使用Hisat2进行序列比对

创建输出数据的文件夹

mkdir cleandata/hisat2_mm10data

因为比对这一过程很耗内存,所以样本多话,计算机内存不够大,需要分批比对,我就先介绍一个样本比对的命令:

hisat2 --dta -t -x /data/RNAseq/mm10/genome -1 /data/RNAseq/cleandata/trim_galoredata/CK-4_1_val_1.fq.gz -2 /data/RNAseq/cleandata/trim_galoredata/CK-4_2_val_2.fq.gz -S cleandata/hisat2_mm10data/CK4.sam 

-1和-2分别表示双端测序的1个文件,后面跟的是文件路径,一定要注意 /data/RNAseq/mm10/genome文件的目录,genome这个不是文件夹,是index文件的前缀,我的mm10文件下并没有这个文件,如果不加genome就会发生如下报错:

## Could not locate a HISAT2 index corresponding to basename "/data/RNAseq/mm10/genome/" ##

另外一定要加上 --dta,后续用Stringtie处理会更容易一些,这是StringTie使用说明里面的一句话:【NOTE: be sure to run HISAT2 with the --dta

option for alignment, or your results will suffer.】后续再介绍。-S后面cleandata/hisat2_mm10data/CK4.sam 是输出文件路径,hisat2_mm10data这个文件夹如果没有建立,需要事先建立。

命令没有问题的话,出现下面提示表示程序正在执行,就去干其他的,等执行结束再往下:

Time loading forward index: 00:02:51
Time loading reference: 00:00:16

我这里,一个样本43分钟。

计算机内存足够大的话,我们可以像前文【转录组分析 | 使用trim-galore去除低质量的reads和adaptor】一样,通过一个脚本一次执行。

创建脚本文件:hisat2_mm10_batch.sh

vim hisat2_mm10_batch.sh

输入下面内容:

#!/bin/bash
# This is for hisat2 batch aligne
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome \
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz \
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz \
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam; done

我这里由于CK-4这个样本前面单独处理了,批量处理就不需要处理。

保存脚本,然后执行脚本。

bash hisat2_mm10_batch.sh

脚本运行时间和计算机配制还有数据量的大小有关。运行结束后我们查看结果:

ll -h ./cleandata/hisat2_mm10data/

原本22G的数据,运行结束后155G,所以自己用本地计算机跑的话,数据多注意一下内存。

四.结果解读

我们比对后得到的是sam文件,关于sam文件是一个什么样的文件,参考文章:生信中常见的数据文件格式


为了快速查看本公众号文章,可阅读文章:公众号文章目录


致谢:本文还参考了以下内容:

【1】https://www.jianshu.com/p/ce3f4afb9b60

【2】https://github.com/ricket-sjtu/bioinformatics/wiki/HiSat2%E5%BA%94%E7%94%A8%E4%BB%8B%E7%BB%8D

对原作者表示感谢!

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-09-17

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 转录组分析 | 使用RSeQC软件对生成的BAM文件进行质控

    RSeQC是发表于2012年的一个RNA-Seq质控工具,属于python包。它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据,比如一些基本...

    DoubleHelix
  • 赶快收藏:R语言常用函数整理

    DoubleHelix
  • R语言数据分析与挖掘(第四章):回归分析(4)——logistic回归

    前面我们介绍的回归方法,一般适用于数值型数据对象,对于分类数据类型就不再适用。对于分类数据对象,我们需要引入广义线性回归方法,比如logistic回归和pois...

    DoubleHelix
  • mongo数据库命令简单学习

    db.getCollection('product').update({status:"offline"},{$set:{status:"online"}},f...

    蓓蕾心晴
  • 光城归来之C语言开发网站

    最近要把防火墙项目做个页面,而底层全部c语言实现,那么就得做个web页面,想了一下,C大法这么厉害,也应该可以的,然后大家就见到了这篇文章。

    公众号guangcity
  • 你不知道的Golang盲点汇总【持续更新】

    根据Go规范, map上的迭代顺序是不确定的,并且在程序运行之间可能会有所不同。实际上,不仅它是不确定的,而且实际上是有意随机化的。这是因为它曾经是可预测的,并...

    sunsky
  • cpuimage 开源之

    前年学习opengl做的一个小东西。 原本计划将gpuimage 的算法一个一个转写成cpu版本 c,c++ 版本。 gpuimage 项目参考: https:...

    cpuimage
  • 停车服务公众号存安全隐患,可能泄漏你的用车行踪

    近些年提供停车服务的微信公众号越来越多,虽然带来很大方便,但也可能存在潜在风险。近期有网友称,某停车服务公众号可随意绑定车牌,据此可以查询到相应车辆的停车记录。...

    FB客服
  • PHP使用了PDO还可能存在sql注入的情况

    “用 PDO 来防止 SQL 注入。”大概学过 PHP 的都听说过这句话。代码中出现了 PDO 就行了吗?答案肯定是否定的。接下来给大家介绍几种使用了 PDO ...

    信安之路
  • MySQL主从复制

    修改mysql配置文件,不同的系统my.cnf路径不同(CentOS中位于/etc/my.cnf)

    十毛

扫码关注云+社区

领取腾讯云代金券