2013-02-04 66 views
1

我想修改其設置這樣的文件:將列添加到基於現有列的文件

chr start ref alt 
chr1 18884 C CAAAA 
chr1 135419 TATACA T 
chr1 332045 T TTG 
chr1 453838 T TAC 
chr1 567652 T TG 
chr1 602541 TTTA T 
chr1 614937 C CTCTCTG 
chr1 654889 C CA 
chr1 736800 AC A 

我想修改它使得: 如果列「參考」是一個字符串> 1(即第2行),那麼我生成2個新的列,其中:

第一新列=起始座標-1 第二新列=起始座標+(參考文獻串的長度)+1

因此,對於第2行輸出將如下所示:

chr1 135419 TATACA T 135418 135426 

或: 如果在 「REF」= 1和列 「ALT」=長度的字符串的字符串的長度> 1(即線1)然後

第一新列=開始座標 第二新列=起始座標+ 2

所以,對於第1行輸出將是:

chr1 18884 C CAAAA 18884 18886 

我在AWK此次嘗試但沒有成功 我的Perl是不存在的,但這是最好的方式嗎?或者也許在R?

+0

如果我讀正確的規格做什麼,你的兩個如果崩潰的一個(如1 + 1 = 2),除非你確實有在案件哪些ALT可能會丟失。如果是這樣(可能會錯過),那麼你可能有一個固定的記錄字段輸入和規格缺失 - 因爲下面的每個解決方案分裂在白色空間。 –

回答

2

Perl解決方案。請注意,您的規格中沒有提到,如果兩個字符串長度爲1

#!/usr/bin/perl 
use warnings; 
use strict; 
use feature qw(say); 

#use Data::Dumper; 
<DATA>; # Skip the header; 
while (<DATA>) { 
    my ($chr, $start, $ref, $alt) = split; 
    my @cols; 
    if (1 < length $ref) { 
      @cols = ($start - 1, $start + 1 + length $ref); 
    } elsif (1 < length $alt) { 
     @cols = ($start, $start + 2); 
    } else { 
     warn "Don't know what to do at $.\n"; 
    } 
    say join "\t", $chr, $start, $ref, $alt, @cols; 
} 


__DATA__ 
chr start ref alt 
chr1 18884 C CAAAA 
chr1 135419 TATACA T 
chr1 332045 T TTG 
chr1 453838 T TAC 
chr1 567652 T TG 
chr1 602541 TTTA T 
chr1 614937 C CTCTCTG 
chr1 654889 C CA 
chr1 736800 AC A 
2

下面是使用awk的一種方法。的script.awk

awk -f script.awk file | column -t 

內容:

NR==1 { 
    next 
} 

length($3)>1 && length($4)==1 { 
    print $0, $2-1, $2+length($3)+1 
    next 
} 

length($3)==1 && length($4)>1 { 
    print $0, $2, $2+2 
    next 
}1 

結果:像運行

chr1 18884 C  CAAAA 18884 18886 
chr1 135419 TATACA T  135418 135426 
chr1 332045 T  TTG  332045 332047 
chr1 453838 T  TAC  453838 453840 
chr1 567652 T  TG  567652 567654 
chr1 602541 TTTA T  602540 602546 
chr1 614937 C  CTCTCTG 614937 614939 
chr1 654889 C  CA  654889 654891 
chr1 736800 AC  A  736799 736803 

另外,這裏是一個班輪:

awk 'NR==1 { next } length($3)>1 && length($4)==1 { print $0, $2-1, $2+length($3)+1; next } length($3)==1 && length($4)>1 { print $0, $2, $2+2; next }1' filem | column -t 

代碼應該很明顯。腳本末尾的1只是啓用每行的默認打印(即'1'返回true)。 HTH。

0

在Perl這樣做實在是小巫見大巫(但這樣是AWK):

#!/usr/bin/perl 
while (<>) { 
    chmop; 
    my ($chr,$start,$ref,$alt)=split(/\s+/,$_); 
    if (len($ref) > 1) { 
print STDOUT 
    "$chr\t$start\t$ref\t$alt\t", 
    $start+len($ref)+1,"\n"; 
    } elsif (len($ref)==1) { 
print STDOUT 
    "$chr\t$start\t$ref\t$alt\t", 
    $start+2,"\n"; 
    } else { 
print STDERR "ERROR: ???\n"; #actually impossible 
    } 
} 

堅持在文件morecols.pl,使用chmod + X morecols.pl,運行更morecols.pl。 (請小心,本代碼/說明中有很多假設)。我有一種感覺,你的實際問題更多的是編程/文本處理,然後是工具或語言。如果是這樣,這個代碼只是一個權宜之計解決方案....

乾杯。