perl/etc/justonce/oyk/FGI/rsCPI.pm

   1 package main;
   2
   3 use strict;
   4 #no strict "subs";
   5 use warnings;
   6 #use Data::Dump qw(ddx);
   7 =head1 NAME && VERSION
   8    cpi.pm
   9    Version:  V1.1
  10 =head1 DESCRIPTION
  11         chosen snp
  12 =head1 AUTHOR && CONTACT
  13     Author :  gaoshengjie huxuesong 2018/05/06
  14     modified at 2019/06/18
  15     Contact:  gaoshengjie@genomics.org.cn
  16 =head1 UPDATE LOG
  17     For missing record, MAF = 0.45, Ref/Alt follow VCF.
  18 =cut
  19
  20 our (%Markers,%MarkerAF);
  21 # %Markers
  22 #   {
  23 #     rs1000381   => ["chr5", 56885878, 0.124304011506483, 0.187151034246758],
  24 #     rs10005624  => ["chr4", 125929227, 0.122668415174707, 0.186323232412646],
  25 #     rs10011941  => ["chr4", 156398437, 0.116374475912992, 0.183033072043504],
  26 #   }
  27 # %MarkerAF
  28 #   {
  29 #     rs9995391   => { A => 0.5391, T => 0.4609 },
  30 #     rs9929923   => { A => 0.4429, C => 0.5522, G => "0.0050" },
  31 #     rs9996524   => { A => "0.5450", G => "0.4550" },
  32 #   }
  33
  34 #print $CHR{"SNP5538"}{"ref"}."\n";
  35 #SNP1110 501     G,A     999     G/G;9996,7      A/A;8,9955      G/A;5895,61;2.58311161233576e-19;G/G;5895,61
  36 our $err=0;
  37 our $special=0;
  38 our $pass=0;
  39 sub getcpiD(@) {
  40         my @a = @_;
  41         if(@a<6){
  42                 $pass++;
  43                 return [1,0,$err,$special,$pass];
  44         }
  45         my @alleles = qw{A C G T};
  46         my @geno=split /\,/,$a[2];
  47         my ($pe);
  48         if (exists $MarkerAF{$a[0]}){
  49                 $pe = $Markers{$a[0]}->[2];
  50                 foreach my $allele (@alleles){
  51                         unless (defined $MarkerAF{$a[0]}{$allele}){
  52                                 $MarkerAF{$a[0]}{$allele} = 0.01;
  53                         }
  54                 }
  55         }else{
  56                 $MarkerAF{$a[0]}{$geno[0]} = 0.55;
  57                 for my $i (1..scalar @geno - 1){
  58                         $MarkerAF{$a[0]}{$geno[$i]} = 0.45;
  59                 }
  60                 $pe = 0.125;
  61         }
  62         my @dad=split /\;/,$a[5];
  63         my @ch=split /\;/,$a[6];
  64         my @genodad=split /\//,$dad[0];
  65         my @genoch=split /\//,$ch[0];
  66         if ($genoch[0] ne $genodad[0] && $genoch[0] ne $genodad[1] && $genoch[1] ne $genodad[0] && $genoch[1] ne $genodad[1]){
  67                 $err++;
  68                 return [1e-4,$pe,$err,$special,$pass];
  69         } elsif ($genoch[0] eq $genoch[1]){
  70                 if ($genodad[0] eq $genodad[1]){
  71                         return [1/$MarkerAF{$a[0]}{$genoch[0]},$pe,$err,$special,$pass];
  72                 }elsif ($genodad[0] ne $genodad[1]){
  73                         return [0.5/$MarkerAF{$a[0]}{$genoch[0]},$pe,$err,$special,$pass];
  74                 }
  75         } elsif ($genoch[0] ne $genoch[1]){
  76                 if ($genodad[0] eq $genodad[1]){
  77                         return [0.5/$MarkerAF{$a[0]}{$genodad[0]},$pe,$err,$special,$pass];
  78                 } elsif ($genodad[0] ne $genodad[1]){
  79                         my @Gchild = sort @genoch;
  80                         my @Gdad = sort @genodad;
  81                         if ($Gchild[0] eq $Gdad[0] && $Gchild[1] eq $Gdad[1]){
  82                                 my $pi = ($MarkerAF{$a[0]}{$Gchild[0]} + $MarkerAF{$a[0]}{$Gchild[1]}) / (4 * $MarkerAF{$a[0]}{$Gchild[0]} * $MarkerAF{$a[0]}{$Gchild[1]});
  83                                 return [$pi,$pe,$err,$special,$pass];
  84                         }else{
  85                                 push @Gchild,@Gdad;
  86                                 my %count;
  87                                 my @need = grep { ++$count{$_} > 1 } @Gchild;
  88                                 return [0.25/$MarkerAF{$a[0]}{$need[0]},$pe,$err,$special,$pass];
  89                         }
  90                 }
  91         }
  92         return [999,$pe,-111,-222,-333];
  93 }
  94 sub getcpiT(@) {
  95         no warnings 'uninitialized';
  96         my @a = @_;
  97         if(@a<6){
  98                 $pass++;
  99                 return [1,0,$err,$special,$pass];
 100         }
 101         my @geno=split /\,/,$a[2];
 102         my $pe;
 103         my @alleles = qw{A C G T};
 104         if (exists $MarkerAF{$a[0]}) {
 105                 $pe = $Markers{$a[0]}->[3];
 106                 foreach my $allele (@alleles){
 107                         unless (defined $MarkerAF{$a[0]}{$allele}){
 108                                 $MarkerAF{$a[0]}{$allele} = 0.01;
 109                         }
 110                 }
 111         } else {
 112                 $MarkerAF{$a[0]}{$geno[0]} = 0.55;
 113                 for my $i (1..scalar @geno - 1){
 114                         $MarkerAF{$a[0]}{$geno[$i]} = 0.45;
 115                 }
 116                 $pe = 0.1875;
 117         }
 118         my @mm=split /\;/,$a[4];
 119         my @dad=split /\;/,$a[5];
 120         my @ch=split /\;/,$a[6];
 121         my @genomm=split /\//,$mm[0];
 122         my @genodad=split /\//,$dad[0];
 123         my @genoch=split /\//,$ch[0];
 124         if($genoch[0] ne $genomm[0] && $genoch[1] ne $genomm[1]){
 125                 $pass++;
 126                 return [1,0,$pe,$err,$special,$pass];
 127         } elsif ($genoch[0] eq $genoch[1]){
 128                 if ($genoch[0] ne $genodad[0] && $genoch[0] ne $genodad[1]){
 129                         $err++;
 130                         return [1e-4,$pe,$err,$special,$pass];
 131                 } elsif ($genodad[0] eq $genodad[1]){
 132                         return [1/$MarkerAF{$a[0]}{$genoch[0]},$pe,$err,$special,$pass];
 133                 } elsif ($genodad[0] ne $genodad[1]){
 134                         return [0.5/$MarkerAF{$a[0]}{$genoch[0]},$pe,$err,$special,$pass];
 135                 }
 136         } elsif ($genoch[0] ne $genoch[1]){
 137                 my $fromF = $genoch[0] eq $genomm[0] ? $genoch[1] : $genoch[0];
 138                 if ($fromF ne $genodad[0] && $fromF ne $genodad[1]){
 139                         $err++;
 140                         return [1e-4,$pe,$err,$special,$pass];
 141                 } elsif ($genodad[0] eq $genodad[1]){
 142                         return [1/$MarkerAF{$a[0]}{$fromF},$pe,$err,$special,$pass];
 143                 } elsif ($genodad[0] ne $genodad[1]){
 144                         return [0.5/$MarkerAF{$a[0]}{$fromF},$pe,$err,$special,$pass];
 145                 }
 146         }
 147         return [999,$pe,-111,-222,-333];
 148 }
 149
 150 1;