BLASTP是一种用于在蛋白质数据库中搜索蛋白质序列相似性的BLAST算法。它可以用于找到与给定蛋白质序列相似的已知蛋白质序列，从而帮助确定其功能、结构等。
不依赖外部模块将fasta文件按>号后的字母排序：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
if len(sys.argv)!= 3:
print("用法: python script.py <输入的FASTA文件路径> <输出的FASTA文件路径>")
sys.exit(1)
# 获取命令行传入的输入文件和输出文件路径
input_fasta_file = sys.argv[1]
output_fasta_file = sys.argv[2]
# 读取FASTA文件
try:
with open(input_fasta_file, "r") as f:
lines = f.readlines()
except FileNotFoundError:
print(f"输入的文件 {input_fasta_file} 不存在，请检查文件路径是否正确。")
sys.exit(1)
# 解析FASTA文件
sequences = {}
current_id = None
current_sequence = ""
for line in lines:
line = line.strip()
if line.startswith(">"):
if current_id:
sequences[current_id] = current_sequence
current_sequence = ""
current_id = line[1:]
else:
current_sequence += line
# 添加最后一个序列
if current_id:
sequences[current_id] = current_sequence
# 按照标识符的字母顺序对序列进行排序
sorted_sequences = sorted(sequences.items(), key=lambda x: x[0])
# 写入排序后的序列到新的FASTA文件
try:
with open(output_fasta_file, "w") as f:
for identifier, sequence in sorted_sequences:
f.write(">" + identifier + "\n")
f.write(sequence + "\n")
except:
print(f"写入输出文件 {output_fasta_file} 时出现错误，请检查相关权限或磁盘空间等情况。")
sys.exit(1)
print(f"FASTA序列已按照标识符的字母顺序排序并写入到 {output_fasta_file} 文件中。")
排序前后
#提取id
grep ">" sorted.fasta |sed s'/>//g' >nip
#blastp参数
cat pepid |while read id
do
blastp -db pep/$id.fa \
-query sorted.fasta \
-out bp/$id.bl \
-evalue 1e-10 \
-best_hit_overhang 0.25 \
# -perc_identity 0.5 \
-perc_identity 0.5 \
-max_target_seqs 5 \
-num_threads 10 -outfmt 6
# -outfmt "6 qacc sacc evalue pident length qstart qend sstart send"
#blastp -db pep/$id.fa -query GHpep -out bl/$id.bl -num_threads 8 -evalue 1e-100
done
blastn是有 -perc_identity 0.5这个参数的，但是blastp没有
提取比对结果
##f1.py
# -*- coding: utf-8 -*-
import sys
# 检查命令行参数数量
if len(sys.argv) != 4:
print("请提供输入文件名、另一个文件名和输出文件名作为命令行参数。")
sys.exit(1)
# 获取命令行参数
input_file_name = sys.argv[1]
another_file_name = sys.argv[2]
output_file_name = sys.argv[3]
# 打开输入文件和输出文件
with open(input_file_name, 'r') as input_file, open(another_file_name, 'r') as another_file, open(output_file_name, 'w') as output_file:
# 遍历输入文件的每一行
for line in input_file:
# 在每一行中查找匹配的字符串
if line.strip(): # 如果行不为空
# 在另一个文件中查找匹配的内容
for another_line in another_file:
if line.strip() in another_line:
# 继续查找并将结果写入到输出文件
for next_line in another_file:
if ">" in next_line:
output_file.write(next_line.strip() + "\n")
break # 找到第一个 '>'后停止继续查找
output_file.write(next_line)
break # 找到匹配后停止继续查找
# 打印完成提示
print("匹配内容已写入到输出文件中。")
cat ../pepid |while read id
do
python f1.py gene $id.bl find/$id
grep ">" find/$id|sed s'/>//g'|sed s'/predict//g' >res/$id
done
NIP里query和target是完全一致的，有43个query，其他的samples里也有43个结果，得到一一匹配的基因
###f2.py
import sys
def find_match(input_file, search_file, output_file):
# 读取输入文件的内容
with open(input_file, 'r') as f:
input_content = f.readlines()
# 读取搜索文件的内容
with open(search_file, 'r') as f:
search_content = f.readlines()
# 是否找到匹配的标志
match_found = False
# 是否已经有相同的匹配的标志
duplicate_match = False
# 查找匹配的内容
for line in input_content:
line = line.strip()
# 判断是否为匹配的行
if line in search_content:
# 如果已经有相同的匹配，则查找第二个">"
if duplicate_match:
match_found = False
duplicate_match = False
continue
match_found = True
continue
# 如果找到匹配，且未出现相同的匹配
if match_found and not duplicate_match:
# 写入输出文件
with open(output_file, 'a') as f:
f.write(line + '\n')
# 如果找到">"，则标记已经有相同的匹配
if line == '>':
duplicate_match = True
if not match_found:
print("未找到匹配的内容")
if __name__ == "__main__":
# 获取命令参数
if len(sys.argv) < 4:
print("请提供输入文件路径、搜索文件路径和输出文件路径作为命令参数")
sys.exit(1)
input_file = sys.argv[1]
search_file = sys.argv[2]
output_file = sys.argv[3]
# 查找匹配的内容并写入输出文件
find_match(input_file, search_file, output_file)
import os
# 获取当前目录下的所有文件名
files = os.listdir()
# 遍历每个文件
for file in files:
# 判断是否为文件，排除文件夹
if os.path.isfile(file):
# 读取文件内容
with open(file, 'r') as f:
content = f.readlines()
# 去除每行开头的空格
modified_content = [line.lstrip() for line in content]
# 将修改后的内容写回文件
with open(file, 'w') as f:
f.writelines(modified_content)
image.png
主要是为了提取出基因组上最符合的匹配，结果和之前手动统计的一样
###，把当前目录下所有的文件合并成一个文件，按列合并，每一列的名字为原文件的名字,注意每一个文件都要有相同的行数
import os
import pandas as pd
# 获取当前目录下的所有文件名
files = os.listdir()
# 用于存储每个文件的内容
file_content = {}
# 遍历每个文件
for file in files:
# 判断是否为文件，排除文件夹
if os.path.isfile(file):
# 获取文件名作为列名
column_name = os.path.splitext(file)[0]
# 读取文件内容并存储到字典中
with open(file, 'r') as f:
content = f.readlines()
file_content[column_name] = content
# 将字典转换为DataFrame
df = pd.DataFrame(file_content)
# 将DataFrame写入合并后的文件
df.to_csv('merged_file.csv', index=False)
得到NIP和其他samples对应的id
通常情况下，蛋白质比对的identity大于等于30%可以用作初步判断两个蛋白质来自同一个基因的依据。这是因为在同一物种中，不同个体的同一个基因序列可能存在一些变异，导致它们的蛋白质序列在某些位置上有一定差异。
然而，需要注意的是，30%的identity只是一个相对较低的阈值，不能作为最终确定两个蛋白质来自同一个基因的绝对标准。为了更准确地判断两个蛋白质是否来自同一个基因，还需要结合其他因素进行综合分析，如比对覆盖度、功能保守性、基因结构一致性等。
此外，对于不同物种之间的蛋白质比对，由于物种的进化差异，所以identity阈值可能需要相应调整。一般来说，物种间的蛋白质identity较高时（如大于70%），可以初步推测这两个蛋白质在进化上可能具有一定的关联性，但最终的判断还需要考虑其他证据和信息。
文章里对稻瘟病鉴定的依据
参考：
[1]林晓涵. 基于水稻泛基因组的抗病等位基因比较分析[D].广西大学,2023.DOI:10.27034/d.cnki.ggxiu.2023.002874.
