blastp结果提取（测试）

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%），可以初步推测这两个蛋白质在进化上可能具有一定的关联性，但最终的判断还需要考虑其他证据和信息。