如何在Linux中进行本地Blast序列比对？

今天分享一篇学习笔记，主要包含本地化blast序列比对和数据提取方法。

首先，需要准备RNA数据和蛋白质数据，本次利用蛋白质数据建立索引库，然后将RNA比对到蛋白质序列。

RNA数据

创建一个目录，导入mRNA序列数据，通常是一个fasta后缀文件。

1. 在工作目录下创建alignment文件夹
2. 将mRNA序列数据文件wheat-test.fasta拷贝到/alignment中

蛋白质数据

新建一个文件夹存放蛋白质序列数据和索引文件（蛋白质序列后缀是fa）

1. 在alignment文件夹内创建database文件夹
2. 将Arabidopsis_protein.fa.gz拷贝到alignment/database中，并解压获得蛋白质序列文件Arabidopsis_protein.fa

建立索引

建立索引需要输入文件为蛋白质序列数据，类型可以选蛋白或者核酸，输出索引需要制定一个名称，最终生成比对索引数据库。

运行makeblastdb软件，对Arabidopsis_protein.fa建立索引文件，以下是该软件的提示信息。

$ makeblastdb -help
USAGE
  makeblastdb [-h] [-help] [-in input_file] [-input_type type]
    -dbtype molecule_type [-title database_title] [-parse_seqids]
    [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
    [-mask_desc mask_algo_descriptions] [-gi_mask]
    [-gi_mask_name gi_based_mask_names] [-out database_name]
    [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
    [-taxid_map TaxIDMapFile] [-version]
-dbtype <String, `nucl', `prot'>

• -in:输入文件（xxx.fa)
• -dbtype:选择类型（prot\nucl)
• out:输出索引文件名（index）

makeblastdb -in input_file -dbtype prot/nucl -out output_prefile

建立索引完成后，可以在当前目录下发现几个新生成的索引文件。

序列比对

运行blast子程序

• blastn：将给定的核酸序列与核酸数据库中的序列进行比对
• blastp：使用蛋白质序列与蛋白质数据库中的序列进行比对
• blastx：将给定的核酸序列按照六种阅读框架将其翻译成蛋白质与蛋白质数据库中的序列进行比对
• tblastn：将给定的氨基酸序列与核酸数据库中的序列（双链）按不同的阅读框进行比对
• tblastx：将DNA被检索的序列和核酸序列数据库中的序列按不同的阅读框全部翻译成蛋白质序列，然后进行蛋白质序列比对

blastn -query input_file -db database -evalue 1e-5 -out wheattest_Aradidopsis.out

实际操作示例代码：

blastx -query ./wheat-test.fasta -db ./database/index -evalue 1e-5 -out wheattest_Aradidopsis.out

query：待查询的序列，可以理解为需要进行比对的序列。

db：本地的数据库，可以理解为由已知信息构建的索引模板，用于对待查询序列进行比对。

evalue：比对的阈值，即判断比对准确性的标准。

最终结果

会自动生产一个包含比对结果的文件，利用less -S打开后可以看到一系列比对信息。

提取满足条件的指定数据

• 利用awk命令

awk '$6 > 5 && sqrt($7*$7) > 1 && $10 < 0.01 {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7 "\t" $10}' ./ matrix.DE_results.txt > name.list

• 利用cat命令

cat matrix.DE_results.txt | awk '$6 > 5 && sqrt($7*$7) > 1 && $10 < 0.01 {print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7 "\t" $10}' > name.list

上面两条命令均表示对matrix.DE_results.txt文件进行筛选，生成新文件name.list

END

© 素材来源于网络，侵权请联系后台删除

往期推荐：

文献丨群体转录组分析锁定关键转录因子

文献丨转录组RNA seq——青年阶段！

笔记丨ggplot2热图入门学习笔记

笔记丨PCA分析基本知识和数学原理