Perl拆分FASTA文件

将FASTA文件按照用户指定的序列条数进行拆分,即每n条序列写到一个文件中,或者按照指定的输出文件数进行拆分,则每个文件中包含的序列数相同

脚本名:splitFasta.pl

#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;
use POSIX;

# 帮助文档
=head1 Description

	This script is used to split fasta file, which is too large with thosands of sequence

=head1 Usage

	$0 -i <input> -o <output_dir> [-n <seq_num_per_file>] [-m <output_file_num>]
	
=head1 Parameters

	-i	[str]	Input raw fasta file
	-o	[str]	Output file to which directory
	-n	[int]	Sequence number per file, alternate chose paramerter "-n" or "-m", if set "-n" and "-m" at the same time, only take "-n" parameter
	-m	[int]	Output file number (default:100)
=cut

my ($input,$output_dir,$seq_num,$file_num);
GetOptions(
	"i:s"=>\$input,
	"o:s"=>\$output_dir,
	"n:i"=>\$seq_num,
	"m:i"=>\$file_num
	);

die `pod2text $0` if ((!$input) or (!$output_dir));

# 设置每个文件的序列条数
if(!defined($seq_num)){
	if(!defined($file_num)){
		$file_num=100;
		my $total_seq_num=`awk 'BEGIN{n=0} /^>/{n++} END{print n}' $input`;
		chomp $total_seq_num;
		$seq_num=ceil($total_seq_num/$file_num);
	}else{
		my $total_seq_num=`awk 'BEGIN{n=0} /^>/{n++} END{print n}' $input`;
		chomp $total_seq_num;
		$seq_num=ceil($total_seq_num/$file_num);
	}
}

open IN,"<$input" or die "Cann't open $input\n";

my $n_seq=0;	# 该变量用于记录当前扫描到的序列数
my $n_file=1;	# 该变量用于记录当前真正写入的文件的计数
my $input_base=`basename $input`;
chomp $input_base;

open OUT,">$output_dir/${input_base}_${n_file}" or die "Cann't create $output_dir/${input_base}_${n_file}\n";

while(<IN>){
	next if (/^\s+$/);	# 跳过空行
	chomp;
	if (/^>/){
		$n_seq++;
		# 判断目前已经扫描到的序列数,若大于设定的split的序列数,则创建新文件
		if ($n_seq>$seq_num){
			$n_seq=1;
			$n_file++;
			close OUT;
			open OUT,">$output_dir/${input_base}_${n_file}" or die "Cann't create $output_dir/${input_base}_${n_file}\n";
			print OUT "$_\n";
		}else{
			print OUT "$_\n";
		}
	}else{
		print OUT "$_\n";
	}
}

close IN;

执行方法:

$ perl splitFasta.pl -i <input> -o <output_dir> [-n <seq_num_per_file>] [-m <output_file_num>]

具体的参数使用说明可以执行 perl splitFasta.pl,查看脚本使用文档

如若转载,请注明出处:https://www.ouq.net/perl%e6%8b%86%e5%88%86fasta%e6%96%87%e4%bb%b6.html

(0)
打赏 微信打赏,为服务器增加100M流量 微信打赏,为服务器增加100M流量 支付宝打赏,为服务器增加100M流量 支付宝打赏,为服务器增加100M流量
上一篇 2022年3月7日 下午11:55
下一篇 2022年3月9日 上午12:12

相关推荐

  • R:使用R连接数据库处理数据

    1.数据库连接 library(DBI) library(dplyr) library(dbplyr) library(odbc) con <- dbConnect(odbc…

    2022年3月6日
    41
  • R:gene symbol与id转换

    利用R将gene symbol与id转换 安装所需包 if (!requireNamespace(“BiocManager”, quietly = TRUE…

    2020年3月16日
    703
  • ggplot2绘图教程

    直方图:geom_histogram 分组直方图的排列方式 position=”dodge” ggplot(small)+geom_histogram(aes(x=price, f…

    2022年3月7日
    83
  • FunRich:功能丰富的分析工具

    FunRich:功能丰富的分析工具http://www.funrich.org/ FunRich 是一个独立的软件工具,主要用于基因和蛋白质的功能富集和相互作用网络分析。此外,分析…

    2022年3月21日
    102
  • CRISPR 设计和Off Target预测:tefor

    1.在线网站 http://crispor.tefor.net/ 2.界面简洁,step1输入KO的序列,step2/3中选择相应的选项即可获得信息。  

    2021年7月15日
    301