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
|
#include "ReadAlign.h"
#include "readLoad.h"
#include "readBarcodeLoad.h"
#include "SequenceFuns.h"
#include "ErrorWarning.h"
#include "GlobalVariables.h"
int ReadAlign::oneRead() {//process one read: load, map, write
//load read name, sequence, quality from the streams into internal arrays
int readStatus[P.readNends];
for (uint32 im=0; im<P.readNends; im++) {
readStatus[im] = readLoad(*(readInStream[im]), P, readLength[im], readLengthOriginal[im], readNameMates[im], Read0[im], Read1[im], Qual0[im], clipMates[im], iReadAll, readFilesIndex, readFilter, readNameExtra[im]);
if (readStatus[im] != readStatus[0]) {//check if the end of file was reached or not for all files
ostringstream errOut;
errOut << "EXITING because of FATAL ERROR: read files are not consistent, reached the end of the one before the other one\n";
errOut << "SOLUTION: Check you your input files: they may be corrupted\n";
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
};
if (readStatus[0]==-1) {//finished with the stream
return -1;
};
if (P.outFilterBySJoutStage != 2) {
for (uint32 im=0; im<P.readNmates; im++) {//not readNends: the barcode quality will be calculated separately
for (uint64 ix=clipMates[im][0].clippedN; ix<readLengthOriginal[im]-clipMates[im][1].clippedN; ix++) {
qualHist[im][(uint8)Qual0[im][ix]]++;
};
};
};
if (P.readNmates==2) {//combine two mates together
Lread=readLength[0]+readLength[1]+1;
readLengthPairOriginal=readLengthOriginal[0]+readLengthOriginal[1]+1;
if (Lread>DEF_readSeqLengthMax) {
ostringstream errOut;
errOut << "EXITING because of FATAL ERROR in reads input: Lread of the pair = " << Lread << " while DEF_readSeqLengthMax=" << DEF_readSeqLengthMax <<endl;
errOut << "Read Name="<<readNameMates[0]<<endl;
errOut << "SOLUTION: increase DEF_readSeqLengthMax in IncludeDefine.h and re-compile STAR"<<endl<<flush;
exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_INPUT_FILES, P);
};
//marker for spacer base
Read1[0][readLength[0]]=MARK_FRAG_SPACER_BASE;
//copy 2nd mate into Read1[0] & reverse-complement
complementSeqNumbers(Read1[1],Read1[0]+readLength[0]+1,readLength[1]);//complement. Here Read1[1] is still the 2nd mate's numeric-sequence. Later Read1[1] will be reverse complement of the combined read.
for (uint ii=0;ii<readLength[1]/2;ii++) {
swap(Read1[0][Lread-ii-1],Read1[0][ii+readLength[0]+1]); //reverse
};
} else {//1 mate
if (readStatus[0]==-1) {//finished with the stream
return -1;
};
Lread=readLength[0];
readLengthPairOriginal=readLengthOriginal[0];
readLength[1]=0;
};
readFileType=readStatus[0];
complementSeqNumbers(Read1[0],Read1[1],Lread); //returns complement of Reads[ii]
for (uint ii=0;ii<Lread;ii++) {//reverse
Read1[2][Lread-ii-1]=Read1[1][ii];
};
statsRA.readN++;
statsRA.readBases += readLength[0]+readLength[1];
//max number of mismatches allowed for this read
outFilterMismatchNmaxTotal=min(P.outFilterMismatchNmax, (uint) (P.outFilterMismatchNoverReadLmax*(readLength[0]+readLength[1])));
//map the read
if (P.pGe.gType==101) {//SpliceGraph
mapOneReadSpliceGraph();
} else {//all other cases - standard alignment algorithm
mapOneRead();
};
peOverlapMergeMap();
multMapSelect();
mappedFilter();
transformGenome();//for now genome transformation happens after multimapper selection, and mapping filter
if (!peOv.yes) {//if the alignment was not mates merged - otherwise the chimeric detection was already done
chimericDetection();
};
if (P.pCh.out.bam && chimRecord) {//chimeric alignment was recorded in main BAM files, and it contains the representative portion, so non-chimeric aligmnent is not output
return 0;
};
waspMap();
#ifdef OFF_BEFORE_OUTPUT
#warning OFF_BEFORE_OUTPUT
return 0;
#endif
//write out alignments
outputAlignments();
return 0;
};
|
