猪代谢相关基因注释

1、参考基因组下载

基因组选择：Genome assembly Sscrofa11.1 （官方参考基因组，杜洛克雌猪）

wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_protein.faa.gz
 
# MD5校验
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/md5checksums.txt
 
md5sum -c md5checksums.txt

2、基因组注释

uniprot注释 （包括下载，建库和注释三部分）

wget -c https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -d uniprot_sprot.fasta.gz

makeblastdb -in ~/0_rawdata/database/uniprot/uniprot_sprot.fasta -dbtype prot -out ~/0_rawdata/database/uniprot/uni
blastdbcmd -info -db  /home/xb/0_rawdata/database/uniprot/uni # 检查数据库构建是否正确

nohup blastx -query /data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_genomic.fna -db ~/0_rawdata/database/uniprot/uni -max_target_seqs 1 -outfmt 6 -evalue 1e-5 > uniprot.out &

# 已放弃，速度慢且达不到想要的结果。
jobs -l 
kill -9 20239

# 重新对faa文件进行注释
diamond makedb --in uniprot_sprot.fasta -d uni

nohup diamond blastp -d ~/0_rawdata/database/uniprot/uni.dmnd -q /data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_protein.faa -o pig.uni.xml -f 5 --sensitive --max-target-seqs 20 -e 1e-5 --id 20 --tmpdir /dev/shm --index-chunks 1 -p 8 &

parsing_blast_result.pl --evalue 1e-5 --HSP-num 1 --out-hit-confidence --suject-annotation pig.uni.xml > pig.uni.tab

根据参考文献方法，发现eggnog注释可能更切合我们的需求，能直接得出带有EC号的蛋白。

软件安装和使用参考：https://www.yunbios.net/eggNOG.html （包括下载安装，数据库下载，比对）

conda create -n eggnog
conda activate eggnog
conda install -c bioconda -y eggnog-mapper

wget -c http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog_proteins.dmnd.gz
gzip -d eggnog_proteins.dmnd.gz   # 测试失败，得用软件提供的脚本下载数据库。

# 首先要手动创建一个data文件夹，脚本默认不创建。
mkdir -p /home/xb/miniconda3/envs/eggnog/lib/python3.12/site-packages/data

download_eggnog_data.py #使用官方脚本下载数据库

nohup emapper.py -i /data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_protein.faa -o pig --cpu 4 --seed_ortholog_evalue 1e-5 --dmnd_db /home/xb/miniconda3/envs/eggnog/lib/python3.12/site-packages/data/eggnog_proteins.dmnd &

# 运行时间大概在30min，下载pig.emapper.annotations至本地，提取注释到含有EC号的蛋白（12359条），并提取序列。
sort -n id.list |uniq > rmdup.list   #去重复，虽然eggnog注释貌似不像blast一样有重复序列

fasta_extract_subseqs_from_list.pl /data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_protein.faa rmdup.list > target.faa

# 使用cd-hit去除重复序列，然后回比到基因组上。
conda install -c bioconda cd-hit
cd-hit -i target.faa -o target.0.95.faa -c 0.95   # 去重后剩余 5270 条

mkdir blastx && cd blastx
makeblastdb -in ../target.0.95.faa -dbtype prot -out ./ec
blastdbcmd -info -db ./ec

nohup blastx -query /data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_genomic.fna -db ./ec  -max_target_seqs 1 -outfmt 5 -evalue 1e-5 -num_threads 8 > pig.out &    # 无结果，转化成tab显示没有比对上的，因为分析时间太长，无法复现

3、根据protein_id提取对应的gene_id（新思路)

查询资料，发现可以直接从gbff中提取gene_id，整理一个脚本测试该功能（我爱gpt！）

from Bio import SeqIO

def extract_gene_ids(gbff_file, protein_ids_file, output_file):
    # 读取蛋白质ID列表
    with open(protein_ids_file, 'r') as f:
        protein_ids = [line.strip() for line in f.readlines()]

    # 创建一个字典来存储蛋白质ID与基因ID的映射
    protein_to_gene = {}

    # 解析GBFF文件
    for record in SeqIO.parse(gbff_file, "genbank"):
        for feature in record.features:
            if feature.type == "CDS":
                # 提取protein_id
                if "protein_id" in feature.qualifiers:
                    protein_id = feature.qualifiers["protein_id"][0]
                    if protein_id in protein_ids:
                        # 提取GeneID
                        gene_ids = [xref.split(":")[1] for xref in feature.qualifiers.get("db_xref", []) if xref.startswith("GeneID")]
                        if gene_ids:
                            protein_to_gene[protein_id] = gene_ids[0]  # 只取第一个GeneID

    # 输出结果到文件
    with open(output_file, 'w') as out_f:
        out_f.write("Protein_ID\tGene_ID\n")
        for protein_id in protein_ids:
            gene_id = protein_to_gene.get(protein_id, "Not found")
            out_f.write(f"{protein_id}\t{gene_id}\n")

# 示例使用
gbff_file = "/data/xb/1_pig_genome/GCF_000003025.6_Sscrofa11.1_genomic.gbff"          # 替换为你的GBFF文件路径
protein_ids_file = "/home/xb/1_results/0906_pig_genome_annot/eggnog/rmdup.list"  # 替换为包含蛋白质ID的txt文件路径
output_file = "/home/xb/1_results/0906_pig_genome_annot/eggnog/extract_from_gbff/gene.id.txt"     # 输出结果的文件路径

extract_gene_ids(gbff_file, protein_ids_file, output_file)

# 提取第二列的gene_id并去重
cat gene.id.txt|cut -f2 > id.list
sort -n id.list |uniq > rmdup.list   #包括 3691 条gene

无效路径，仅做记录。

# 再整理一个脚本，根据gene_id先从NCBI中批量提取对应的gene 序列文件（需要先获取API接口）。
# gpt回答的无法实现

# 要根据GeneID提取gene的fasta序列需要分为两步：将GeneID转换为相应的nucleotide ID或者RefSeq ID。GeneID本身通常在 `gene` 数据库中用于注释和搜索，但实际的序列数据在 `nucleotide` 数据库中。
# gpt回复任然无效

参考一个帖子找到灵感：如何从NCBI上的Gene数据库批量下载基因序列数据_ncbi批量下载基因序列-CSDN博客 [官方指导文件](Supported programming languages (nih.gov))

实操过程如下：第一步，通过OpenAPI java libraries建立Build Python NCBI Datasets API v2alpha library

#!/usr/bin/env bash  编写一个bash脚本，直接运行即可（安装路径为：/home/xb/3_opt/ncbi_api）。
OUTPUT_DIR="python_lib"
 
wget https://www.ncbi.nlm.nih.gov/datasets/docs/v2/openapi3/openapi3.docs.yaml
 
wget https://repo1.maven.org/maven2/org/openapitools/openapi-generator-cli/7.2.0/openapi-generator-cli-7.2.0.jar -O openapi-generator-cli.jar
 
java -jar openapi-generator-cli.jar generate -g python -i openapi3.docs.yaml --package-name "ncbi.datasets.openapi" --additional-properties=pythonAttrNoneIfUnset=true,projectName="ncbi-datasets-pylib"

这一步之后要输入： pip install . （非常重要，官方文档没有提到，但是不运行此步骤后面会报错）

第二步，编写提取脚本gene_get_info.py

import sys
import csv
import io
import os
from typing import List
from zipfile import ZipFile
 
from ncbi.datasets.openapi import ApiClient as DatasetsApiClient
from ncbi.datasets.openapi import ApiException as DatasetsApiException
from ncbi.datasets.openapi import GeneApi as DatasetsGeneApi
 
zipfile_name = "gene.zip"  #自定义下载的压缩包名称
output_file_name = "protein.faa"  # 自定义输出的文件名称

 
# 从CSV文件中读取基因ID并存储在列表中
gene_ids_from_csv = []
csv_file_path = "/home/xb/1_results/0906_pig_genome_annot/eggnog/extract_from_gbff/rmdup.list"  # 基因ID的CSV文件路径
 
with open(csv_file_path, newline='') as csvfile:
    reader = csv.reader(csvfile)
    for row in reader:
        gene_id = int(row[0])  # 基因ID在CSV文件的第二列
        gene_ids_from_csv.append(gene_id)
 
# 将基因ID列表转换为字符串格式，形如 "[1, 2, 3, ...]"
gene_ids_string = str(gene_ids_from_csv)
 
# 将字符串中的单引号替换为双引号，使其成为合法的Python列表表示形式
gene_ids_string = gene_ids_string.replace("'", '"')
 
# 将字符串转换为Python列表
gene_ids_list = eval(gene_ids_string)
 
 
with DatasetsApiClient() as api_client:
    gene_ids: List[int] = gene_ids_list
    gene_api = DatasetsGeneApi(api_client)
    try:
        gene_dataset_download = gene_api.download_gene_package_without_preload_content(
            gene_ids,
            include_annotation_type=["FASTA_GENE", "FASTA_PROTEIN"],  #选择下载的数据格式
        )
 
        with open(zipfile_name, "wb") as f:
            f.write(gene_dataset_download.data)
    except DatasetsApiException as e:
        sys.exit(f"Exception when calling GeneApi: {e}\n")
 
 
try:
    dataset_zip = ZipFile(zipfile_name)
    zinfo = dataset_zip.getinfo(os.path.join("ncbi_dataset/data", "protein.faa"))
    with io.TextIOWrapper(dataset_zip.open(zinfo), encoding="utf8") as fh:
        with open(output_file_name, "w", encoding="utf8") as output_file:
            output_file.write(fh.read())
except KeyError as e:
    sys.exit(f"File {output_file_name} not found in zipfile: {e}")

第三步，运行程序。

source ~/3_opt/cobra/bin/activate
pip install python_lib

chmod 755 gene_get_info.py
python3 gene_get_info.py

第四步，性能测试，受限于ncbi，一次性无法下载3691条gene，经测试发现单次最多只能下载700条（循环下载的话只能600），因此需要先分割gene_id_list，然后再提取gene 序列，脚本改进如下：

# 分割脚本，按照500个切割表格，形成多个小表格供后续使用。
import csv
import os

def split_list(lst, n):
    """将列表分割成多个小列表，每个列表最多包含n个元素。"""
    for i in range(0, len(lst), n):
        yield lst[i:i + n]

def save_split_files(gene_ids, batch_size, output_dir):
    """将基因ID列表分割并保存到多个CSV文件中。"""
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    for idx, batch in enumerate(split_list(gene_ids, batch_size), start=1):
        file_path = os.path.join(output_dir, f'batch_{idx}.csv')
        with open(file_path, 'w', newline='') as csvfile:
            writer = csv.writer(csvfile)
            for gene_id in batch:
                writer.writerow([gene_id])

if __name__ == "__main__":
    csv_file_path = "/home/xb/1_results/0906_pig_genome_annot/eggnog/extract_from_gbff/rmdup.list"  # 输入CSV文件路径
    output_dir = "/home/xb/1_results/0906_pig_genome_annot/eggnog/extract_from_gbff/split_files"  # 输出目录
    batch_size = 600  # 每个文件的基因ID数量

    gene_ids_from_csv = []
    with open(csv_file_path, newline='') as csvfile:
        reader = csv.reader(csvfile)
        for row in reader:
            if len(row) > 0:
                gene_id = int(row[0])  # 基因ID在CSV文件的第一列
                gene_ids_from_csv.append(gene_id)

    save_split_files(gene_ids_from_csv, batch_size, output_dir)

然后再访问NCBI API批量下载数据（压缩包形式，按gene_id命名）

# 下载脚本，根据id下载gene序列，有的gene id在性染色体上，对应能查出两条基因序列。
import sys
import csv
import io
import os
from typing import List
from zipfile import ZipFile
from ncbi.datasets.openapi import ApiClient as DatasetsApiClient
from ncbi.datasets.openapi import ApiException as DatasetsApiException
from ncbi.datasets.openapi import GeneApi as DatasetsGeneApi

def download_and_extract_genes(gene_ids: List[int], zipfile_name: str, output_file_name: str):
    """下载基因数据并提取到文件中。"""
    with DatasetsApiClient() as api_client:
        gene_api = DatasetsGeneApi(api_client)
        try:
            gene_dataset_download = gene_api.download_gene_package_without_preload_content(
                gene_ids,
                include_annotation_type=["FASTA_GENE", "FASTA_PROTEIN"],  # 选择下载的数据格式
            )

            with open(zipfile_name, "wb") as f:
                f.write(gene_dataset_download.data)
        except DatasetsApiException as e:
            sys.exit(f"Exception when calling GeneApi: {e}\n")

    try:
        with ZipFile(zipfile_name) as dataset_zip:
            zinfo = dataset_zip.getinfo(os.path.join("ncbi_dataset/data", "protein.faa"))
            with io.TextIOWrapper(dataset_zip.open(zinfo), encoding="utf8") as fh:
                with open(output_file_name, "a", encoding="utf8") as output_file:
                    output_file.write(fh.read())
    except KeyError as e:
        sys.exit(f"File {output_file_name} not found in zipfile: {e}")

def process_csv_files(input_dir: str, output_file_name: str):
    """处理目录中的所有CSV文件并提取基因数据。"""
    csv_files = [f for f in os.listdir(input_dir) if f.endswith('.csv')]
    if not csv_files:
        sys.exit(f"No CSV files found in directory: {input_dir}")

    # 确保输出文件存在
    open(output_file_name, 'w').close()

    for csv_file in csv_files:
        csv_file_path = os.path.join(input_dir, csv_file)
        with open(csv_file_path, newline='') as csvfile:
            reader = csv.reader(csvfile)
            gene_ids_from_csv = [int(row[0]) for row in reader if len(row) > 0]
        
        # 使用 CSV 文件名作为唯一标识的一部分
        zipfile_name = f"gene_{os.path.splitext(csv_file)[0]}.zip"
        download_and_extract_genes(gene_ids_from_csv, zipfile_name, output_file_name)

if __name__ == "__main__":
    input_dir = "/home/xb/1_results/0906_pig_genome_annot/eggnog/extract_from_gbff/split_files"  # CSV文件所在目录
    output_file_name = "combined_protein.faa"  # 合并后的输出文件

    process_csv_files(input_dir, output_file_name)

解压、重命名及合并fna文件

#!/bin/bash   解压&重命名

# 确保目标文件夹存在
mkdir -p gene_output

# 遍历所有以 gene_batch_ 开头的 zip 文件
for zip_file in gene_batch_*.zip; do
    # 提取文件名（去掉扩展名）
    base_name=$(basename "$zip_file" .zip)
    
    # 解压指定文件并重命名
    unzip -j "$zip_file" ncbi_dataset/data/gene.fna -d ./gene_output/ &&
    mv ./gene_output/gene.fna ./gene_output/${base_name}.fna
done

cd gene_output/

cat *.fna > output.fna
seqtk seq -A output.fna | sort -u > output.rmdup.fna

grep ">" output.rmdup.eggnog.fna |wc -l   # 最后得到3694条gene,结果文件为output.rmdup.eggnog.fna