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，条件解释：

1. 第6列大于5
2. 第7列绝对值大于1，第10列小于0.01

对满足上面条件的数据输出第1、2、3、4、5、6、7、10列，保存为新文件。

---