SamTags.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "SamTags.h"
00019 #include "BaseUtilities.h"
00020
00021 const char* SamTags::BQ_TAG = "BQ";
00022 const char SamTags::BQ_TAG_TYPE = 'Z';
00023 const char* SamTags::MD_TAG = "MD";
00024 const char SamTags::MD_TAG_TYPE = 'Z';
00025 const char* SamTags::ORIG_POS_TAG = "OP";
00026 const char SamTags::ORIG_POS_TAG_TYPE = 'i';
00027 const char* SamTags::ORIG_CIGAR_TAG = "OC";
00028 const char SamTags::ORIG_CIGAR_TAG_TYPE = 'Z';
00029 const char* SamTags::ORIG_QUAL_TAG = "OQ";
00030 const char SamTags::ORIG_QUAL_TAG_TYPE = 'Z';
00031
00032
00033
00034 bool SamTags::createMDTag(String& outputMDtag, SamRecord& inputRec,
00035 GenomeSequence& genome)
00036 {
00037 outputMDtag.Clear();
00038
00039 Cigar* cigarInfo = inputRec.getCigarInfo();
00040 if(cigarInfo == NULL)
00041 {
00042 throw(std::runtime_error("Cannot createMDTag - failed to get the cigar"));
00043 return(false);
00044 }
00045 int32_t queryIndex = Cigar::INDEX_NA;
00046
00047
00048 uint32_t startOfReadOnRefIndex =
00049 genome.getGenomePosition(inputRec.getReferenceName());
00050 if(startOfReadOnRefIndex == (uint32_t)INVALID_CHROMOSOME_INDEX)
00051 {
00052
00053 return(false);
00054 }
00055 startOfReadOnRefIndex += inputRec.get0BasedPosition();
00056
00057
00058 int32_t matchCount = 0;
00059
00060
00061 bool currentDeletion = false;
00062
00063
00064 for(int refOffset = 0;
00065 refOffset < cigarInfo->getExpectedReferenceBaseCount();
00066 ++refOffset)
00067 {
00068
00069 queryIndex = cigarInfo->getQueryIndex(refOffset);
00070
00071 char refBase = genome[startOfReadOnRefIndex + refOffset];
00072
00073 if(queryIndex != Cigar::INDEX_NA)
00074 {
00075
00076 char readBase = inputRec.getSequence(queryIndex);
00077 currentDeletion = false;
00078
00079
00080
00081 if(!BaseUtilities::isAmbiguous(readBase) &&
00082 !BaseUtilities::isAmbiguous(refBase) &&
00083 (BaseUtilities::areEqual(readBase, refBase)))
00084 {
00085
00086 ++matchCount;
00087 }
00088 else
00089 {
00090
00091 if(matchCount != 0)
00092 {
00093 outputMDtag += matchCount;
00094 matchCount = 0;
00095 }
00096 outputMDtag += refBase;
00097 }
00098 }
00099 else
00100 {
00101
00102
00103 if(matchCount != 0)
00104 {
00105 outputMDtag += matchCount;
00106 matchCount = 0;
00107 }
00108
00109 if(!currentDeletion)
00110 {
00111
00112 outputMDtag += '^';
00113 }
00114
00115 outputMDtag += refBase;
00116 currentDeletion = true;
00117 }
00118 }
00119
00120
00121 outputMDtag += matchCount;
00122 return(true);
00123 }
00124
00125
00126 bool SamTags::isMDTagCorrect(SamRecord& inputRec, GenomeSequence& genome)
00127 {
00128 String calcMDtag;
00129 if(!createMDTag(calcMDtag, inputRec, genome))
00130 {
00131
00132 return(false);
00133 }
00134
00135 String* origMDtag = inputRec.getStringTag(MD_TAG);
00136
00137 if(origMDtag == NULL)
00138 {
00139
00140
00141
00142
00143
00144 return(calcMDtag.IsEmpty());
00145 }
00146 else
00147 {
00148
00149 return(calcMDtag == *origMDtag);
00150 }
00151 }
00152
00153
00154
00155 bool SamTags::updateMDTag(SamRecord& inputRec, GenomeSequence& genome)
00156 {
00157
00158 String calcMDtag;
00159 createMDTag(calcMDtag, inputRec, genome);
00160
00161
00162
00163
00164 return(inputRec.addTag(MD_TAG, MD_TAG_TYPE, calcMDtag.c_str()));
00165 }