2017-04-10 60 views
2

我有大的製表符分隔的文件,如下面的例子:如何優化獨特行的搜索?

scaffold1443 182629 182998 chr1.1.1.241051.241420 367  99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462  99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330  99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322  99.40 
scaffold526  2870014 2870523 chr1.1.5.488372.488881 507  99.90 
scaffold526  2865956 2866314 chr1.1.6.490869.491234 357  98.10 
scaffold526  2867666 2868024 chr1.1.6.490869.491234 357  98.10 
scaffold526  2485557 2485867 chr1.1.7.610677.610987 310  100.00 

我想在一個新的文件只有行的第4列是唯一的打印。 在前面的示例中,除了第4列中包含「chr1.1.6.490869.491234」的兩行外,應打印所有行。

我編寫的以下腳本(它是較大管道的一部分)完成這項工作,但速度非常慢,尤其是當輸入文件非常大時。

#!/usr/bin/perl 
use strict; 
use warnings; 
#This script takes the best hits output and finds the unique elements that up to only one scaffold. 

my $target = $ARGV[0]; 
my $chromosome = $ARGV[1]; 
my @mykeys = `cat OUTPUT_$target/psl_score_byname_$target/$chromosome.table| awk '{print \$4}'| sort -u`; 

foreach (@mykeys) 
{ 
    my $key = $_; 
    chomp($key); 
    my $command = "cat OUTPUT_$target/psl_score_byname_$target/$chromosome.table|grep -w $key"; 
    my @belongs= `$command`; 
    chomp(@belongs); 
    my $count = scalar(@belongs); 
    if ($count == 1) 
    { 
      open FILE, ">>OUTPUT_$target/unique_hces_$target/$chromosome.txt" or die $!; 
      print FILE "@belongs\n"; 

      @belongs =(); 
    } 
    else { 

      @belongs =(); 
    } 
} 

有沒有更智能,更快捷的方法來做到這一點? 非常感謝您提前。

+0

重複拍攝哪個項目有重要嗎? –

+0

因爲您必須掃描整個文件,所以在這裏排序文件似乎並不需要,您可以選擇將第一個或最後一個項目放入一組重複項中。 –

+0

不,在這一點上,我想避免所有的重複。在前面的示例中,我不想保留包含chr1.1.6.490869.491234 – Vasilis

回答

1

鑑於您不想打印具有重複項目的行,您需要在打印之前查看整個文件,以首先找到具有重複項的行。然後返回並打印其他人。

這可以通過將整個文件與輔助數據結構一起保存在內存中,或者通過兩次傳遞來完成。由於該文件是「很大」這裏是更少的內存,使勁方式

use warnings; 
use strict; 

my $file = 'skip.txt'; 
open my $fh, '<', $file or die "Can't open $file: $!"; 

my (%seen, %dupe); 
while (<$fh>) 
{ 
    my $patt = (split)[3]; 

    # Record line numbers if the 4th field has been seen 
    if (exists $seen{$patt}) { 
     $dupe{ $seen{$patt} }++; # num of line with it seen first, with count 
     $dupe{$.} = 1;    # this line's number as well 
    } 
    else { $seen{$patt} = $. }  # first time this 4th field is seen 
} 

# Now we know all lines which carry duplicate fourth field 
my $outfile = 'filtered_' . $file; 
open my $fh_out, '>', $outfile or die "Can't open $outfile: $!"; 

seek $fh, 0, 0; # rewind to the beginning 
$. = 0;   # seek doesn't reset $. 
while (<$fh>) { 
    print $fh_out $_ if not exists $dupe{$.} 
} 
close $fh_out; 

重複發現其原線也需要被記錄,$dupe{$seen{$patt}}++,在該分支的第一次。這個需求只需要完成一次,而我們可以檢查(是否已經被記錄),我們可能會選擇一個可能有用的重複計數。

我已經添加了一些更多的副本(一些超過兩倍)到您的發佈示例,這會產生正確的輸出。


評論上張貼的代碼

張貼的代碼檢查關於對整個文件中的每一行中的第四字段中,從而處理該文件作爲有線多次。這是很多工作,需要花時間,特別是對於大文件。

此外,沒有理由爲該作業使用外部程序。

+0

的行,我只是將它們應用於我的管道,看起來它工作正常。還有一些測試要做,並讓你知道。非常感謝你。 – Vasilis

1

由於oneliner:

perl -F"\t" -lanE 'push @l,[@F];$s{$F[3]}++}{say join"\t",@$_ for grep{$s{$_->[3]}==1}@l' <<EOF 
scaffold1443 182629 182998 chr1.1.1.241051.241420 367 99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462 99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330 99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322 99.40 
scaffold526 2870014 2870523 chr1.1.5.488372.488881 507 99.90 
scaffold526 2865956 2866314 chr1.1.6.490869.491234 357 98.10 
scaffold526 2867666 2868024 chr1.1.6.490869.491234 357 98.10 
scaffold526 2485557 2485867 chr1.1.7.610677.610987 310 100.00 
EOF 

輸出

scaffold1443 182629 182998 chr1.1.1.241051.241420 367 99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462 99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330 99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322 99.40 
scaffold526 2870014 2870523 chr1.1.5.488372.488881 507 99.90 
scaffold526 2485557 2485867 chr1.1.7.610677.610987 310 100.00 

更具可讀性:

perl -F"\t" -lanE ' 
    push @lines, [ @F ]; $seen{ $F[3] }++; 
    END { 
     say join("\t",@$_) for grep { $seen{ $_->[3] } == 1 } @lines 
    } 
' 

你可以,如果願意,我創造了這個作爲oneliner因爲你說其翻譯爲完整的腳本:它是更大的管道的一部分。

另外需要注意的是,上面的內容首先將整個文件讀入內存 - 所以非常大的文件可能會導致問題。

+0

Downvoter,我該如何/應該改進答案? – jm666

+2

不知道。 。 。 。 。 – ikegami

+0

我會嘗試在完整的腳本中翻譯它,因爲它會更適合我的管道。謝謝。 – Vasilis

1

簡單的方法包括使用關聯數組來識別重複項。

perl -F'\t' -lane' 
    push @{ $h{ $F[3] } }, $_; 
    END { 
     for (values(%h)) { 
      print(@$_) if @$_ == 1; 
     } 
    } 
' file.tsv 

上述方法需要儘可能多的內存,因爲文件很大。如果你的文件真的很大,那麼這是一個不行。


如果你有真正的大文件,簡單的辦法是使用sort命令行實用程序對文件進行排序(這是相當快,能夠處理任意大的文件)。通過首先重新排列文件以使重複文件彼此相鄰,我們可以輕鬆過濾出重複文件,而不用擔心內存問題。

sort -t$'\t' -k 4,4 file.tsv | perl -F'\t' -lane' 
    if ($key ne $F[3]) { 
     print(@buf) if @buf == 1; 
     @buf =(); 
    } 

    $key = $F[3]; 
    push @buf, $_; 

    END { print(@buf) if @buf == 1; } 
' 

如果你有真正的大文件,另一種比較簡單的辦法是加載在數據庫(例如一個sqlite3的數據庫)中的數據。您可以使用這種方法輕鬆維護原始訂單。