使用HISAT2对转录组数据进行比对 编辑

HISAT2(Hical Indexing for Spliced Alignment of Transcripts)是由约翰霍普金斯大学开发,能够将RNA-Seq的reads与基因组进行快速比对的一款程序。这项成果发表在2015年3月9日的《Nature Methods》上。HISAT利用大量FM索引,以覆盖整个基因组。以人类基因组为例,它需要48, 000个索引,每个索引代表~64, 000bp的基因组区域。这些小的索引结合几种比对策略,实现了RNA-Seq读取的高效比对,特别是那些跨越多个外显子的读取。尽管它利用大量索引,但HISAT只需要4.3GB的内存,而运行速度比Tophat快了一个数量级。

该软件比 tophat 比对速度快很多倍。

HISAT2相关链接 编辑

HISAT2的下载和安装 编辑

软件和基因组索引下载地址

# installing HISAT2
wget https://cloud.biohpc.swmed.edu/index.php/s/oTtGWbWjaxsQ2Ho/download -P ~/software/
unzip ~/software/hisat2-2.2.1-Linux_x86_64.zip -d /opt/biosoft/
echo 'PATH=$PATH:/opt/biosoft/hisat2-2.2.1-Linux_x86_64' >> ~/.bashrc
source ~/.bashrc
# extract_splice_sites.py程序的运行需要python2.7
tar zxf ~/software/Python-2.7.11.tgz
./configure --prefix=/opt/sysoft/Python-2.7.11
make -j 4
echo 'PATH=/opt/sysoft/Python-2.7.11/bin/:$PATH' >> ~/.bashrc
source ~/.bashrc
cd ..
rm Python-2.7.11 -rf

HISAT2的使用 编辑

Hisat2的使用方式和命令参数与Bowtie2非常类似。

首先,构建参考序列的索引:

hisat2-build genome.fasta genome

hisat2-build 常用参数:

-p <int> default:1
		使用多线程运行。当基因组较大的时候,构建索引数据库非常耗时间。该参数能加快速度。
--snp <path>
		输入一个包含SNP信息的文件。该文件有5列以tab分割的数据:SNPID、参考序列ID、SNP类型 (single, deletion或insertion)、SNP位点(以第一个碱基位点为0计算)、变异喊基信息。例如:“rs58784443 single 13 18447947 T"。
--haplotype <path>
		输入一个单倍型信息文件,表明--snp参数指定的某些变异位点在要分析的样品中是单倍型的,和其变异位点碱基信息一致。该文件有5列以tab分割的数据:Haplotype ID、参考序列ID、起始位点(以第一个碱基位点为0计算)、结束位点(以第一个碱基位点为0计算)、逗号分隔的多个SNP ID。 例如:ht35 13 18446877 18446945 rs12381094, rs12381056, rs192016659, rs538569910O 输入该参数能减小索引数据库的大小。	'
--ss <path>
		输入一个包含有剪接位点(Splicing Site)信息的文件。该文件可以利用HISAT2软件自带的hisat2_extract_splice_sites.py程序对编码蛋白基因结构注释GTF文件转换获得。
--exon <path>
		输入一个含有exon信息的文件。利用HISAT2软件自带的hisat2_extract_exons.py程序对编码蛋白基因结构注释GTF文件转换获得该文件。

然后,使用HISAT2对转录组数据进行比对。 常用示例:

hisat2 -x genome -u 1000000 -p 24 -I 0 -X 500 —fr --min-intronlen 20 --max-intronlen 4000 -1 reads.1.fastq -2 reads.2.fastq -U single.fastq -S result.sam

Hisat2的常用参数:

#一般参数
--version
打印程序版本并退出
-h/--help
打印用法信息并推出
#主要参数:
-x <hisat-idx>
索引数据文件的前缀。
-1 <m1>
双末端测序结果的第一个文件。若有多组数据,使用逗号将文件分隔。Reads的长度可以不一致。
-2 <m2>
双末端测序结果的第二个文件。若有多组数据,使用逗号将文件分隔,并且文件顺序要和-1参数的对应。Reads的长度可以不一致。
-U <r>
单端数据文件。若有多组数据,使用逗号将文件分隔。可以和-1,-2参数同时使用。Reads的长度可以不一致。
--sra-acc <SRA accession number>
输入SRA登录号,比如SRR353653,SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载数据并识别数据类型,进行比对。该参数的正常使用需要安装NCBI-NGS toolkit。
-S <hit>
设置输出文件名。默认情况下输出到标准输出。

#输入参数:
-q
输入的Reads是fastq格式。默认情况下开启了该参数。
-qseq
输入的文件为QSEQ格式文件。
--qc-filter
滤除QSEQ fileter filed为非0的reads. 仅当有—qseq选项时有效. Default: off. 
--seed <int>
使用<int>作为随机数产生的种子. Default: 0.
-f
输入的Reads是fasta格式。
-r
输入的文件中,每一行代表一条序列,没有序列名和测序质量等。选择此项时,表示--ignore-quals也被选择了。
-c
此参数后直接是比对的reads序列,而不是包含序列的文件名。序列间用逗号隔开。选择此项时,表示--ignore-quals也被选择了。
-s | --skip <int>
跳过前<int>个read pairs,再进行比对。
-u | --qupto <int>
仅对前<int>个read pairs进行比对。
-5 | --trim5 <int>
将reads 5‘端碱基去除<int>个后,再进行比对。
-3 | --trim3 <int>
将reads 3‘端碱基去除<int>个后,再进行比对。
--phred33
输入的fastq文件的碱基质量为phred33。默认开启了此参数。
--phred64
输入的fastq文件的碱基质量为phred64。
--solexa-quals
将Solexa的碱基质量转换为Phred。在老的GA Pipeline版本中得以运用。Default: off.
--int-quals
输入文件中的碱基质量为用“ ”分隔的数值,而不是ASCII码。比如40 40 30 40...。 Default: off.

#比对参数:
--n-ceil <func>
设定read中允许含有不确定碱基(非GTAC,通常为N)的最大数目。Default: L,0,0.15。计算公式为:f(x) = 0 + 0.15 * x,表示长度为100的read最多运行存在15个不确定碱基。一旦不确定碱基数超过15,则该条read会被过滤掉。
--ignore-quals
计算错配罚分的时候不考虑碱基质量.当输入序列的模式为-f, -r或者-c的时候,该设置自动成为默认设置.
--nofw/--norc
--nofw设置reads不和前导链(forward reference strand)进行比对;--norc设置reads和后随链(reverse-complement reference strand)进行比对。Default: both strands enabled。

#得分参数:
--ma <int>
设定匹配得分.--local模式下每个read上碱基和参考序列上碱基匹配,则加<int>分.在--end-to-end 模式中无效.Default: 2.
--mp MX,MN
设定错配罚分。其中MX为所罚最高分,MN为所罚最低分。默认设置下罚分与碱基质量相关。罚分遵循的公式为:MN + floor( (MX-MN)(MIN(Q, 40.0)/40.0))。其中Q为碱基的质量值。如果设置了--ignore-qual参数,则错配总是罚最高分。Default: MX = 6, MN = 2。
--sp <int>,<int> default: 1,2
设定soft-clipping罚分,罚分与碱基质量值有关,与参数--mp的罚分算法一致。Read首尾碱基若和参考基因组不匹配的时候,程序将read首尾序列进行软截断(即匹配信息结果中没有计算被截断的序列,同时又保留了整条read信息)。该参数的罚分较低,表示允许Read首尾碱基不匹配。Hisat2默认设置下该参数生效,其SAM文件结果中,第4列参考序列上的比对起始位点对应着截断后read的起始位点,第6列中的S表示发生了soft-clipping。
--no-softclip
加入该参数后,则禁止soft-clipping功能。有些表达量或intron分析软件不能正确处理有soft-clipping的比对结果,可以加入该参数禁用soft-clipping功能,从而提高SAM文件的兼容性。
--np <int>
当匹配位点中read或reference上有不确定碱基(比如N)时所设定的罚分值。Default: 1。
--rdg <int1>,<int2>
设置在read上打开gap罚分<int1>,延长gap罚分<int2>。Default: 5, 3。
--rfg <int1>,<int2>
设置在reference上打开gap罚分 <int1>,延长 gap 罚分<int2>。Default: 5, 3。
--score-min <func>
设定成为有效比对的最小分值。Default: L,0,-0.2。表示对于100bp的read,允许的最小得分是-20。

#剪接比对参数:
--pen-cansplice <int>
当剪接位点为常规位点(比如GT/AG)时的罚分。Default: 0。
--pen-noncansplice <int>
当剪接位点为非常规位点(比如non-GT/AG)时的罚分。Default: 3。一般情况下,有少部分剪接位点为GC/AG和AT/AC。
--min-intronlen <int>
设置最小的intron长度。Default: 20。
--max-intronlen <int>
设置最大的 intron 长度。Default: 500000。
--known-splicesite-infile <path>
提供一个包含剪接位点信息的文本文件。该文件可以使用extract_splice_sites.py命令对GTF文件分析得到。
--nove1-sp1icesite-outfile <path>
HISAT2报告新剪接位点的信息文件。
--novel-splicesite-infile <path>
用于输入上一个参数得到的文件。
--no-spliced-alignment
不进行剪接性比对。
--rna-strandness <string>
设置链特异性参数。若是单端测序数据,设置为F表示reads都比对到transcripts序列正链上, 设置为R表示reads都比对到transcripts负链上。若是双末端测序数据,则设置为RF或FR。其 中,RF表示双末端测序中第一个Fastq文件中的reads都比对到transcripts序列负链上,且第 二个Fastq文件的中reads都比对到transcripts正链上。FR表示双末端测序中第一个Fastq文件中的reads都比对到transcripts序列正链上,且第二个Fastq文件的中reads都比对到 transcripts负链上。一般链特异性测序中,该参数的值应该是RF。该参数和Tophat软件中的参 数--library-type 类似。
--tmo | --transcriptome-mapping-only
开启该参数后则仅报告比对到transcripts序列上的比对结果。
--dta | --downstream-transcriptome-assembly
开启该参数后,能报告和转录子组装软件StringTie兼容的比对结果。开启该参数后,对于重头预 测的剪接位点,HISAT2需要更长的锚定长度。
--dta-cufflinks
开启该参数后,能报告与cufflinks软件兼容的比对结果。开启该参数,相当于开启了--dta参数, 并且HISAT2寻找新的剪接位点时,则必须满足是(GT/AG, GC/AG和AT/AC),同时比对结果中多了 XS:A:[+-]的标签信息。若非链特异性测序结果不加此参数,则HISAT2的结果进行Cufflinks分 析的时候会报错。

#报告参数:
-k <int>
设置报告的最优比对结果数目。Hisat会寻找最优比对结果,当超过此数目时,则停止寻找。 Default: 5 (HFM) or 10 (HGFM)。当基因组重复序列较多,该参数设置较大,则运行速度慢。

#Paired-end 参数:
-I/--minins <int> Default: 0
设置插入片段长度的最小值。插入片段长度为:两端序列长度+ Gap长度。-I和-X参数设置 的范围越大,则运行速度越慢。若其值分别设置为2的和400,则Hisat运行比较高效。
-X/—maxins <int> Default: 500
设置插入片段长度的最大值。插入片段长度为:两端序列长度+ Gap长度。
--fr/--rf/--ff
	设置双末端测序数据的方向。Default: --fr。
	--fr: (==>> <<==)匹配时,mate1在5’端上游,和前导链一致,mate2在3’端下游,和前导链反向互补;或者mate2在上游,和前导链一致,mate1在下游反向互补。
	--rf: (<<== ==>>)匹配时,mate1在5’端上游,和前导链反向互补,mate2在3’端下游,和前导链一致。
	--ff: (=>> ==>>) matel  mate2 都和前导链一致。
--no-mixed
默认设置下,一对reads不能成对比对到参考序列上,则单独对每个read进行比对。该选项则阻止此行为。
--no-discordant
默认设置下,一对reads不能和谐比对(concordant alignment,即满足-I、-X和--fr/--rf/--ff的条件)到参考序列上,则搜寻其不和谐比对(disconcordant alignment,即两条reads都能独一无二地比对到参考序列上,但是不满足-I、-X和--fr/--rf/--ff的条件)。该选项阻止此行为。

#SAM参数:
--rg-id <text>
设定read group ID为text。在SAM文件的头中增加一行@RG,在输出的SAM文件中添加Tag ”RG:Z:text”。
--rg <text>
使用text作为@RG的一列,比如”SM:sample1、”PL:illumina”等。在@RG中加入多列,则多次使用该参数即可。

#性能参数:
-p/--threads NTHREADS
设置程序运行的线程数。

HISAT2结果统计 编辑

Bowtie2输出的比对结果文件为SAM格式,同时在比对完毕时,将比对的统计结果输 出到标准错误输出。更为详细的统计结果可以自行编写程序从SAM文件分析。

$ RNA_Seq_Sam_Statistics.pl hisat2.sam > hisat2.stats

hisat2.stats文件内容示例:
Total reads			 95202044		100.00%
Total Basepairs		 13822693238	100.00%
Total Mapped reads 	 93152492 		 97.85%
Fully match			 71067818		 74.65%
<=2bp mismatch		 19578773		 20.57%
Unique match		 73750114		 77.47%
Multi-position match 19402378		 20.38%
One Best Score match 88880927		 93.36%
Total unmatched reads 2049552		  2.15%

程序RNA_Seq_Sam_Statisties.pl根据HISAT2的SAM文件中的标签来识别reads的匹配情况:
NM:i:<N>
该标签表示不匹配的碱基数(read和参考序列差异的最小碱基个数)ZS:i:<N>
该标签考示read的次优比对得分。若是双末端测序,单端read的次优比对得分可能高于最优比对得分。若存在ZS标签,则表明read一定能比对到基因组多个位置。若不存在ZS标签,read也有可能比对到基因组上多个位置(这多个位置的比对得分相同)。
NH:i:<N>
该标签表示read或read pair比对到基因组上的位置数目。