1 #include "ReadAlign.h" 2 #include "SequenceFuns.h" 3 outputTranscriptCIGARp(Transcript const & trOut)4string ReadAlign::outputTranscriptCIGARp(Transcript const &trOut) {//generates CIGARp string for the transcript 5 //p is a special CIGAR operation to encode gap between mates. This gap is negative for overlapping mates 6 7 string CIGAR; 8 samStreamCIGAR.str(std::string()); 9 10 uint leftMate=0; 11 if (P.readFilesIn.size()>1) leftMate=trOut.Str; 12 13 uint trimL=trOut.exons[0][EX_R] - (trOut.exons[0][EX_R]<readLengthOriginal[leftMate] ? 0 : readLengthOriginal[leftMate]+1); 14 if (trimL>0) { 15 samStreamCIGAR << trimL << "S"; //initial trimming 16 }; 17 18 for (uint ii=0;ii<trOut.nExons;ii++) {//cycle over all exons, record CIGAR 19 if (ii>0) {//record gaps 20 uint gapG=trOut.exons[ii][EX_G]-(trOut.exons[ii-1][EX_G]+trOut.exons[ii-1][EX_L]); 21 22 if (trOut.exons[ii][EX_G] >= (trOut.exons[ii-1][EX_G]+trOut.exons[ii-1][EX_L]) ) {// 23 if (trOut.canonSJ[ii-1]==-3) {//gap between mates 24 //soft clipping of the second mate 25 uint s1=readLengthOriginal[leftMate]-(trOut.exons[ii-1][EX_R]+trOut.exons[ii-1][EX_L]); 26 uint s2=trOut.exons[ii][EX_R]-(readLengthOriginal[leftMate]+1); 27 if (s1>0){ 28 samStreamCIGAR << s1 << "S"; 29 }; 30 samStreamCIGAR << gapG << "p"; 31 if (s2>0){ 32 samStreamCIGAR << s2 << "S"; 33 }; 34 35 } else { 36 //it's possible to have a D or N and I for at the same time 37 uint gapR=trOut.exons[ii][EX_R]-trOut.exons[ii-1][EX_R]-trOut.exons[ii-1][EX_L]; //gapR>0 always 38 if (gapR>0){ 39 samStreamCIGAR << gapR << "I"; 40 }; 41 if (trOut.canonSJ[ii-1]>=0 || trOut.sjAnnot[ii-1]==1) {//junction: N 42 samStreamCIGAR << gapG << "N"; 43 } else if (gapG>0) {//deletion 44 samStreamCIGAR << gapG << "D"; 45 }; 46 }; 47 } else {//mates overlap 48 samStreamCIGAR << "-" << (trOut.exons[ii-1][EX_G]+trOut.exons[ii-1][EX_L]) - trOut.exons[ii][EX_G] << "p"; 49 }; 50 }; 51 samStreamCIGAR << trOut.exons[ii][EX_L] << "M"; 52 }; 53 54 55 trimL=(trOut.exons[trOut.nExons-1][EX_R]<readLengthOriginal[leftMate] ? readLengthOriginal[leftMate] : readLengthPairOriginal) - trOut.exons[trOut.nExons-1][EX_R]-trOut.exons[trOut.nExons-1][EX_L]; 56 if ( trimL > 0 ) { 57 samStreamCIGAR << trimL << "S"; //final trimming 58 }; 59 CIGAR=samStreamCIGAR.str(); 60 61 return CIGAR; 62 63 64 65 }; 66