求替換snp位點(diǎn)的腳本

yang7473453

現(xiàn)有兩個(gè)文件，一個(gè)是fasta文件，一個(gè)snp位點(diǎn)信息文件，如下
文件一：
>gh2222
TTTTTTAAAAAAAAAAGGGGGGCCCCCCCCCAAAAAA
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
>gh2242
TTTTTTTTTTTTTTTTTAAAAAAAAAAAAAAGGGGGGG
CCCCCCCCCCCCCTT

文件二：
基因號(hào)       位置      snp堿基
gh2222         5              A
gh2222         22            A
gh2222         9              G
gh2242         15            C

需要根據(jù)文件二中的信息將文件中相應(yīng)的snp位點(diǎn)替換掉，生成新的替換后的文件，求perl腳本！

huang6894

substr......

yang7473453

大神，可以具體一點(diǎn)嗎？回復(fù) 2# huang6894


   

ntwarren

位置是從1開始算還是從0開始算, 文件中:
TTTTTTAAAAAAAAAAGGGGGGCCCCCCCCCAAAAAA
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
只要替換第一行的還是兩行都要?文件的具體格式是不是都是兩行組成的?

問題描述不清,則無法回答.
你要給個(gè)簡(jiǎn)單的例子說明你要想的結(jié)果.

yang7473453

位置從0開始算，序列不管有幾行都把它當(dāng)連續(xù)的一段來看回復(fù) 4# ntwarren


   

yang7473453

給你一個(gè)具體一點(diǎn)的，文件一：
>gg22
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>gg23
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
>gg25
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

文件二：      位置信息是按同一個(gè)序列的所有堿基計(jì)算的
序列名       位置       堿基
gg22             5             T
gg22             15           T
gg25             8             A

需要生成的文件：   文件三
>gg22
AAAATAAAAAAAAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>gg23
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
>gg25
GGGGGGGAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
            
文件一中每個(gè)序列的行數(shù)不一樣，根據(jù)文件二中的信息將文件一中相應(yīng)的位置替換堿基，
沒有替換過的序列就按文件一中的再次輸出一遍，得到替換堿基后的文件三。不知道腳本要怎么寫？

回復(fù) 4# ntwarren


   

zxy877298415

1. awk -v RS=">" -v OFS="\n"  'FNR==NR{for (i=4;i<=NF;i+=3) if($i~/[0-9a-Z]/)a[$i]=$(i+1)"#"$(i+2)"#"a[$i];next}NR>FNR&&FNR>1{
   
2. > if(a[$1]){split(a[$1],b,"#");split($2,c,"");for(j=1;j<=length(c);j++){for(s=1;s<=length(b);s+=2) if(j==b[s]) c[j]=b[s+1]}
   
3. > for (k=1;k<=length(c);k++){t=t c[k];$2=t}t=""}}{if($0!~/^$/) print ">"$0}' file2 file1
   
4. >gg22
   
5. AAAATAAAAAAAAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
   
6. AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
   
7. >gg23
   
8. TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
   
9. TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
   
10. TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
    
11. 
    
12. >gg25
    
13. GGGGGGGAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

復(fù)制代碼

回復(fù) 1# yang7473453 perl不會(huì)，上個(gè)awk的吧！


   

yang7473453

很感謝您的幫助，只是不懂a(chǎn)wk，很抱歉回復(fù) 7# zxy877298415


   

huang6894

1. #!/usr/bin/perl -w
   
2. my($fa, $snp) = @ARGV;
   
3. my %hash;
   
4. open SNP,"< $snp" or die "$!";
   
5. while(<SNP>){
   
6.         chomp;
   
7.         my($chr, $pos, $mut) = split;
   
8.         push @{ $hash{$chr} },{ P => $pos, M => $mut };
   
9. }
   
10. close SNP;
    
11. 
    
12. $/ = '>';
    
13. open FA,"< $fa" or die "$!";
    
14. <FA>;
    
15. while(<FA>){
    
16.         chomp;
    
17.         my ( $CHR, @SEQ ) = split;
    
18.         unless(exists $hash{$CHR}){
    
19.                 print ">".$_;
    
20.                 next;
    
21.         }
    
22.         my $SEQ = join '', @SEQ;
    
23.         for my $C ( @{ $hash{$CHR} } ) {
    
24.                 substr( $SEQ, $C->{P}, 1 ) = $C->{M};
    
25.         }
    
26.         print ">$CHR\n$SEQ\n";
    
27. }
    
28. close FA;
    
29. $/ = '\n';
    
30. 
    

復(fù)制代碼

ntwarren

用這個(gè)吧, 能保留你每一小段的換行風(fēng)格不變.

use strict;
use warnings;

open FILE1, "<", "file1.txt";
open FILE2, "<", "file2.txt";
open FILE3, ">", "file3.txt";

my @file2 = map {chomp; [split /\s+/] } <FILE2>;

our $name;
our $content;

while(<FILE1>){
    if(/^>(\w+)/){
        if(defined($name)){
            for my $a (@file2){
                if($a->[0] eq $name){
                    my @totalLineFeed = substr($content,0,$a->[1]) =~ /\n/g;
                    substr($content,$a->[1]-1 + @totalLineFeed,1)=$a->[2];
                }
            }
            print FILE3 ">".$name."\n".$content;
            
        }
        $name = $1;
        $content = "";
    } else {
        $content .= $_;
    }
}