range_fdr.cpp

降低精度的快速FDR算法

// Range_FDR.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

long double Range_FDR(wchar_t*	negLogP_filename, long double totalCount, double threshold);

int _tmain(int argc, _TCHAR* argv[])
{
	long double	FDR_pValue;
	ofstream	outfile;

	FDR_pValue = Range_FDR(L"Data/SNPvsPS_negLogP.rangeCount_LM2.txt", 1797942727, 0.05);
	cout<<"FDR pValue: "<<FDR_pValue<<endl;
	outfile.open(L"Result/SNPvsPS_LM2_FDR.txt");
	outfile<<FDR_pValue<<endl;
	outfile.close();

	system("pause");
	return 0;
}

long double Range_FDR(wchar_t*	negLogP_filename, long double totalCount, double threshold)
{
	ifstream	infile;
	string		buffer, pRange;
	int			index, pCount, pPos;
	long double		pRangeStart, pRangeStop, pValue, qValue;

	infile.open(negLogP_filename);
	// read title
	getline(infile, buffer);
	pPos = 0;
	// read lines
	while (!infile.eof())
	{
		// 1st string process
		getline(infile, buffer);
		index = buffer.find('	');
		pRange = buffer.substr(0,index);
		pCount = atoi(buffer.substr(index+1, buffer.length()-index-1).c_str());
		// 2nd string process
		index = pRange.find('-');
		pRangeStart = atof(pRange.substr(0,index).c_str());
		pRangeStop = atof(pRange.substr(index+1, pRange.length()-index-1).c_str());
		// calculate FDR cut
		pValue = pow((long double)10.0, -pRangeStart);
		qValue = (long double)(totalCount-pPos)*threshold/(long double)totalCount;
		//cout<<totalCount<<"-"<<pPos<<endl;
		//cout<<qValue<<"-"<<pValue<<"="<<qValue-pValue<<endl;
		if ( pValue < qValue )
			return pValue;
		pPos += pCount;
	}
	cout<<pPos<<endl;
	infile.close();
	return -1;
}