.abi / .ab1 reader

#pragma rtGlobals=1		// Use modern global access method.


// staden.sourceforge.net/scf-rfc.html


Structure SCFHeader
	uint32 magic_number		// Number of elements in Samples matrix
	uint32 samples				// Byte offset from start of file
	uint32 samples_offset		// Number of bases in Bases matrix
	uint32 bases				// OBSOLETE: No. bases in left clip (vector)
	uint32 bases_left_clip		// OBSOLETE: No. bases in right clip (qual)
	uint32 bases_right_clip	// Byte offset from start of file
	uint32 bases_offset			// Number of bytes in Comment section
	uint32 comments_size		// Byte offset from start of file
	uint32 comments_offset		// "version.revision", eg '3' '.' '0' '0'
	char version[4]
	uint32 sample_size
	uint32 code_set
	uint32 private_size
	uint32 private_offset
	uint32 spare[18]
EndStructure

Structure SequenceEntry
	uint32 peak_index
	uchar prob_A
	uchar prob_C
	uchar prob_G
	uchar prob_T
	char   base
	uchar prob_sub
	uchar prob_ins
	uchar prob_del
EndStructure

Function ReadSCFHeader()
	
	variable refNum, i
	
	Open /D/R /F="Chromatograms:.scf" refNum
	if (strlen(S_fileName) == 0)
		return 0
	endif
	
	Open /R refNum as S_fileName
	
	Struct SCFHeader header
	
	string magicString = "1234"
	variable bigEndian
	
	FBinRead refNum, magicString
	if (StringMatch(magicString, ".scf") == 1)
		// the format is big endian
		bigEndian = 1
	else
		bigEndian = 0
	endif
	
	FSetPos refNum, 0
	
	if (bigEndian == 1)	// big-endian
		FBinRead /B=2 refNum, header
	else
		FBinRead /B=3 refNum, header
	endif
	
	if (header.sample_size != 2)
		Close refNum
		Abort "The sample size is not 2"
	endif
	
	Make /W/O/N=(header.samples) ChromatogramA, ChromatogramC, ChromatogramG, ChromatogramT
	
	FSetPos refNum, header.samples_offset
	
	// read the raw chromatogram data
	if (bigEndian == 1)
		FBinRead /B=2 refNum, ChromatogramA
		FBinRead /B=2 refNum, ChromatogramC
		FBinRead /B=2 refNum, ChromatogramG
		FBinRead /B=2 refNum, ChromatogramT
	else
		FBinRead /B=3 refNum, ChromatogramA
		FBinRead /B=3 refNum, ChromatogramC
		FBinRead /B=3 refNum, ChromatogramG
		FBinRead /B=3 refNum, ChromatogramT
	endif
	
	ConvertSamplesToAbsoluteValues(ChromatogramA)
	ConvertSamplesToAbsoluteValues(ChromatogramC)
	ConvertSamplesToAbsoluteValues(ChromatogramG)
	ConvertSamplesToAbsoluteValues(ChromatogramT)
	
	// read the base assignments
	Make /T/O/N=(header.bases) BaseAssignments
	Make /O/N=(header.bases) BasePositions
	
	Struct SequenceEntry sequence
	
	FSetPos refNum, header.bases_offset
	
	for (i = 0; i < header.bases; i+=1)
		if (bigEndian == 1)
			FBinRead /B=2 refNum, sequence
		else
			FBinRead /B=3 refNum, sequence
		endif
		BaseAssignments[i] = GetMostLikelyBase(sequence)
		BasePositions[i] = sequence.peak_index
	endfor
	
	Close refNum
	
	// Display the chromatogram
	DoWindow /F ChromatogramViewer
	if (V_flag != 1)
		Display /N=ChromatogramViewer as "Chromatogram"
		AppendToGraph ChromatogramA, ChromatogramC, ChromatogramT, ChromatogramG
		ModifyGraph rgb(ChromatogramA)=(0,65535,0),rgb(ChromatogramC)=(0,0,65535);DelayUpdate
		ModifyGraph rgb(ChromatogramG)=(0,0,0)
	endif
	
End

Function ConvertSamplesToAbsoluteValues(samples)
	wave samples
	
	variable i, p_sample = 0
	variable nSamples = DimSize(samples, 0)
	
	for (i=0; i < nSamples ;i+=1)
		samples[i] = samples[i] + p_sample
		p_sample = samples[i]
	endfor
	p_sample = 0
	for (i=0; i < nSamples ;i+=1)
		samples[i] = samples[i] + p_sample
		p_sample = samples[i]
	endfor
End
	
Function /S GetMostLikelyBase(sequence)
	Struct SequenceEntry &sequence
	
	variable probA, probC, probT, probG
	
	probA = sequence.prob_A
	probT = sequence.prob_T
	probC = sequence.prob_C
	probG = sequence.prob_G
	
	variable highestProbability = max(max(probA, probT), max(probC, probG))
	
	if (highestProbability == probA)
		return "A"
	elseif (highestProbability == probT)
		return "T"
	elseif (highestProbability == probC)
		return "C"
	elseif (highestProbability == probG)
		return "G"
	endif
End