2014-10-22 50 views
1

我正在去除格式化爲這樣的2個文件匹配的數據:多列比較perl的

scaffold_1 21786 . A G 198 . DP=44;VDB=0.0402;AF1=1 

我所試圖做的是檢查這2個文件相互爲1列5(在這種情況下爲scaffold_1G)。如果所有這些匹配,那麼我想跳過這些行,但是如果我在這5個內部有不匹配的地方,那麼他們會去他們自己的寫入文件。

#!/usr/bin/perl 
#bothfixed.pl 
use strict; use warnings; 

die "usage <file1> <file2> <write1> <write2>\n" unless @ARGV ==4; 

open (my $file1, "<", "$ARGV[0]") or die "Can't open $file1:$!"; 
open (my $file2, "<", "$ARGV[1]") or die "Can't open $file1:$!"; 
open (my $write1, ">", "$ARGV[2]"); 
open (my $write2, ">", "$ARGV[3]"); 


my $line1; 
my $line2; 
my @array1; 
my @array2; 
while (not eof $file1 && not eof $file2) { 
    $line1 = <$file1> ;$line2 = <$file2>; #print $line1; 

    @array1 = split (/\t/, $line1); 
    my @slice1 = @array1[0..4]; 

    @array2 = split (/\t/, $line2); 
    my @slice2 = @array2[0..4]; 

    my $string1 = join ('',@slice1); 
    my $string2 = join ('',@slice2); 

    if ($string1 eq $string2) {next} 
    else {print $write1 "$line1" ; print $write2 "$line2"} 
} 

我只看到一條線在此過程中去除,而且我認爲,這意味着什麼毛病我while循環,但這就是好猜測,因爲我有。

任何想法?

我已經添加了一些測試數據。從文件1.

scaffold_1 721 . T C 222 . DP=67;VDB=0.0411;AF1=1;AC1=2;DP4=0,0,30,33;MQ=60;FQ=-217 GT:PL:DP:GQ 1/1:255,190,0:63:99 

scaffold_1 1282 . T G 36 . DP=67;VDB=0.0396;AF1=0.5;AC1=1;DP4=23,23,15,5;MQ=26;FQ=39;PV4=0.1,0.097,0.039,1 GT:PL:DP:GQ 0/1:66,0,255:66:69 

scaffold_1 15917 . AATATATATATATATATATATATATATATATATATAT AATATATATATATATATATATATATATATATATATATAT 127 . INDEL;DP=71;VDB=0.0387;AF1=0.5;AC1=1;DP4=3,13,5,6;MQ=57;FQ=130;PV4=0.21,1,0.21,1 GT:PL:DP:GQ 0/1:165,0,255:27:99 

scaffold_1 19183 . TAC TACAC 217 . INDEL;DP=83;VDB=0.0408;AF1=0.5;AC1=1;DP4=24,18,16,19;MQ=60;FQ=217;PV4=0.36,1,1,0.0074 GT:PL:DP:GQ 0/1:255,0,255:77:99 

scaffold_1 21786 . A G 198 . DP=44;VDB=0.0402;AF1=1;AC1=2;DP4=0,0,21,22;MQ=60;FQ=-156 GT:PL:DP:GQ 1/1:231,129,0:43:99 

scaffold_1 26031 . G A 169 . DP=83;VDB=0.0263;AF1=0.5;AC1=1;DP4=23,14,21,24;MQ=60;FQ=172;PV4=0.19,6.7e-28,1,1 GT:PL:DP:GQ 0/1:199,0,255:82:99 

scaffold_1 33033 . A T 206 . DP=61;VDB=0.0411;AF1=0.5;AC1=1;DP4=17,22,13,8;MQ=60;FQ=209;PV4=0.28,3.4e-05,1,1 GT:PL:DP:GQ 0/1:236,0,255:60:99 

scaffold_1 33799 . C A 146 . DP=56;VDB=0.0394;AF1=0.5;AC1=1;DP4=13,14,13,14;MQ=60;FQ=149;PV4=1,6.2e-28,0.16,0.14 GT:PL:DP:GQ 0/1:176,0,255:54:99 

scaffold_1 35051 . CAAAAAAAAAA CAAAAAAAAAAA 32.5 . INDEL;DP=51;VDB=0.0447;AF1=1;AC1=2;DP4=1,1,22,13;MQ=60;FQ=-118;PV4=1,1,1,1 GT:PL:DP:GQ 1/1:73,83,0:37:99 

線從文件2

scaffold_1 721 . T C 221 . DP=57;VDB=0.0407;AF1=1;AC1=2;DP4=0,0,23,32;MQ=60;FQ=-193 GT:PL:DP:GQ 1/1:254,166,0:55:99 

scaffold_1 1282 . T G 80 . DP=82;VDB=0.0383;AF1=0.5;AC1=1;DP4=29,30,13,9;MQ=26;FQ=83;PV4=0.46,0.19,0.5,1 GT:PL:DP:GQ 0/1:110,0,238:81:99 

scaffold_1 10472 . A C 23 . DP=44;VDB=0.0402;AF1=0.5;AC1=1;DP4=6,8,1,11;MQ=60;FQ=26;PV4=0.081,9.7e-12,1,1 GT:PL:DP:GQ 0/1:53,0,246:26:56 

scaffold_1 15917 . AATATATATATATATATATATATATATATATATATAT AATATATATATATATATATATATATATATATATATATAT 186 . INDEL;DP=39;VDB=0.0416;AF1=0.5;AC1=1;DP4=5,6,5,5;MQ=60;FQ=189;PV4=1,1,1,0.43 GT:PL:DP:GQ 0/1:224,0,237:21:99 

scaffold_1 19183 . TAC TACAC 217 . INDEL;DP=76;VDB=0.0383;AF1=0.5;AC1=1;DP4=16,12,20,23;MQ=60;FQ=217;PV4=0.47,1,1,0.11 GT:PL:DP:GQ 0/1:255,0,255:71:99 
scaffold_1 21786 . A G 196 . DP=58;VDB=0.0365;AF1=1;AC1=2;DP4=0,0,33,24;MQ=60;FQ=-199 GT:PL:DP:GQ 1/1:229,172,0:57:99 

scaffold_1 26031 . G A 169 . DP=70;VDB=0.0407;AF1=0.5;AC1=1;DP4=13,12,22,23;MQ=60;FQ=172;PV4=1,6.1e-26,1,1 GT:PL:DP:GQ 0/1:199,0,255:70:99 

scaffold_1 33033 . A T 225 . DP=41;VDB=0.0404;AF1=0.5;AC1=1;DP4=8,9,13,10;MQ=60;FQ=225;PV4=0.75,0.0052,1,1 GT:PL:DP:GQ 0/1:255,0,255:40:99 

scaffold_1 33799 . C A 116 . DP=61;VDB=0.0410;AF1=0.5;AC1=1;DP4=18,20,11,12;MQ=60;FQ=119;PV4=1,3.6e-28,1,0.46 GT:PL:DP:GQ 0/1:146,0,255:61:99 
scaffold_1 35051 . CAAAAAAAAAA CAAAAAAAAAAA,CAAAAAAAAAAAA 47.5 . INDEL;DP=40;VDB=0.0384;AF1=1;AC1=2;DP4=0,0,16,17;MQ=60;FQ=-128 GT:PL:DP:GQ 1/1:88,93,0,82,60,76:33:99 
+0

如果所有5次比較必須爲真,每次檢查後都無法打印,則必須立即進行所有測試。如果你正在尋找字符串相等,你可能會捕獲一個字符串中的前5列並進行比較。 – TLP 2014-10-22 09:45:44

+0

我編輯了代碼來匹配你所描述的(我認爲),但我仍然遇到同樣的問題。 – aupadhyaya 2014-10-22 09:56:00

+0

爲什麼你放棄其他每一行? (在最後一段) – choroba 2014-10-22 10:15:57

回答

0

這應該做你想要什麼。你需要將第一個文件寫入內存並將其存儲在一個用$string1比較字符串索引的散列中,然後可以檢查第二個文件中的每一行是否我們曾經見過。這段代碼會將第一個文件中的輸出行混洗;如果您需要保持它們的順序,您可以添加另一個散列來保存行號,並使用它對最終輸出循環進行排序。

#!/usr/bin/perl 
#bothfixed.pl 
use strict; use warnings; 

die "usage <file1> <file2> <write1> <write2>\n" unless @ARGV ==4; 

open (my $file1, "<", "$ARGV[0]") or die "Can't open $file1:$!"; 
open (my $file2, "<", "$ARGV[1]") or die "Can't open $file1:$!"; 
open (my $write1, ">", "$ARGV[2]"); 
open (my $write2, ">", "$ARGV[3]"); 


my $line1; 
my $line2; 
my @array1; 
my @array2; 
my %file1; 
my %match; 

# Slurp first file into memory 
while ($line1 = <file1>) { 
    @array1 = split (/\t/, $line1); 
    my @slice1 = @array1[0..4]; 
    my $string1 = join ('',@slice1); 
    $file1{$string1} = $line1; 
} 

# Run through second file 
while ($line2 = <$file2>) { 
    @array2 = split (/\t/, $line2); 
    my @slice2 = @array2[0..4]; 
    my $string2 = join ('',@slice2); 

    if ($file1{$string2}) { 
     # If we've matched then we wont want to print it out at the end 
     $match{$string2}++; 
    } else { 
     print $write2 "$line2"; 
    } 
} 

foreach my $string1 (keys %file1) { 
    next if $match{$string1}; 
    print $write1 "$line1"; 
}