1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
|
#include "ReadAlign.h"
#include "SequenceFuns.h"
#include "ErrorWarning.h"
uint ReadAlign::outputTranscriptSAM(Transcript const &trOut, uint nTrOut, uint iTrOut, uint mateChr, uint mateStart, char mateStrand, int unmapType, bool *mateMap, ostream *outStream) {
if (P.outSAMmode=="None") return 0; //no SAM output
uint outStreamPos0=(uint)outStream->tellp();
if (unmapType>=0)
{//unmapped reads: SAM
for (uint imate=0;imate<P.readNmates;imate++) {//cycle over mates //not readNends: this is alignment
if (!mateMap[imate]) {
uint16 samFLAG=0x4;
if (P.readNmates==2) {//paired read //not readNends: this is alignment
samFLAG|=0x1 + (imate==0 ? 0x40 : 0x80);
if (mateMap[1-imate]) {//mate mapped
if ( trOut.Str != (1-imate) ) {
samFLAG|=0x20;//mate strand reverted
};
} else {//mate unmapped
samFLAG|=0x8;
};
};
if (readFilter=='Y')
samFLAG|=0x200; //not passing quality control
if (mateMap[1-imate] && !trOut.primaryFlag && P.outSAMunmapped.keepPairs) {//mapped mate is not primary, keep unmapped mate for each pair, hence need to mark some as not primary
samFLAG|=0x100;
};
*outStream << readName+1 <<"\t"<< samFLAG \
<<"\t"<< '*' <<"\t"<< '0' <<"\t"<< '0' <<"\t"<< '*';
if (mateMap[1-imate]) {//mate is mapped
*outStream <<"\t"<< genOut.chrName[trOut.Chr] <<"\t"<< trOut.exons[0][EX_G] + 1 - genOut.chrStart[trOut.Chr];
} else {
*outStream <<"\t"<< '*' <<"\t"<< '0';
};
*outStream <<"\t"<< '0' <<"\t"<< Read0[imate] <<"\t"<< (readFileType==2 ? Qual0[imate]:"*") \
<<"\tNH:i:0" <<"\tHI:i:0" <<"\tAS:i:"<<trOut.maxScore <<"\tnM:i:"<<trOut.nMM<<"\tuT:A:" <<unmapType;
if (!P.outSAMattrRG.empty()) *outStream<< "\tRG:Z:" <<P.outSAMattrRG.at(readFilesIndex);
if (P.readFilesTypeN==10 && !readNameExtra[imate].empty()) {//SAM files as input - output extra attributes
*outStream << "\t" <<readNameExtra[imate];
};
*outStream <<"\n";
};
};
return (uint)outStream->tellp()-outStreamPos0;
};//if (unmapType>=0 && outStream != NULL) //unmapped reads: SAM
bool flagPaired = P.readNmates==2; //not readNends: this is alignment
string CIGAR;
//for SAM output need to split mates
uint iExMate; //last exon of the first mate
uint nMates=1;
for (iExMate=0;iExMate<trOut.nExons-1;iExMate++) {
if (trOut.canonSJ[iExMate]==-3){
nMates=2;
break;
};
};
uint samFlagCommon=0;//FLAG common for both mates
if (flagPaired)
{//paired reads
samFlagCommon=0x0001;
if (iExMate==trOut.nExons-1)
{//single mate
if (mateChr>genOut.nChrReal) samFlagCommon+=0x0008; //not mapped as pair
} else
{//paired align
if (P.alignEndsProtrude.concordantPair || \
( (trOut.exons[0][EX_G] <= trOut.exons[iExMate+1][EX_G]+trOut.exons[0][EX_R]) && \
(trOut.exons[iExMate][EX_G]+trOut.exons[iExMate][EX_L] <= trOut.exons[trOut.nExons-1][EX_G]+Lread-trOut.exons[trOut.nExons-1][EX_R]) ) )
{//properly paired
samFlagCommon+=0x0002;
};
};
} else
{//single end
samFlagCommon=0;
};
if (readFilter=='Y') samFlagCommon+=0x200; //not passing quality control
uint Str= trOut.Str;//note that Strand = the mate on the left
uint leftMate=0; //the mate (0 or 1) which is on the left
if (flagPaired) {
leftMate=Str;
};
if (P.outSAMattrPresent.MC) {
calcCIGAR(trOut, nMates, iExMate, leftMate);
};
uint samFLAG;
for (uint imate=0;imate<nMates;imate++) {
samFLAG=samFlagCommon;
uint iEx1 = (imate==0 ? 0 : iExMate+1);
uint iEx2 = (imate==0 ? iExMate : trOut.nExons-1);
uint Mate=trOut.exons[iEx1][EX_iFrag];
if (Mate==0) {
samFLAG|= Str*0x10;
if (nMates==2) samFLAG|= (1-Str)*0x20;
} else {//second mate strand need to be reverted
samFLAG|= (1-Str)*0x10;
if (nMates==2) samFLAG|= Str*0x20;
};
if (flagPaired) {
samFLAG|= (Mate==0 ? 0x0040 : 0x0080);
if (flagPaired && nMates==1 && mateStrand==1) samFLAG|=0x20;//revert strand using inout value of mateStrand (e.g. for chimeric aligns)
};
//not primary align?
if (!trOut.primaryFlag) samFLAG|=0x100;
//empty streams
samStreamCIGAR.str(std::string());
samStreamSJmotif.str(std::string());
samStreamSJintron.str(std::string());
// samStreamSJannot.str(std::string());
uint trimL;
if (Str==0 && Mate==0) {
trimL=clipMates[Mate][0].clippedN;
} else if (Str==0 && Mate==1) {
trimL=clipMates[Mate][1].clippedN;
} else if (Str==1 && Mate==0) {
trimL=clipMates[Mate][1].clippedN;
} else {
trimL=clipMates[Mate][0].clippedN;
};
uint trimL1 = trimL + trOut.exons[iEx1][EX_R] - (trOut.exons[iEx1][EX_R]<readLength[leftMate] ? 0 : readLength[leftMate]+1);
if (trimL1>0) {
samStreamCIGAR << trimL1 << "S"; //initial trimming
};
for (uint ii=iEx1;ii<=iEx2;ii++) {
if (ii>iEx1) {//record gaps
uint gapG=trOut.exons[ii][EX_G]-(trOut.exons[ii-1][EX_G]+trOut.exons[ii-1][EX_L]);
uint gapR=trOut.exons[ii][EX_R]-trOut.exons[ii-1][EX_R]-trOut.exons[ii-1][EX_L];
//it's possible to have a D or N and I at the same time
if (gapR>0){
samStreamCIGAR << gapR;
samStreamCIGAR << "I";
};
if (trOut.canonSJ[ii-1]>=0 || trOut.sjAnnot[ii-1]==1) {//junction: N
samStreamCIGAR << gapG;
samStreamCIGAR << "N";
samStreamSJmotif <<','<< trOut.canonSJ[ii-1] + (trOut.sjAnnot[ii-1]==0 ? 0 : SJ_SAM_AnnotatedMotifShift); //record junction type
// samStreamSJannot <<','<< (int) trOut.sjAnnot[ii-1]; //record annotation type
samStreamSJintron <<','<< trOut.exons[ii-1][EX_G] + trOut.exons[ii-1][EX_L] + 1 - genOut.chrStart[trOut.Chr] <<','\
<< trOut.exons[ii][EX_G] - genOut.chrStart[trOut.Chr]; //record intron loci
} else if (gapG>0) {//deletion: N
samStreamCIGAR << gapG;
samStreamCIGAR << "D";
};
};
samStreamCIGAR << trOut.exons[ii][EX_L] << "M";
};
string SJmotif = samStreamSJmotif.str();
string SJintron = samStreamSJintron.str();
// string SJannot = samStreamSJannot.str();
if (SJmotif.length()==0) {//no junctions recorded, mark with -1
SJmotif=",-1";
SJintron=",-1";
// SJannot=",-1";
};
uint trimR1=(trOut.exons[iEx1][EX_R]<readLength[leftMate] ? \
readLengthOriginal[leftMate] : readLength[leftMate]+1+readLengthOriginal[Mate]) \
- trOut.exons[iEx2][EX_R]-trOut.exons[iEx2][EX_L] - trimL;
if ( trimR1 > 0 ) {
samStreamCIGAR << trimR1 << "S"; //final trimming
};
CIGAR=samStreamCIGAR.str();
char seqMate[DEF_readSeqLengthMax+1], qualMate[DEF_readSeqLengthMax+1];
char *seqOut=NULL, *qualOut=NULL;
if ( Mate==Str ) {//seq strand is correct
seqOut=Read0[Mate];
qualOut=Qual0[Mate];
} else {
revComplementNucleotides(Read0[Mate], seqMate, readLengthOriginal[Mate]);
seqMate[readLengthOriginal[Mate]]=0;
for (uint ii=0;ii<readLengthOriginal[Mate]; ii++) qualMate[ii]=Qual0[Mate][readLengthOriginal[Mate]-1-ii];
qualMate[readLengthOriginal[Mate]]=0;
seqOut=&seqMate[0];
qualOut=&qualMate[0];
};
// return;
int MAPQ=P.outSAMmapqUnique;
if (nTrOut>=5) {
MAPQ=0;
} else if (nTrOut>=3) {
MAPQ=1;
} else if (nTrOut==2) {
MAPQ=3;
};
*outStream << readName+1 <<"\t"<< ((samFLAG & P.outSAMflagAND) | P.outSAMflagOR) <<"\t"<< genOut.chrName[trOut.Chr] <<"\t"<< trOut.exons[iEx1][EX_G] + 1 - genOut.chrStart[trOut.Chr]
<<"\t"<< MAPQ <<"\t"<< CIGAR;
if (nMates>1) {
*outStream <<"\t"<< "=" <<"\t"<< trOut.exons[(imate==0 ? iExMate+1 : 0)][EX_G]+ 1 - genOut.chrStart[trOut.Chr]
<<"\t"<< (imate==0? "":"-") << trOut.exons[trOut.nExons-1][EX_G]+trOut.exons[trOut.nExons-1][EX_L]-trOut.exons[0][EX_G];
} else if (mateChr<genOut.nChrReal){//mateChr is given in the function parameters
*outStream <<"\t"<< genOut.chrName[mateChr] <<"\t"<< mateStart+1-genOut.chrStart[mateChr] <<"\t"<< 0;
} else {
*outStream <<"\t"<< "*" <<"\t"<< 0 <<"\t"<< 0;
};
*outStream <<"\t"<< seqOut;
if (readFileType==2 && P.outSAMmode != "NoQS") {//fastq
*outStream <<"\t"<< qualOut ;
} else {
*outStream <<"\t"<< "*";
};
// vector<string> customAttr(outSAMattrN,"");
uint tagNM=0;
string tagMD("");
if (P.outSAMattrPresent.NM || P.outSAMattrPresent.MD) {
char* R=Read1[trOut.roStr==0 ? 0:2];
uint matchN=0;
for (uint iex=iEx1;iex<=iEx2;iex++) {
for (uint ii=0;ii<trOut.exons[iex][EX_L];ii++) {
char r1 = R[ii+trOut.exons[iex][EX_R]];
char g1 = genOut.G[ii+trOut.exons[iex][EX_G]];
if ( r1!=g1 || r1==4 || g1==4) {
++tagNM;
// if (matchN>0 || (ii==0 && iex>0 && trOut.canonSJ[iex]==-1) ) {
tagMD+=to_string(matchN);
// };
tagMD+=P.genomeNumToNT[(uint8) g1];
matchN=0;
} else {
matchN++;
};
};
if (iex<iEx2) {
if (trOut.canonSJ[iex]==-1) {//deletion
tagNM+=trOut.exons[iex+1][EX_G]-(trOut.exons[iex][EX_G]+trOut.exons[iex][EX_L]);
tagMD+=to_string(matchN) + "^";
for (uint ii=trOut.exons[iex][EX_G]+trOut.exons[iex][EX_L];ii<trOut.exons[iex+1][EX_G];ii++) {
tagMD+=P.genomeNumToNT[(uint8) genOut.G[ii]];
};
matchN=0;
} else if (trOut.canonSJ[iex]==-2) {//insertion
tagNM+=trOut.exons[iex+1][EX_R]-trOut.exons[iex][EX_R]-trOut.exons[iex][EX_L];
};
};
};
tagMD+=to_string(matchN);
};
for (uint ii=0;ii<P.outSAMattrOrder.size();ii++) {
switch (P.outSAMattrOrder[ii]) {
case ATTR_NH:
*outStream <<"\tNH:i:" << nTrOut;
break;
case ATTR_HI:
*outStream <<"\tHI:i:"<<iTrOut+P.outSAMattrIHstart;
break;
case ATTR_AS:
*outStream<<"\tAS:i:"<<trOut.maxScore;
break;
case ATTR_nM:
*outStream<<"\tnM:i:"<<trOut.nMM;
break;
case ATTR_jM:
*outStream<<"\tjM:B:c"<< SJmotif;
break;
case ATTR_jI:
*outStream<<"\tjI:B:i"<< SJintron;
break;
case ATTR_XS:
if (trOut.sjMotifStrand==1) {
*outStream<<"\tXS:A:+";
} else if (trOut.sjMotifStrand==2) {
*outStream<<"\tXS:A:-";
};
break;
case ATTR_NM:
*outStream<< "\tNM:i:" <<tagNM;
break;
case ATTR_MD:
*outStream<< "\tMD:Z:" <<tagMD;
break;
case ATTR_RG:
*outStream<< "\tRG:Z:" <<P.outSAMattrRG.at(readFilesIndex);
break;
case ATTR_MC:
if (nMates>1) {
*outStream<< "\tMC:Z:" <<matesCIGAR[1-imate];
};
break;
case ATTR_ha:
if (mapGen.pGe.transform.type==2) {
*outStream<< "\tha:i:" <<trOut.haploType;
};
break;
//do nothing - this attributes only work for BAM output
case ATTR_ch:
case ATTR_CR:
case ATTR_CY:
case ATTR_UR:
case ATTR_UY:
case ATTR_CB:
case ATTR_UB:
case ATTR_sM:
case ATTR_sS:
case ATTR_sQ:
case ATTR_rB:
case ATTR_vG:
case ATTR_vA:
case ATTR_vW:
case ATTR_GX:
case ATTR_GN:
break;
default:
ostringstream errOut;
errOut <<"EXITING because of FATAL error: unknown/unimplemented SAM atrribute (tag): "<<P.outSAMattrOrder[ii] <<"\n";
errOut <<"SOLUTION: contact Alex Dobin at dobin@cshl.edu\n";
exitWithError(errOut.str(), std::cerr, P.inOut->logMain, EXIT_CODE_PARAMETER, P);
};
};
if (P.readFilesTypeN==10 && !readNameExtra[imate].empty()) {//SAM files as input - output extra attributes
*outStream << "\t" << readNameExtra.at(imate);
};
*outStream << "\n"; //done with one SAM line
};//for (uint imate=0;imate<nMates;imate++)
return (uint)outStream->tellp()-outStreamPos0;
};
|
