Code covered by the BSD License  

Highlights from
CODIS

image thumbnail

CODIS

by

 

14 May 2009 (Updated )

Combined DNA Index System: a GUI tool for forensic paternity lawsuit

PI=CODIS3REV(handles)
function PI=CODIS3REV(handles) %#ok<INUSD>
%CODIS3 is a CODIS ancillary routine.
%This function is recalled by CODIS to compute the Paternity Index in the cases
%of reverse parentage (i.e. when the suspect is that the child was exchanged) 
%
% Syntax: 	CODIS3REV
%      
%     Inputs:
%           all needed variables are global, so CODIS3 directly access to the
%           RAM
%     Outputs:
%           - ordered loci, PI, and all needed flags
%
%      Example: 
%           run CODIS and select the DEMO menu->Trio attribution or Trio exclusion
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it

%each variable is inheredited by CODIS
global M C F I CODISpop EP RMNE matF matM Amb pmatF loci NL

%M, C and F are the matrices of mother's, child's and alleged father's STRs data;
%sort M C and F vectors
M=sort(M,2); C=sort(C,2); F=sort(F,2);
%HM, HC and HF are the flag matrices of heterozygosity (HM - mother; HC - child; 
%HF - alleged father) HX(I) is 0 if the I-th locus is homozygous, 1 otherwise
HM=~(M(:,1)==M(:,2)); HC=~(C(:,1)==C(:,2)); HF=~(F(:,1)==F(:,2));
%Set the match-flag vectors
matF=zeros(18,1); matM=matF; 
%pmatF: there are some particular cases in which there is a probable match: i.e.
%M=XP C=PQ F=PP. In this case, the alleged father doesn't bring the Q allele so
%there is a mismatch; anyway it is possible the alleged father gave the P allele
%and the mother gave X or P allele and it was mutated in Q.
pmatF=matF;
%Amb is the ambigous match-flag: if M=PQ C=PQ F=PQ there is a match, but
%it is impossible to extablish which allele is inheredited allele from the
%mother and which from the father, so Amb(I)=2; otherwise Amb(I)=1;
Amb=matF;
%set the Parentage Index
PI=matF;

%Start to compute the reverse parentage index for each locus
for I=1:18
    if (all(C(I,:)) && all(F(I,:)) && all(M(I,:))) 
        eval(sprintf('loci{I}=get(handles.text%i,''String'');',I))
        TF=sum(strcmp(loci{I},CODISpop(:,7)).*(1:1:NL)');
        %set the frequencies of the STRs
        x=CODISpop{TF,1}; y=CODISpop{TF,2};
        %Set the mutation rate and the STRs label array
        Fmut=CODISpop{TF,3}; Mmut=CODISpop{TF,4}; Z=CODISpop{TF,5};
        %update the STR frequencies to avoid that PI=inf if you find a child
        %with an allele with a null frequency
        z=round(x.*y);
        z(C(I,:)-1)=z(C(I,:)-1)+1;
        y=y+2; x=z./y;
        %EP is the exclusion probability, the fraction of couples that can be
        %excluded as parents
        %the EP is (1-FreqC1)^2 or (1-FreqC1-FreqC2)^2 for each parent.
        %So the couples that can be excluded are:
        %1) Both parents are not compatible with child EP^2
        %2) Only one parent can be excluded 2*EP*(1-EP)
        %summing: EP^2+2*EP*(1-EP)=EP*(2-EP)
        EPtmp=(1-sum(x(C(I,:)-1))/(1+(~HC(I))))^2;
        EP(I)=EPtmp*(2-EPtmp);
        al=prod(x(C(I,:)-1));
        clear y x
        if sum(ismember(M(I,:),C(I,:)))==0 %there is not a match between mother and child....
            if sum(ismember(F(I,:),C(I,:)))==0 %...and there is not a match between father and child
                if HC(I)==0 %the child locus is homozygous
                     c=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,1)-1))-1); %c is always computed
                     d=HM(I)*(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,1)-1))-1); %d is only computed when the mother is heterozygous
                     e=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,1)-1))-1); %e is always computed
                     f=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,1)-1))-1); %f is only computed when the father id heterozygous
                     PI=PItmp*Mmut*(c+d)*Fmut*(e+f);
                else %the child locus is heterozygous
                     a=(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,1)-1))-1);
                     b=(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,2)-1))-1);
                     c=HM(I)*(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,1)-1))-1);
                     d=HM(I)*(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,2)-1))-1);
                     e=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,1)-1))-1);
                     f=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1);
                     g=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,1)-1))-1);
                     h=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,2)-1))-1);
                     PI=PItmp*Mmut*(a+b+c+d)*Fmut*(e+f+g+h);
                end
            else %...but there is a match between father and child
                matF(I)=1; %switch on the flag
                %tidy up the loci considering only the father and the child
                %Check if F=PQ & C=PQ
                eqflag=(HF(I) && HC(I) && isequal(F(I,:),C(I,:)));
                if HC(I)==0 %the child's locus is homozygous
                    %the father's locus is heterozygous and it needs to be swapped
                    if HF(I)==1 && F(I,2)==C(I,1) %C=PP F=QP
                        CODISswap(0,0,1) %C=PP F=PQ
                    end
                else %the child's locus is heterozygous
                    if HF(I)==0 %the father's locus is homozygous
                        %and the child's locus need to be swapped
                        if C(I,1)==F(I,1) %C=PQ F=PP
                            CODISswap(0,1,0) %C=QP F=PP
                        end
                    else %the father's locus is heterozygous
                        if ~eqflag %if we are not in the situation C=PQ F=PQ
                            if C(I,2)==F(I,2)
                                CODISswap(0,0,1) %swap only the father's locus
                            elseif C(I,1)==F(I,2)
                                CODISswap(0,1,1) %swap both loci
                            elseif C(I,1)==F(I,1)
                                CODISswap(0,1,0) %swap only the child's locus
                            end
                        end
                    end
                end
            end
            %now compute a temporary PI
            %The PI is a likelihood ratio.
            %Numerator is computed on the basis of the Punnet square.
            %It is the probability to observe the genotype of the child given the
            %genotypes of the alleged parents. It could be 0,1,1/2 and 1/4.
            %Denominator depens by child genotype under Hardy-Weinberg law: p^2 or 2pq (allele frequencies)
            %Using a little trick....
            PItmp=1/(al*2^(HM(I)+HC(I)+HF(I)));
            %now check the mother
            if ~eqflag %if we are not in the situation C=PQ F=PQ
                c=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,1)-1))-1); %c is always computed
                d=HM(I)*(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,1)-1))-1); %d is only computed when the mother is heterozygous
                PI(I)=PItmp*Mmut*(c+d);
            else
                if HM(I)==0 %if the mother locus is homozygous
                    c=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,1)-1))-1); 
                    d=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,2)-1))-1); 
                    PI(I)=PItmp*Mmut*(c+d);
                else
                    c=(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,1)-1))-1);
                    d=(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,2)-1))-1);
                    e=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,1)-1))-1);
                    f=(1/2)*(1/10)^(abs(Z(M(I,2)-1)-Z(C(I,2)-1))-1);
                    PI(I)=PItmp*Mmut*(c+d+e+f);
                    Amb(I)=1;
                end
            end
        else %if there is the match between mother and child...
            matM(I)=1; %switch on the flag
            %tidy up the loci considering only the mother and the child
            %Check if M=PQ & C=PQ
            eqflag=(HM(I) && HC(I) && isequal(M(I,:),C(I,:)));
            if HC(I)==0 %the child's locus is homozygous
                %the mother's locus is heterozygous and it needs to be swapped
                if HM(I)==1 && M(I,1)==C(I,1) %M=PQ C=PP
                    CODISswap(1,0,0) %M=QP C=PP
                end
            else %the child's locus is heterozygous
                if HM(I)==0 %the mother's locus is homozygous
                    %and the child's locus need to be swapped
                    if C(I,2)==M(I,2) %M=PP C=QP
                        CODISswap(0,1,0) %M=PP C=PQ
                    end
                else %the mother's locus is heterozygous
                    if ~eqflag %if we are not in the situation M=PQ C=PQ
                        if M(I,1)==C(I,1)
                            CODISswap(1,0,0) %swap only the mother's locus
                        elseif M(I,1)==C(I,2)
                            CODISswap(1,1,0) %swap both loci
                        elseif M(I,2)==C(I,2)
                            CODISswap(0,1,0) %swap only the child's locus
                        end
                    end
                end
            end
            %now compute a temporary PI
            %compute the PI.
            %The PI is a likelihood ratio.
            %Numerator is computed on the basis of the Punnet square.
            %It is the probability to observe the genotype of the child given the
            %genotypes of the alleged parents. It could be 0,1,1/2 and 1/4.
            %Denominator depens by child genotype under Hardy-Weinberg law: p^2 or 2pq (allele frequencies)
            %Using a little trick....
            PItmp=1/(al*2^(HM(I)+HC(I)+HF(I)));
            clear y z al 
            %now check the father
            if ~eqflag %if we are not in the case M=PQ C=PQ...
                %.. the alleged father can transmit only the C(I,2) allele
                if ismember(C(I,2),F(I,:)) %there is a match between child's and alleged father's loci
                    %The PI doesn't need correction and the match-flag is
                    %switched on
                    PI(I)=PItmp; matF(I)=1;
                    %check if a swap of the alleged father's locus is needed
                    if C(I,2)==F(I,2)
                        CODISswap(0,0,1)
                    end
                else  %there isn't a match between child's and alleged father's loci
                    if HC(I)==0 %the child's locus is homozygous
                        c=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1); %c is always computed
                        d=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,2)-1))-1); %d is only computed when the father is heterozygous
                        PI(I)=PItmp*Fmut*(c+d);
                    else %the child's locus is heterozygous
                        %check if could be a possible match between child's allele 1 and the alleged father's alleles
                        pmatch=ismember(F(I,:),C(I,1));
                        if any(pmatch) %mutation in the mother? Consider it!
                            pmatF(I)=1; %swith on the possible match flag
                            %check is a swap is needed
                            if HF(I) && find(pmatch)==2
                                CODISswap(0,0,1)
                            end
                            c=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1); %c is always computed
                            d=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,2)-1))-1); %d is only computed when the father is heterozygous
                            e=HM(I)*(1/2)*(1/10)^(abs(Z(M(I,1)-1)-Z(C(I,2)-1))-1); %d is only computed when the mother is heterozygous
                            PI(I)=PItmp*(Fmut*(c+d)+Mmut*(c+e));
                        else %there isn't any match between child and alleged father
                            c=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1); %c is always computed
                            d=HF(I)*(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,2)-1))-1); %d is only computed when the father is heterozygous
                            PI(I)=PItmp*Fmut*(c+d);
                        end
                    end
                end
            else %M=PQ C=PQ
                Amb(I)=(HM(I) && HC(I) && HF(I)) && isequal(M(I,:),C(I,:),F(I,:));
                if Amb(I) %M=PQ C=PQ F=PQ
                    PI(I)=PItmp*2; matF(I)=1;
                else %M=PQ C=PQ F=...
                    switch HF(I) %check the alleged father's locus
                        case 0 %The alleged father's locus is homozygous
                            if any(ismember(C(I,:),F(I,1))) %there is a match between child's and alleged father's loci
                                PI(I)=PItmp; matF(I)=1;
                                %check if a swap is needed
                                if C(I,1)==F(I,1)
                                    CODISswap(0,1,0) %swap the child's locus
                                else
                                    CODISswap(1,0,0) %swap the mother's locus
                                end
                            else %there isn't a match between child's and alleged father's loci
                                c=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,1)-1))-1);
                                d=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1);
                                PI(I)=PItmp*Fmut*(c+d);
                                Amb(I)=1;
                            end
                        case 1 %The alleged father's locus is heterozygous
                            if any(ismember(C(I,:),F(I,:))) %there is a match between child's and alleged father's loci
                                PI(I)=PItmp; matF(I)=1;
                                %check if a swap is needed
                                if C(I,1)==F(I,1)
                                    CODISswap(0,1,0)
                                elseif C(I,1)==F(I,2)
                                    CODISswap(0,1,1)
                                elseif C(I,2)==F(I,1)
                                    CODISswap(1,0,0)
                                elseif C(I,2)==F(I,2)
                                    CODISswap(1,0,1)
                                end
                            else %there isn't a match between child's and alleged father's loci
                                c=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,1)-1))-1);
                                d=(1/2)*(1/10)^(abs(Z(F(I,1)-1)-Z(C(I,2)-1))-1);
                                e=(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,1)-1))-1);
                                f=(1/2)*(1/10)^(abs(Z(F(I,2)-1)-Z(C(I,2)-1))-1);
                                PI(I)=PItmp*Fmut*(c+d+e+f);
                                Amb(I)=1;
                            end
                    end
                end
            end
        end
    end
end
%OEP is the overall exclusion probability
OEP=1-prod(EP(EP~=0));
%compute the Random Couple Not Excluded percentage (RCNE)
%The RMNE is an estimation of the percentage of random men that can be
%compatible with paternity. Pratically, if the child's locus is PP a man that
%brings the P allele cannot be excluded as the father of the child (RM=PP or
%RM=PX); if the child's locus is PQ all the men that bring the P and the Q
%alleles cannot be excluded (RM=PP, QQ, PX, QX, PQ)
%The same for the mother.
RMNE=1-OEP;

Contact us