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 <num_mismatch> <probe_file> <trans_file> <output_file> [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
网友评论