1 #include "Transcriptome.h"
2 #include "ReadAlign.h"
3 #include "Transcript.h"
4 #include "serviceFuns.cpp"
5 #include <random>
6
quantTranscriptome(Transcriptome * Tr,uint nAlignG,Transcript ** alignG,Transcript * alignT)7 uint ReadAlign::quantTranscriptome (Transcriptome *Tr, uint nAlignG, Transcript **alignG, Transcript *alignT) {
8 uint nAlignT=0;
9 for (uint iag=0; iag<nAlignG; iag++) {//transform all alignments
10
11 Transcript *align1=alignG[iag];
12
13 if (!P.quant.trSAM.indel && (align1->nDel>0 || align1->nIns>0) ) {
14 //prevent indels if requested
15 continue;
16 };
17 if (!P.quant.trSAM.singleEnd && (P.readNmates==2 && align1->exons[0][EX_iFrag]==align1->exons[align1->nExons-1][EX_iFrag]) ) {//not readNends: this is alignment
18 //prevent single end alignments
19 continue;
20 };
21
22 if (!P.quant.trSAM.softClip) {
23 //soft clipping not allowed, extend them if possible
24 uint nMM1=0;
25 char* R=Read1[align1->roStr==0 ? 0:2];
26 Transcript align2=*align1; //copy this transcript to avoid changing the original one
27
28 for (uint32 iab=0; iab<align2.nExons; iab++) {
29 uint left1=0,right1=0;//how many bases to move left or right
30 if (iab==0) {
31 left1=align2.exons[iab][EX_R];
32 } else if (align2.canonSJ[iab-1]==-3) {
33 left1=align2.exons[iab][EX_R]-readLength[align2.exons[iab-1][EX_iFrag]]-1;
34 };
35 if (iab==align2.nExons-1) {//last block of left mates
36 right1=Lread-align2.exons[iab][EX_R]-align2.exons[iab][EX_L];
37
38 } else if (align2.canonSJ[iab]==-3) {//last block of the right mate (i.e. whole read)
39 right1=readLength[align2.exons[iab][EX_iFrag]]-align2.exons[iab][EX_R]-align2.exons[iab][EX_L];
40 };
41
42 for (uint b=1; b<=left1 ; b++) {//extend to the left
43 char r1=R[align2.exons[iab][EX_R]-b];
44 char g1=mapGen.G[align2.exons[iab][EX_G]-b];
45 if ( r1!=g1 && r1<4 && g1<4) ++nMM1;
46 };
47 for (uint b=0; b<right1 ; b++) {//extend to the left
48 char r1=R[align2.exons[iab][EX_R]+align2.exons[iab][EX_L]+b];
49 char g1=mapGen.G[align2.exons[iab][EX_G]+align2.exons[iab][EX_L]+b];
50 if ( r1!=g1 && r1<4 && g1<4) ++nMM1;
51 };
52 align2.exons[iab][EX_R] -= left1;
53 align2.exons[iab][EX_G] -= left1;
54 align2.exons[iab][EX_L] += left1+right1;
55 };
56
57 if ( (align2.nMM + nMM1) > min(outFilterMismatchNmaxTotal, (uint) (P.outFilterMismatchNoverLmax*(Lread-1)) ) ) {
58 //extension of soft clips yielded too many mismatches, no output
59 continue;
60 };
61
62 align1 = &align2;
63 };
64
65 nAlignT += Tr->quantAlign(*align1,alignT+nAlignT);
66 };
67
68 if (P.quant.trSAM.bamYes) {//output Aligned.toTranscriptome.bam
69 alignT[int(rngUniformReal0to1(rngMultOrder)*nAlignT)].primaryFlag=true;
70
71 for (uint iatr=0;iatr<nAlignT;iatr++) {//write all transcripts
72 alignBAM(alignT[iatr], nAlignT, iatr, 0, (uint) -1, (uint) -1, 0, -1, NULL, P.outSAMattrOrderQuant, outBAMoneAlign, outBAMoneAlignNbytes);
73 for (uint imate=0; imate<P.readNmates; imate++) {//output each mate //not readNends: this is alignment
74 outBAMquant->unsortedOneAlign(outBAMoneAlign[imate], outBAMoneAlignNbytes[imate], (imate>0 || iatr>0) ? 0 : (outBAMoneAlignNbytes[0]+outBAMoneAlignNbytes[1])*2*nAlignT);
75 };
76 };
77 };
78
79 //not used anymore per Colin Dewey's request
80 // if (nAlignT==0 && P.outSAMunmapped=="Within") {//read could be mapped to genome, but not transcriptome - output as unmapped
81 // uint unmapType=5;
82 // bool mateMapped[2]={false,false};
83 // alignBAM(*alignG[0], 0, 0, mapGen.chrStart[alignG[0]->Chr], (uint) -1, (uint) -1, 0, unmapType, mateMapped, P.outSAMattrOrder);
84 // for (uint imate=0; imate<P.readNmates; imate++) {//output each mate //not readNends: this is alignment
85 // outBAMquant->unsortedOneAlign(outBAMoneAlign[imate], outBAMoneAlignNbytes[imate], imate>0 ? 0 : outBAMoneAlignNbytes[0]+outBAMoneAlignNbytes[1]);
86 // };
87 //
88 // };
89
90 return nAlignT;
91 };
92