求替换snp位点的脚本

yang7473453

现有两个文件，一个是fasta文件，一个snp位点信息文件，如下
文件一：
>gh2222
TTTTTTAAAAAAAAAAGGGGGGCCCCCCCCCAAAAAA
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
>gh2242
TTTTTTTTTTTTTTTTTAAAAAAAAAAAAAAGGGGGGG
CCCCCCCCCCCCCTT

文件二：
基因号       位置      snp碱基
gh2222         5              A
gh2222         22            A
gh2222         9              G
gh2242         15            C

需要根据文件二中的信息将文件中相应的snp位点替换掉，生成新的替换后的文件，求perl脚本！

huang6894

substr......

yang7473453

大神，可以具体一点吗？回复 2# huang6894


   

ntwarren

位置是从1开始算还是从0开始算, 文件中:
TTTTTTAAAAAAAAAAGGGGGGCCCCCCCCCAAAAAA
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
只要替换第一行的还是两行都要?文件的具体格式是不是都是两行组成的?

问题描述不清,则无法回答.
你要给个简单的例子说明你要想的结果.

yang7473453

位置从0开始算，序列不管有几行都把它当连续的一段来看回复 4# ntwarren


   

yang7473453

给你一个具体一点的，文件一：
>gg22
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>gg23
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
>gg25
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

文件二：      位置信息是按同一个序列的所有碱基计算的
序列名       位置       碱基
gg22             5             T
gg22             15           T
gg25             8             A

需要生成的文件：   文件三
>gg22
AAAATAAAAAAAAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>gg23
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
>gg25
GGGGGGGAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
            
文件一中每个序列的行数不一样，根据文件二中的信息将文件一中相应的位置替换碱基，
没有替换过的序列就按文件一中的再次输出一遍，得到替换碱基后的文件三。不知道脚本要怎么写？

回复 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

复制代码

回复 1# yang7473453 perl不会，上个awk的吧！


   

yang7473453

很感谢您的帮助，只是不懂awk，很抱歉回复 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. 
    

复制代码

ntwarren

用这个吧, 能保留你每一小段的换行风格不变.

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 .= $_;
    }
}