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
|
#include "ReadAlignChunk.h"
#include "ThreadControl.h"
#include "ErrorWarning.h"
#include "SequenceFuns.h"
#include "GlobalVariables.h"
inline uint64 fastqReadOneLine(ifstream &streamIn, char *arrIn);
inline void removeStringEndControl(string &str);
void ReadAlignChunk::processChunks() {//read-map-write chunks
noReadsLeft=false; //true if there no more reads left in the file
bool newFile=false; //new file marker in the input stream
while (!noReadsLeft) {//continue until the input EOF
//////////////read a chunk from input files and store in memory
if (P.outFilterBySJoutStage<2) {//read chunks from input file
if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexInRead);
chunkInSizeBytesTotal={0,0};
while (chunkInSizeBytesTotal[0] < P.chunkInSizeBytes && chunkInSizeBytesTotal[1] < P.chunkInSizeBytes && P.inOut->readIn[0].good() && P.inOut->readIn[1].good()) {
char nextChar=P.inOut->readIn[0].peek();
if (P.iReadAll==P.readMapNumber) {//do not read any more reads
break;
///////////////////////////////////////////////////////////////////////////////////// SAM
} else if (P.readFilesTypeN==10 && P.inOut->readIn[0].good() && P.outFilterBySJoutStage!=2) {//SAM input && not eof && not 2nd stage
if (nextChar=='@') {//with SAM input linest that start with @ are headers
P.inOut->readIn[0].ignore(DEF_readNameSeqLengthMax,'\n'); //read line and skip it
continue;
};
string str1;
P.inOut->readIn[0] >> str1;
if (str1=="FILE") {
newFile=true;
} else {
P.iReadAll++; //increment read number
uint64 flag1;
P.inOut->readIn[0] >> flag1;
uint imate1=0;
for (uint imate=0;imate<P.readNmates;imate++) {//not readNends: this is SAM input
if (imate>0) {
string str2;
uint64 flag2;
P.inOut->readIn[0] >> str2; //for imate=0 str1 was already read
P.inOut->readIn[0] >> flag2; //read name and flag
if ( str1 != str2 ) {
ostringstream errOut;
errOut << ERROR_OUT <<" EXITING because of FATAL ERROR in input BAM file: the consecutive lines in paired-end BAM have different read IDs:\n"
<< str1 <<" vs "<< str2 << '\n'
<< "\n SOLUTION: fix BAM file formatting. Paired-end reads should be always consecutive lines, with exactly 2 lines per paired-end read" ;
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
if (! ( ((flag1 & 0x40) && (flag2 & 0x80)) || ((flag2 & 0x40) && (flag1 & 0x80)) ) ) {
ostringstream errOut;
errOut << ERROR_OUT <<" EXITING because of FATAL ERROR in input BAM file: the consecutive lines in paired-end BAM have wrong mate FLAG bits:\n"
<< str1 <<" "<< flag1 <<" vs "<< str2 <<" "<< flag2 << '\n'
<< "\n SOLUTION: fix BAM file formatting. Paired-end reads should be always consecutive lines, with exactly 2 lines per paired-end read."
<< " Mate1 should have 0x40 bit set in the FLAG, Mate2 should have 0x80 bit set in the FLAG";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
str1 = str2; //used below for both mates
flag1 = flag2; //used below for both mates
};
char passFilterIllumina=(flag1 & 0x800 ? 'Y' : 'N');
if (imate==1) {//2nd line is always opposite of the 1st one
imate1=1-imate1;
} else if (P.readNmates==2 && (flag1 & 0x80)) {//not readNends: this is SAM input
imate1=1;
} else {
imate1=0;
};
//read ID or number
if (P.outSAMreadID=="Number") {
chunkInSizeBytesTotal[imate1] += sprintf(chunkIn[imate1] + chunkInSizeBytesTotal[imate1], "@%llu", P.iReadAll);
} else {
chunkInSizeBytesTotal[imate1] += sprintf(chunkIn[imate1] + chunkInSizeBytesTotal[imate1], "@%s", str1.c_str());
};
//iReadAll, passFilterIllumina, passFilterIllumina
chunkInSizeBytesTotal[imate1] += sprintf(chunkIn[imate1] + chunkInSizeBytesTotal[imate1], " %llu %c %i", P.iReadAll, passFilterIllumina, P.readFilesIndex);
string dummy;
for (int ii=3; ii<=9; ii++)
P.inOut->readIn[0] >> dummy; //skip fields until sequence
string seq1,qual1;
P.inOut->readIn[0] >> seq1 >> qual1;
if (flag1 & 0x10) {//sequence reverse-coomplemented
revComplementNucleotides(seq1);
reverse(qual1.begin(),qual1.end());
};
string attrs;
getline(P.inOut->readIn[0], attrs); //rest of the SAM line: str1 is now all SAM attributes - it's added to the read ID line (1st "fastq" line)
chunkInSizeBytesTotal[imate1] += sprintf(chunkIn[imate1] + chunkInSizeBytesTotal[imate1], "%s\n%s\n+\n%s\n", attrs.c_str(), seq1.c_str(), qual1.c_str());
};
};
///////////////////////////////////////////////////////////////////////////////////// FASTQ
} else if (nextChar=='@') {//fastq, not multi-line
P.iReadAll++; //increment read number
if (P.outFilterBySJoutStage!=2) {//not the 2nd stage of the 2-stage mapping, read ID from the 1st read
string readID;
P.inOut->readIn[0] >> readID;
removeStringEndControl(readID);
if (P.outSAMreadIDnumber) {
readID="@"+to_string(P.iReadAll);
};
//read the second field of the read name line
char passFilterIllumina='N';
if (P.inOut->readIn[0].peek()!='\n') {//2nd field exists
string field2;
P.inOut->readIn[0] >> field2;
if (field2.length()>=3 && field2[1]==':' && field2[2]=='Y' && field2[3]==':' )
passFilterIllumina='Y';
};
//add extra information to readID line
readID += ' '+ to_string(P.iReadAll)+' '+passFilterIllumina+' '+to_string(P.readFilesIndex);
//ignore the rest of the read name for both mates
for (uint imate=0; imate<P.readNends; imate++)
P.inOut->readIn[imate].ignore(DEF_readNameSeqLengthMax,'\n');
//copy the same readID to both mates
for (uint imate=0; imate<P.readNends; imate++) {
chunkInSizeBytesTotal[imate] += 1 + readID.copy(chunkIn[imate] + chunkInSizeBytesTotal[imate], readID.size(),0);
chunkIn[imate][chunkInSizeBytesTotal[imate]-1]='\n';
};
};
//copy 3 (4 for stage 2) lines: sequence, dummy, quality
for (uint imate=0; imate<P.readNends; imate++) {
// read 1st line for 2nd stage only
if (P.outFilterBySJoutStage == 2)
chunkInSizeBytesTotal[imate] += fastqReadOneLine(P.inOut->readIn[imate], chunkIn[imate] + chunkInSizeBytesTotal[imate]);
//sequence
chunkInSizeBytesTotal[imate] += fastqReadOneLine(P.inOut->readIn[imate], chunkIn[imate] + chunkInSizeBytesTotal[imate]);
//skip 3rd line, record '+'
P.inOut->readIn[imate].ignore(DEF_readNameSeqLengthMax, '\n');
chunkIn[imate][chunkInSizeBytesTotal[imate]] = '+';
chunkIn[imate][chunkInSizeBytesTotal[imate]+1] = '\n';
chunkInSizeBytesTotal[imate] += 2;
//quality
uint64 lenIn = fastqReadOneLine(P.inOut->readIn[imate], chunkIn[imate] + chunkInSizeBytesTotal[imate]);
chunkInSizeBytesTotal[imate] += lenIn;
};
} else if (nextChar=='>') {//fasta, can be multiline, which is converted to single line
P.iReadAll++; //increment read number
for (uint imate=0; imate<P.readNends; imate++) {
if (P.outFilterBySJoutStage!=2) {//not the 2nd stage of the 2-stage mapping
if (P.outSAMreadID=="Number") {
chunkInSizeBytesTotal[imate] += sprintf(chunkIn[imate] + chunkInSizeBytesTotal[imate], ">%llu", P.iReadAll);
} else {
P.inOut->readIn[imate] >> (chunkIn[imate] + chunkInSizeBytesTotal[imate]);
chunkInSizeBytesTotal[imate] += strlen(chunkIn[imate] + chunkInSizeBytesTotal[imate]);
};
P.inOut->readIn[imate].ignore(DEF_readNameSeqLengthMax,'\n');
chunkInSizeBytesTotal[imate] += sprintf(chunkIn[imate] + chunkInSizeBytesTotal[imate], " %llu %c %i \n", P.iReadAll, 'N', P.readFilesIndex);
};
//read multi-line fasta
nextChar=P.inOut->readIn[imate].peek();
while (nextChar!='@' && nextChar!='>' && nextChar!=' ' && nextChar!='\n' && P.inOut->readIn[imate].good()) {
P.inOut->readIn[imate].getline(chunkIn[imate] + chunkInSizeBytesTotal[imate], DEF_readSeqLengthMax + 1 );
if (P.inOut->readIn[imate].gcount()<2)
break; //no more input
chunkInSizeBytesTotal[imate] += P.inOut->readIn[imate].gcount()-1; //-1 because \n was counted, bu wee need to remove it
if ( int(chunkIn[imate][chunkInSizeBytesTotal[imate]-1]) < 33 ) {//remove control char at the end if present
chunkInSizeBytesTotal[imate]--;
};
nextChar=P.inOut->readIn[imate].peek();
};
chunkIn[imate][chunkInSizeBytesTotal[imate]]='\n';
chunkInSizeBytesTotal[imate] ++;
};
} else if (nextChar==' ' || nextChar=='\n' || !P.inOut->readIn[0].good()) {//end of stream
P.inOut->logMain << "Thread #" <<iThread <<" end of input stream, nextChar="<<int(nextChar) <<endl;
break;
} else {
string word1;
P.inOut->readIn[0] >> word1;
if (word1=="FILE") {//new file marker
newFile=true;
} else {//error
ostringstream errOut;
errOut << ERROR_OUT <<" EXITING because of FATAL ERROR in input reads: unknown file format: the read ID should start with @ or > \n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
};
if (newFile) {
P.inOut->readIn[0] >> P.readFilesIndex;
pthread_mutex_lock(&g_threadChunks.mutexLogMain);
P.inOut->logMain << "Starting to map file # " << P.readFilesIndex<<"\n";
for (uint imate=0; imate<P.readFilesNames.size(); imate++) {
P.inOut->logMain << "mate " <<imate+1 <<": "<<P.readFilesNames.at(imate).at(P.readFilesIndex) <<"\n";
P.inOut->readIn[imate].ignore(numeric_limits<streamsize>::max(),'\n');
};
P.inOut->logMain<<flush;
pthread_mutex_unlock(&g_threadChunks.mutexLogMain);
newFile=false;
};
};
//TODO: check here that both mates are zero or non-zero
if (chunkInSizeBytesTotal[0]==0) {
noReadsLeft=true; //true if there no more reads left in the file
iChunkIn=g_threadChunks.chunkInN;//to keep things consistent
g_threadChunks.chunkInN++;
} else {
noReadsLeft=false;
iChunkIn=g_threadChunks.chunkInN;
g_threadChunks.chunkInN++;
};
for (uint imate=0; imate<P.readNends; imate++)
chunkIn[imate][chunkInSizeBytesTotal[imate]]='\n';//extra empty line at the end of the chunks
if (P.runThreadN>1) pthread_mutex_unlock(&g_threadChunks.mutexInRead);
} else {//read from one file per thread
noReadsLeft=true;
for (uint imate=0; imate<P.readNends; imate++) {
RA->chunkOutFilterBySJoutFiles[imate].flush();
RA->chunkOutFilterBySJoutFiles[imate].seekg(0,ios::beg);
RA->readInStream[imate]=& RA->chunkOutFilterBySJoutFiles[imate];
};
};
mapChunk();
if (iThread==0 && P.runThreadN>1 && P.outSAMorder=="PairedKeepInputOrder") {//concatenate Aligned.* files
chunkFilesCat(P.inOut->outSAM, P.outFileTmp + "/Aligned.out.sam.chunk", g_threadChunks.chunkOutN);
};
};//cycle over input chunks
if (P.outFilterBySJoutStage!=1 && RA->iRead>0) {//not the first stage of the 2-stage mapping
if (P.outBAMunsorted) chunkOutBAMunsorted->unsortedFlush();
if (P.outBAMcoord) chunkOutBAMcoord->coordFlush();
if (chunkOutBAMquant!=NULL) chunkOutBAMquant->unsortedFlush();
//the thread is finished mapping reads, concatenate the temp files into output files
if (P.pCh.segmentMin>0) {
chunkFstreamCat (RA->chunkOutChimSAM, P.inOut->outChimSAM, P.runThreadN>1, g_threadChunks.mutexOutChimSAM);
chunkFstreamCat (*RA->chunkOutChimJunction, P.inOut->outChimJunction, P.runThreadN>1, g_threadChunks.mutexOutChimJunction);
};
if (P.outReadsUnmapped=="Fastx" ) {
if (P.runThreadN>1)
pthread_mutex_lock(&g_threadChunks.mutexOutUnmappedFastx);
for (uint ii=0;ii<P.readNends;ii++) {
chunkFstreamCat (RA->chunkOutUnmappedReadsStream[ii],P.inOut->outUnmappedReadsStream[ii], false, g_threadChunks.mutexOutUnmappedFastx);
};
if (P.runThreadN>1)
pthread_mutex_unlock(&g_threadChunks.mutexOutUnmappedFastx);
};
};
if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexLogMain);
P.inOut->logMain << "Completed: thread #" <<iThread <<endl;
if (P.runThreadN>1) pthread_mutex_unlock(&g_threadChunks.mutexLogMain);
};
inline uint64 fastqReadOneLine(ifstream &streamIn, char *arrIn)
{
uint64 lenIn;
streamIn.getline(arrIn, DEF_readNameSeqLengthMax+1 );
lenIn = streamIn.gcount(); //=seqLength+1: includes \0 but not \n. We will replace \0 with \n
if ( int(arrIn[lenIn-2]) < 33 ) {//remove control char at the end if present
--lenIn;
};
arrIn[lenIn-1]='\n'; //replace \0 with \n
return lenIn; //lenIn contains \n at the end
};
inline void removeStringEndControl(string &str)
{//removes control character (including space) from the end of the string
if (int(str.back())<33)
str.pop_back();
};
|
