使用eggNOG-mapper获取gene的COG

安装eggnog-mapper

安装eggnog-mapper，参照其他教程，配置python环境，下载数据库文件

人类gene DNA序列提取

使用seqkit进行提取

注意，seqkit提取出来的序列id使用的其gene id，导致提取出来的序列文件有多个相同的gene id，因此需要指定seqkit所提取的序列以及对应的id。

./seqkit subseq --feature gene --gtf-tag transcript_id --gtf data/human/Homo_sapiens.GRCh38.107.chr.gtf -o data/human/extracted_dna_seqs.fa data/human/Homo_sapiens.GRCh38.dna.primary_assembly.fa

将fa文件中序列多行转为单行

# fa_process.py
file=open('hm_gene.fa','r')
file_out=open('hm_gene_out.fa','a')
seq=''
for line in file.readlines():
    if '>' in line:
        if seq != '':
            file_out.write(seq+'\n')
        file_out.write(line)
        seq=''
    else:
        seq+=line.replace('\n','')
file_out.write(seq)
file.close()
file_out.close()

同时注意fa文件中>后为gene id

使用eggNOG-mapper进行注释

这里将输入格式设置为CDS（–itype CDS），且不设置–translate，这里的gene序列实际上是ORF （参照github相关提问 提问1 提问2）
e.g.

python emapper.py -i gene_sequence/Zm-B73-REFERENCE-NAM-5.0_Zm00001eb.1.gene.fa --output output/maize_cds_out --tax_scope root --itype CDS --cpu 24

cog结果处理

import numpy as np

# 在输出结果中整理gene对应的COG
gene_cog={} # key gene_id value cog_id
file = open('maize_cds_out.emapper.annotations','r')
for line in file.readlines():
    if line[0]=='#':
        continue
    data=line.split()
    gene_id=data[0]
    cogs=data[4].split(',')
    if gene_id not in gene_cog.keys():
        gene_cog[gene_id]=[]
    cog=cogs[-1]
    cog_id = cog[:cog.index('@')]
    gene_cog[gene_id].append(cog_id)
file.close()
np.save('maize_gene_cog',gene_cog)

# 读取
# gene_cog=dict(np.load('maize_gene_cog.npy',allow_pickle=True).item())
# print(1)

以下内容为针对论文实验的笔记

出现的问题
human中所使用的cds序列中有部分gene不存在于gtf中，因此将这部分序列单独提取出来，进行上述分析，得到COG结果。

首先将cds fasta多余换行去掉

# process cds fasta
file=open('Homo_sapiens.GRCh38.cds.all.fa','r')
file_out=open('hm_cds.fa','a')

seq=''
for line in file.readlines():
    if '>' in line:
        if seq != '':
            file_out.write(seq+'\n')
        geneid = line.split()[3][5:]
        if '.' in geneid:
            geneid = geneid[:geneid.index('.')]
        file_out.write('>' + geneid + '\n')
        # file_out.write(line)
        seq=''
    else:
        seq+=line.replace('\n','')
file_out.write(seq)

file.close()
file_out.close()

提取缺失基因序列

# missing gene id in human dataset
sel_gene_id=np.load('sel_gene_id.npy')
gene_id_list=[]
file = open('human_cds_out.emapper.annotations','r')
for line in file.readlines():
    if line[0]=='#':
        continue
    data=line.split()
    gene_id=data[0]
    gene_id_list.append(gene_id)
file.close()

missing_gene_ids=[]
for geneid in sel_gene_id:
    if geneid not in gene_id_list:
        missing_gene_ids.append(geneid)
print(missing_gene_ids)
print(len(missing_gene_ids))

file=open('hm_cds.fa','r')
file_out=open('human_missing_genes.fa','a')

cds_genes=[]
lines = file.readlines()
for i, line in enumerate(lines):
    if '>' in line:
        geneid=line[1:-1]
        print(geneid)
        cds_genes.append(geneid)
        if geneid in missing_gene_ids:
            file_out.write('>'+geneid+'\n')
            file_out.write(lines[i+1])

file.close()
file_out.close()