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

PICC=CODISCHECK(handles,A,B,HA,HB,flag)
function PICC=CODISCHECK(handles,A,B,HA,HB,flag) %#ok<INUSL>
global I CODISpop NL
if flag
    KC=[2; ...%Parent-child
        1/2; ... %half-siblings or Grandparent-Grandchild or Uncle-nephew
        1/4; ... %First cousins
        1/16;]./4; %Second cousins
else
    KC=0.5;
end
L=length(KC);

PICCM=ones(18,L);
for I=1:18
    if (all(A(I,:)) && all(B(I,:)))
        %set the frequencies of the STRs
        eval(sprintf('locus=get(handles.text%i,''String'');',I))
        TF=sum(strcmp(locus,CODISpop(:,7)).*(1:1:NL)');
        %set the frequencies of the STRs
        x=CODISpop{TF,1}; y=CODISpop{TF,2};
        if flag
            Fmut=CODISpop{TF,3}; 
        else
            Fmut=CODISpop{TF,4}; 
        end
        Z=CODISpop{TF,5};
        for J=1:L
            PICCM(I,J)=Check(KC(J),x,y,Z,Fmut,I,A,B,HA,HB);
        end
    end
end
PICC=prod(PICCM);

function PI2=Check(KC,x,y,Z,Fmut,I,A,B,HA,HB)
z=round(x.*y); y=y+2; z(B(I,:)-1)=z(B(I,:)-1)+1+HB(I); x=z./y;
if (~HA(I) && ~HB(I)) && isequal(A(I,:),B(I,:)) %C=PQ & F=PQ
    z=round(x.*y); y=y+1; z(A(I,:)-1)=z(A(I,:)-1)+0.5; x=z./y;    
    PItmp=sum(x(A(I,:)-1))/prod(x(A(I,:)-1))/4;
    if KC==0.5
        PI2=PItmp;
    else
        PI2=PItmp*2*KC+(1-2*KC); %eventually transform PI in AI (Avuncular index)
    end
else % we are not in the condition C=PQ & F=PQ
    if any(ismember(A(I,:),B(I,:))) %there is a match        
        %check if swaps are required
        if ~HA(I) %the child's locus is heterozygous...
            if HB(I) %If the alleged father's locus is homozygous...
                 %swap the child's locus only if C=PQ & F=PP
                if A(I,1)==B(I,1)
                    A(I,:)=fliplr(A(I,:)); %C=QP & F=PP
                end
            else %the alleged father's locus is heterozygous...
                if A(I,1)==B(I,1) %if C=PQ & F=PX
                    A(I,:)=fliplr(A(I,:)); %C=QP & F=PX
                elseif A(I,1)==B(I,2) %if C=PQ & F=XP
                    A(I,:)=fliplr(A(I,:)); %C=QP & F=PX
                end
            end
        end
        %Update the frequencies
        %In this case, we know which allele is transmitted by the alleged father
        %to the child and so we can add the untransmitted allele
        z=round(x.*y); y=y+1; z(A(I,1)-1)=z(A(I,1)-1)+1; x=z./y;
        PItmp=1/((2-HA(I))*(2-HB(I))*x(A(I,2)-1)); %there is a match and so the correction factor is not required
        if KC==0.5
            PI2=PItmp;
        else
            PI2=PItmp*2*KC+(1-2*KC); %eventually transform PI in AI (Avuncular index)
        end
    else %there is not a match between the child's and the alleged father's allele        
        if KC==0.5
            %Now we must compute the correction factor
            if HA(I) %child's locus is homozygous C=PP
                %If the alleged father's locus is homozygous, F=RR, the quantity
                %|R-P|-1 will be summed twice, but the denominator will be 2 and so
                %the proper quantity will be computed.
                %On the contrary, if the locus is heterozygous, F=RS, the quantities
                %|R-P|-1 and |S-P|-1 will be both computed and the denominator will
                %be 1
                if HB(I)
                    corr=Fmut*(1/2)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,1)-1))-1);
                else
                    corr=Fmut*(1/2)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,1)-1))-1)+Fmut*(1/2)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,2)-1))-1);
                end
                PItmp=corr/Z(A(I,1)-1)^2; 
            else %child's locus is heterozygous C=PQ
                %also in this case the denominator will insure that the quantities
                %will not be summed twice
                if HB(I)
                    corr=Fmut*(1/2)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,1)-1))-1)+Fmut*(1/2)*(1/10)^(abs(Z(A(I,2)-1)-Z(B(I,1)-1))-1);
                else
                    corr=Fmut*(1/4)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,1)-1))-1)+Fmut*(1/4)*(1/10)^(abs(Z(A(I,2)-1)-Z(B(I,1)-1))-1)+...
                         Fmut*(1/4)*(1/10)^(abs(Z(A(I,1)-1)-Z(B(I,2)-1))-1)+Fmut*(1/4)*(1/10)^(abs(Z(A(I,2)-1)-Z(B(I,2)-1))-1);
                end
                PItmp=corr/(2*Z(A(I,1)-1)*Z(A(I,2)-1)); 
            end
        else
            PItmp=0;
        end
        %Finally, compute the PI and eventually transform PI in AI (Avuncular
        %index)
        PI2=PItmp*2*KC+(1-2*KC);
    end
end

Contact us