跳至內容

生物信息學/hisat

新增連結
維基教科書，自由的教學讀本

使用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比对到基因组上的位置数目。
取自「https://zh.wikibooks.org/w/index.php?title=生物信息学/hisat&oldid=142064」