#gene的ID会和mRNA一致,mRNA的ID和Parent信息也会相同,可能会出错
Scaffold_1 GLEAN gene 13403 73001 0.716366 - . ID=BTA000001.1;
Scaffold_1 GLEAN mRNA 13403 73001 0.716366 - . ID=BTA000001.1;Parent=BTA000001.1;
#正确
Scaffold_1 GLEAN gene 13403 73001 0.716366 - . ID=BTA000001;
Scaffold_1 GLEAN mRNA 13403 73001 0.716366 - . ID=BTA000001.1;Parent=BTA000001;
import re
# 读入文件
with open('input.gff', 'r') as f:
lines = f.readlines()
# 初始化基因ID
gene_id = ''
# 遍历每一行
for line in lines:
if line.startswith('##'):
# 忽略以##开头的注释行
continue
elif line.startswith('#'):
# 对于以#开头的注释行,直接输出
print(line.strip())
else:
# 对于数据行
fields = line.strip().split('\t')
# 如果是mRNA行
if fields[2] == 'mRNA':
# 获取mRNA ID
mrna_id = re.findall(r'ID=([^;]+)', fields[8])[0]
# 生成新的基因ID
gene_id = re.sub(r'\.\d+', '', mrna_id)
# 添加gene行
gene_line = [fields[0], fields[1], 'gene', fields[3], fields[4], fields[5], fields[6], fields[7], f'ID={gene_id};']
print('\t'.join(gene_line))
# 修改mRNA行
mrna_line = fields
mrna_line[2] = 'mRNA'
mrna_line[8] = f'ID={mrna_id};Parent={gene_id};'
mrna_line[8] = re.sub(r';source_id=[^;]+', '', mrna_line[8])
print('\t'.join(mrna_line))
else:
# 对于其他行直接输出
print(line.strip())
对于gene行,直接输出;
对于mRNA行,根据ID属性值生成gene行的ID属性值,并输出gene行和修改后的mRNA行;
对于其他行(CDS/exon等),根据Parent属性值生成修改后的行,并输出;
3.将处理后的结果写入输出文件。
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!