附件:Easyseq源码([7]mismatch)

#!/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);
    if($list_elm[0]=~/Target_F/){
      $ORFF=$list_elm[1];
      $OFlen=$list_elm[1]=~s/[A-Z]/1/g;
    }
    if($list_elm[0]=~/Target_P/){
      $ORFP=$list_elm[1];
      $OPlen=$list_elm[1]=~s/[A-Z]/1/g;
    }
    if($list_elm[0]=~/Target_R/){
      $rev_seq= reverse ($list_elm[1]);
      $rev_seq=~tr/[ATCGatcg]/[TAGCtagc]/;
      $ORFR=$rev_seq;
      $ORlen=$rev_seq[1]=~s/[A-Z]/1/g;
    }
    if($list_elm[0]=~/Ic_F/){
      $NF=$list_elm[1];
      $NFlen=$list_elm[1]=~s/[A-Z]/1/g;
    }
    if($list_elm[0]=~/Ic_P/){
      $NP=$list_elm[1];
      $NPlen=$list_elm[1]=~s/[A-Z]/1/g;
    }
    if($list_elm[0]=~/Ic_R/){
      $rev_seq= reverse ($list_elm[1]);
      $rev_seq=~tr/[ATCGatcg]/[TAGCtagc]/;
      $NR=$rev_seq;
      $NRlen=$rev_seq=~s/[A-Z]/1/g;
    }
     if($list_elm[0]=~/Merger_acquisition_bases/){
      shift @list_elm;
      @Merger_acquisition_bases=@list_elm;
    } 
    }
    close DATA;
    } 

join_mismatch();
mutation_location();
sub join_mismatch{
    $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    @file = glob("$location_dir\\mismatch\\*.mismatch");
    $dirpath="$location_dir\\mismatch";
    mkdir ($dirpath);
    foreach $file_txt(@file){
    open(DATA,"<$file_txt") or die "Can't open $file_txt";
    print " successful import $file_txt\n";
    @list = <DATA>;
    foreach $list (@list){
    chomp $list;
    @lines=split("    ",$list);
    if($file_txt=~/[A-Za-z]_ORF1ab/){
    if($lines[3] eq None){
        $ORFF2=$ORFF;
    }else{
    $ORFF2=$lines[3];
    }
    if($lines[7] eq None){
        $ORFP2=$ORFP;
    }else{
    $ORFP2=$lines[7];
    }
    if($lines[11] eq None){
        $ORFR2=$ORFR;
    }else{
    $ORFR2=$lines[11];
    }
    $ORF1ab2="$ORFF2"."$ORFP2"."$ORFR2";
    open WH,">>$file_txt.hit";
    print WH ">$lines[0]\n";
    print WH "$ORF1ab2\n";
    close WH;
    }elsif($file_txt=~/[A-Za-z]_N/){
    if($lines[3] eq None){
        $NF2=$NF;
    }else{
    $NF2=$lines[3];
    }
    if($lines[7] eq None){
        $NP2=$NP;
    }else{
    $NP2=$lines[7];
    }
    if($lines[11] eq None){
        $NR2=$NR;
    }else{
    $NR2=$lines[11];
    }
    $N2="$NF2"."$NP2"."$NR2";
    
    open WH,">>$file_txt.hit";
    print WH ">$lines[0]\n";
    print WH "$N2\n";
    close WH;
    }else{}
    }
    close DATA;
    }
}
sub mutation_location{
    $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    #print "$location_dir\n";
    @file = glob("$location_dir\\mismatch\\*.hit");
    my $dirpath2="$location_dir\\mismatch";
    mkdir ($dirpath2);
    open WH,">>$dirpath2\\Results.txt";
    print WH "Variant_Strain Gene Mutation_Type ID\n";
    close WH;
    foreach $file(@file){
    print "$file\n";
    $n=0;
    $ORF1ab="$ORFF"."$ORFP"."$ORFR";
    $N="$NF"."$NP"."$NR";
    if($file=~/ORF1ab/){
    $dir=$&;
    $IN=$ORF1ab;
    $name=O;
    }
    if($file=~/N/){
    $dir=$&;
    $IN=$N;
    $name=N;
    }
    $num=0;
    $files=$file;
    $file2=$file;
    $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    $files=~s/\://g;
        #print "$files\n";
    $files=~s/\\/_/g;
    $num=$file2=~s/\\//g;
        #print "$num\n";
    @files=split("_",$files);
        #print "$location_dir\n";
        #$results="$location_dir\\results";
        #mkdir ($results);
        #print "$results\n";
    $dirpath="$location_dir\\mismatch\\$files[$num]";
    mkdir ($dirpath); 
        #print "$dirpath\n";
    $dirpathwhole="$dirpath\\$dir";
    mkdir ($dirpathwhole);
        #print "$dirpathwhole\n";
        #print "loading $file successiful\n";
    $fa=Bio::SeqIO->new(-file=>$file,-format=>'fasta');
    LOOP:do{
    while($seq_obj=$fa->next_seq){
    @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/g){
                print "$id contains Merger acquisition bases: $Merger_acquisition_bases\n";
                goto LOOP;
            }
    }
            $seq_length=@seq;
            $primer_length=@IN;
            while($c1<=$primer_length-1){
                if($seq[$c1] eq $IN[$c1]){
                    $c1=$c1+1;}
                else{
        
                if($file=~/ORF1ab/){
                    if($c1<=$OFlen-1){
                        $location=F;
                        $c2=$c1+1;
                    }elsif($c1<=$OFlen+$OPlen-1){
                        $location=P;
                        $c2=$c1-$OFlen+1;
                    }else{
                        $location=R;
                        $c2=$c1-$OFlen-$OPlen+1;
                    }
                }
                if($file=~/N/){
                    if($c1<=$NFlen-1){
                        $location=F;
                        $c2=$c1+1;
                    }elsif($c1<=$NFlen+$NPlen-1){
                        $location=P;
                        $c2=$c1-$NFlen+1;
                    }else{
                        $location=R;
                        $c2=$c1-$NFlen-$NPlen+1;
                    }
                }
                $mu="$location$c2$IN[$c1]\-$seq[$c1]";
                push(@location, $mu);
                $c1=$c1+1;                     
                }
            }
                $mu_loc=join("",@location);
                #print "$mu_loc\n";   
     
    open WH,">>$dirpath2\\Results.txt";
    print WH "$files[$num] $dir $mu_loc $id\n";
    close WH;
    open WH,">>$dirpathwhole\\$mu_loc.txt";
    print WH "$id\n";
    close WH;
    }
    }
    }
}
 

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值