infer codon flag from input

zqfang · Aug 22, 2024 · dc0e2cc · dc0e2cc
1 parent c3a9ea4
commit dc0e2cc
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 9 deletions.
diff --git a/haplomap/src/ghmap.cpp b/haplomap/src/ghmap.cpp
@@ -14,6 +14,7 @@ std::unordered_map<std::string, GeneSummary *> geneTable; // for gene-oriented i
 std::vector<BlockSummary *> blocks;                       // global vector of all blocks.
 std::vector<std::string> 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<std::string, int> PRIOR // from constants.h;
 // std::unordered_map<std::string, int> 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
@@ -73,9 +73,8 @@ GhmapOptions *parseGhmapOptions(int argc, char **argv)
                         "                                Output gene-summaried results by default.\n"
                         "\nOptional arguments:\n"
                         "    -n, --name             <name of phenotype dataset> \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         <genetic relation file> \n" 
                         "                                <.rel> file for population structure analysis.\n"
@@ -85,9 +84,9 @@ GhmapOptions *parseGhmapOptions(int argc, char **argv)
                         ///"    -t, --goterms          <name of file>\n"
                         ///"    -i, --goterms_include  <name of file> \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"