function writescf(filename,dat,primeridstring)
% WRITESCF writes SCF formatted files.
%
% SCF format files are used to store data from DNA sequencing instruments.
% This version of WRITESCF only adds the raw traces to the SCF files, it
% does not write the bases or confidence values.
%
% WRITESCF(FILENAME,DAT,PRIMERIDSTRING)
%
% FILENAME - is a char string with the file name
% DAT - is a [nx4] matrix with the raw traces ordered as GACT
% PRIMERID - is a scalar from 1 to 84 or any of the char strings in
% the following list
%
% 1 "DP6%25Ac{-21M13}" primer rhodamine ABI_373_377
% 2 "DP6%Ac{-21M13}" primer rhodamine ABI_373_377
% 3 "DP6%25Ac{M13Rev}" primer rhodamine ABI_373_377
% 4 "DP6%Ac{M13Rev}" primer rhodamine ABI_373_377
% 5 "DyePrimer{-21m13}" primer rhodamine ABI_373_377
% 6 "DyePrimer{KS}" primer rhodamine ABI_373_377
% 7 "DyePrimer{M13RP1}" primer rhodamine ABI_373_377
% 8 "DyePrimer{SK}" primer rhodamine ABI_373_377
% 9 "DyePrimer{SP6}" primer rhodamine ABI_373_377
% 10 "DyePrimer{T3}" primer rhodamine ABI_373_377
% 11 "DyePrimer{T7}" primer rhodamine ABI_373_377
% 12 "DyePrimer{-21M13LR}" primer rhodamine ABI_373_377
% 13 "DP5%LR{HS}" primer rhodamine ABI_373_377
% 14 "DP4%25Ac{-21M13}" primer rhodamine ABI_373_377
% 15 "DP4%Ac{-21M13}" primer rhodamine ABI_373_377
% 16 "DP4%25Acv2{M13Rev}" primer rhodamine ABI_373_377
% 17 "DP4%Acv2{M13Rev}" primer rhodamine ABI_373_377
% 18 "DP4%Ac{KS}" primer rhodamine ABI_373_377
% 19 "DP4%Ac{M13Rev}" primer rhodamine ABI_373_377
% 20 "DP4%Ac{SK}" primer rhodamine ABI_373_377
% 21 "DP4%Ac{SP6}" primer rhodamine ABI_373_377
% 22 "DP4%Ac{T3}" primer rhodamine ABI_373_377
% 23 "DP4%Ac{T7}" primer rhodamine ABI_373_377
% 24 "DyePrimer48cm{-21M13}" primer rhodamine ABI_373_377
% 25 "DP5%CEHV(KS)" primer rhodamine ABI_373_377
% 26 "DP5%CEHV(SK)" primer rhodamine ABI_373_377
% 27 "DP5%LR(KS)" primer rhodamine ABI_373_377
% 28 "DP5%LR(SK)" primer rhodamine ABI_373_377
% 29 "DP6%Ac{SP6}" primer rhodamine ABI_373_377
% 30 "DyePrimer{RP1LR}" primer rhodamine ABI_373_377
% 31 "P+DyePrimer{-21m13}" primer rhodamine ABI_373_377
% 32 "mob.REGA+ET" primer energy-transfer ABI_373_377
% 33 "DP4%Ac{ETPrimer}" primer energy-transfer ABI_373_377
% 34 "DyePrimer48cm{RP1}" primer energy-transfer ABI_373_377
% 35 "ET{-28m13rev}960221" primer energy-transfer ABI_373_377
% 36 "ET{-28m13rev}" primer energy-transfer ABI_373_377
% 37 "ET{-40m13}" primer energy-transfer ABI_373_377
% 38 "ET{-40m13F}" primer energy-transfer ABI_373_377
% 39 "ET{-40m13}960208" primer energy-transfer ABI_373_377
% 40 "ET{-40m13}longreader" primer energy-transfer ABI_373_377
% 41 "ET{-21m13fwd}longreader" primer energy-transfer ABI_373_377
% 42 "ET{-21m13fwd}960909" primer energy-transfer ABI_373_377
% 43 "ET{-21m13fwd}" primer energy-transfer ABI_373_377
% 44 "ET{-21m13}fwd longreader" primer energy-transfer ABI_373_377
% 45 "mob377longrngrxlET" primer energy-transfer ABI_373_377
% 46 "mob.REGA+ET" primer energy-transfer ABI_373_377
% 47 "ET{40m13}960208" primer energy-transfer ABI_373_377
% 48 "DP5%LR{BD_M13_FWD_&_REV}" primer big-dye ABI_373_377
% 49 "DP5%LR{BD M13 FWD & REV}" primer big-dye ABI_373_377
% 50 "ET Primers" primer energy-transfer MolDyn_MegaBACE
% 51 "DT6%Ac{A Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 52 "DT6%Ac{A_Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 53 "DT6%25Ac{A_Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 54 "DyeTerm{T7}-Set B" terminator rhodamine ABI_373_377
% 55 "DyeTerm{T7}-Set_B" terminator rhodamine ABI_373_377
% 56 "FS_Dye_Terminator{AnyPrimer}" terminator rhodamine ABI_373_377
% 57 "FS Dye Terminator{AnyPrimer}" terminator rhodamine ABI_373_377
% 58 "DT4%25Ac{A Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 59 "DT4%Ac{A Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 60 "DT4%25Ac{A_Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 61 "DT4%Ac{A_Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 62 "DyeTerminator{AnyPrimer}" terminator rhodamine ABI_373_377
% 63 "DT4%25Ac{B Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 64 "DT4%Ac{B Set-AnyPrimer}" terminator rhodamine ABI_373_377
% 65 "DT_{dR_Set_Any-Primer}" terminator d-rhodamine ABI_373_377
% 66 "DT {dR Set Any-Primer}" terminator d-rhodamine ABI_373_377
% 67 "DT_{dRhod_Set_Any-Primer}" terminator d-rhodamine ABI_373_377
% 68 "DT DSP{dR Set-AnyPrimer}" terminator d-rhodamine ABI_373_377
% 69 "DT_{BD_Set_Any-Primer}" terminator big-dye ABI_373_377
% 70 "DT{BD_Set_Any-Primer}" terminator big-dye ABI_373_377
% 71 "DT {BD Set Any-Primer}" terminator big-dye ABI_373_377
% 72 "DT373{BD Set-AnyPrimer}" terminator big-dye ABI_373_377
% 73 "DT_{BigDye_Set_Any-Primer}" terminator big-dye ABI_373_377
% 74 "DT_BigDye_Set_Any-Primer" terminator big-dye ABI_373_377
% 75 "373 BDT" terminator big-dye ABI_373_377
% 76 "373_BDT" terminator big-dye ABI_373_377
% 77 "ET Terminators" terminator energy-transfer MolDyn_MegaBACE
% 78 "MegaBACE Mobility File" unknown energy-transfer MolDyn_MegaBACE
% 79 "DT3700POP5{BD}v2.mob" terminator big-dye ABI_3700
% 80 "DP3700POP6{BD-M13Rev}v1.mob" primer big-dye ABI_3700
% 81 "DT3700POP6{BD}v0_3.mob" terminator big-dye ABI_3700
% 82 "DT3700POP6{BD-MIT}v1.mob" terminator big-dye ABI_3700
% 83 "DT3700POP6{BD}v4.mob" terminator big-dye ABI_3700
% 84 "DT3700PO.MOB" terminator big-dye ABI_3700
% by Lucio A. Cetto (last modification April 2005)
validids ={...
'DP6%25Ac{-21M13}', 'DP6%Ac{-21M13}', 'DP6%25Ac{M13Rev}', 'DP6%Ac{M13Rev}',...
'DyePrimer{-21m13}', 'DyePrimer{KS}', 'DyePrimer{M13RP1}', 'DyePrimer{SK}',...
'DyePrimer{SP6}', 'DyePrimer{T3}', 'DyePrimer{T7}', 'DyePrimer{-21M13LR}',...
'DP5%LR{HS}', 'DP4%25Ac{-21M13}','DP4%Ac{-21M13}', 'DP4%25Acv2{M13Rev}',...
'DP4%Acv2{M13Rev}', 'DP4%Ac{KS}', 'DP4%Ac{M13Rev}', 'DP4%Ac{SK}',...
'DP4%Ac{SP6}', 'DP4%Ac{T3}', 'DP4%Ac{T7}', 'DyePrimer48cm{-21M13}',...
'DP5%CEHV(KS)', 'DP5%CEHV(SK)', 'DP5%LR(KS)', 'DP5%LR(SK)',...
'DP6%Ac{SP6}', 'DyePrimer{RP1LR}','P+DyePrimer{-21m13}','mob.REGA+ET',...
'DP4%Ac{ETPrimer}', 'DyePrimer48cm{RP1}','ET{-28m13rev}960221','ET{-28m13rev}',...
'ET{-40m13}', 'ET{-40m13F}', 'ET{-40m13}960208', 'ET{-40m13}longreader',...
'ET{-21m13fwd}longreader','ET{-21m13fwd}960909','ET{-21m13fwd}','ET{-21m13}fwd longreader',...
'mob377longrngrxlET', 'mob.REGA+ET', 'ET{40m13}960208', 'DP5%LR{BD_M13_FWD_&_REV}',...
'DP5%LR{BD M13 FWD & REV}','ET Primers', 'DT6%Ac{A Set-AnyPrimer}','DT6%Ac{A_Set-AnyPrimer}',...
'DT6%25Ac{A_Set-AnyPrimer}','DyeTerm{T7}-Set B','DyeTerm{T7}-Set_B','FS_Dye_Terminator{AnyPrimer}',...
'FS Dye Terminator{AnyPrimer}','DT4%25Ac{A Set-AnyPrimer}','DT4%Ac{A Set-AnyPrimer}','DT4%25Ac{A_Set-AnyPrimer}',...
'DT4%Ac{A_Set-AnyPrimer}','DyeTerminator{AnyPrimer}','DT4%25Ac{B Set-AnyPrimer}','DT4%Ac{B Set-AnyPrimer}',...
'DT_{dR_Set_Any-Primer}','DT {dR Set Any-Primer}','DT_{dRhod_Set_Any-Primer}','DT DSP{dR Set-AnyPrimer}',...
'DT_{BD_Set_Any-Primer}','DT{BD_Set_Any-Primer}','DT {BD Set Any-Primer}','DT373{BD Set-AnyPrimer}',...
'DT_{BigDye_Set_Any-Primer}','DT_BigDye_Set_Any-Primer','373 BDT','373_BDT',...
'ET Terminators','MegaBACE Mobility File','DT3700POP5{BD}v2.mob','DP3700POP6{BD-M13Rev}v1.mob',...
'DT3700POP6{BD}v0_3.mob','DT3700POP6{BD-MIT}v1.mob','DT3700POP6{BD}v4.mob','DT3700PO.MOB'};
if ischar(primeridstring)
h = strmatch(str,validids);
if numel(h)>1 error('Ambiguous primer ID string'); end
if numel(h)==0 h=0; end
else
h = primeridstring;
end
if ~isnumeric(h) | h<1 | h>84
error('Can not identify the primer ID string')
end
comments = ['DYEP=' validids{h} 10];
disp(['Primer ID set to ' comments(1:end-1)])
disp(['Writting ' num2str(length(dat)) ' samples to ' filename])
magic_number = (('.'*256+'s')*256+'c')*256+'f'; % magic number
num_samples = length(dat); % Number of elements in Samples matrix
sample_size = 2; % Size of samples in bytes 1=8bits, 2=16bits
samples_offset = 128; % Header = 128
num_bases = 0; % Number of bases in Bases matrix
bases_left_clip = 0; % OBSOLETE: No. bases in left clip (vector)
bases_right_clip = 0; % OBSOLETE: No. bases in right clip (qual)
bases_offset = 4 * sample_size * num_samples + samples_offset;
% Byte offset from start of file
comments_size = length(comments); % Number of bytes in Comment section
comments_offset = bases_offset + 12 * num_bases; % Byte offset from start of file
version = '3.00'; % "version.revision", eg '3' '.' '0' '0'
code_set = 0; % code set used (but ignored!)
private_size = 0; % No. of bytes of Private data, 0 if none
private_offset = comments_offset + private_size; % Byte offset from start of file
spare = zeros(18,1); % Unused
fid = fopen(filename,'w','b'); % force ieee big endian (is what phred reads)
fwrite(fid,magic_number,'uint32');
fwrite(fid,num_samples,'uint32');
fwrite(fid,samples_offset,'uint32');
fwrite(fid,num_bases,'uint32');
fwrite(fid,bases_left_clip,'uint32');
fwrite(fid,bases_right_clip,'uint32');
fwrite(fid,bases_offset,'uint32');
fwrite(fid,comments_size,'uint32');
fwrite(fid,comments_offset,'uint32');
fwrite(fid,version,'char');
fwrite(fid,sample_size,'uint32');
fwrite(fid,code_set,'uint32');
fwrite(fid,private_size,'uint32');
fwrite(fid,private_offset,'uint32');
fwrite(fid,spare,'uint32');
%remap dat to the scf order (i.e. ACGT)
% we use for dnadata GACT
dat = dat(:, [2 3 1 4]);
dat = min(10,dat); % saturates at 10 times the normalized value, so phred does not crash
dat = max(0,dat);
ma = max(max(dat))
% prepare data to write
dat(:,1) = dat(:,1) / ma;
dat(:,2) = dat(:,2) / ma;
dat(:,3) = dat(:,3) / ma;
dat(:,4) = dat(:,4) / ma;
dat = round(dat * 32000);
deltadelta=[0 0 0 0;diff([0 0 0 0;diff(dat)])];
deltadelta = max(-32768,deltadelta);
deltadelta = min(32767,deltadelta);
% write sample data
fwrite(fid,deltadelta,'int16'); % sample_size = 2;
% write comments
fwrite(fid,comments,'char');
fclose(fid);