perl - 使用 snp 位置修改序列并在同一文件中输出

Question

我有两个文件，一个是位置信息，另一个是序列信息。现在我需要读取位置并在位置获取snps，然后用序列中的snp信息替换该位置基数并将其写入snp信息文件中..例如

SNP 文件包含

10 A C A/C

序列文件包含

ATCGAACTCTACATTAC

这里第 10 个元素是 T 所以我将用 [A/C] 替换 T 所以最终输出应该是

10 A C A/C ATCGAACTC[A/C]ACATTAC

示例文件是

Snp 文件

SNP Ref Alt
10  A   C
19  G   C
30  C   T
42  A   G

序列 :

sequence1 CTAGAATCAAAGCAAGAATACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT

最终输出:

SNP Ref Alt Output
10  A   C   CTAGAATCA[A/C]AGCAAGAANACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT
19  G   C   CTAGAATCANAGCAAGAA[C/G]ACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT
30  C   T   CTAGAATCAAAGCAAGAATACACTCTTTT[T/C]TTTGGAAAAGAATATCTCATGTTTGCTCTT
42  A   G   CTAGAATCAAAGCAAGAATACACTCTTTTTTTTGGAAAAGA[A/G]TATCTCATGTTTGCTCTT

在从 Ref 和 A​​lt 列替换这里的 snps 时，我们需要考虑 {A,T,C,G} 的顺序，如 [Ref/Alt] 总是第一个碱基应该是 A 或 T 或 C，然后是命令。

另一件事是，如果我们取snp 位置，如果有10 个碱基差异的snp，我们需要用“N”替换那个snp 位置。在上面的示例中，前两个位置的差为 9，我们将另一个元素替换为“N”。

我已经编写了代码，它按顺序打印位置，并用该 snp 位置替换序列，但无法读取附近的位置并用 N 替换。

我的方法可能完全错误，因为我是编码初学者。我认为通过使用散列，我们可以轻松实现这一点，但我对散列不太熟悉..请帮忙提供一些建议..我不必只坚持使用 perl ,

my $input_file = $ARGV[0];
my $snp_file = $ARGV[1];
my $output_file = $ARGV[2];

%sequence_hash = ();

open SNP, $snp_file || die $!;
$indel_count = 0;
$snp_count = 0;
$total_count = 0;

#### hashes and array
@id_array = ();

while ($line = <SNP>) {

    next if $line =~ /^No/;
    $line =~ s/\r|\n//g;


   # if ($line =~ /\tINDEL/) {

    #    $indel_count++;
     #   $snp_type = "INDEL";

    #} else {
     #   $snp_count++;
      #  $snp_type = "SNP";
    #}

    @parts = split(/\t/,$line);

    $id = $parts[0];
    $pos = $parts[1];
    #$ref_base = $parts[3];
    @temp_ref = split(",",$parts[2]);
    $ref_base = $temp_ref[0];
    @alt = split(",",$parts[3]);
    $alt_base = $alt[0];
    $snpformat = '';

    if ($ref_base eq "A" || $alt_base eq "A")
    {

        if ($ref_base eq "A"){
            $snpformat = "[".join("/",$ref_base,$alt_base)."]";}
            #$snpformat = $ref_base/$alt_base;}
            #print "refbase is A $ref_base \t $alt_base \t $snpformat \n"; }

            else 
            {$snpformat = "[".join("/",$alt_base,$ref_base)."]";}
            #print "Altbase is A $ref_base \t $alt_base \t $snpformat \n";}
    }

    elsif ($ref_base eq "T" || $alt_base eq "T")
    {

        if ($ref_base eq "T"){
            $snpformat = "[".join("/",$ref_base,$alt_base)."]";}
            #$snpformat = $ref_base/$alt_base;}
            #print "refbase is A $ref_base \t $alt_base \t $snpformat \n"; }

            else 
            {$snpformat = "[".join("/",$alt_base,$ref_base)."]";}
            #print "Altbase is A $ref_base \t $alt_base \t $snpformat \n";}
    }

    elsif ($ref_base eq "C" || $alt_base eq "C")
    {

        if ($ref_base eq "C"){
            $snpformat = "[".join("/",$ref_base,$alt_base)."]";}
            #$snpformat = $ref_base/$alt_base;}
            #print "refbase is A $ref_base \t $alt_base \t $snpformat \n"; }

            else 
            {$snpformat = "[".join("/",$alt_base,$ref_base)."]";}
            #print "Altbase is A $ref_base \t $alt_base \t $snpformat \n";}
    }


    else 
    {$snpformat = "-/-" ;}
    print " $id \t $pos \t $ref_base \t $alt_base \t $snpformat \n  ";
}

open SEQ, $input_file ||die $!;

$header = '';
$sequence = '';
$num_sequences = 0;

while ($line = <SEQ>) {

    $line =~ s/\r|\n//g;
    next if $line =~ //;

    if ($line =~ /^>(.+)/) {
        if ($header eq '') {

            $header = $1;
            $sequence = '';
            next;
        } else {

            $sequence_len = length($sequence);

            $sequence_hash{$header} = $sequence;
            push (@headers,$header);
            #print $header."\t".$sequence_len."\n";
            #print $sequence."\n";
            $num_sequences++;

            $header = $1;
            $sequence = '';

        }


    } else {

        $sequence .= $line;

    }

}
$sequence_len = length($sequence);
$sequence_hash{$header} = $sequence;
push (@headers,$header);
#print $header."\t".$sequence_len."\n";

$num_sequences++;

close (SEQ);

$pos = '4';
substr($sequence,$pos,1) = "[A/G]";
print $sequence."\n";   
print "$pos \n";

Answer 1

此 awk脚本可能可以帮助您获得所需的结果。

awk '
BEGIN {
print "SNP\tRef\tAlt\tOutput"
}
NR==FNR { 
    a[++i]=$0
    next 
} 
FNR>1 { 
    x=substr(a[i],1,($1-1))
    z=substr(a[i],($1+1))
    if ($2=="A") {
        y="["$2"/"$3"]"
    } 
    else if ($2=="T" && $3=="A") {
        y="["$3"/"$2"]"
    }
    else if ($2=="C" && ($2=="A" || $2=="T")) {
        y="["$3"/"$2"]"
    }
    else if ($2=="G" && ($2=="A" || $2=="T" || $2=="C")) {
        y="["$3"/"$2"]"
    }
    else 
        y="["$3"/"$2"]"
    print $1"\t"$2"\t"$3"\t"x""y""z
}' sequence snp

测试:

[jaypal:~/temp] cat sequence
CTAGAATCAAAGCAAGAATACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT

[jaypal:~/temp] cat snp
SNP Ref Alt
10  A   C
19  G   C
30  C   T
42  A   G

[jaypal:~/temp] awk '
BEGIN {
print "SNP\tRef\tAlt\tOutput"
}
NR==FNR { 
    a[++i]=$0
    next 
} 
FNR>1 { 
    x=substr(a[i],1,($1-1))
    z=substr(a[i],($1+1))
    if ($2=="A") {
        y="["$2"/"$3"]"
    } 
    else if ($2=="T" && $3=="A") {
        y="["$3"/"$2"]"
    }
    else if ($2=="C" && ($2=="A" || $2=="T")) {
        y="["$3"/"$2"]"
    }
    else if ($2=="G" && ($2=="A" || $2=="T" || $2=="C")) {
        y="["$3"/"$2"]"
    }
    else 
        y="["$3"/"$2"]"
        print $1"\t"$2"\t"$3"\t"x""y""z
}' sequence snp
SNP Ref Alt Output
10  A   C   CTAGAATCA[A/C]AGCAAGAATACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT
19  G   C   CTAGAATCAAAGCAAGAA[C/G]ACACTCTTTTTTTTGGAAAAGAATATCTCATGTTTGCTCTT
30  C   T   CTAGAATCAAAGCAAGAATACACTCTTTT[T/C]TTTGGAAAAGAATATCTCATGTTTGCTCTT
42  A   G   CTAGAATCAAAGCAAGAATACACTCTTTTTTTTGGAAAAGA[A/G]TATCTCATGTTTGCTCTT