perl应用：SNP的提取（2）：从对比序列中找到SNP位点并输出 a.pl

我们从第一步中得到的是一个如下的文件：

##maf version=1 scoring=lastz.v1.02.00
# lastz.v1.02.00 --ambiguous=iupac --notransition --step=20 --nogapped --format=maf 
#
# hsp_threshold      = 3000
# gapped_threshold   = 3000
# x_drop             = 910
# y_drop             = 9400
# gap_open_penalty   = 400
# gap_extend_penalty = 30
#        A    C    G    T
#   A   91 -114  -31 -123
#   C -114  100 -125  -31
#   G  -31 -125  100 -114
#   T -123  -31 -114   91
a score=4855
s Chr1 0 115 + 30427671 CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAATCTTTAAATCCTACATCCATGAATCCCTAAATACCTAATTCCCTAAACCCGAAACC
s Chr1 0 115 + 29952646 CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAATCCCAAAATCCCTAAATCACTTAATCCTTAAATCCCTAAATCCCTATACCCTAAACCCTAAACCC

我们来观察一下：

第二行显示了输出这个文件，我们所用到lastz的命令。然后一些其他的东西。

后面跟着一个类似表格的东西，表示的是如果ref序列中的A能够匹配到sample序列的A，那么我们就在这一个位点加91分。如果在ref中是A，在smaple中能够匹配到的序列中是C那么就要减去114分，其他的分数同理。最后

a score=4855这一行的分数就是这么来的。

然后下面就是配匹的内容。

Chr1：表示的是在一号染色体的配匹

0：表示的就是从染色体的第0个位置开始匹配。

115表示的是配匹了115个碱基

+代表的正链方向配匹

30427671：表示的是ref序列一共有这么多的碱基

后面就是序列的具体内容。

我们在这里要做的就是把序列中不同的碱基给找出来。第一行的（也就是ref序列）和第二行的(也就是sample)序列在其他的地方都一样，都能够匹配，那么我们就可以认为这个点两条序列不同的原因是发生了突变那我们就把这个位点给统计出来。

具体的程序如下：

#!/usr/bin/perl
#从.maf文件中提取SNP位点

use strict;
use warnings;

my $dna_filename;#{{{
$/="\n";
my @mout='';
my $z;
my @score;
my @info1;
my @info2;
my @sequ1;
my @sequ2;
my $len1;
my $cout;#}}}
my $pos1;
my $pos2;
my $snp;


print "please input the file input\n";
chomp($dna_filename=<STDIN>);
open (DNA,$dna_filename)||die("can not open the file!");

open(RESULT,">snp$dna_filename.txt")||die("can not open the file!");

$z=<DNA>;#{{{
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;
$z=<DNA>;#}}}

$/="\n\n";      #把行定义为空白行分割，也就是一次读入一个匹配单元，包括score，序列1，序列2
while(<DNA>)
{
	@mout   = split/\n/,$_;
	@score  = split/=/,$mout[0]; #$mout[0]是score的信息，[1]是序列1的信息，[2]是序列2的信息  
	unless ($score[1]<90000)     #判断lastz匹配出来的序列的值，只有当值满足一定条件时才进行操作
	{
		@info1 = split/\s+/,$mout[1]; #分别将序列1和序列2进行拆分，为了得到独立出序列的信息
		@info2 = split/\s+/,$mout[2];

		@sequ1 = split//  ,$info1[6]; #将序列再拆分成一个个的碱基，为了方便两个序列之间的比对
		@sequ2 = split//  ,$info2[6];

		for($cout=0;$cout<$info1[3];$cout++) #info1[3]是匹配的序列长度#这里错用了逗号代替分号导致很久的排错
		{
			if($sequ1[$cout] eq  $sequ2[$cout])
			{
				next;
			}
			else
			{
				$pos1=$info1[2]+$cout+1;
				$pos2=$info2[2]+$cout+1;
				print RESULT "$pos1 $sequ1[$cout] $pos2  $sequ2[$cout]\n";
			}

		}

	}
}

我们最后得到的结果如下：

989 C 1109  T
992 T 1112  G
996 T 1116  C
1001 G 1121  A
5266 C 5386  M
13677 C 13797  T
13678 T 13798  A
13679 T 13799  C
13681 A 13801  T

添加一个修改过的程序，可以一次处理多个文件：

#!/usr/bin/perl

use strict;
use warnings;

my @mout;
my @score;
my @info1;
my @info2;
my @seque1;
my @seque2;
my %snp;
my $z;
my $cout;
my $pos1;
my $pos2;
my $key;
my $sample;
my $filename;
my @sample=qw/TAIR_vs_bur.maf TAIR_vs_can.maf TAIR_vs_ct.maf   TAIR_vs_edi.maf
              TAIR_vs_hi.maf  TAIR_vs_kn.maf  TAIR_vs_ler.maf  TAIR_vs_mt.maf 
			  TAIR_vs_no.maf  TAIR_vs_oy.maf  TAIR_vs_po.maf   TAIR_vs_rsch.maf 
			  TAIR_vs_sf.maf  TAIR_vs_tsu.maf TAIR_vs_wil.maf  TAIR_vs_ws.maf 
			  TAIR_vs_wu.maf  TAIR_vs_zu.maf/; 
foreach $sample(@sample)
{
    $filename=$sample;
    open(DNA,"$filename") || die("can not open!");
    
    $/="\n";
    $z=<DNA>;#{{{
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;
    $z=<DNA>;#}}}
    
    $/="\n\n";
    while(<DNA>)
    {#{{{
	    @mout  = split/\n/ ,        $_;
	    @score = split/=/  ,  $mout[0];
	    @info1 = split/\s+/,  $mout[1];
	    @info2 = split/\s+/,  $mout[2];
	    @seque1 = split//   , $info1[6];
	    @seque2 = split//   , $info2[6];
    	
	    for($cout=0;$cout<$info1[3];$cout++)
	    {
		    if($seque1[$cout] eq $seque2[$cout])
		    {
			    next;
		    }
		    else
		    {
			    $pos1 = $info1[2]+$cout+1;
			    $pos2 = $info2[2]+$cout+1;
			    if(exists $snp{$pos1})
			    {
				    $key    = $snp{$pos1};
				    if($key>=$score[1])
				    {
					    next;
				    }
				    else
				    {
				        $snp{$pos1}{$score[1]}{$seque1[$cout]} = $seque2[$cout];
				    }
			    }
			    else
			    {
				    $snp{$pos1}{$score[1]}{$seque1[$cout]} = $seque2[$cout];
			    }
    
		    }
	    }
    
    }#}}}
     
    my $key1;
    my $key2;
    my $key3;
    my $value;
    
    open(OUT,">snp_$filename.txt") || die("can not open!");
    foreach $key1 (sort keys %snp)
    {#{{{
	    print OUT "$key1 ";
	    foreach $key2 (sort keys $snp{$key1})
	    {
		    foreach $key3 (sort keys $snp{$key1}{$key2})
		    {
			    $value = $snp{$key1}{$key2}{$key3};
			    print OUT "$key3 $value\n";
		    }
	    }
    }#}}}
	close DNA;
	close OUT;
	%snp=();
}