linux序列比对程序,序列比对软件简单使用教程

linux可以使用的序列比对的工具有三个。blast、blat、seqmap。这三个软件都需要把待blast的序列做成fa格式

构建fa格式的序列

如果有个待比对的序列是含有两列，其中包括第一列(ID),第二列(sequence)。如果需要形成fa格式的话，可以使用下面的linux代码

awk '{print">"$1"\n"$2}' file

blast

linux 的blast软件分为基本上分为两个个步骤：

构建参考数据库

###下载软件

conda install blast

##下载genecode的参考基因组的fa

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.transcripts.fa.gz

##解压文件

gunzip gencode.v29.transcripts.fa.gz

##构建基因组的离线数据库

makeblastdb -in gencode.v29.transcripts -dbtype nucl -out humanGenome

构建离线数据库的参数中dbtype含有两种参数：nucl，prot分别代表核苷酸和蛋白

构建完成的数据库包括三个以out参数为开头的文件。比如示例的三个文件分别为：humanGenome.nhr humanGenome.nin humanGenome.nsq

选择blast的工具(blastn/blastp)对序列进行blast

blast可以分为很多的工具，

具体工具的选择看下表

img

blast数据库参数详解

blast软件详细的参数信息可以参见，官网上的描述。

-db 格式化了的数据库路径及数据库名

-query: 检索文件

-query_loc : 指定检索的位置

-strand： 搜索正义链还是反义链，还是都要

out : 输出文件

-remote： 可以用NCBI的远程数据库， 一般与 -db nr

-evalue 科学计数法，比如说1e3，定义期望值阈值。E值表明在随机的情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性

-outfmt: 输出的格式。有18个选项。其中6,7,8为自定义选型。6为正常的blast m8格式。

-num_descriptions：tabular格式输出结果的条数

-num_threads：线程数

-task:比对的时候的选项。有四个选项。1.)megablast，用于非常相似的序列(例如，测序错误)，2. dc-megablast，通常用于种间比较，3. blastn，用于种间的传统程序 比较，4. blastn-short，针对小于30个核苷酸的序列进行了优化。

blastn -query seq.fasta -out seq.blast -db dbname -outfmt 6 -task blastn-short -evalue 1e-5 -num_descriptions 10 -num_threads 8

blat软件使用

blat是UCSC用来比对序列序列的方式。网页版也是可以使用的。这里介绍linux版的。

###安装软件

conda install blat

blat软件参数详解

软件的基本格式：blat database query [-ooc=11.ooc] output.psl

软件的具体参数可以参见官方网站。这里介绍一下常见参数

-t=type: 参考数据库的数据类型。接受三个选项。1.dna(默认选项) ;2.prot；3.dnax(DNA sequence translated in six frames to protein)

-q=type：想要blat的数据类型。接受五个选项。1.dna - DNA sequence；2.rna - RNA sequence；3. prot - protein sequence；4.dnax - DNA sequence translated in six frames to protein；5. rnax - DNA sequence translated in three frames to protein

-out=type: 输出的格式。接受9中参数。1.psl - (Default) tab-separated format, no sequence；2. pslx - tab separated format with sequence；3.axt - blastz-associated axt format；4.maf - multiz-associated maf format；5. sim4 - similar to sim4 format；6. wublast - similar to wublast format；7.blast - similar to NCBI blast format；8. blast8 - NCBI blast tabular format；9.blast9 - NCBI blast tabular format with comments

blat常规设置

表达序列标签(EST)是cDNA序列的短子序列。

Mapping expressed sequence tag (EST) to the genome within the same species: -ooc=11.ooc

Mapping full length mRNAs to the genome in the same species: -ooc=11.ooc -fine -q=rna

Mapping ESTs to the genome across species: -q=dnax -t=dnax

Mapping mRNA to the genome across species: -q=rnax -t=dnax

Mapping proteins to the genome: -q=prot -t=dnax

Mapping DNA to DNA in the same species: -ooc=11.ooc -fastMap

Mapping DNA from one species to another species: -q=dnax -t=dnax

##比对芯片序列到基因组上且输出为blast格式

blat GCF.fa test_R1.fasta -out=blast8 -ooc=11.ooc

seqmap

seqmap是用于短序列比对特别快的工具。但是它出来的结果没有blast和blat多。如果要对芯片的序列进行重注释。是很好的一个工具

软件的安装

conda install seqmap

seqmap常规参数

软件的基本格式为：seqmap [options]

1.输入格式中参考基因组和比对的基因组必须是fa格式

2.num_mismatch代表比对的时候不匹配的个数

3.输出文件的格式分为两种。其中默认的是：Eland格式。另外一种是我们可以看得比较清楚的。用来显示所有匹配结果的格式：/output_all_matches

seqmap 0 GPL.fasta gencode.v29.transcripts.fa seqmap_gene.tmp /output_all_matches

在使用seqmap的时候。这个顺序不能错

上述的显示结果为

trans_id trans_coord target_seq probe_id probe_seq num_mismatch

1 313902 AACTCCGGGAGGGCCGCTTTGTATG 509644 AACTCCGGGAGTGCCGCTTTGTAGG 2

1 423680 TTTCACAATCAATGGATCAGGCCGC 129326 TTTCACAATCATTGGATCAGGCCAC 2

1 537816 CTTGAATTCAGTAAATAGTTTAACG 330515 CTTGAATTTAGTAAATAGTTTACCG 2

2 297292 CGTCAAATTTCGTCCTTTTCGCTGT 636826 CGTCAATTTTCGTCCTTTTCGGTGT 2

2 326279 CGTAGGACCATTCAGGCCGTTAAGC 986424 CGTAGGAGCATTCAGGCCGTTATGC 2

2 870729 GTTAACCTGTGGTAAGTAACGTAGT 433048 GTTAACCTGGGGTAAGTAACGTATT 2

3 204747 TAGCTCATTAACAGGGGATCTTAGG 917614 TAGCTCATTAATAGCGGATCTTAGG 2

3 601827 GTCGTTTTATTCCGCCTGGAGAGGT 321632 GTCGTCTGATTCCGCCTGGAGAGGT 2

3 674797 TCGCACTTGGGGCTAAATGGGCATC 336321 TCGCACTTCGGGCTAAATGGGAATC 2

3 927627 CAGCCAAAGATACGCAGCTCAGTCT 619563 GAGGCAAAGATACGCAGCTCAGTCT 2

4 305440 GACGGAAATCCATATAAGGTAGGGA 80583 GACGGAAATCGAGATAAGGTAGGGA 2