附件:Easyseq源码([4]mutation)

#!/usr/bin/perl
use Bio::SeqIO;
use Bio::Seq;
use Tie::File;
use Cwd;
use File::Path;
    $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    @file = glob("$location_dir\\*.filter");
foreach $file(@file){
    open(DATA,"<$file") or die "Can't open $file";
    @list = <DATA>;
  foreach $list (@list){
     @list_elm=split(" ",$list);
     $pairs=1;
     if($list_elm[0]=~/primer_pairs_count/){
     $primer_pairs_count=$list_elm[1];
     }
     while($pairs<$primer_pairs_count+1){
      if($list_elm[0]=~/F_$pairs/){
          $ORFF[$pairs]=$list_elm[1];
          $ORFFlength[$pairs]=$list_elm[1]=~s/[A-Za-z]/1/g;
      }
      if($list_elm[0]=~/P_$pairs/){
          $ORFP[$pairs]=$list_elm[1];
          $ORFPlength[$pairs]=$list_elm[1]=~s/[A-Za-z]/1/g;
      }
      if($list_elm[0]=~/R_$pairs/){
          $rev_seq= reverse ($list_elm[1]);
          $rev_seq=~tr/[ATCGatcg]/[TAGCtagc]/;
          $ORFR[$pairs]=$rev_seq;
          $ORFRlength[$pairs]=$list_elm[1]=~s/[A-Za-z]/1/g;
      }
      $pairs=$pairs+1;
     }
     if($list_elm[0]=~/Seq_N/){
      $Seq_N=$list_elm[1];
     }
     if($list_elm[0]=~/Seq_min_length/){
      $Seq_min_length=$list_elm[1];
     }
     if($list_elm[0]=~/Seq_max_length/){
      $Seq_max_length=$list_elm[1];
     }
     if($list_elm[0]=~/Species_names/){
      shift @list_elm;
      @Species_names=@list_elm;
     }
     if($list_elm[0]=~/Merger_acquisition_bases/){
      shift @list_elm;
      @Merger_acquisition_bases=@list_elm;
     } 
  }
        close DATA;
}
    $pai=1;
    $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    @file = glob("$location_dir\\mapped\\*.hit");
        open WH,">$location_dir\\mutation\\Results.txt";
        print WH "Variant_Strain;Primer;Mutation_Type;ID\n";
        close WH;
while($pai<=$primer_pairs_count){
 foreach $file(@file){
        $sequence[$pai]=$ORFF[$pai].$ORFP[$pai].$ORFR[$pai];
        $files=$file;
        $file2=$file;
        $location_dir=getcwd;
        $location_dir=~s/\//\\/;
        $files=~s/\://g;
        $files=~s/\\/_/g;
        $num=$file2=~s/\\//g;
        @files=split("_",$files);
        $files[$num+1]=~s/\.hit//;
        $LL=$files[$num+1];
    if($LL eq Primer.$pai){
        $IN=$sequence[$pai];
        $name=Primer.$pai;
        $num=0;
        $files=$file;
        $file2=$file;
        $location_dir=getcwd;
        $location_dir=~s/\//\\/;
        $files=~s/\://g;
        $files=~s/\\/_/g;
        $num=$file2=~s/\\//g;
        @files=split("_",$files);
        $dirpath="$location_dir\\mutation";
            mkdir ($dirpath); 
        $dirpath="$dirpath\\$files[$num]";
            mkdir ($dirpath); 
        $dirpathwhole="$dirpath\\$name";
            mkdir ($dirpathwhole);
        $fa=Bio::SeqIO->new(-file=>$file,-format=>'fasta');
            print " successful import $file\n";
            LOOP:do{
             while($seq_obj=$fa->next_seq){
                $name=Primer.$pai;
                @location=(); 
                $id=$seq_obj->id;
                $seq=$seq_obj->seq;
                @IN=split("",$IN);
                $c1=0;
                $c2=0;
                %seq=();
                @seq=split('', $seq);    
                foreach $Merger_acquisition_bases(@Merger_acquisition_bases){
                    if($seq=~/$Merger_acquisition_bases/){
                    goto LOOP;
                    }
                }
                    $seq_length=@seq;
                    $primer_length=@IN;
                while($c1<=$primer_length-1){
                    if($seq[$c1] eq $IN[$c1]){
                        $c1=$c1+1;
                    }else{   
                        
                            if($c1<=$ORFFlength[$pai]-1){
                                $location=F;
                                $c2=$c1+1;
                            }elsif($c1<=$ORFFlength[$pai]+$ORFPlength[$pai]-1){
                                $location=P;
                                $c2=$c1-$ORFFlength[$pai]+1;
                            }else{
                                $location=R;
                                $c2=$c1-($ORFFlength[$pai]+$ORFPlength[$pai]-1);
                            }
                            $mu="$location$c2\_$IN[$c1]\-$seq[$c1]";
                            push(@location, $mu);
                        
                        $c1=$c1+1;                         
                    }
                }
                $mu_loc=join(" ", @location);
                if($mu_loc){
                open WH,">>$location_dir\\mutation\\Results.txt";
                print WH "$files[$num];$LL;$mu_loc;$id\n";
                close WH;
                open WH,">>$dirpathwhole\\$mu_loc.txt";
                print WH "$id\n";
                close WH; 
                }else{
                $mu_loc="Matched";   
                open WH,">>$location_dir\\mutation\\Results.txt";
                print WH "$files[$num];$LL;$mu_loc;$id\n";
                close WH;
                open WH,">>$dirpathwhole\\$mu_loc.txt";
                print WH "$id\n";
                close WH; 

                }  
             } 
            }
    }
 }
    $pai=$pai+1;
}

  • 15
    点赞
  • 13
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值