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
|
#include "ReadAlignChunk.h"
#include "GlobalVariables.h"
#include "ThreadControl.h"
#include "ErrorWarning.h"
#include SAMTOOLS_BGZF_H
void ReadAlignChunk::mapChunk() {//map one chunk. Input reads stream has to be setup in RA->readInStream[ii]
for (uint32 im=0; im<1; im++) {//hardcoded mate 1 5p onyl for now
RA->clipMates[im][0].clipChunk(chunkIn[im], chunkInSizeBytesTotal[im]);
};
RA->statsRA.resetN();
for (uint ii=0;ii<P.readNends;ii++) {//clear eof and rewind the input streams
RA->readInStream[ii]->clear();
RA->readInStream[ii]->seekg(0,ios::beg);
};
if ( P.outSAMorder == "PairedKeepInputOrder" && P.runThreadN>1 ) {//open chunk file
ostringstream name1("");
name1 << P.outFileTmp + "/Aligned.tmp.sam.chunk"<<iChunkIn;
chunkOutBAMfileName = name1.str();
chunkOutBAMfile.open(chunkOutBAMfileName.c_str());
};
int readStatus=0;
while (readStatus==0) {//main cycle over all reads
readStatus=RA->oneRead(); //map one read
if (readStatus==0) {//there was a read processed
RA->iRead++;
// chunkOutBAMtotal=(uint) RA->outSAMstream->tellp();
chunkOutBAMtotal+=RA->outBAMbytes;
// uint ddd=(uint) RA->outSAMstream->tellp();
};
//write SAM aligns to chunk buffer
if (P.outSAMbool) {
if ( chunkOutBAMtotal > P.chunkOutBAMsizeBytes ) {//this should not happen!
ostringstream errOut;
errOut <<"EXITING because of fatal error: buffer size for SAM/BAM output is too small\n";
errOut <<"Solution: increase input parameter --limitOutSAMoneReadBytes\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
} else if ( chunkOutBAMtotal + P.limitOutSAMoneReadBytes > P.chunkOutBAMsizeBytes || (readStatus==-1 && noReadsLeft) ) {//write buffer to disk because it's almost full, or all reads are mapped
if ( P.outSAMorder == "PairedKeepInputOrder" && P.runThreadN>1 ) {//output chunks into separate files
chunkOutBAMfile.write(chunkOutBAM,chunkOutBAMtotal);
chunkOutBAMfile.clear(); //in case 0 bytes were written which could set fail bit
} else {//standard way, directly into Aligned.out.sam file
//SAM output
if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexOutSAM);
P.inOut->outSAM->write(chunkOutBAM,chunkOutBAMtotal);
P.inOut->outSAM->clear();//in case 0 bytes were written which could set fail bit
if (P.runThreadN>1) pthread_mutex_unlock(&g_threadChunks.mutexOutSAM);
};
RA->outSAMstream->seekp(0,ios::beg); //rewind the chunk storage
chunkOutBAMtotal=0;
};
};
//collapse SJ buffer if needed
if ( chunkOutSJ->N > P.limitOutSJcollapsed ) {//this means the number of collapsed junctions is larger than the chunks size
ostringstream errOut;
errOut <<"EXITING because of fatal error: buffer size for SJ output is too small\n";
errOut <<"Solution: increase input parameter --limitOutSJoneRead\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
} else if ( chunkOutSJ->N + P.limitOutSJoneRead > P.limitOutSJcollapsed || (readStatus==-1 && noReadsLeft) ) {//write buffer to disk because it's almost full, or all reads are mapped
chunkOutSJ->collapseSJ();
if ( chunkOutSJ->N + 2*P.limitOutSJoneRead > P.limitOutSJcollapsed ) {
ostringstream errOut;
errOut <<"EXITING because of fatal error: buffer size for SJ output is too small\n";
errOut <<"Solution: increase input parameter --limitOutSJcollapsed\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
};
//collapse SJ1 buffer if needed
if ( chunkOutSJ1->N > P.limitOutSJcollapsed ) {//this means the number of collapsed junctions is larger than the chunks size
ostringstream errOut;
errOut <<"EXITING because of fatal error: buffer size for SJ output is too small\n";
errOut <<"Solution: increase input parameter --limitOutSJoneRead\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
} else if ( chunkOutSJ1->N + P.limitOutSJoneRead > P.limitOutSJcollapsed || (readStatus==-1 && noReadsLeft) ) {//write buffer to disk because it's almost full, or all reads are mapped
chunkOutSJ1->collapseSJ();
if ( chunkOutSJ1->N + 2*P.limitOutSJoneRead > P.limitOutSJcollapsed ) {
ostringstream errOut;
errOut <<"EXITING because of fatal error: buffer size for SJ output is too small\n";
errOut <<"Solution: increase input parameter --limitOutSJcollapsed\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
};
}; //reads cycle
if ( P.outSAMbool && P.outSAMorder == "PairedKeepInputOrder" && P.runThreadN>1 ) {//write the remaining part of the buffer, close and rename chunk files
chunkOutBAMfile.write(chunkOutBAM,chunkOutBAMtotal);
chunkOutBAMfile.clear(); //in case 0 bytes were written which could set fail bit
chunkOutBAMfile.close();
RA->outSAMstream->seekp(0,ios::beg); //rewind the chunk storage
chunkOutBAMtotal=0;
ostringstream name2("");
name2 << P.outFileTmp + "/Aligned.out.sam.chunk"<<iChunkIn;
rename(chunkOutBAMfileName.c_str(),name2.str().c_str());//marks files as completedly written
};
//add stats, write progress if needed
if (P.runThreadN>1) pthread_mutex_lock(&g_threadChunks.mutexStats);
g_statsAll.addStats(RA->statsRA);
g_statsAll.progressReport(P.inOut->logProgress);
if (P.runThreadN>1) pthread_mutex_unlock(&g_threadChunks.mutexStats);
};
|
