BLASTP是一种用于在蛋白质数据库中搜索蛋白质序列相似性的BLAST算法。它可以用于找到与给定蛋白质序列相似的已知蛋白质序列，从而帮助确定其功能、结构等。
不依赖外部模块将fasta文件按>号后的字母排序：
# -*- coding: utf-8 -*-
# 读取FASTA文件
fasta_file = "nip"
# 读取文件内容
with open(fasta_file, "r") as f:
lines = f.readlines()
# 解析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文件
output_file = "sorted.fasta"
with open(output_file, "w") as f:
for identifier, sequence in sorted_sequences:
f.write(">" + identifier + "\n")
f.write(sequence + "\n")
print("FASTA序列已按照标识符的字母顺序排序并写入到sorted.fasta文件中。")
排序前后
#提取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没有
提取比对结果
# -*- 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("匹配内容已写入到输出文件中。")
#以下是一个Python脚本的示例，用于读取当前目录下的所有文件，并去除文件开头的空格：
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)
通常情况下，蛋白质比对的identity大于等于30%可以用作初步判断两个蛋白质来自同一个基因的依据。这是因为在同一物种中，不同个体的同一个基因序列可能存在一些变异，导致它们的蛋白质序列在某些位置上有一定差异。
然而，需要注意的是，30%的identity只是一个相对较低的阈值，不能作为最终确定两个蛋白质来自同一个基因的绝对标准。为了更准确地判断两个蛋白质是否来自同一个基因，还需要结合其他因素进行综合分析，如比对覆盖度、功能保守性、基因结构一致性等。
此外，对于不同物种之间的蛋白质比对，由于物种的进化差异，所以identity阈值可能需要相应调整。一般来说，物种间的蛋白质identity较高时（如大于70%），可以初步推测这两个蛋白质在进化上可能具有一定的关联性，但最终的判断还需要考虑其他证据和信息。
网友评论