1 #include "samHeaders.h"
2 #include <iostream>
3 #include "BAMfunctions.h"
4
samHeaders(Parameters & P,Genome & genomeOut,Transcriptome & transcriptomeMain)5 void samHeaders(Parameters &P, Genome &genomeOut, Transcriptome &transcriptomeMain)
6 {
7 /////////////////////////////////////////////////////////////////////////////////// transcriptome BAM header
8 if ( P.quant.trSAM.bamYes ) {//header for transcriptome BAM
9 ostringstream samHeaderStream;
10 vector <uint> trlength;
11 for (uint32 ii=0;ii<transcriptomeMain.trID.size();ii++) {
12 uint32 iex1=transcriptomeMain.trExI[ii]+transcriptomeMain.trExN[ii]-1; //last exon of the transcript
13 trlength.push_back(transcriptomeMain.exLenCum[iex1]+transcriptomeMain.exSE[2*iex1+1]-transcriptomeMain.exSE[2*iex1]+1);
14 samHeaderStream << "@SQ\tSN:"<< transcriptomeMain.trID.at(ii) <<"\tLN:"<<trlength.back()<<"\n";
15 };
16 for (uint32 ii=0;ii<P.outSAMattrRGlineSplit.size();ii++) {//@RG lines
17 samHeaderStream << "@RG\t" << P.outSAMattrRGlineSplit.at(ii) <<"\n";
18 };
19 outBAMwriteHeader(P.inOut->outQuantBAMfile,samHeaderStream.str(),transcriptomeMain.trID,trlength);
20 };
21
22 //////////////////////////////////////////////////////////////////////////////// main headers
23 if (P.outSAMmode == "None" || P.outSAMtype[0] == "None") {//no SAM output
24 return;
25 };
26
27 ostringstream samHeaderStream;
28
29 for (uint ii=0;ii<genomeOut.nChrReal;ii++) {
30 samHeaderStream << "@SQ\tSN:"<< genomeOut.chrName.at(ii) <<"\tLN:"<<genomeOut.chrLength[ii]<<"\n";
31 };
32
33 genomeOut.chrNameAll=genomeOut.chrName;
34 genomeOut.chrLengthAll=genomeOut.chrLength;
35 {//add exra references
36 ifstream extrastream (P.pGe.gDir + "/extraReferences.txt");
37 while (extrastream.good()) {
38 string line1;
39 getline(extrastream,line1);
40 istringstream stream1 (line1);
41 string field1;
42 stream1 >> field1;//should check for @SQ
43
44 if (field1!="") {//skip blank lines
45 samHeaderStream << line1 <<"\n";
46
47 stream1 >> field1;
48 genomeOut.chrNameAll.push_back(field1.substr(3));
49 stream1 >> field1;
50 genomeOut.chrLengthAll.push_back((uint) stoll(field1.substr(3)));
51 };
52 };
53 extrastream.close();
54 };
55
56 if (P.outSAMheaderPG.at(0)!="-") {
57 samHeaderStream << P.outSAMheaderPG.at(0);
58 for (uint ii=1;ii<P.outSAMheaderPG.size(); ii++) {
59 samHeaderStream << "\t" << P.outSAMheaderPG.at(ii);
60 };
61 samHeaderStream << "\n";
62 };
63
64 samHeaderStream << "@PG\tID:STAR\tPN:STAR\tVN:" << STAR_VERSION <<"\tCL:" << P.commandLineFull <<"\n";
65
66 if (P.outSAMheaderCommentFile!="-") {
67 ifstream comstream (P.outSAMheaderCommentFile);
68 while (comstream.good()) {
69 string line1;
70 getline(comstream,line1);
71 if (line1.find_first_not_of(" \t\n\v\f\r")!=std::string::npos) {//skip blank lines
72 samHeaderStream << line1 <<"\n";
73 };
74 };
75 comstream.close();
76 };
77
78
79 for (uint32 ii=0;ii<P.outSAMattrRGlineSplit.size();ii++) {//@RG lines
80 samHeaderStream << "@RG\t" << P.outSAMattrRGlineSplit.at(ii) <<"\n";
81 };
82
83
84 samHeaderStream << "@CO\t" <<"user command line: " << P.commandLine <<"\n";
85
86 samHeaderStream << P.samHeaderExtra;
87
88 if (P.outSAMheaderHD.at(0)!="-") {
89 P.samHeaderHD = P.outSAMheaderHD.at(0);
90 for (uint ii=1;ii<P.outSAMheaderHD.size(); ii++) {
91 P.samHeaderHD +="\t" + P.outSAMheaderHD.at(ii);
92 };
93 } else {
94 P.samHeaderHD = "@HD\tVN:1.4";
95 };
96
97
98 P.samHeader=P.samHeaderHD+"\n"+samHeaderStream.str();
99 //for the sorted BAM, need to add SO:cooridnate to the header line
100 P.samHeaderSortedCoord=P.samHeaderHD + (P.outSAMheaderHD.size()==0 ? "" : "\tSO:coordinate") + "\n" + samHeaderStream.str();
101
102 if (P.outSAMbool) {//
103 *P.inOut->outSAM << P.samHeader;
104 };
105 if (P.outBAMunsorted){
106 outBAMwriteHeader(P.inOut->outBAMfileUnsorted,P.samHeader,genomeOut.chrNameAll,genomeOut.chrLengthAll);
107 };
108 };