序列比对-BLAST

一.BWA

BWA主要是将reads比对到大型基因组上，主要功能是：序列比对。首先通过BWT(Burrows-Wheeler Transformation，BWT压缩算法)为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。由三种类似算法组成：

BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。

三种算法的使用范围

BWA-backtrack：reads长度<70bp时，推荐本算法，建议输入reads长度 < 100bp。

BWA-SW：在reads具有频繁的gap时，比对更敏感，推荐本算法。reads长度一般为70bp-1Mbp，支持long-reads，split alignment。

BWA-MEM(推荐)：在reads长度在70bp-1Mbp范围时，推荐本算法(除了上面两种情况)。支持long-reads，split alignment。

bwa主要功能与参数

1.建立索引Index

ndex Usage：
      bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>
OPTIONS: 
      -p STR   输出数据库的前缀；【默认和输入的文件名一致，输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀。】
      -a [is|bwtsw]   构建index的算法，有两个算法： is 是默认的算法，虽然相对较快，但是需要较大的内存，当构建的数据库大于
               2GB的时候就不能正常工作了。 bwtsw 对于短的参考序列式不工作的，必须要大于等于10MB, 但能用于较大的基因组数据，比如人的全基因组。

#根据reference genome data(e.g. ref.fa) 建立 Index File例子：
$ bwa index ref.fa -p genome###可以不加-p genome，这样建立索引都是以ref.fa为前缀

2.mem比对

mem Usage: bwa mem [options] ref.fa reads.fq [mates.fq]

常用参数如下：

-t   INT 线程数，默认是1。
-M   将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
-p   若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R   STR 完整的read group的头部，可以用 '\t' 作为分隔符， 在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部。
-T   INT   当比对的分值比 INT 小时，不输出该比对结果，这个参数只影响输出的结果，不影响比对的过程。-a 将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。
#例子：
$ bwa mem ref.fa reads.fq > mem-se.sam
$ bwa mem ref.fa read1.fq read2.fq > mem-pe.sam
$ bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" ref.fa read1.fastq read2.fastq > mem-pe.sam 2> ./mem-pe.log

二.HISAT

简介

HISAT2是TopHat2/Bowti2的继任者，使用改进的BWT算法，实现了更快的速度和更少的资源占用，作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。 官网

建立索引

建立基因组索引

hisat2-build –p 4 genome.fa genome

建立基因组+转录组+SNP索引： bowtie2的索引只有基因组序列信息，tophat2比对时，转录组信息通过-G参数指定。HISAT2建立索引时，就应该把转录组信息加进去。 HISAT2提供两个Python脚本将GTF文件转换成hisat2-build能使用的文件：

extract_exons.py Homo_sapiens.GRCh38.83.chr.gtf > genome.exon
extract_splice_sites.py Homo_sapiens.GRCh38.83.chr.gtf > genome.ss

此外，HISAT2还支持将SNP信息加入到索引中，这样比对的时候就可以考虑SNP的情况。这仍然需要将SNP文件转换成hisat2-build能使用的文件：

extract_snps.py snp142Common.txt > genome.snp

最后，将基因组、转录组、SNP建立索引：

hisat2-build -p 4 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran

官网提供了人和小鼠的索引文件下载，压缩包有make_grch38_tran.sh文件，详细记录了创建索引的过程。

比对与主要参数

运行HISAT2

hisat2 -p 16 -x ./grch38_tran/genome_tran -1 SRR534293_1.fastq -2 SRR534293_2.fastq –S SRR534293.sam

-x 指定基因组索引

-1 指定第一个fastq文件

-2 指定第二个fastq文件

-S 指定输出的SAM文件

更多参数请查看HISAT2的操作手册

官方操作手册简要版

用法：

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]

主要参数： -x 参考基因组索引文件的前缀。

-1 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。

-2 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。

-U 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。

–sra-acc 输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。

-S 指定输出的SAM文件。

输入选项： -q 输入文件为FASTQ格式。FASTQ格式为默认参数。

-qseq 输入文件为QSEQ格式。 -f 输入文件为FASTA格式。

-r 输入文件中，每一行代表一条序列，没有序列名和测序质量等。选择此项时，–ignore-quals参数也会被选择。

-c 此参数后是直接比对的序列，而不是包含序列的文件名。序列间用逗号隔开。选择此项时，–ignore-quals参数也会被选择。

-s/–skip 跳过输入文件中前条序列进行比对。

-u/–qupto 只使用输入文件中前条序列进行比对，默认是没有限制。

-5/–trim5 比对前去除每条序列5’端个碱基

-3/–trim3 比对前去除每条序列3’端个碱基

–phred33 输入的FASTQ文件碱基质量值编码标准为phred33，phred33为默认参数。

–phred64 输入的FASTQ文件碱基质量值编码标准为phred64。

–solexa-quals 将Solexa的碱基质量值编码标准转换为phred。

–int-quals 输入文件中的碱基质量值为用空格分隔的数值，而不是ASCII码，例如40 30 30 40。

 

三.BLAST

BLAST: basic local alignment search tool,是目前最常用的序列比对软件，它能够对生物不同蛋白质的氨基酸序列或不同基因的DNA序列进行比对，并从相应数据库中找到相同或相似的序列。

BLAST比对的要点是片段对的概念，它是指两个给定序列中的一对子序列，它们的长度相等，且可以形成无空格的完全匹配。

BLAST比对的基本过程是它首先找出查询序列和目标序列间所有匹配程度超过一定阈值的片段对，然后对片段对根据给定的相似性阈值进行延伸，得到一定长度的相似性片段，最后给出高分值片段对。

BLAST的在线使用

BLAST软件提供了在线的比对工具，它不仅可以进行核酸序列之间的比对，也可以进行蛋白质与蛋白质以及蛋白质与核酸之间的比对。

这里，我们以核酸序列比对为例。首先我们可以选取我们想要比对的序列，这里我从NCBI上选取了家鼠的Tlr4基因和Tnf基因的部分序列进行比对，并分别获得它们fasta格式：

之后将选取的DNA序列粘贴到BLAST的两个比对框里，勾选相应选项（注意！由于选取的两个序列相似性程度较低，所以在这里勾选了第三个选项）

 

最后点击BLAST即得到相应结果从结果中我们可以看出，比对得分仅有44.1，因此这两个序列相似性程度很低。