diff --git a/haplomap/src/ghmap.cpp b/haplomap/src/ghmap.cpp index 4fbbdf1..bf344da 100644 --- a/haplomap/src/ghmap.cpp +++ b/haplomap/src/ghmap.cpp @@ -14,6 +14,7 @@ std::unordered_map geneTable; // for gene-oriented i std::vector blocks; // global vector of all blocks. std::vector geneExprHeader; int traceFStat = false; +int codonFlag = 0; // 0: combined mode (snp+indel+sv), 1: impact score only, 2: SNP codon score only // std::unordered_map PRIOR // from constants.h; // std::unordered_map CSQs // from constants.h @@ -143,6 +144,8 @@ int BlockSummary::updateCodonScore(std::string str) if (pos == std::string::npos) break; endpos = str.find(">", pos + 1); + // mark and update codon type + codonFlag = 2; // means this annotation is SNP int aa1 = str[pos - 1] - 'A'; int aa2 = str[endpos + 5] - 'A'; @@ -205,6 +208,8 @@ int BlockSummary::updateCodonScore(std::string str) tokstart = dpos + 1; // don't care if it overflows } } + // annotation type is only impact score + codonFlag = std::max(codonFlag, 1); return flag_max; } @@ -375,7 +380,7 @@ GhmapWriter::GhmapWriter(char *outputFileName, char *datasetName, bool categoric os = std::ofstream(outputFileName); if (!os.is_open()) { - std::cout << "Open of file \"" << outputFileName << "\" failed: "; + std::cerr << "Open of file \"" << outputFileName << "\" failed: "; std::exit(1); } os << "##" << _dataset_name << "\t" @@ -384,11 +389,11 @@ GhmapWriter::GhmapWriter(char *outputFileName, char *datasetName, bool categoric std::string dn(_dataset_name); upcase(dn); std::string flag; - if ((dn.find("SNP") != std::string::npos) || (dn.find("_SNV") != std::string::npos)) + if ((dn.find("SNP") != std::string::npos) || (dn.find("_SNV") != std::string::npos) || codonFlag == 2) { // Non-Coding -> (INTRONIC,intergenic,5PRIME_UTR,3PRIME_UTR) flag = "##CodonFlag\t3:Stop\t2:Splicing\t1:Non-Synonymous\t0:Synonymous\t-1:Non-Coding"; - } else if ((dn.find("INDEL") != std::string::npos) || (dn.find("_SV") != std::string::npos)) + } else if ((dn.find("INDEL") != std::string::npos) || (dn.find("_SV") != std::string::npos) || codonFlag == 1 ) { flag = "##CodonFlag:\t2:High\t1:Moderate\t0:Low\t-1:Modifier"; } diff --git a/haplomap/src/quantTraitMap.cpp b/haplomap/src/quantTraitMap.cpp index bb6e549..f105ad4 100644 --- a/haplomap/src/quantTraitMap.cpp +++ b/haplomap/src/quantTraitMap.cpp @@ -73,9 +73,8 @@ GhmapOptions *parseGhmapOptions(int argc, char **argv) " Output gene-summaried results by default.\n" "\nOptional arguments:\n" " -n, --name \n" - " To show SNP CodonFlag explicity, please add `_SNP`\n" - " as a suffix to the name, e.g. MPD123_SNP. \n" - " Default: show impact score and SNP CodonFlag.\n" + " NOTE:: Add suffix `_SNP`|`_INDEL`|`_SV`(e.g. MPD123_SNP)\n" + " can help select correct CodonFlag in the output.\n" " -c, --categorical phenotype (-p) is categorical\n" " -r, --relation \n" " <.rel> file for population structure analysis.\n" @@ -85,9 +84,9 @@ GhmapOptions *parseGhmapOptions(int argc, char **argv) ///" -t, --goterms \n" ///" -i, --goterms_include \n" ///" Output only genes with these terms\n" - " -f, --filter_coding Filter out non-coding blocks\n" - " -g, --gene Output gene-summaried results. Default.\n" - " NOTE:: Only write the overlapped halpoblock with best pvalue/Fstat, representing all overlapped blocks. \n" + " -f, --filter_coding Filter out non-coding blocks\n" + " -g, --gene Output gene-summaried results. Default.\n" + " NOTE:: Only write the overlapped halpoblock with best pvalue/Fstat, representing all overlapped blocks. \n" " The CodonFlag is an aggregated indicator showing that a gene has blocks with coding change. \n" " The best block itself might not contain any coding changes. \n" " Run the ghmap with -a/k/m tag will give you all overlapped blocks with correct CodonFlag. \n"