我有数百个文件,每个文件都有不同数量的条目(> xxxx),并且希望单独保留所有文件之间的共享条目。我不确定什么是最好的方法,也许是Perl!我使用了排序,bash的uniq,但我没有得到正确的答案。ID的格式以>>开头,并在所有文件中遵循4个字符。
1.fa
>abcd
CTGAATGCC
2.fa
>abcd
AAATGCGCG
>efgh
CGTAC
3.fa
>abcd
ATGCAATA
>efgh
TAACGTAA
>ijkl
TGCAA
最终结果,此示例的结果是:
1.fa
>abcd
CTGAATGCC
2.fa
>abcd
AAATGCGCG
3.fa
>abcd
ATGCAATA
此Perl程序将按照您的要求进行。它使用perl的内置编辑 功能,并将原始文件重命名为 1.fa.bak
等。只要序列始终立即在一行上,数据中的空白行就不会有问题遵循ID
use strict;
use warnings 'all';
my @files = glob '*.fa';
printf "Processing %d file%sn", scalar @files, @files == 1 ? "" : "s";
exit if @files < 2;
my %ids;
{
local @ARGV = @files;
while ( <> ) {
++$ids{$1} if /^>(S+)/;
}
}
# remove keys that aren't in all files
delete @ids{ grep { $ids{$_} < @files } keys %ids };
my $n = keys %ids;
printf "%d ID%s common to all filesn", $n, $n == 1 ? '' : "s";
exit unless $n;
{
local @ARGV = @files;
local $^I = '.bak';
while ( <> ) {
next unless /^>(S+)/ and $ids{$1};
print;
print scalar <>;
}
}
这是Perl解决方案,可能会对您有所帮助:
use feature qw(say);
use strict;
use warnings;
my $file_dir = 'files';
chdir $file_dir;
my @files = <*.fa>;
my $num_files = scalar @files;
my %ids;
for my $file (@files) {
open ( my $fh, '<', $file) or die "Could not open file '$file': $!";
while (my $id = <$fh>) {
chomp $id;
chomp (my $sequence = <$fh>);
$ids{$id}++;
}
close $fh;
}
for my $file (@files) {
open ( my $fh, '<', $file) or die "Could not open file '$file': $!";
my $new_name = $file . '.new';
open ( my $fh_write, '>', $new_name ) or die "Could not open file '$new_name': $!";
while (my $id = <$fh>) {
chomp $id;
chomp (my $sequence = <$fh>);
if ( $ids{$id} == $num_files ) {
say $fh_write $id;
say $fh_write $sequence;
}
}
close $fh_write;
close $fh;
}
假设所有.fa
文件都位于名为$file_dir
的目录中,并且将新序列写入同一目录中的新文件。新文件名获取.new
扩展名。