Was this page helpful?

filterSNP

    Table of contents
    No headers
     
    # Main program
    #
    #-------------------------------------------------------------------------------
    #!/usr/bin/perl
    # ddbi.pl#
    #-------------------------------------------------------------------------------------------
    #
    # name: Parse_Filter_SNP.pl
    #
    # author: Jie Yin     
    #
    # date: 11-27-13
    # 1. Pasrse SNP files (.vcf)
    # 2. Implement SNP filters
    #     a.Biallelic SNP only
    #     b.QUAL (Confidence of SNP) >=30
    #     c.Reads support ALT >= 10% of all informative reads
    # 3. Generate a filtered SNP file (.vcf).
    #
    #-------------------------------------------------------------------------------
     
    #example usage:
     
    # perl filterSNP.pl /Users/jieyin/Desktop/Trancriptomics/NewAssebly/SangerdataforSNP/Sanger_var_filtered.vcf  /Users/jieyin/Desktop/Trancriptomics/NewAssebly/SangerdataforSNP/Sanger_SH_var_filtered.vcf  
     
    use DBI;
    use strict;
    use Data::Dumper;
     
    #Main part
     
    my %filter;
    my %sh;
    my @header;   # all lines starting with "#".
     
     
    #functions
    ParseFilter($ARGV[0]);
    SHSNP($ARGV[1]);
    Print();
     
     
     
     
    sub ParseFilter {
     
     
    my ($allSNP) = @_;
    my $indx;
     
    my @header;
            my $inner;
    my %inner;
     
     
     
    open my $read_line, "<", $allSNP or die " unable to open file $allSNP\n";
     
    while( my $line = <$read_line> ) { 
     
    chomp $line;
            
            # $indx, $chr, $pos, $ID, $ref, $alt, $indel, $qua.
    if ( $line =~ /^scaffold/ ) {
                   $indx ++;
    ( my $sca, my $pos, my $ID, my $ref, my $alt, my $qua, my $filter, my $other ) = split "\t", $line, 8;
    my $scanum=substr($sca,9);
                    my $uniq_ID = join('_',$scanum,$pos);          
                           
                if ( ($alt !~ /\,/) && ($filter eq "PASS" )   ){
        
        $filter{$uniq_ID}->{sca}=$sca;
          
            $filter{$uniq_ID}->{scanum}=$scanum;
        
         $filter{$uniq_ID}->{pos}=$pos;
     
      $filter{$uniq_ID}->{ref}=$ref;
      
        $filter{$uniq_ID}->{alt}=$alt;
       
          
     
     
     
    # get rid of the effect of PCR error.
       }
     
       else{
             
       }
             }
            
            else{
                push @header, $line;
            }
      }
     
    }
     
     
    sub SHSNP {
    my ($shSNP) = @_;
    my $indx;
     
    my @header;
            my $inner;
    my %inner;
     
     
     
    open my $read_line, "<", $shSNP or die " unable to open file $shSNP\n";
     
    while( my $line = <$read_line> ) { 
     
    chomp $line;
                    if ( $line =~ /^scaffold/ ) {
                         my  $indx ++;
    ( my $sca, my $pos,  my $ref, my $alt ) = split "\t", $line, 4;
    my $scanum=substr($sca,9);
     
                    my $uniq_ID = join('_',$scanum,$pos);
     
        $sh{$uniq_ID}->{sca}=$sca;
        
         $sh{$uniq_ID}->{scanum}=$scanum;
        
          $sh{$uniq_ID}->{pos}=$pos;
            
      $sh{$uniq_ID}->{ref}=$ref;
      
        $sh{$uniq_ID}->{alt}=$alt;
     
     
             }
            
            
      }
     
    #print Dumper (\%filter);
    }
    print Dumper (\%sh);
     
      
    sub Print {
     
    my %inner1;
     
     
     
        my $filename1 = "/Users/jieyin/Desktop/Trancriptomics/NewAssebly/SangerdataforSNP/Sanger_YKvsSH_SNP.txt";
        open ( OUTFILE, ">$filename1" ) || die "Unable to open file $filename1\n";
        
     
         foreach my $key (  sort {  { $filter{$a}->{scanum} <=> $filter{$b}->{scanum} }  &&  { $filter{$a}->{pos} <=> $filter{$b}->{pos} }   } keys %filter) {
     
    if (  !defined $sh{$key} ) {
     
            %inner1=%{$filter{$key}};
     
              printf OUTFILE "$inner1{sca}\t$inner1{pos}\t$inner1{ref}\t$inner1{alt}\n";  
    }
       }
                
           close OUTFILE;
        
     
    }
     
     
    Was this page helpful?
    Tag page (Edit tags)
    • No tags
    You must login to post a comment.