perl脚本将embl格式文件转换成fasta及gff3格式

已于 2023-11-24 12:01:52 修改
阅读量186 收藏 1
点赞数 1
分类专栏： 生物信息学 perl 学习笔记 文章标签： perl 生物信息学 数据库
于 2023-11-03 20:10:45 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Alter_X/article/details/134209832
版权

 最近找基因组注释的时候发现一个基因组只上传了embl格式，为了下游分析想把它分成fasta和gff3，但是该embl格式有一点奇怪，常用的转格式软件转不成功，那就自己写个脚本好了。
 该embl文件大概长这样：

ID   XXX; XXX; linear; other RNA; XXX; XXX; 177118352 BP.
XX
AC   XXX; 
XX
AC * _scaf_1
XX
PR   Project:PRJEB51672;
XX
DT   04-NOV-2022 (Rel. 133, Created)
XX
DE   XXX
XX
KW   .
XX
OS   Poa pratensis
XX
OC   cellular organisms; Eukaryota; Viridiplantae; Streptophyta; Streptophytina;
OC   Embryophyta; Tracheophyta; Euphyllophyta; Spermatophyta; Magnoliopsida;
OC   Mesangiospermae; Liliopsida; Petrosaviidae; commelinids; Poales; Poaceae;
OC   BOP clade; Pooideae; Poodae; Poeae; Poeae Chloroplast Group 2 (Poeae type);
OC   Poodinae; Poinae; Poa.
XX
RN   [1]
RP   1-177118352
RG   XXX
RT   ;
RL   Submitted (04-NOV-2022) to the INSDC.
XX
FH   Key             Location/Qualifiers
FH
FT   source          1..177118352
FT                   /mol_type="other RNA"
FT                   /organism="Poa pratensis"
FT   gene            complement(11455981..11457619)
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462"
FT                   /note="source:PanAnd"
FT   mRNA            complement(11455981..11457619)
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462_T001"
FT                   /note="source:PanAnd"
FT   exon            complement(11455981..11457619)
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462_T001.exon.1"
FT                   /note="source:PanAnd"
FT   3'UTR           complement(11455981..11456845)
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462_T001.three_prime_utr.4"
FT                   /note="source:PanAnd"
FT   CDS             complement(11456846..11457412)
FT                   /codon_start=1
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462_P001"
FT                   /note="source:PanAnd"
FT                   /transl_table=1
FT   5'UTR           complement(11457413..11457619)
FT                   /locus_tag="PPR_LOCUS256280"
FT                   /note="ID:Ppr00001aa254462_T001.five_prime_utr.3"
FT                   /note="source:PanAnd"
FT   gene            complement(11455981..11461780)
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461"
FT                   /note="source:PanAnd"
FT   mRNA            complement(join(11455981..11456778,11461090..11461780))
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_T001"
FT                   /note="source:PanAnd"
FT   exon            complement(11455981..11456778)
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_T001.exon.1"
FT                   /note="source:PanAnd"
FT   exon            complement(11461090..11461780)
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_T001.exon.2"
FT                   /note="source:PanAnd"
FT   3'UTR           complement(join(11455981..11456778,11461090..11461145))
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_T001.three_prime_utr.5"
FT                   /note="ID:Ppr00001aa254461_T001.three_prime_utr.6"
FT                   /note="source:PanAnd"
FT   CDS             complement(11461146..11461712)
FT                   /codon_start=1
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_P001"
FT                   /note="source:PanAnd"
FT                   /transl_table=1
FT   5'UTR           complement(11461713..11461780)
FT                   /locus_tag="PPR_LOCUS256281"
FT                   /note="ID:Ppr00001aa254461_T001.five_prime_utr.4"
FT                   /note="source:PanAnd"
//

 好多东西暂时没用，比如UTR之类的，就先不提取了。

#注释过几天再补！
#!/usr/bin/perl
use strict; use warnings; use Getopt::Std;

my $usage = "embl2other.pl -e <embl file> -g <output gff3 file> or -f <output fasta file>";

my %opt;
getopts('e:g:f:', \%opt);

my $AC;

#如果需要转成fasta
if($opt{f}){
    my $SQ;

	#创建或覆盖目标文件
	open(FASTA,">","$opt{f}");
	print "";
	close(FASTA);

	#打开文件
	open(FASTA,">>","$opt{f}");
	open(EMBL,"$opt{e}") or die $!;
	while(<EMBL>){
		#提取AC，即染色体/contig编号，源文件中有两个AC，抛弃内容只有XXX的一个，再存入文件
	    if(m/^AC/){
		    next if (split(/\s+/))[1] eq "XXX;";
		    $AC = (split(/\s+/))[2];
		    print FASTA ">".$AC."\n";
		}

		#从SQ行到\\行之间提取序列并存入目标文件
		if(m/^SQ/ .. m/^\/\//){
			next if (m/^SQ/);
			next if (m/^\/\//);
			s/[^A-Z]//g;
			$SQ = $_;
			print FASTA $SQ."\n";
		}
	}
	#关闭文件
	close(FASTA);
	close(EMBL);
}

#如果需要转成gff3
if($opt{g}){
    my $feature;
    my @start;
    my @end;
    my $strand;
    my $LOC;
    my @ID;
    my $source;
    my $gene_id;
    my @mRNA_id;
	my $cds_num;

	#创建或覆盖目标文件
	open(GFF,">","$opt{g}");
        print "";
        close(GFF);

	#打开文件
	open(GFF,">>","$opt{g}");
	open(EMBL,"$opt{e}") or die $!;
	while(<EMBL>){
		#提取AC，抛弃只含XXX的
		if(m/^ID/ .. m/^\/\//){
			if(m/^AC/){
			    next if (split(/\s+/))[1] eq "XXX;";
			    $AC = (split(/\s+/))[2];
			}

			#当到达FT注释行时
			if(m/^FT/){
				#跳过source注释和gap注释
				next if (m/ source|\/mol_type|\/organism|gap|\/estimated_length/);
				#按空格分割每行
				my @FT_line = split(/\s+/);
				#空格分隔后有三列时，提取feature、start、end、strand
				if(@FT_line == 3){
					$feature = $FT_line[1];
					#调用自定义函数getcoor
					my ($start_res,$end_res) = getcoor($_);
					@start = @$start_res;
					@end = @$end_res;
					#包含complement的认为是反向元素
					$strand = ($FT_line[2] =~ /complement/ ? "-" : "+");
					#初始化ID
					@ID = ();
				}

				#空格分隔后有两列时，提取LOC、ID、source
				if(@FT_line == 2){
					#跳过3'UTR、5'UTR的注释以及CDS注释的最后一行
					next if ($feature =~ /3'UTR|5'UTR/);
					next if (m/\/transl_table=/);
					$LOC = $1 if(m/\/locus_tag="(.*)"/);
					push(@ID,$1) if(m/\/note="ID:(.*)"/);
					$source = $1 if(m/\/note="source:(\w+)"/);
				}
	
				#当读取到feature为gene，且读到source，即每个feature的最后一行时将对应gene元素存入文件
				if($feature eq "gene" and $source){
					my @gene_output = ($AC,$source,$feature,$start[0],$end[0],"\.",$strand,"\.","ID=$ID[0];locus_tag=$LOC");
					print GFF join("\t",@gene_output)."\n";
					#初始化source、mRNA_id、cds_num，保存gene_id
					$source = undef;
					$gene_id = $ID[0];
					@mRNA_id = ();
					$cds_num = 0;
				}

				#当读取到feature为mRNA，且读到source，即每个feature的最后一行时将对应mRNA、exon元素存入文件
				if($feature eq "mRNA" and $source){
					#存入mRNA元素
					my @mRNA_output = ($AC,$source,$feature,$start[0],$end[-1],"\.",$strand,"\.","ID=$ID[0];Parent=$gene_id");
					print GFF join("\t",@mRNA_output)."\n";
					push(@mRNA_id,$ID[0]);

					#存入exon元素
					for my $e_num (0 .. @start-1){
						my @exon_output = ($AC,$source,"exon",$start[$e_num],$end[$e_num],"\.",$strand,"\.","ID=$ID[0].exon.".($e_num+1).";Parent=$ID[0]");
						print GFF join("\t",@exon_output)."\n";
						###如果源文件中exon命名不是从1开始连续的增加，那这里会有和源文件对不上的问题，暂时不知道该怎么优化
					}

					#初始化source
					$source = undef;
				}

				#原本想在这里提取exon保证exon命名统一，但是不同mRNA的共享exon该embl文件中只有1个，
				#只靠exon这个feature无法得到与mRNA的对应关系，如果通过坐标判断就太繁琐了，
				#所以改在处理mRNA的时候一并提取exon
				#这里就相当于跳过exon的feature段
				if($feature eq "exon" and $source){
					$source = undef;
					next;
				}

				#当读取到feature为CDS，且读到source，即每个feature的最后一行时将对应CDS元素存入文件
				if($feature eq "CDS" and $source){
					#定义phase，按坐标来说，正向的第一个CDS和反向的最后一个CDS的phase一般一直都是0
					my @phase = (0);
					#当strand为正向，且CDS数量大于一个时
					if($strand eq "+" and @start > 1){
						#计算第二个到最后一个CDS的phase
						for my $num (0 .. @start-2){
							$phase[$num+1] = (3-($end[$num]-$start[$num]+1+$phase[$num])%3);
							$phase[$num+1] = 0 if ($phase[$num+1] == 3);
						}
						#从第一个CDS开始写入GFF
						for my $i (0 .. @start-1){
							my @CDS_output = ($AC,$source,$feature,$start[$i],$end[$i],"\.",$strand,$phase[$i],"ID=$ID[0].cds.".($i+1).";Parent=$mRNA_id[$cds_num]");
							print GFF join("\t",@CDS_output)."\n";
						}
						
					#当strand为反向，且CDS数量大于1个时
					}elsif($strand eq "-" and @start > 1){
						#计算倒数第二个到第一个CDS的phase 
						for my $num (0 .. @start-2){
							$phase[$num+1] = (3-($end[@start-2-$num+1]-$start[@start-2-$num+1]+1+$phase[$num])%3);
							$phase[$num+1] = 0 if ($phase[$num+1] == 3);
						}
						#从第一个CDS开始写入GFF
						for my $i (0 .. @start-1){
							my @CDS_output = ($AC,$source,$feature,$start[$i],$end[$i],"\.",$strand,$phase[@start-1-$i],"ID=$ID[0].cds.".($i+1).";Parent=$mRNA_id[$cds_num]");
							print GFF join("\t",@CDS_output)."\n";
						}
					
					#当只存在一个CDS的时候，直接按phase=0写入
					}elsif(@start == 1){
						my @CDS_output = ($AC,$source,$feature,$start[0],$end[0],"\.",$strand,"0","ID=$ID[0].cds."."1".";Parent=$mRNA_id[$cds_num]");
                                                print GFF join("\t",@CDS_output)."\n";
					}
					#重置source变量，准备下一轮判断
					$source = undef;
					#每完成一轮CDS的输出，计数加1
					$cds_num += 1;
					#该格式文件中mRNA数量应该和cds输出轮数相等，如果cds轮数大于mRNA数量说明解析出错
					#其实小于应该也要die掉，但判断就不能放在这里了，还得修改一下
					die("Number of CDS and Number of mRNA could not be paired!") if ($cds_num > scalar(@mRNA_id));
				}
			}
		}
	}
	close(EMBL);
	close(GFF);
}

#自定义函数，作用是给一个feature的起始行进去，成对输出若干个元素的起始位点和终止位点
sub getcoor{
	my ($line) = @_;
	my @start;
	my @end;
	#去掉行尾换行符
	chomp $line;
	#当行末尾存在逗号时
	if($line =~ /,$/){
		#去掉行尾换行符
		chomp $line;
		#按空格分割取第三列，即保存坐标的列
		$line = (split(/\s+/,$line))[2];
		#读取该行的下一行
		my $next_line = <EMBL>;
		#当下一行末尾也存在逗号时
		while($next_line =~ /,$/){
			#去掉行尾的换行符
			chomp $next_line;
			#按空格分割取第二列，即保存坐标的列
			my $next_line_split = (split(/\s+/,$next_line))[1];
			#将该行添加到上一行末尾存入line变量
			$line .= $next_line_split;
			#提取再下一行继续while判断
			$next_line = <EMBL>;
		}
		#去掉不满足while条件的行末尾的换行符
		#其实这个地方有点小bug，如果没进入上面的循环的话next_line应该处于未声明状态，可能会导致warning
		chomp $next_line;
		#按空格分隔取该行第二列，即保存坐标的列
		my $next_line_split = (split(/\s+/,$next_line))[1];
		#将该行添加到上一行末尾存入line变量
		$line .= $next_line_split;
	}
	#按逗号分隔line并按组存入group中
	my @group = split(",",$line);
	#遍历group分别按顺序提取起始位点和终止位点存入两个数组中
	foreach my $group (@group){
		if($group =~ /(\d+)\.\.(\d+)/){
			push(@start,$1);
			push(@end,$2);
		}
	}
	#返回两个数组
	return(\@start,\@end);
}

写完之后感觉还是面向对象编程的能力不是很够，这个脚本还是特解，想快速写出能获得通解的脚本还是能力不够，继续锻炼吧。（如果有大佬路过请多指教，谢谢！）

Alter X
关注
  • 1
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
perl 小练习
Cassiel60的博客
03-12 977
1.以scaffold.fasta作为输入文件，计算GC含量以及N50和N90 2.根据给定的基因组scaffold.fasta文件和相对用的基因注释gff文件提取基因的cds区域，并以每行60个碱基的格式输出到cds.fasta文件中 3.以cds.fasta作为输出文件，将其翻译成蛋白质序列并以每行60个氨基酸的格式输出到pep.fasta文件 ...
【MATLAB生信分析】MATLAB生物信息分析工具箱（一）
W_LAILAI的博客
11-10 4721
这里给出 MATLAB_R2017a 的生物信息学工具箱中 自带的样例 一览： 如果是默认安装，则相关目录在： C:\Program Files\MATLAB\R2017a\examples\bioinfo 按字母排序如下： AlignMultipleSequencesExample.m AlignQuerySequenceToProfileUsingHMMModelAlignmentExa...
perl脚本根据gff3文件的提取CDS最长的转录本
最新发布
Alter_X的博客
12-05 884
perl脚本根据Coge的gff3文件提取CDS最长的转录本
生信自学笔记（三）分子数据库
羊城迷鹿的博客
07-27 4303
分类 基因组数据库：GDB、GenBank、Ensembl 核酸数据库：ENA（EMBL）、GenBank、DDBJ 蛋白质数据库：SWISS-PROF、PIR、PDB 一个数据库记录由两部分组成：原始序列数据和描述这些数据生物学信息的注释。 记录格式 FASTA 在生物信息学中，FASTA 格式是一种用于记录核酸序列或肽序列的文本格式，其中的核酸或氨基酸均以单...
从gff3文件获取fasta序列(2)
xxxxx
09-06 4010
get fasta sequence from gff3#taken from https://www.biostars.org/p/46281/ and modified as needed 10/15/15 use strict; use warnings; use Bio::Seq; use Bio::SeqIO; use Bio::DB::Fasta;#add a help message
EMBLmyGFF3:将gff3批注文件转换为准备提交的EMBL格式的有效方法
05-07
EMBLmyGFF3GFF3到EMBL的转换工具 EMBLmyGFF3将FASTA格式的程序集以及GFF3格式的关联注释转换为，这是将带注释的程序集提交给所需的格式。 [类似地，准备将数据提交给NCBI时，请使用 。 ] ！ NCBI和ENA是一部分，它们...
x2fasta:将常见序列格式重新格式化为Fasta的工具-开源
04-29
一组将常见的生物序列格式（例如EMBL，SWISS-PROT，UniProtKB，GenBank和RefSeq）转换为fasta序列格式的工具。
EMBLRead:EMBLRead 的修改版本，用于更快地将 EMBL 文件读入 Matlab 结构-matlab开发
06-01
这是基于官方版本的 EMBLRead 函数的修改版本，可以将 EMBL 文件读入 Matlab 结构。 修改原因有2个： 1. 正式版读取大文件（>200k，例如“U18997”）的时间长得不可原谅（728.6秒）。 在这个修改版本中，它被极大地...
embl-adapter：包含用于将EMBL内容连接到GBIF的适配器
02-19
最初，此存储库将包含一个EMBL API搜寻器，该搜寻器可以构建适合于提取到GBIF的DwC-A文件。 结果将替换当前的。 该代码库将基于Java，并且可以在GBIF或EMBL基础结构上运行（尚待确定，并且可能会随时间变化）。 ...
pyEntrezId:将Ensembl，Uniprot和HGNC ID转换为Entrez基因ID
05-02
PyEntrezId Python软件包可在各种基因或蛋白质ID之间进行转换 ...将Ensemble转录本基因ID转换为Entrez基因ID from PyEntrezId import Conversion EnsemblId = 'ENST00000407559' # include your email a
根据位置从gff中提取gene
热门推荐
GRC
06-29 1万+
小程序-根据位置从gff中提取gene Table of Contents 1 需要的文件格式2 程序 1 需要的文件格式 1: 染色体 位置 2 程序 1: use strict; 2: use warnings; 3: 4: my (@information,$chr,%hash,$key1,$key2,); 5
gff数据处理
ruby912的博客
06-04 3899
Perl脚本练习 通过数组和哈希引用读取存储复杂数据 要求 gff文件，记录每个基因的名称、起始终止位置、染色体、转录本信息 数据 #注释行 分割：\t 9列：gff格式 思路 构建一个哈希，结构 my%gene = ( geneID => { 'location' => [chrs, start, end] 'transcripts' => [ { ...
Mysql中的join、cross join、inner join是等效的，某些场景下left join被转成inner join可以有更好的性能提升
Saintyyu的博客
08-31 3699
今天在Mysql官网（参考博客1）看到一段话： At the parser stage, queries with right outer join operations are converted to equivalent queries containing only left join operations. In the general case, the conversion is ...
gtf与gff3文件【格式】【转换】
GeekFocus
08-31 7383
GFF3 官方 General Feature Format Version 3 存储序列结构信息的一种数据格式。序列结构就是一个scaffold或者染色体上面每个位置都是什么序列元件。 GFF每一行代表一个序列元件（以#为开头的注释行除外），一行9列9个属性，必须tab分割，属性为空用“.”代替。 1. seqid - scaffold或者chromosome的名称说明 2. source - 产生一个序列元件的软件的名称或者数据源（数据库名称或者项目名称） 3. type - 序列元件的类型，例如：
gff3格式转换
xxxxx
11-19 7094
使用gmap软件输出的gff3文件并不能直接用来作为gffcompare的输入文件 故此，需要进行转换，转换也简单就是讲target字段删除即可。#!/usr/bin/env python # -*- coding: utf-8 -*- __author__ = 'shengwei ma' __author_email__ = 'shengweima@icloud.com'with open('
gff文件_gb格式注释文件转换成gff3注释文件格式
weixin_39788382的博客
11-30 4019
今天在NCBI下载了酵母的参考基因组，没有找到gff格式的基因组注释文件，只找到了genbank格式的基因组注释文件。应该会有现成的工具来实现常用的基因组注释文件不同格式之间的相互转换。比如gtf、gff、和genbank之间的相互转换。经过搜索找到三款工具可以把gb格式文件转换成gff格式注释文件。第一个是 EMBOSS工具中的seqret命令参考https://www.biosta...
亚马逊AWS S3云下载大文件断点续传及重下错误片段的方法
Alter_X的博客
02-17 2116
通过Amazon S3 MultiThread Resume Upload Tool - Single Node实现下载aws文件断点续传，结合aws cli和dd命令处理下载错误的小块片段
非R包的ballgown提取FPKM的方法及R中dplyr包和Reduce函数的使用实例
Alter_X的博客
03-28 1269
R中利用dplyr包和Reduce函数提取stringtie生成的.ctab文件中的FPKM值。
从meme的motif查询结果中获取含有某个motif的基因名称列表（python解析xml）
Alter_X的博客
04-20 1264
python中利用ElementTree包对xml文件的解析
iOS系统有可以编辑fasta格式文件的软件吗
05-31
Sequin支持多种格式的序列文件，包括fasta、Genbank、EMBL等。 2. BioEdit: 这是一个免费的生物信息学软件，可以在Windows和Mac OS X上使用。BioEdit也有一个iOS版本，可以在iPhone和iPad上使用。BioEdit可以编辑...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值