%%%This function generates a PgSNP file format for UCSC Genome Browser
%%%visualization
function []=GeneratePgSNPFile(chrz,SaNum,path1)
tail='.txt';
head='SA';
SaNumStr=num2str(SaNum);
filename=strcat(head,SaNumStr,chrz,tail);
fid = fopen(filename);
C =textscan(fid, '%s %f %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s');
fclose(fid);
%create a new cell array
Chrome{1,1}=C{1,1};       %Chromosome number
Chrome{1,2}=C{1,2};       %Start of SNP
Chrome{1,3}=C{1,2}+1;     %End of SNP
Chrome{1,4}=C{1,4};       %SNP
Chrome{1,5}=ones(length(Chrome{1,1}),1);   % Number of alleles listed in the  SNP field 
Chrome{1,6}=C{1,6};  %Reads for each SNP
Chrome{1,7}=C{1,5}; %Score

% Change the IUPAC code for the SNP to a single nucleotide change
for i=1:length(Chrome{1,4})
    if  (strcmp(Chrome{1,4}(i), 'R')==1) % R=A/G
        if (strcmp(C{1,3}(i), 'A')==1)
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'R','G');
        else
         Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'R','A')
        end
    end
    if  (strcmp(Chrome{1,4}(i), 'Y')==1) % Y=C/T
        if (strcmp(C{1,3}(i), 'C')==1)
       Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'Y','T')
        else
         Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'Y','C')
        end
    end
if  (strcmp(Chrome{1,4}(i), 'S')==1) % S=G/C
        if (strcmp(C{1,3}(i), 'G')==1)
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'S','C')
        else
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'S','G')
        end
end
    if  (strcmp(Chrome{1,4}(i), 'W')==1) % W=A/T
        if (strcmp(C{1,3}(i), 'A')==1)
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'W','T')
        else
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'W','A')
        end
    end
    if  (strcmp(Chrome{1,4}(i), 'K')==1) % K=G/T
        if (strcmp(C{1,3}(i), 'G')==1)
       Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'K','T')
        else
        Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'K','G')  
        end
    end
    if  (strcmp(Chrome{1,4}(i), 'M')==1) % M=A/C
        if (strcmp(C{1,3}(i), 'A')==1)
      Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'M','C')
        else
         Chrome{1,4}(i)=strrep(Chrome{1,4}(i),'M','A') 
        end
    end
end


for i=1:length(Chrome{1,4})
    %Check if there is information for Reads
    if not (strcmp(Chrome{1,6}(i), '0')==1) %if there is information for Reads
        Chrome{1,5}(i)=2;% Then there are 2 alleles
        a=C{1,3}(i); %Consensus nt
        b=Chrome{1,4}(i); %SNP nt
        c=Chrome{1,7}(i); %Score for the SNP
       Chrome{1,4}(i)= strcat(a,'/',b);  %Output in format "Consensus/SNP"
       Chrome{1,7}(i)= strcat('0',',',c); % Output in format " Score for Consensus,Score for SNP"
      d= Chrome{1,6}(i); %Reads for SNP
      e=C{1,7}(i);  %Reads for the consensus

   %Outputs the consensus and SNP reads in format: "Consensus Reads, SNP Reads". 
      read1=d{1}(3);
      read2=e{1}(3);
      if length(d{1})>3
      for m=4:(length(d{1})) 
          read1=strcat(read1,d{1}(m));
      end
      end
      if length(e{1})>3
      for m=4:(length(e{1}))
          read2=strcat(read2,e{1}(m));
      end
      end
   Chrome{1,6}(i)= strrep(Chrome{1,6}(i), d, strcat(read2,',',read1)); 

    end
end

%checking for the min and max positions of SNPs in the defined chromosome to assess  the range for UCSC visualiation
max_position=num2str(max(Chrome{1,2})+1000);
min_position=num2str(min(Chrome{1,2})-1000);

    
% creating a new text file in the PgSNP format which is compatible with UCSC
% custom track and can be easily uploaded and viewed
filename = strcat(head,SaNumStr,'PgSNP',chrz,tail);
fid2 = fopen(strcat(path1,'\', filename),'w');
title1=strcat('track type=pgSnp visibility=3 db=hg19 name="Genome Gems Track" description="Genome Gems-Personal Genome SNP file" '); 
title2=('browser position');
title3=strcat(chrz,':',min_position,'-',max_position);
    chrname=char(Chrome{1,1}(1));
    snploc1(1)=Chrome{1,2}(1);
    snploc2(1)=Chrome{1,3}(1);
    snp=char(Chrome{1,4}(1));
    numofsnp(1)=Chrome{1,5}(1);
    reads=char(Chrome{1,6}(1));
    score=char(Chrome{1,7}(1));
    fprintf(fid2,'%s\n%s\t%s\n%s %d %d %s %d %s %s \n',title1, title2, title3,chrname, snploc1(1), snploc2(1), snp, numofsnp(1), reads, score);
for i=2:length(Chrome{1,1})
    chrname=char(Chrome{1,1}(i));
    snploc1(i)=Chrome{1,2}(i);
    snploc2(i)=Chrome{1,3}(i);
    snp=char(Chrome{1,4}(i));
    numofsnp(i)=Chrome{1,5}(i);
    reads=char(Chrome{1,6}(i));
    score=char(Chrome{1,7}(i));
    fprintf(fid2,'%s %d %d %s %d %s %s \n',chrname, snploc1(i), snploc2(i), snp, numofsnp(i), reads, score);
end