（24）放射盘石耳基因与南极石耳基因相对【本地blast大法】

背景：因为放射盘石耳的全基因组测序和注释预测结果较好，因此在last比对的时候是利用了放射盘石耳作为参考基因组（orthofinder同源比对的时候咱是拿了南极石耳做的），在进行PAML正选择分析的时候虽然基因gapfiltered目录下也都是放射盘石耳的基因序列，但是会选择南极石耳作为前景枝，又介于咱们最终想分析南极石耳对于极寒适应有哪些关键基因，所以要做一步放射盘石耳到南极石耳的对应。

（一）借助orthofinder、MEGA图形比对法

1.根据Orthogroups.tsv中查找同源基因，也可以去Orthogroups.Genecount中查一下每个物种中数量，看是否单拷贝）

 Orthogroups.tsv，可以看谁和谁是同源的

Orthogroups.Genecount，可以看基因在每个物种中的数量 

2. 找到XX和XX对应之后，去gapfiltered目录下找到相应的UmXXcXXXX序列，less一下

可以看到是24个物种的比对情况，其中就有比对的南极石耳的相关序列，但是可能会有空缺，不全等，不要紧

3.将gapfitered中找到的比对序列和查到的同源南极石耳序列放到一起，用MEGA比对一下，看序列是否一致

如Um18c0050和U.anFUNXX一致，则把这两条放一起用MEGA看一下

优点：傻瓜，三五个的还好弄一些，多了不行

缺点：比较麻烦，有的时候比不准

（二）利用代码、本地blast进行比对

比如说，想找一系列放射盘石耳基因对应的南极石耳基因是哪些

基础知识：Blastn核酸与核酸比对

Blastx核苷酸翻译成蛋白，在蛋白层面进行比对

Blastp蛋白与蛋白库比对

Tblastn蛋白与核酸库进行比对

（就是对于未知的蛋白序列，蛋白库中可能找不到对应的，可以在核酸库中进行寻找，把核酸库翻译成蛋白再进行比对）

1.先把这些放射盘石耳基因记到grep.txt文件中去，每行一个基因文件名字

（2）创建新目录，把这些文件先从gapfiltered中提取出来，注意代码运行时路径引用规范

mkdir new_directory  # 创建新目录
while IFS= read -r file; do cp "/0926bidui/$file" new_directory/; done < grep.txt #提取文件
如：
mkdir branchsite_ok55/
while IFS= read -r file; do cp "$file" ../branchsite_ok55/; done < ../grep.txt

（3）对于提取得到的文件中每一个序列后面加上文件名

（因为本身都是物种比对，都一样，加上文件名好区分，像下面这样）

脚本：add_name.py 存放在~/Qxy/qxyjiaoben/fromzyr下了

python3 ~/Qxy/qxyjiaoben/fromzyr/add_name.py

注意脚本中file_name.replace和文件路径可能要修改，修改之后文件名字变成XX_modified

import os

def add_filename_to_sequences(file_path):
    # 获取文件名和扩展名
    file_name, file_ext = os.path.splitext(file_path)

    with open(file_path, 'r') as file:
        lines = file.readlines()

    new_lines = []
    i = 0
    while i < len(lines):
        line = lines[i]
        if line.startswith('>'):
            # 获取原始序列名字
            sequence_name = line.strip()
            # 提取文件名部分并替换特定字符串
            file_name_part = file_name.replace("/ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0926bidui/getsquence/", "")
            # 添加文件名到序列名字后面，并用_下划线连接
            new_sequence_name = f"{sequence_name}_{file_name_part}"
            new_lines.append(new_sequence_name + '\n')
        else:
            new_lines.append(line)
        i += 1

    new_file_path = f"{file_name}_modified{file_ext}"

    with open(new_file_path, 'w') as new_file:
        new_file.writelines(new_lines)

    print(f"已生成修改后的文件：{new_file_path}")


# 替换为实际的目录路径
directory_path = "/ifs1/User/wuqi/Qxy/knowngenome/gongxianxing/0926bidui/getsquence"

# 遍历目录下所有.fa文件
for file_name in os.listdir(directory_path):
    if file_name.endswith("_modified.fa"):
        file_path = os.path.join(directory_path, file_name)
        add_filename_to_sequences(file_path)

 （4）从得到的修改名字的序列中，单独把南极石耳比对序列挑出来

for file in ./*_modified.fa; do awk '/^>Umbilicaria_antarctica_genomic_/{print; getline; print}' "$file"; done > Umbilicaria_antarctica_branchesite.fa

这个代码的意思是对文件夹下的每一个XX_modified文件中U.an序列挑出来形成一个新的文件

（5）本地blast，其实就俩步，一建库，二比对

首先需要找到F:\NCBI\blast-BLAST_VERSION+\bin目录

然后在路径那里输入cmd变成命令行模式

输入下面代码建库和比对

makeblastdb.exe -in ALL_Usnea_nanji_prot.fasta -hash_index -dbtype prot  #in后面需要改，in后面跟建库的文件名称
#建立本地的blast数据库，用in 后面的fasta文件建库，建立索引index，数据库类型是氨基酸就填prot，DNA就填nucl

makeblastdb.exe -in Umbilicaria_antarctica_cds.fasta -hash_index -dbtype nucl

blastn.exe -db Umbilicaria_muehlenbergii_genomic.fasta -query Umbilicaria_muehlenbergii.cds-transcripts.fasta -out Umbilicaria_muehlenbergii_genomic.txt -outfmt=6 –evalue 10e-3 -max_hsps 1 -num_alignments 1
#本地blast2.0版本，添加了设置 -outfmt=6 输出格式6的样式，以table形式显示，–evalue E-value值，改小了，缩小比对限制  -max_hsps 1 -num_alignments 1 设置输出结果为只输出一个最好的

blastn.exe -db Umbilicaria_antarctica_cds.fasta -query branchsite_ok55.fa -out branchsite_ok55.txt -outfmt=6 –evalue 10e-3 -max_hsps 1 -num_alignments 1

比对速度很快的~

结果从左往右依次是

query id

subject id

%identity

alignment length

mismatches

gap opens

q.start

q.end

s.start

s.end

evalue

bit score

附：blast比对原始代码1.0，结果很冗长，不太方便阅读

blastn.exe -db Umbilicaria_antarctica_cds.fasta -query branchsite_ok55.fa -out branchsite_ok55_2.txt

整理之后的结果

当然还会有一些比对不到，再议再议~