[med-svn] [Git][med-team/rna-star][upstream] New upstream version 2.7.6a+dfsg
Andreas Tille
gitlab at salsa.debian.org
Fri Oct 2 13:54:40 BST 2020
Andreas Tille pushed to branch upstream at Debian Med / rna-star
Commits:
ababb9e7 by Andreas Tille at 2020-10-02T14:49:26+02:00
New upstream version 2.7.6a+dfsg
- - - - -
30 changed files:
- CHANGES.md
- README.md
- RELEASEnotes.md
- doc/STARmanual.pdf
- extras/doc-latex/STARmanual.tex
- extras/docker/Dockerfile
- source/BAMfunctions.cpp
- source/ChimericAlign.cpp
- source/ChimericAlign.h
- + source/ChimericAlign_chimericBAMoutput.cpp
- source/ChimericDetection.h
- source/ChimericDetection_chimericDetectionMult.cpp
- source/Makefile
- source/Parameters.cpp
- source/Parameters.h
- source/ParametersChimeric_initialize.cpp
- + source/Parameters_samAttributes.cpp
- source/ReadAlign.h
- source/ReadAlignChunk.cpp
- source/ReadAlign_alignBAM.cpp
- source/ReadAlign_chimericDetection.cpp
- source/ReadAlign_chimericDetectionOldOutput.cpp
- source/ReadAlign_chimericDetectionPEmerged.cpp
- source/ReadAlign_outputAlignments.cpp
- source/ReadAlign_outputTranscriptSAM.cpp
- source/ReadAlign_peOverlapMergeMap.cpp
- source/SoloFeature_cellFiltering.cpp
- source/Transcript.h
- source/Transcript_alignScore.cpp
- source/VERSION
Changes:
=====================================
CHANGES.md
=====================================
@@ -1,3 +1,17 @@
+STAR 2.7.6a --- 2020/09/19
+==========================
+**Major new feature:**
+Output multimapping chimeric alignments in BAM format using
+```--chimMultimapNmax N>1 --chimOutType WithinBAM --outSAMtype BAM Unsorted [and/or] SortedByCoordinate```
+Many thanks to Sebastian @suhrig who implemented this feature!
+More detailed description from Sebastian in PR #802.
+
+**Minor features and bug fixes:**
+* Issue #1008: fixed the problem with Unmapped.out.mate? output for --soloType CB_samTagOut output.
+* PR # 1012: fixed the bug with --soloCellFiltering TopCells option.
+* Issue #786: fixed the bug causing the *Different SJ motifs problem* for overlapping mates.
+* Issue #945: GX/GN can be output for all --soloType, as well as for non-solo runs.
+
STAR 2.7.5c --- 2020/08/16
==========================
Bug-fix release.
=====================================
README.md
=====================================
@@ -35,9 +35,9 @@ Download the latest [release from](https://github.com/alexdobin/STAR/releases) a
```bash
# Get latest STAR source from releases
-wget https://github.com/alexdobin/STAR/archive/2.7.5c.tar.gz
-tar -xzf 2.7.5c.tar.gz
-cd STAR-2.7.5c
+wget https://github.com/alexdobin/STAR/archive/2.7.6a.tar.gz
+tar -xzf 2.7.6a.tar.gz
+cd STAR-2.7.6a
# Alternatively, get STAR source using git
git clone https://github.com/alexdobin/STAR.git
=====================================
RELEASEnotes.md
=====================================
@@ -1,3 +1,11 @@
+STAR 2.7.6a --- 2020/09/19
+==========================
+**Major new feature:**
+Output multimapping chimeric alignments in BAM format using
+```--chimMultimapNmax N>1 --chimOutType WithinBAM --outSAMtype BAM Unsorted [and/or] SortedByCoordinate```
+Many thanks to Sebastian @suhrig who implemented this feature!
+More detailed description from Sebastian in PR #802.
+
STAR 2.7.5a 2020/06/16
======================
**Major new features:
=====================================
doc/STARmanual.pdf
=====================================
Binary files a/doc/STARmanual.pdf and b/doc/STARmanual.pdf differ
=====================================
extras/doc-latex/STARmanual.tex
=====================================
@@ -34,7 +34,7 @@
\newcommand{\sechyperref}[1]{\hyperref[#1]{Section \ref{#1}. \nameref{#1}}}
-\title{STAR manual 2.7.5c}
+\title{STAR manual 2.7.6a}
\author{Alexander Dobin\\
dobin at cshl.edu}
\maketitle
=====================================
extras/docker/Dockerfile
=====================================
@@ -2,7 +2,7 @@ FROM debian:stable-slim
MAINTAINER dobin at cshl.edu
-ARG STAR_VERSION=2.7.5c
+ARG STAR_VERSION=2.7.6a
ENV PACKAGES gcc g++ make wget zlib1g-dev unzip
=====================================
source/BAMfunctions.cpp
=====================================
@@ -59,9 +59,10 @@ int bam_read1_fromArray(char *bamChar, bam1_t *b) //modified from samtools bam_r
if (b->m_data < b->l_data) {
b->m_data = b->l_data;
kroundup32(b->m_data);
- b->data = (uint8_t*)realloc(b->data, b->m_data);
- if (!b->data)
- return -4;
+// no need to realloc b->data, because it is overwritten later with bamChar
+// b->data = (uint8_t*)realloc(b->data, b->m_data);
+// if (!b->data)
+// return -4;
}
// // if (bgzf_read(fp, b->data, b->l_data) != b->l_data) return -4;
// // //b->l_aux = b->l_data - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;
=====================================
source/ChimericAlign.cpp
=====================================
@@ -1,7 +1,7 @@
#include "ChimericAlign.h"
ChimericAlign::ChimericAlign(ChimericSegment &seg1in, ChimericSegment &seg2in, int chimScoreIn, Genome &genomeIn, ReadAlign *RAin)
- : seg1(seg1in), seg2(seg2in),chimScore(chimScoreIn), P(seg1in.P), pCh(P.pCh), mapGen(genomeIn), RA(RAin) {
+ : seg1(seg1in), seg2(seg2in),chimScore(chimScoreIn), RA(RAin), P(seg1in.P), pCh(P.pCh), mapGen(genomeIn) {
stitchingDone=false;
al1=&seg1.align;
=====================================
source/ChimericAlign.h
=====================================
@@ -24,16 +24,17 @@ class ChimericAlign
ChimericAlign(ChimericSegment &seg1in, ChimericSegment &seg2in, int chimScoreIn, Genome &genomeIn, ReadAlign *RAin); //allocate
void chimericJunctionOutput(fstream &outStream, uint chimN, int maxNonChimAlignScore, bool PEmerged_flag, int chimScoreBest, int maxPossibleAlignScore);
+ static void chimericBAMoutput(Transcript *al1, Transcript *al2, ReadAlign *RA, const uint iTr, const uint chimN, const bool isBestChimAlign, const Parameters& P);
void chimericStitching(char *genSeq, char **Read1);
bool chimericCheck();
bool stitchingDone;
+ ReadAlign *RA;
private:
Parameters &P;
ParametersChimeric &pCh;
Genome &mapGen;
- ReadAlign *RA;
};
=====================================
source/ChimericAlign_chimericBAMoutput.cpp
=====================================
@@ -0,0 +1,105 @@
+#include "ChimericAlign.h"
+#include "ReadAlign.h"
+#include "BAMfunctions.h"
+
+#include <vector>
+
+void ChimericAlign::chimericBAMoutput(Transcript *al1, Transcript *al2, ReadAlign *RA, const uint iTr, const uint chimN, const bool isBestChimAlign, const Parameters& P)
+{
+ vector<Transcript*> trChim(2);
+ trChim[0] = al1;
+ trChim[1] = al2;
+
+ int chimRepresent=-999, chimType=0;
+ if (trChim[0]->exons[0][EX_iFrag]!=trChim[0]->exons[trChim[0]->nExons-1][EX_iFrag]) {//tr0 has both mates
+ chimRepresent = 0;
+ chimType = 1;
+ } else if (trChim[1]->exons[0][EX_iFrag]!=trChim[1]->exons[trChim[1]->nExons-1][EX_iFrag]) {//tr1 has both mates
+ chimRepresent = 1;
+ chimType = 1;
+ } else if (trChim[0]->exons[0][EX_iFrag]!=trChim[1]->exons[0][EX_iFrag]) {//tr0 and tr1 are single different mates
+ chimRepresent = -1;
+ chimType = 2;
+ } else {//two chimeric segments are on the same mate - this can only happen for single-end reads
+ chimRepresent = (trChim[0]->maxScore > trChim[1]->maxScore) ? 0 : 1;
+ chimType = 3;
+ };
+
+ int alignType, bamN=0, bamIsuppl=-1, bamIrepr=-1;
+ uint bamBytesTotal=0;//estimate of the total size of all bam records, for output buffering
+ uint mateChr,mateStart;
+ uint8_t mateStrand;
+ for (uint itr=0;itr<trChim.size();itr++) {//generate bam for all chimeric pieces
+ trChim[itr]->primaryFlag=isBestChimAlign;
+ if (chimType==2) {//PE, encompassing
+ mateChr=trChim[1-itr]->Chr;
+ mateStart=trChim[1-itr]->exons[0][EX_G];
+ mateStrand=(uint8_t) (trChim[1-itr]->Str!=trChim[1-itr]->exons[0][EX_iFrag]);
+ alignType=-10;
+ } else {//spanning chimeric alignment, could be PE or SE
+ mateChr=-1;mateStart=-1;mateStrand=0;//no need fot mate info unless this is the supplementary alignment
+ if (chimRepresent==(int)itr) {
+ alignType=-10; //this is representative part of chimeric alignment, record is as normal; if encompassing chimeric junction, both are recorded as normal
+ bamIrepr=bamN;
+ if (trChim[itr]->exons[0][EX_iFrag]!=trChim[1-itr]->exons[0][EX_iFrag]) {//the next mate is chimerically split
+ ++bamIrepr;
+ };
+//TODO check flags of SE split read
+// if (chimType==3) {
+// bamIrepr=bamN;
+// } else if (trChim[itr]->exons[0][EX_iFrag]==trChim[1-itr]->exons[0][EX_iFrag]) {
+//
+// };
+// bamIrepr=( (itr%2)==(trChim[itr]->Str) && chimType!=3) ? bamN+1 : bamN;//this is the mate that is chimerically split
+ } else {//"supplementary" chimeric segment
+ alignType=P.pCh.out.bamHardClip ? ( ( itr%2==trChim[itr]->Str ) ? -12 : -11) : -13 ; //right:left chimeric junction
+ bamIsuppl=bamN;
+ if (chimType==1) {//PE alignment, need mate info for the suppl
+ uint iex=0;
+ for (;iex<trChim[chimRepresent]->nExons-1;iex++) {
+ if (trChim[chimRepresent]->exons[iex][EX_iFrag]!=trChim[itr]->exons[0][EX_iFrag]) {
+ break;
+ };
+ };
+ mateChr=trChim[chimRepresent]->Chr;
+ mateStart=trChim[chimRepresent]->exons[iex][EX_G];
+ mateStrand=(uint8_t) (trChim[chimRepresent]->Str!=trChim[chimRepresent]->exons[iex][EX_iFrag]);
+ };
+ };
+
+ };
+
+ bamN+=RA->alignBAM(*trChim[itr], chimN, iTr, RA->mapGen.chrStart[trChim[itr]->Chr], mateChr, \
+ mateStart-RA->mapGen.chrStart[(mateChr<RA->mapGen.nChrReal ? mateChr : 0)], mateStrand, alignType, \
+ NULL, P.outSAMattrOrder, RA->outBAMoneAlign+bamN, RA->outBAMoneAlignNbytes+bamN);
+ bamBytesTotal+=RA->outBAMoneAlignNbytes[0]+RA->outBAMoneAlignNbytes[1];//outBAMoneAlignNbytes[1] = 0 if SE is recorded
+ };
+
+ //write all bam lines
+ for (int ii=0; ii<bamN; ii++) {//output all pieces
+ int tagI=-1;
+ if (ii==bamIrepr) {
+ tagI=bamIsuppl;
+ } else if (ii==bamIsuppl) {
+ tagI=bamIrepr;
+ };
+ if (tagI>=0) {
+ bam1_t *b;
+ b=bam_init1();
+ bam_read1_fromArray(RA->outBAMoneAlign[tagI], b);
+ uint8_t* auxp=bam_aux_get(b,"NM");
+ uint32_t auxv=bam_aux2i(auxp);
+ string tagSA1="SAZ"+RA->mapGen.chrName[b->core.tid]+','+to_string((uint)b->core.pos+1) +',' + ( (b->core.flag&0x10)==0 ? '+':'-') + \
+ ',' + bam_cigarString(b) + ',' + to_string((uint)b->core.qual) + ',' + to_string((uint)auxv) + ';' ;
+
+ memcpy( (void*) (RA->outBAMoneAlign[ii]+RA->outBAMoneAlignNbytes[ii]), tagSA1.c_str(), tagSA1.size()+1);//copy string including \0 at the end
+ RA->outBAMoneAlignNbytes[ii]+=tagSA1.size()+1;
+ * ( (uint32*) RA->outBAMoneAlign[ii] ) = RA->outBAMoneAlignNbytes[ii]-sizeof(uint32);
+ free(b); // don't use bam_destroy1(), because bam_read1_fromArray does not allocate memory for b->data
+ };
+
+ if (P.outBAMunsorted) RA->outBAMunsorted->unsortedOneAlign(RA->outBAMoneAlign[ii], RA->outBAMoneAlignNbytes[ii], ii>0 ? 0 : bamBytesTotal);
+ if (P.outBAMcoord) RA->outBAMcoord->coordOneAlign(RA->outBAMoneAlign[ii], RA->outBAMoneAlignNbytes[ii], (RA->iReadAll<<32) );
+ };
+
+};
=====================================
source/ChimericDetection.h
=====================================
@@ -17,16 +17,12 @@ class ChimericDetection {
uint nW, *nWinTr;
char** Read1;
Genome &outGen;
- uint *readLength;
public:
- uint chimN;
vector <ChimericAlign> chimAligns;
- bool chimRecord;
- int chimScoreBest;
ChimericDetection(Parameters &Pin, Transcript ***trAll, uint *nWinTr, char** Read1in, Genome &genomeIn, fstream *ostreamChimJunctionIn, ReadAlign *RA);
- bool chimericDetectionMult(uint nWin, uint *readLengthIn, int maxNonChimAlignScore, bool PEmerged_flag);
+ bool chimericDetectionMult(uint nWin, uint *readLength, int maxNonChimAlignScore, ReadAlign *PEunmergedRA);
fstream *ostreamChimJunction;
};
=====================================
source/ChimericDetection_chimericDetectionMult.cpp
=====================================
@@ -1,6 +1,7 @@
//#include "blocksOverlap.h"
#include "ChimericDetection.h"
#include "ChimericSegment.h"
+#include "ReadAlign.h"
int chimericAlignScore (ChimericSegment & seg1, ChimericSegment & seg2)
{
@@ -19,7 +20,7 @@ int chimericAlignScore (ChimericSegment & seg1, ChimericSegment & seg2)
};
/////////////////////////////////////////////////////////////
-bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int maxNonChimAlignScore, bool PEmerged_flag) {
+bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int maxNonChimAlignScore, ReadAlign *PEunmergedRA) {
// for (uint ii=0;ii<chimAligns.size();ii++) {//deallocate aligns
// if (chimAligns.at(ii).stitchingDone) {//al1,al2 were allocated
@@ -36,7 +37,8 @@ bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int max
};
chimAligns.clear();
- chimScoreBest=0;
+ int chimScoreBest=0;
+ std::size_t bestChimAlign=0; // points to element of chimAligns with highest chimScoreBest
int maxPossibleAlignScore = (int)(readLength[0]+readLength[1]);
int minScoreToConsider = P.pCh.scoreMin;
@@ -81,10 +83,11 @@ bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int max
if (chimAligns.back().chimScore > chimScoreBest) {
chimScoreBest=chimAligns.back().chimScore;
+ bestChimAlign = chimAligns.size()-1;
if ((chimScoreBest - (int)P.pCh.multimapScoreRange) > minScoreToConsider)
// best score increased, so subsequent alignment candidates must score higher
minScoreToConsider = chimScoreBest - (int)P.pCh.multimapScoreRange;
- }
+ };
} // endif stitched chimera survived.
else {
// al1, al2 allocated during stitching
@@ -101,7 +104,7 @@ bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int max
if (chimScoreBest==0)
return false;
- chimN=0;
+ uint chimN=0;
for (auto cAit=chimAligns.begin(); cAit<chimAligns.end(); cAit++) {
//scan all chimeras, find the number within score range
if (cAit->chimScore >= minScoreToConsider)
@@ -111,9 +114,24 @@ bool ChimericDetection::chimericDetectionMult(uint nW, uint *readLength, int max
if (chimN > P.pCh.multimapNmax) //too many loci
return false;
- for (auto cAit=chimAligns.begin(); cAit<chimAligns.end(); cAit++) {//output chimeras within score range
- if (cAit->chimScore >= minScoreToConsider)
- cAit->chimericJunctionOutput(*ostreamChimJunction, chimN, maxNonChimAlignScore, PEmerged_flag, chimScoreBest, maxPossibleAlignScore);
+ uint iTr = 0; //keep track of "HI" SAM attribute
+ for (std::size_t i = 0; i < chimAligns.size(); i++) {//output chimeras within score range
+ if (chimAligns[i].chimScore >= minScoreToConsider) {
+
+ if (P.pCh.out.junctions)
+ chimAligns[i].chimericJunctionOutput(*ostreamChimJunction, chimN, maxNonChimAlignScore, PEunmergedRA != NULL, chimScoreBest, maxPossibleAlignScore);
+
+ if (P.pCh.out.bam) {
+ // convert merged SE chimera to PE chimera if this is a merged chimera
+ if (PEunmergedRA != NULL) {
+ chimAligns[i].RA = PEunmergedRA;
+ chimAligns[i].RA->peOverlapChimericSEtoPE(chimAligns[i].al1, chimAligns[i].al2, chimAligns[i].al1, chimAligns[i].al2);
+ };
+ chimAligns[i].chimericBAMoutput(chimAligns[i].al1, chimAligns[i].al2, chimAligns[i].RA, iTr, chimN, i == bestChimAlign, P);
+ };
+ iTr++;
+
+ };
};
return chimN > 0;
=====================================
source/Makefile
=====================================
@@ -45,7 +45,7 @@ OBJECTS = soloInputFeatureUMI.o SoloFeature_countSmartSeq.o SoloFeature_redistri
GTF.o GTF_transcriptGeneSJ.o GTF_superTranscript.o SuperTranscriptome.o \
ReadAlign_outputAlignments.o \
ReadAlign.o STAR.o \
- SharedMemory.o PackedArray.o SuffixArrayFuns.o Parameters.o InOutStreams.o SequenceFuns.o Genome.o Stats.o \
+ SharedMemory.o PackedArray.o SuffixArrayFuns.o Parameters.o Parameters_samAttributes.o InOutStreams.o SequenceFuns.o Genome.o Stats.o \
Transcript.o Transcript_alignScore.o Transcript_generateCigarP.o Chain.o \
Transcript_variationAdjust.o Variation.o ReadAlign_waspMap.o \
ReadAlign_storeAligns.o ReadAlign_stitchPieces.o ReadAlign_multMapSelect.o ReadAlign_mapOneRead.o readLoad.o \
@@ -58,8 +58,8 @@ OBJECTS = soloInputFeatureUMI.o SoloFeature_countSmartSeq.o SoloFeature_redistri
ReadAlign_peOverlapMergeMap.o ReadAlign_mappedFilter.o \
ParametersChimeric_initialize.o ReadAlign_chimericDetection.o ReadAlign_chimericDetectionOld.o ReadAlign_chimericDetectionOldOutput.o\
ChimericDetection.o ChimericDetection_chimericDetectionMult.o ReadAlign_chimericDetectionPEmerged.o \
- ChimericSegment.cpp ChimericAlign.cpp ChimericAlign_chimericJunctionOutput.o ChimericAlign_chimericStitching.o \
stitchWindowAligns.o extendAlign.o stitchAlignToTranscript.o \
+ ChimericSegment.cpp ChimericAlign.cpp ChimericAlign_chimericJunctionOutput.o ChimericAlign_chimericBAMoutput.o ChimericAlign_chimericStitching.o \
Genome_genomeGenerate.o genomeParametersWrite.o genomeScanFastaFiles.o genomeSAindex.o \
Genome_insertSequences.o Genome_consensusSequence.o \
insertSeqSA.o funCompareUintAndSuffixes.o funCompareUintAndSuffixesMemcmp.o \
=====================================
source/Parameters.cpp
=====================================
@@ -1003,208 +1003,13 @@ void Parameters::inputParameters (int argInN, char* argIn[]) {//input parameters
exitWithError(errOut.str(),std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
};
- //outSAMattributes
- outSAMattrPresent.NH=false;//TODO re-write as class with constructor?
- outSAMattrPresent.HI=false;
- outSAMattrPresent.AS=false;
- outSAMattrPresent.NM=false;
- outSAMattrPresent.MD=false;
- outSAMattrPresent.nM=false;
- outSAMattrPresent.jM=false;
- outSAMattrPresent.jI=false;
- outSAMattrPresent.RG=false;
- outSAMattrPresent.MC=false;
- outSAMattrPresent.XS=false;
- outSAMattrPresent.ch=false;
-
- outSAMattrPresent.vA=false;
- outSAMattrPresent.vG=false;
- outSAMattrPresent.vW=false;
- outSAMattrPresent.rB=false;
- outSAMattrPresent.ha=false;
-
- outSAMattrPresent.CR=false;
- outSAMattrPresent.CY=false;
- outSAMattrPresent.UR=false;
- outSAMattrPresent.UY=false;
- outSAMattrPresent.CB=false;
- outSAMattrPresent.UB=false;
- outSAMattrPresent.GX=false;
- outSAMattrPresent.GN=false;
- outSAMattrPresent.sM=false;
- outSAMattrPresent.sS=false;
- outSAMattrPresent.sQ=false;
+ //SAM attributes
+ samAttributes();
- //for quant SAM output only NH and HI flags
- outSAMattrPresentQuant=outSAMattrPresent;
- outSAMattrPresentQuant.NH=true;
- outSAMattrPresentQuant.HI=true;
- outSAMattrOrderQuant.push_back(ATTR_NH);
- outSAMattrOrderQuant.push_back(ATTR_HI);
-
- vector<string> vAttr1;
- if (outSAMattributes.at(0)=="None") {
- } else if (outSAMattributes.at(0)=="All"){
- vAttr1={"NH","HI","AS","nM","NM","MD","jM","jI","MC","ch"};
- } else if (outSAMattributes.at(0)=="Standard"){
- vAttr1={"NH","HI","AS","nM"};
- } else {
- vAttr1=outSAMattributes;
- };
-
- for (uint ii=0;ii<vAttr1.size();ii++) {
- if (vAttr1.at(ii)== "NH") {
- outSAMattrOrder.push_back(ATTR_NH);
- outSAMattrPresent.NH=true;
- } else if (vAttr1.at(ii)== "HI") {
- outSAMattrOrder.push_back(ATTR_HI);
- outSAMattrPresent.HI=true;
- } else if (vAttr1.at(ii)== "AS") {
- outSAMattrOrder.push_back(ATTR_AS);
- outSAMattrPresent.AS=true;
- } else if (vAttr1.at(ii)== "NM") {
- outSAMattrOrder.push_back(ATTR_NM);
- outSAMattrPresent.NM=true;
- } else if (vAttr1.at(ii)== "MD") {
- outSAMattrOrder.push_back(ATTR_MD);
- outSAMattrPresent.MD=true;
- } else if (vAttr1.at(ii)== "nM") {
- outSAMattrOrder.push_back(ATTR_nM);
- outSAMattrPresent.nM=true;
- } else if (vAttr1.at(ii)== "jM") {
- outSAMattrOrder.push_back(ATTR_jM);
- outSAMattrPresent.jM=true;
- } else if (vAttr1.at(ii)== "jI") {
- outSAMattrOrder.push_back(ATTR_jI);
- outSAMattrPresent.jI=true;
- } else if (vAttr1.at(ii)== "vA") {
- outSAMattrOrder.push_back(ATTR_vA);
- outSAMattrPresent.vA=true;
- } else if (vAttr1.at(ii)== "vG") {
- outSAMattrOrder.push_back(ATTR_vG);
- outSAMattrPresent.vG=true;
- } else if (vAttr1.at(ii)== "vW") {
- outSAMattrOrder.push_back(ATTR_vW);
- outSAMattrPresent.vW=true;
- } else if (vAttr1.at(ii)== "ha") {
- outSAMattrOrder.push_back(ATTR_ha);
- outSAMattrPresent.ha=true;
- } else if (vAttr1.at(ii)== "RG") {
- outSAMattrOrder.push_back(ATTR_RG);
- outSAMattrOrderQuant.push_back(ATTR_RG);
- outSAMattrPresent.RG=true;
- } else if (vAttr1.at(ii)== "rB") {
- outSAMattrOrder.push_back(ATTR_rB);
- outSAMattrOrderQuant.push_back(ATTR_rB);
- outSAMattrPresent.rB=true;
- } else if (vAttr1.at(ii)== "ch") {
- outSAMattrOrder.push_back(ATTR_ch);
- outSAMattrOrderQuant.push_back(ATTR_ch);
- outSAMattrPresent.ch=true;
- } else if (vAttr1.at(ii)== "MC") {
- outSAMattrOrder.push_back(ATTR_MC);
- outSAMattrOrderQuant.push_back(ATTR_MC);
- outSAMattrPresent.MC=true;
- } else if (vAttr1.at(ii)== "CR") {
- outSAMattrOrder.push_back(ATTR_CR);
- outSAMattrOrderQuant.push_back(ATTR_CR);
- outSAMattrPresent.CR=true;
- } else if (vAttr1.at(ii)== "CY") {
- outSAMattrOrder.push_back(ATTR_CY);
- outSAMattrOrderQuant.push_back(ATTR_CY);
- outSAMattrPresent.CY=true;
- } else if (vAttr1.at(ii)== "UR") {
- outSAMattrOrder.push_back(ATTR_UR);
- outSAMattrOrderQuant.push_back(ATTR_UR);
- outSAMattrPresent.UR=true;
- } else if (vAttr1.at(ii)== "UY") {
- outSAMattrOrder.push_back(ATTR_UY);
- outSAMattrOrderQuant.push_back(ATTR_UY);
- outSAMattrPresent.UY=true;
- } else if (vAttr1.at(ii)== "CB") {
- outSAMattrOrder.push_back(ATTR_CB);
- outSAMattrOrderQuant.push_back(ATTR_CB);
- outSAMattrPresent.CB=true;
- } else if (vAttr1.at(ii)== "UB") {
- outSAMattrOrder.push_back(ATTR_UB);
- outSAMattrOrderQuant.push_back(ATTR_UB);
- outSAMattrPresent.UB=true;
- } else if (vAttr1.at(ii)== "GX") {
- outSAMattrOrder.push_back(ATTR_GX);
- outSAMattrOrderQuant.push_back(ATTR_GX);
- outSAMattrPresent.GX=true;
- } else if (vAttr1.at(ii)== "GN") {
- outSAMattrOrder.push_back(ATTR_GN);
- outSAMattrOrderQuant.push_back(ATTR_GN);
- outSAMattrPresent.GN=true;
- } else if (vAttr1.at(ii)== "sM") {
- outSAMattrOrder.push_back(ATTR_sM);
- outSAMattrOrderQuant.push_back(ATTR_sM);
- outSAMattrPresent.sM=true;
- } else if (vAttr1.at(ii)== "sS") {
- outSAMattrOrder.push_back(ATTR_sS);
- outSAMattrOrderQuant.push_back(ATTR_sS);
- outSAMattrPresent.sS=true;
- } else if (vAttr1.at(ii)== "sQ") {
- outSAMattrOrder.push_back(ATTR_sQ);
- outSAMattrOrderQuant.push_back(ATTR_sQ);
- outSAMattrPresent.sQ=true;
- } else if (vAttr1.at(ii)== "XS") {
- outSAMattrOrder.push_back(ATTR_XS);
- outSAMattrPresent.XS=true;
- if (outSAMstrandField.type!=1) {
- inOut->logMain << "WARNING --outSAMattributes contains XS, therefore STAR will use --outSAMstrandField intronMotif" <<endl;
- outSAMstrandField.type=1;
- };
- } else {
- ostringstream errOut;
- errOut <<"EXITING because of FATAL INPUT ERROR: unknown/unimplemented SAM atrribute (tag): "<<vAttr1.at(ii) <<"\n";
- errOut <<"SOLUTION: re-run STAR with --outSAMattributes that contains only implemented attributes\n";
- exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
- };
- };
-
- if (!var.yes && (outSAMattrPresent.vA | outSAMattrPresent.vG)) {
- ostringstream errOut;
- errOut <<"EXITING because of fatal PARAMETER error: --outSAMattributes contains vA and/or vG tag(s), but --varVCFfile is not set\n";
- errOut <<"SOLUTION: re-run STAR with a --varVCFfile option, or without vA/vG tags in --outSAMattributes\n";
- exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
- };
-
- if (!wasp.yes && outSAMattrPresent.vW) {
- ostringstream errOut;
- errOut <<"EXITING because of fatal PARAMETER error: --outSAMattributes contains vW tag, but --waspOutputMode is not set\n";
- errOut <<"SOLUTION: re-run STAR with a --waspOutputMode option, or without vW tags in --outSAMattributes\n";
- exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
- };
-
- if (outSAMattrRGline[0]!="-" && !outSAMattrPresent.RG) {
- outSAMattrOrder.push_back(ATTR_RG);
- outSAMattrOrderQuant.push_back(ATTR_RG);
- outSAMattrPresent.RG=true;
- inOut->logMain << "WARNING --outSAMattrRG defines a read group, therefore STAR will output RG attribute" <<endl;
- } else if (outSAMattrRG.size()==0 && outSAMattrPresent.RG) {
- ostringstream errOut;
- errOut <<"EXITING because of FATAL INPUT ERROR: --outSAMattributes contains RG tag, but --outSAMattrRGline is not set\n";
- errOut <<"SOLUTION: re-run STAR with a valid read group parameter --outSAMattrRGline\n";
- exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
- };
-
- if (outSAMstrandField.type==1 && !outSAMattrPresent.XS) {
- outSAMattrOrder.push_back(ATTR_XS);
- inOut->logMain << "WARNING --outSAMstrandField=intronMotif, therefore STAR will output XS attribute" <<endl;
- };
-
- if (wasp.yes && !outSAMattrPresent.vW) {
- outSAMattrOrder.push_back(ATTR_vW);
- outSAMattrOrderQuant.push_back(ATTR_vW);
- outSAMattrPresent.vW=true;
- inOut->logMain << "WARNING --waspOutputMode is set, therefore STAR will output vW attribute" <<endl;
- };
-
//solo
pSolo.initialize(this);
+ //alignEnds
alignEndsType.ext[0][0]=false;
alignEndsType.ext[0][1]=false;
alignEndsType.ext[1][0]=false;
@@ -1392,8 +1197,6 @@ void Parameters::inputParameters (int argInN, char* argIn[]) {//input parameters
inOut->logMain << "Finished loading and checking parameters\n" <<flush;
};
-
-
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Parameters::scanAllLines (istream &streamIn, int inputLevel, int inputLevelRequested) {//scan
// istringstream stringInStream (stringIn);
@@ -1468,6 +1271,3 @@ int Parameters::scanOneLine (string &lineIn, int inputLevel, int inputLevelReque
};
return 0;
};
-
-
-
=====================================
source/Parameters.h
=====================================
@@ -348,6 +348,7 @@ class Parameters {
void readFilesInit();
void closeReadsFiles();
void readSAMheader(const string readFilesCommandString, const vector<string> readFilesNames);
-
+ void samAttributes();
+ void samAttrRequiresBAM(bool attrYes, string attrTag);
};
#endif // Parameters.h
=====================================
source/ParametersChimeric_initialize.cpp
=====================================
@@ -80,10 +80,10 @@ void ParametersChimeric::initialize(Parameters *pPin)
exitWithError(errOut.str(), std::cerr, pP->inOut->logMain, EXIT_CODE_PARAMETER, *pP);
};
- if (multimapNmax>0 && (out.bam || out.samOld)) {
+ if (multimapNmax>0 && out.samOld) {
ostringstream errOut;
- errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax > 0 (new chimeric detection) presently only works with --chimOutType Junctions\n";
- errOut <<"SOLUTION: re-run with --chimOutType Junctions\n";
+ errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax > 0 (new chimeric detection) presently only works with --chimOutType Junctions/WithinBAM\n";
+ errOut <<"SOLUTION: re-run with --chimOutType Junctions/WithinBAM\n";
exitWithError(errOut.str(), std::cerr, pP->inOut->logMain, EXIT_CODE_PARAMETER, *pP);
};
=====================================
source/Parameters_samAttributes.cpp
=====================================
@@ -0,0 +1,240 @@
+#include "Parameters.h"
+#include "ErrorWarning.h"
+
+void Parameters::samAttributes(){//everything related to SAM attributes
+ //outSAMattributes
+ outSAMattrPresent.NH=false;//TODO re-write as class with constructor?
+ outSAMattrPresent.HI=false;
+ outSAMattrPresent.AS=false;
+ outSAMattrPresent.NM=false;
+ outSAMattrPresent.MD=false;
+ outSAMattrPresent.nM=false;
+ outSAMattrPresent.jM=false;
+ outSAMattrPresent.jI=false;
+ outSAMattrPresent.RG=false;
+ outSAMattrPresent.MC=false;
+ outSAMattrPresent.XS=false;
+ outSAMattrPresent.ch=false;
+
+ outSAMattrPresent.vA=false;
+ outSAMattrPresent.vG=false;
+ outSAMattrPresent.vW=false;
+ outSAMattrPresent.rB=false;
+ outSAMattrPresent.ha=false;
+
+ outSAMattrPresent.CR=false;
+ outSAMattrPresent.CY=false;
+ outSAMattrPresent.UR=false;
+ outSAMattrPresent.UY=false;
+ outSAMattrPresent.CB=false;
+ outSAMattrPresent.UB=false;
+ outSAMattrPresent.GX=false;
+ outSAMattrPresent.GN=false;
+ outSAMattrPresent.sM=false;
+ outSAMattrPresent.sS=false;
+ outSAMattrPresent.sQ=false;
+
+ //for quant SAM output only NH and HI flags
+ outSAMattrPresentQuant=outSAMattrPresent;
+ outSAMattrPresentQuant.NH=true;
+ outSAMattrPresentQuant.HI=true;
+ outSAMattrOrderQuant.push_back(ATTR_NH);
+ outSAMattrOrderQuant.push_back(ATTR_HI);
+
+ vector<string> vAttr1;
+ if (outSAMattributes.at(0)=="None") {
+ } else if (outSAMattributes.at(0)=="All"){
+ vAttr1={"NH","HI","AS","nM","NM","MD","jM","jI","MC","ch"};
+ } else if (outSAMattributes.at(0)=="Standard"){
+ vAttr1={"NH","HI","AS","nM"};
+ } else {
+ vAttr1=outSAMattributes;
+ };
+
+ for (uint ii=0;ii<vAttr1.size();ii++) {
+ if (vAttr1.at(ii)== "NH") {
+ outSAMattrOrder.push_back(ATTR_NH);
+ outSAMattrPresent.NH=true;
+ } else if (vAttr1.at(ii)== "HI") {
+ outSAMattrOrder.push_back(ATTR_HI);
+ outSAMattrPresent.HI=true;
+ } else if (vAttr1.at(ii)== "AS") {
+ outSAMattrOrder.push_back(ATTR_AS);
+ outSAMattrPresent.AS=true;
+ } else if (vAttr1.at(ii)== "NM") {
+ outSAMattrOrder.push_back(ATTR_NM);
+ outSAMattrPresent.NM=true;
+ } else if (vAttr1.at(ii)== "MD") {
+ outSAMattrOrder.push_back(ATTR_MD);
+ outSAMattrPresent.MD=true;
+ } else if (vAttr1.at(ii)== "nM") {
+ outSAMattrOrder.push_back(ATTR_nM);
+ outSAMattrPresent.nM=true;
+ } else if (vAttr1.at(ii)== "jM") {
+ outSAMattrOrder.push_back(ATTR_jM);
+ outSAMattrPresent.jM=true;
+ } else if (vAttr1.at(ii)== "jI") {
+ outSAMattrOrder.push_back(ATTR_jI);
+ outSAMattrPresent.jI=true;
+ } else if (vAttr1.at(ii)== "vA") {
+ outSAMattrOrder.push_back(ATTR_vA);
+ outSAMattrPresent.vA=true;
+ } else if (vAttr1.at(ii)== "vG") {
+ outSAMattrOrder.push_back(ATTR_vG);
+ outSAMattrPresent.vG=true;
+ } else if (vAttr1.at(ii)== "vW") {
+ outSAMattrOrder.push_back(ATTR_vW);
+ outSAMattrPresent.vW=true;
+ } else if (vAttr1.at(ii)== "ha") {
+ outSAMattrOrder.push_back(ATTR_ha);
+ outSAMattrPresent.ha=true;
+ } else if (vAttr1.at(ii)== "RG") {
+ outSAMattrOrder.push_back(ATTR_RG);
+ outSAMattrOrderQuant.push_back(ATTR_RG);
+ outSAMattrPresent.RG=true;
+ } else if (vAttr1.at(ii)== "rB") {
+ outSAMattrOrder.push_back(ATTR_rB);
+ outSAMattrOrderQuant.push_back(ATTR_rB);
+ outSAMattrPresent.rB=true;
+ } else if (vAttr1.at(ii)== "ch") {
+ outSAMattrOrder.push_back(ATTR_ch);
+ outSAMattrOrderQuant.push_back(ATTR_ch);
+ outSAMattrPresent.ch=true;
+ } else if (vAttr1.at(ii)== "MC") {
+ outSAMattrOrder.push_back(ATTR_MC);
+ outSAMattrOrderQuant.push_back(ATTR_MC);
+ outSAMattrPresent.MC=true;
+ } else if (vAttr1.at(ii)== "CR") {
+ outSAMattrOrder.push_back(ATTR_CR);
+ outSAMattrOrderQuant.push_back(ATTR_CR);
+ outSAMattrPresent.CR=true;
+ } else if (vAttr1.at(ii)== "CY") {
+ outSAMattrOrder.push_back(ATTR_CY);
+ outSAMattrOrderQuant.push_back(ATTR_CY);
+ outSAMattrPresent.CY=true;
+ } else if (vAttr1.at(ii)== "UR") {
+ outSAMattrOrder.push_back(ATTR_UR);
+ outSAMattrOrderQuant.push_back(ATTR_UR);
+ outSAMattrPresent.UR=true;
+ } else if (vAttr1.at(ii)== "UY") {
+ outSAMattrOrder.push_back(ATTR_UY);
+ outSAMattrOrderQuant.push_back(ATTR_UY);
+ outSAMattrPresent.UY=true;
+ } else if (vAttr1.at(ii)== "CB") {
+ outSAMattrOrder.push_back(ATTR_CB);
+ outSAMattrOrderQuant.push_back(ATTR_CB);
+ outSAMattrPresent.CB=true;
+ } else if (vAttr1.at(ii)== "UB") {
+ outSAMattrOrder.push_back(ATTR_UB);
+ outSAMattrOrderQuant.push_back(ATTR_UB);
+ outSAMattrPresent.UB=true;
+ } else if (vAttr1.at(ii)== "GX") {
+ outSAMattrOrder.push_back(ATTR_GX);
+ outSAMattrOrderQuant.push_back(ATTR_GX);
+ outSAMattrPresent.GX=true;
+ } else if (vAttr1.at(ii)== "GN") {
+ outSAMattrOrder.push_back(ATTR_GN);
+ outSAMattrOrderQuant.push_back(ATTR_GN);
+ outSAMattrPresent.GN=true;
+ } else if (vAttr1.at(ii)== "sM") {
+ outSAMattrOrder.push_back(ATTR_sM);
+ outSAMattrOrderQuant.push_back(ATTR_sM);
+ outSAMattrPresent.sM=true;
+ } else if (vAttr1.at(ii)== "sS") {
+ outSAMattrOrder.push_back(ATTR_sS);
+ outSAMattrOrderQuant.push_back(ATTR_sS);
+ outSAMattrPresent.sS=true;
+ } else if (vAttr1.at(ii)== "sQ") {
+ outSAMattrOrder.push_back(ATTR_sQ);
+ outSAMattrOrderQuant.push_back(ATTR_sQ);
+ outSAMattrPresent.sQ=true;
+ } else if (vAttr1.at(ii)== "XS") {
+ outSAMattrOrder.push_back(ATTR_XS);
+ outSAMattrPresent.XS=true;
+ if (outSAMstrandField.type!=1) {
+ inOut->logMain << "WARNING --outSAMattributes contains XS, therefore STAR will use --outSAMstrandField intronMotif" <<endl;
+ outSAMstrandField.type=1;
+ };
+ } else {
+ ostringstream errOut;
+ errOut <<"EXITING because of FATAL INPUT ERROR: unknown/unimplemented SAM atrribute (tag): "<<vAttr1.at(ii) <<"\n";
+ errOut <<"SOLUTION: re-run STAR with --outSAMattributes that contains only implemented attributes\n";
+ exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+ };
+ };
+
+ if (!var.yes && (outSAMattrPresent.vA | outSAMattrPresent.vG)) {
+ ostringstream errOut;
+ errOut <<"EXITING because of fatal PARAMETER error: --outSAMattributes contains vA and/or vG tag(s), but --varVCFfile is not set\n";
+ errOut <<"SOLUTION: re-run STAR with a --varVCFfile option, or without vA/vG tags in --outSAMattributes\n";
+ exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+ };
+
+ if (!wasp.yes && outSAMattrPresent.vW) {
+ ostringstream errOut;
+ errOut <<"EXITING because of fatal PARAMETER error: --outSAMattributes contains vW tag, but --waspOutputMode is not set\n";
+ errOut <<"SOLUTION: re-run STAR with a --waspOutputMode option, or without vW tags in --outSAMattributes\n";
+ exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+ };
+
+ if (outSAMattrRGline[0]!="-" && !outSAMattrPresent.RG) {
+ outSAMattrOrder.push_back(ATTR_RG);
+ outSAMattrOrderQuant.push_back(ATTR_RG);
+ outSAMattrPresent.RG=true;
+ inOut->logMain << "WARNING --outSAMattrRG defines a read group, therefore STAR will output RG attribute" <<endl;
+ } else if (outSAMattrRG.size()==0 && outSAMattrPresent.RG) {
+ ostringstream errOut;
+ errOut <<"EXITING because of FATAL INPUT ERROR: --outSAMattributes contains RG tag, but --outSAMattrRGline is not set\n";
+ errOut <<"SOLUTION: re-run STAR with a valid read group parameter --outSAMattrRGline\n";
+ exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+ };
+
+ if (outSAMstrandField.type==1 && !outSAMattrPresent.XS) {
+ outSAMattrOrder.push_back(ATTR_XS);
+ inOut->logMain << "WARNING --outSAMstrandField=intronMotif, therefore STAR will output XS attribute" <<endl;
+ };
+
+ if (wasp.yes && !outSAMattrPresent.vW) {
+ outSAMattrOrder.push_back(ATTR_vW);
+ outSAMattrOrderQuant.push_back(ATTR_vW);
+ outSAMattrPresent.vW=true;
+ inOut->logMain << "WARNING --waspOutputMode is set, therefore STAR will output vW attribute" <<endl;
+ };
+
+ //TODO: BAM-only flag should all have values > 100
+ samAttrRequiresBAM(outSAMattrPresent.ch, "ch");
+ samAttrRequiresBAM(outSAMattrPresent.CR, "CR");
+ samAttrRequiresBAM(outSAMattrPresent.CY, "CY");
+ samAttrRequiresBAM(outSAMattrPresent.UR, "UR");
+ samAttrRequiresBAM(outSAMattrPresent.UY, "UY");
+ samAttrRequiresBAM(outSAMattrPresent.CB, "CB");
+ samAttrRequiresBAM(outSAMattrPresent.UB, "UB");
+ samAttrRequiresBAM(outSAMattrPresent.sM, "sM");
+ samAttrRequiresBAM(outSAMattrPresent.sS, "sS");
+ samAttrRequiresBAM(outSAMattrPresent.sQ, "sQ");
+ samAttrRequiresBAM(outSAMattrPresent.rB, "rB");
+ samAttrRequiresBAM(outSAMattrPresent.vG, "vG");
+ samAttrRequiresBAM(outSAMattrPresent.vA, "vA");
+ samAttrRequiresBAM(outSAMattrPresent.vW, "vW");
+ samAttrRequiresBAM(outSAMattrPresent.GX, "GX");
+ samAttrRequiresBAM(outSAMattrPresent.GN, "GN");
+
+ if (outSAMattrPresent.GX || outSAMattrPresent.GN) {//turn on quantification
+ quant.gene.yes=true;
+ quant.yes = true;
+ };
+};
+
+void Parameters:: samAttrRequiresBAM(bool attrYes, string attrTag) {
+ if (!attrYes)
+ return; // nothing to do
+ if (!outBAMunsorted && !outBAMcoord) {
+ ostringstream errOut;
+ errOut <<"EXITING because of fatal PARAMETER error: --outSAMattributes contains "<< attrTag <<" tag, which requires BAM output.\n";
+ errOut <<"SOLUTION: re-run STAR with --outSAMtype BAM Unsorted (and/or) SortedByCoordinate option, or without "<< attrTag <<" tag in --outSAMattributes\n";
+ exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+ };
+
+ if (outSAMbool)
+ inOut->logMain << "WARNING: --outSAMattributes contains "<< attrTag <<" tag. It will be output into BAM file(s), but not SAM file.\n";
+};
\ No newline at end of file
=====================================
source/ReadAlign.h
=====================================
@@ -62,11 +62,17 @@ class ReadAlign {
ReadAlign *peMergeRA; //ReadAlign for merged PE mates
ChimericDetection *chimDet;
+ void peOverlapChimericSEtoPE(const Transcript *seTrIn1, const Transcript *seTrIn2, Transcript *peTrOut1, Transcript *peTrOut2);
SoloRead *soloRead; //counts reads per CB per and outputs CB/UMI/gene into file, per thread
SpliceGraph *splGraph;
+ //input,output
+ char** outBAMoneAlign;
+ uint* outBAMoneAlignNbytes;
+ int alignBAM(Transcript const &trOut, uint nTrOut, uint iTrOut, uint trChrStart, uint mateChr, uint mateStart, char mateStrand, int unmapType, bool *mateMapped, vector<int> outSAMattrOrder, char** outBAMarray, uint* outBAMarrayN);
+
private:
Parameters& P; //pointer to the parameters, will be initialized on construction
@@ -82,9 +88,6 @@ class ReadAlign {
//input,output
- char** outBAMoneAlign;
- uint* outBAMoneAlignNbytes;
-
ostringstream samStreamCIGAR, samStreamSJmotif, samStreamSJintron;
vector <string> matesCIGAR;
@@ -165,10 +168,11 @@ class ReadAlign {
void storeAligns (uint iDir, uint Shift, uint Nrep, uint L, uint indStartEnd[2], uint iFrag);
bool outputTranscript(Transcript *trOut, uint nTrOut, ofstream *outBED);
+
uint64 outputTranscriptSAM(Transcript const &trOut, uint nTrOut, uint iTrOut, uint mateChr, uint mateStart, char mateStrand, int unmapType, bool *mateMapped, ostream *outStream);
+
uint64 outputSpliceGraphSAM(Transcript const &trOut, uint nTrOut, uint iTrOut, ostream *outStream);
- int alignBAM(Transcript const &trOut, uint nTrOut, uint iTrOut, uint trChrStart, uint mateChr, uint mateStart, char mateStrand, int unmapType, bool *mateMapped, vector<int> outSAMattrOrder, char** outBAMarray, uint* outBAMarrayN);
void samAttrNM_MD (Transcript const &trOut, uint iEx1, uint iEx2, uint &tagNM, string &tagMD);
void outputTranscriptSJ(Transcript const &trOut, uint nTrOut, OutSJ *outStream, uint sjReadStartN );
=====================================
source/ReadAlignChunk.cpp
=====================================
@@ -77,13 +77,15 @@ ReadAlignChunk::ReadAlignChunk(Parameters& Pin, Genome &genomeIn, Transcriptome
chunkFstreamOpen(P.outFileTmp + "/Chimeric.out.junction.thread", iChunk, *RA->chunkOutChimJunction);
};
};
+
if (P.outReadsUnmapped=="Fastx" ) {
- chunkFstreamOpen(P.outFileTmp + "/Unmapped.out.mate1.thread",iChunk, RA->chunkOutUnmappedReadsStream[0]);
- if (P.readNmatesIn==2) chunkFstreamOpen(P.outFileTmp + "/Unmapped.out.mate2.thread",iChunk, RA->chunkOutUnmappedReadsStream[1]);
+ for (uint32 imate=0; imate < P.readNmatesIn; imate++)
+ chunkFstreamOpen(P.outFileTmp + "/Unmapped.out.mate"+ to_string(imate) +".thread",iChunk, RA->chunkOutUnmappedReadsStream[imate]);
};
+
if (P.outFilterType=="BySJout") {
chunkFstreamOpen(P.outFileTmp + "/FilterBySJoutFiles.mate1.thread",iChunk, RA->chunkOutFilterBySJoutFiles[0]);
- if (P.readNmates==2) chunkFstreamOpen(P.outFileTmp + "/FilterBySJoutFiles.mate2.thread",iChunk, RA->chunkOutFilterBySJoutFiles[1]);
+ if (P.readNmates==2) chunkFstreamOpen(P.outFileTmp + "/FilterBySJoutFiles.mate2.thread",iChunk, RA->chunkOutFilterBySJoutFiles[1]); //here we do not output barcode read
};
if (P.wasp.yes) {
@@ -94,6 +96,8 @@ ReadAlignChunk::ReadAlignChunk(Parameters& Pin, Genome &genomeIn, Transcriptome
delete RA->peMergeRA->chunkOutChimJunction;
RA->peMergeRA->chunkOutChimJunction=RA->chunkOutChimJunction;//point to the same out-stream
RA->peMergeRA->chimDet->ostreamChimJunction=RA->peMergeRA->chunkOutChimJunction;
+ RA->peMergeRA->outBAMunsorted=RA->outBAMunsorted;
+ RA->peMergeRA->outBAMcoord=RA->outBAMcoord;
};
};
=====================================
source/ReadAlign_alignBAM.cpp
=====================================
@@ -184,11 +184,8 @@ int ReadAlign::alignBAM(Transcript const &trOut, uint nTrOut, uint iTrOut, uint
if (readFilter=='Y') samFLAG|=0x200; //not passing quality control
- if (alignType==-11 || alignType==-12 || alignType==-13) {
- samFLAG|=0x800; //mark chimeric alignments
- } else {//only non-chimeric alignments will be marked as non-primary, since chimeric are already marked with 0x800
- if (!trOut.primaryFlag) samFLAG|=0x100;//mark not primary align
- };
+ if (alignType==-11 || alignType==-12 || alignType==-13) samFLAG|=0x800; //mark chimeric alignments
+ if (!trOut.primaryFlag) samFLAG|=0x100; //mark not primary align
iEx1 = (imate==0 ? 0 : iExMate+1);
iEx2 = (imate==0 ? iExMate : trOut.nExons-1);
=====================================
source/ReadAlign_chimericDetection.cpp
=====================================
@@ -46,7 +46,7 @@ void ReadAlign::chimericDetection() {
chimRecord=chimericDetectionOld();
chimericDetectionOldOutput();
} else if (trBest->maxScore <= (int) (readLength[0]+readLength[1]) - (int) P.pCh.nonchimScoreDropMin) {//require big enough drop in the best score
- chimRecord=chimDet->chimericDetectionMult(nW, readLength, trBest->maxScore, false);
+ chimRecord=chimDet->chimericDetectionMult(nW, readLength, trBest->maxScore, NULL);
};
if ( chimRecord ) {
=====================================
source/ReadAlign_chimericDetectionOldOutput.cpp
=====================================
@@ -1,5 +1,6 @@
#include "ReadAlign.h"
#include "BAMfunctions.h"
+#include "ChimericAlign.h"
void ReadAlign::chimericDetectionOldOutput() {
@@ -7,110 +8,31 @@ void ReadAlign::chimericDetectionOldOutput() {
return;
};
- chimN=2; //this is hard-coded for now
//re-calculate the score for chimeric transcripts
trChim[0].alignScore(Read1, mapGen.G, P);
trChim[1].alignScore(Read1, mapGen.G, P);
- int chimRepresent=-999, chimType=0;
- if (trChim[0].exons[0][EX_iFrag]!=trChim[0].exons[trChim[0].nExons-1][EX_iFrag]) {//tr0 has both mates
- chimRepresent = 0;
- chimType = 1;
- trChim[0].primaryFlag=true;//paired portion is primary
- trChim[1].primaryFlag=false;
- } else if (trChim[1].exons[0][EX_iFrag]!=trChim[1].exons[trChim[1].nExons-1][EX_iFrag]) {//tr1 has both mates
- chimRepresent = 1;
- chimType = 1;
- trChim[1].primaryFlag=true;//paired portion is primary
- trChim[0].primaryFlag=false;
- } else if (trChim[0].exons[0][EX_iFrag]!=trChim[1].exons[0][EX_iFrag]) {//tr0 and tr1 are single different mates
- chimRepresent = -1;
- chimType = 2;
- trChim[0].primaryFlag=true;
- trChim[1].primaryFlag=true;
- } else {//two chimeric segments are on the same mate - this can only happen for single-end reads
- chimRepresent = (trChim[0].maxScore > trChim[1].maxScore) ? 0 : 1;
- chimType = 3;
- trChim[chimRepresent].primaryFlag=true;
- trChim[1-chimRepresent].primaryFlag=false;
- };
-
- if (P.pCh.out.bam) {//BAM output
- int alignType, bamN=0, bamIsuppl=-1, bamIrepr=-1;
- uint bamBytesTotal=0;//estimate of the total size of all bam records, for output buffering
- uint mateChr,mateStart;
- uint8_t mateStrand;
- for (uint itr=0;itr<chimN;itr++) {//generate bam for all chimeric pieces
- if (chimType==2) {//PE, encompassing
- mateChr=trChim[1-itr].Chr;
- mateStart=trChim[1-itr].exons[0][EX_G];
- mateStrand=(uint8_t) (trChim[1-itr].Str!=trChim[1-itr].exons[0][EX_iFrag]);
- alignType=-10;
- } else {//spanning chimeric alignment, could be PE or SE
- mateChr=-1;mateStart=-1;mateStrand=0;//no need fot mate info unless this is the supplementary alignment
- if (chimRepresent==(int)itr) {
- alignType=-10; //this is representative part of chimeric alignment, record is as normal; if encompassing chimeric junction, both are recorded as normal
- bamIrepr=bamN;
- if (trChim[itr].exons[0][EX_iFrag]!=trChim[1-itr].exons[0][EX_iFrag]) {//the next mate is chimerically split
- ++bamIrepr;
- };
-// if (chimType==3) {
-// bamIrepr=bamN;
-// } else if (trChim[itr].exons[0][EX_iFrag]==trChim[1-itr].exons[0][EX_iFrag]) {
-//
-// };
-// bamIrepr=( (itr%2)==(trChim[itr].Str) && chimType!=3) ? bamN+1 : bamN;//this is the mate that is chimerically split
- } else {//"supplementary" chimeric segment
- alignType=P.pCh.out.bamHardClip ? ( ( itr%2==trChim[itr].Str ) ? -12 : -11) : -13 ; //right:left chimeric junction
- bamIsuppl=bamN;
- if (chimType==1) {//PE alignment, need mate info for the suppl
- uint iex=0;
- for (;iex<trChim[chimRepresent].nExons-1;iex++) {
- if (trChim[chimRepresent].exons[iex][EX_iFrag]!=trChim[itr].exons[0][EX_iFrag]) {
- break;
- };
- };
- mateChr=trChim[chimRepresent].Chr;
- mateStart=trChim[chimRepresent].exons[iex][EX_G];
- mateStrand=(uint8_t) (trChim[chimRepresent].Str!=trChim[chimRepresent].exons[iex][EX_iFrag]);
- };
- };
-
- };
+ if (P.pCh.out.bam) //BAM output
+ ChimericAlign::chimericBAMoutput(&trChim[0], &trChim[1], this, 0, 1, true, P);
- bamN+=alignBAM(trChim[itr], 1, 0, mapGen.chrStart[trChim[itr].Chr], mateChr, mateStart-mapGen.chrStart[(mateChr<mapGen.nChrReal ? mateChr : 0)], mateStrand, \
- alignType, NULL, P.outSAMattrOrder, outBAMoneAlign+bamN, outBAMoneAlignNbytes+bamN);
- bamBytesTotal+=outBAMoneAlignNbytes[0]+outBAMoneAlignNbytes[1];//outBAMoneAlignNbytes[1] = 0 if SE is recorded
- };
-
- //write all bam lines
- for (int ii=0; ii<bamN; ii++) {//output all pieces
- int tagI=-1;
- if (ii==bamIrepr) {
- tagI=bamIsuppl;
- } else if (ii==bamIsuppl) {
- tagI=bamIrepr;
- };
- if (tagI>=0) {
- bam1_t *b;
- b=bam_init1();
- bam_read1_fromArray(outBAMoneAlign[tagI], b);
- uint8_t* auxp=bam_aux_get(b,"NM");
- uint32_t auxv=bam_aux2i(auxp);
- string tagSA1="SAZ"+mapGen.chrName[b->core.tid]+','+to_string((uint)b->core.pos+1) +',' + ( (b->core.flag&0x10)==0 ? '+':'-') + \
- ',' + bam_cigarString(b) + ',' + to_string((uint)b->core.qual) + ',' + to_string((uint)auxv) + ';' ;
-
- memcpy( (void*) (outBAMoneAlign[ii]+outBAMoneAlignNbytes[ii]), tagSA1.c_str(), tagSA1.size()+1);//copy string including \0 at the end
- outBAMoneAlignNbytes[ii]+=tagSA1.size()+1;
- * ( (uint32*) outBAMoneAlign[ii] ) = outBAMoneAlignNbytes[ii]-sizeof(uint32);
- };
+ if (P.pCh.out.samOld) {
- if (P.outBAMunsorted) outBAMunsorted->unsortedOneAlign(outBAMoneAlign[ii], outBAMoneAlignNbytes[ii], ii>0 ? 0 : bamBytesTotal);
- if (P.outBAMcoord) outBAMcoord->coordOneAlign(outBAMoneAlign[ii], outBAMoneAlignNbytes[ii], (iReadAll<<32) );
+ chimN=2; //this is hard-coded for now
+ if (trChim[0].exons[0][EX_iFrag]!=trChim[0].exons[trChim[0].nExons-1][EX_iFrag]) {//tr0 has both mates
+ trChim[0].primaryFlag=true;//paired portion is primary
+ trChim[1].primaryFlag=false;
+ } else if (trChim[1].exons[0][EX_iFrag]!=trChim[1].exons[trChim[1].nExons-1][EX_iFrag]) {//tr1 has both mates
+ trChim[1].primaryFlag=true;//paired portion is primary
+ trChim[0].primaryFlag=false;
+ } else if (trChim[0].exons[0][EX_iFrag]!=trChim[1].exons[0][EX_iFrag]) {//tr0 and tr1 are single different mates
+ trChim[0].primaryFlag=true;
+ trChim[1].primaryFlag=true;
+ } else {//two chimeric segments are on the same mate - this can only happen for single-end reads
+ int chimRepresent = (trChim[0].maxScore > trChim[1].maxScore) ? 0 : 1;
+ trChim[chimRepresent].primaryFlag=true;
+ trChim[1-chimRepresent].primaryFlag=false;
};
- };
- if (P.pCh.out.samOld) {
for (uint iTr=0;iTr<chimN;iTr++)
{//write all chimeric pieces to Chimeric.out.sam/junction
if (P.readNmates==2) {//PE: need mate info
@@ -149,4 +71,4 @@ void ReadAlign::chimericDetectionOldOutput() {
};
return;
-};
\ No newline at end of file
+};
=====================================
source/ReadAlign_chimericDetectionPEmerged.cpp
=====================================
@@ -1,6 +1,7 @@
#include "ReadAlign.h"
#include "BAMfunctions.h"
+
void ReadAlign::chimericDetectionPEmerged(ReadAlign &seRA) {
chimRecord=false;
@@ -20,59 +21,14 @@ void ReadAlign::chimericDetectionPEmerged(ReadAlign &seRA) {
return;
};
- //convert merged into PE
- for (uint ii=0; ii<2; ii++) {
- trChim[ii]=*trInit;
- trChim[ii].peOverlapSEtoPE(peOv.mateStart,seRA.trChim[ii]);
- };
-
- uint segLen[2][2]; //segment length [trChim][mate]
- uint segEx[2];//last exon of the mate0 [trChim]
- uint segLmin=-1LLU, i1=0,i2=0;
- for (uint ii=0; ii<2; ii++) {
- segLen[ii][0]=0;
- segLen[ii][1]=0;
- for (uint iex=0; iex<trChim[ii].nExons; iex++) {
- if (trChim[ii].exons[iex][EX_iFrag]==trChim[ii].exons[0][EX_iFrag]) {
- segLen[ii][0]+=trChim[ii].exons[iex][EX_L];
- segEx[ii]=iex;
- } else {
- segLen[ii][1]+=trChim[ii].exons[iex][EX_L];
- };
- };
- for (uint jj=0; jj<2; jj++) {
- if (segLen[ii][jj]<segLmin || (segLen[ii][jj]==segLmin && trChim[ii].exons[0][EX_G]>trChim[ii-1].exons[0][EX_G])) {
- segLmin=segLen[ii][jj];
- i1=ii;//trChim of the shortest segment length
- i2=jj;//mate of the shortest segment length
- };
- };
- };
-
- if (i2==1) {//eliminate mate1: simply cut the exons that belong to mate1
- trChim[i1].nExons=segEx[i1]+1;
- } else {//eliminate mate 0: shift mate1 exon to the beginning
- for (uint iex=0; iex<trChim[i1].nExons; iex++) {
- uint iex1=iex+segEx[i1]+1;
- for (uint ii=0; ii<EX_SIZE; ii++) {
- trChim[i1].exons[iex][ii]=trChim[i1].exons[iex1][ii];
- };
- trChim[i1].canonSJ[iex]=trChim[i1].canonSJ[iex1];
- trChim[i1].sjAnnot[iex]=trChim[i1].sjAnnot[iex1];
- trChim[i1].sjStr[iex]=trChim[i1].sjStr[iex1];
- trChim[i1].shiftSJ[iex][0]=trChim[i1].shiftSJ[iex1][0];
- trChim[i1].shiftSJ[iex][1]=trChim[i1].shiftSJ[iex1][1];
- };
- trChim[i1].nExons=trChim[i1].nExons-segEx[i1]-1;
- };
-
+ peOverlapChimericSEtoPE(&seRA.trChim[0], &seRA.trChim[1], &trChim[0], &trChim[1]);
chimericDetectionOldOutput();
} else if (trBest->maxScore <= (int) (readLength[0]+readLength[1]) - (int) P.pCh.nonchimScoreDropMin) {//require big enough drop in the best score
// new chimeric detection routine
- chimRecord=seRA.chimDet->chimericDetectionMult(seRA.nW, seRA.readLength, seRA.trBest->maxScore, true);
+ chimRecord=seRA.chimDet->chimericDetectionMult(seRA.nW, seRA.readLength, seRA.trBest->maxScore, this);
};
if ( chimRecord ) {
=====================================
source/ReadAlign_outputAlignments.cpp
=====================================
@@ -250,11 +250,10 @@ void ReadAlign::outputAlignments() {
chunkOutUnmappedReadsStream[im] << Qual0[im] <<"\n";
};
};
- if (P.pSolo.type>0) {//need to output 2nd (barcode) read
- chunkOutUnmappedReadsStream[1] << readNameMates[0] <<"\n";
- uint32 qualStart = readNameExtra[0].find(' ');
- chunkOutUnmappedReadsStream[1] << readNameExtra[0].substr(0,qualStart) <<"\n+\n";
- chunkOutUnmappedReadsStream[1] << readNameExtra[0].substr(qualStart+1) <<"\n";
+ if (P.readNmatesIn > P.readNmates) {//need to output 2nd (barcode) read, FASTQ only
+ chunkOutUnmappedReadsStream[P.readNmatesIn-1] << readNameMates[0] <<"\n";
+ chunkOutUnmappedReadsStream[P.readNmatesIn-1] << soloRead->readBar->bSeq <<"\n+\n";
+ chunkOutUnmappedReadsStream[P.readNmatesIn-1] << soloRead->readBar->bQual <<"\n";
};
};
};
=====================================
source/ReadAlign_outputTranscriptSAM.cpp
=====================================
@@ -344,6 +344,8 @@ uint ReadAlign::outputTranscriptSAM(Transcript const &trOut, uint nTrOut, uint i
case ATTR_vG:
case ATTR_vA:
case ATTR_vW:
+ case ATTR_GX:
+ case ATTR_GN:
break;
default:
ostringstream errOut;
=====================================
source/ReadAlign_peOverlapMergeMap.cpp
=====================================
@@ -141,7 +141,7 @@ void ReadAlign::peMergeMates() {
return;
};
-void Transcript::peOverlapSEtoPE(uint* mateStart, Transcript &t) {//convert alignment from merged-SE to PE
+void Transcript::peOverlapSEtoPE(uint* mateStart, const Transcript &t) {//convert alignment from merged-SE to PE
uint mLen[2];
mLen[0]=readLength[t.Str];
@@ -222,6 +222,11 @@ void Transcript::peOverlapSEtoPE(uint* mateStart, Transcript &t) {//convert alig
};
++nExons;
+ if (nExons>MAX_N_EXONS) {//cannot transform this alignment to PE, too many exons
+ maxScore=0;
+ nExons=0;
+ return;
+ };
};
canonSJ[nExons-1]=-3; //marks "junction" between mates
sjAnnot[nExons-1]=0;
@@ -268,28 +273,97 @@ void Transcript::peOverlapSEtoPE(uint* mateStart, Transcript &t) {//convert alig
void ReadAlign::peOverlapSEtoPE(ReadAlign &seRA) {//ReAdAlign: convert SE to PE and copy
- nW=seRA.nW;
- memcpy((void*) nWinTr, (void*) seRA.nWinTr, nW*sizeof(*nWinTr));
+ //nW=seRA.nW;
+ //memcpy((void*) nWinTr, (void*) seRA.nWinTr, nW*sizeof(*nWinTr));
uint trNtotal=0;
intScore bestScore=-10*Lread;
trBest=trArray;//just to initialize - to the 0th spot in the trArray
- for (uint iW=0; iW<nW; iW++) {//scan windows
- trAll[iW]=trArrayPointer+trNtotal;
- for (uint iTr=0; iTr<nWinTr[iW]; iTr++) {//scan transcripts
- ++trNtotal;
- *trAll[iW][iTr]=*trInit;
- trAll[iW][iTr]->peOverlapSEtoPE(peOv.mateStart, *seRA.trAll[iW][iTr]);
- trAll[iW][iTr]->alignScore(Read1,mapGen.G,P);
- if (trAll[iW][iTr]->maxScore > trAll[iW][0]->maxScore) {
- swap(trAll[iW][iTr],trAll[iW][0]);
+
+ uint64 iW1=0;
+ for (uint iW=0; iW<seRA.nW; iW++) {//scan windows
+ trAll[iW1]=trArrayPointer+trNtotal;
+ uint64 iTr1=0;
+ for (uint iTr=0; iTr<seRA.nWinTr[iW]; iTr++) {//scan transcripts
+ *trAll[iW1][iTr1]=*trInit;
+
+ trAll[iW1][iTr1]->peOverlapSEtoPE(peOv.mateStart, *seRA.trAll[iW][iTr]);
+ if (trAll[iW1][iTr1]->nExons==0)
+ continue; //conversion did not work
+
+ trAll[iW1][iTr1]->alignScore(Read1,mapGen.G,P);
+ if (trAll[iW1][iTr1]->maxScore > trAll[iW1][0]->maxScore) {
+ swap(trAll[iW][iTr1],trAll[iW][0]);
};
+
+ ++iTr1;
+ ++trNtotal;
};
- if (trAll[iW][0]->maxScore>bestScore) {
- trBest=trAll[iW][0];
- bestScore=trBest->maxScore;
+ if (iTr1>0) {//if conversion worked for at least one align
+ nWinTr[iW1]=iTr1;
+ if (trAll[iW1][0]->maxScore > bestScore) {
+ trBest=trAll[iW1][0];
+ bestScore=trBest->maxScore;
+ };
+ ++iW1;
};
};
+ nW=iW1;
return;
};
+
+void ReadAlign::peOverlapChimericSEtoPE(const Transcript *seTrIn1, const Transcript *seTrIn2, Transcript *peTrOut1, Transcript *peTrOut2) {
+
+ //convert merged into PE
+ Transcript tempTrChim[2];
+ tempTrChim[0]=*trInit;
+ tempTrChim[1]=*trInit;
+ tempTrChim[0].peOverlapSEtoPE(peOv.mateStart,*seTrIn1);
+ tempTrChim[1].peOverlapSEtoPE(peOv.mateStart,*seTrIn2);
+
+ uint segLen[2][2]; //segment length [tempTrChim][mate]
+ uint segEx[2];//last exon of the mate0 [tempTrChim]
+ uint segLmin=-1LLU, i1=0,i2=0;
+ for (uint ii=0; ii<2; ii++) {
+ segLen[ii][0]=0;
+ segLen[ii][1]=0;
+ for (uint iex=0; iex<tempTrChim[ii].nExons; iex++) {
+ if (tempTrChim[ii].exons[iex][EX_iFrag]==tempTrChim[ii].exons[0][EX_iFrag]) {
+ segLen[ii][0]+=tempTrChim[ii].exons[iex][EX_L];
+ segEx[ii]=iex;
+ } else {
+ segLen[ii][1]+=tempTrChim[ii].exons[iex][EX_L];
+ };
+ };
+ for (uint jj=0; jj<2; jj++) {
+ if (segLen[ii][jj]<segLmin || (segLen[ii][jj]==segLmin && tempTrChim[ii].exons[0][EX_G]>tempTrChim[ii-1].exons[0][EX_G])) {
+ segLmin=segLen[ii][jj];
+ i1=ii;//tempTrChim of the shortest segment length
+ i2=jj;//mate of the shortest segment length
+ };
+ };
+ };
+
+ if (i2==1) {//eliminate mate1: simply cut the exons that belong to mate1
+ tempTrChim[i1].nExons=segEx[i1]+1;
+ } else {//eliminate mate 0: shift mate1 exon to the beginning
+ for (uint iex=0; iex<tempTrChim[i1].nExons; iex++) {
+ uint iex1=iex+segEx[i1]+1;
+ for (uint ii=0; ii<EX_SIZE; ii++) {
+ tempTrChim[i1].exons[iex][ii]=tempTrChim[i1].exons[iex1][ii];
+ };
+ tempTrChim[i1].canonSJ[iex]=tempTrChim[i1].canonSJ[iex1];
+ tempTrChim[i1].sjAnnot[iex]=tempTrChim[i1].sjAnnot[iex1];
+ tempTrChim[i1].sjStr[iex]=tempTrChim[i1].sjStr[iex1];
+ tempTrChim[i1].shiftSJ[iex][0]=tempTrChim[i1].shiftSJ[iex1][0];
+ tempTrChim[i1].shiftSJ[iex][1]=tempTrChim[i1].shiftSJ[iex1][1];
+ };
+ tempTrChim[i1].nExons=tempTrChim[i1].nExons-segEx[i1]-1;
+ };
+
+ *peTrOut1=tempTrChim[0];
+ *peTrOut2=tempTrChim[1];
+ return;
+};
+
=====================================
source/SoloFeature_cellFiltering.cpp
=====================================
@@ -21,7 +21,7 @@ void SoloFeature::cellFiltering()
nUMImax = nUMIperCBsorted[min(nCB-1,pSolo.cellFilter.cr2maxCellInd)];//robust estimate of the max UMI
nUMImin = int(nUMImax/pSolo.cellFilter.cr2maxMinRatio+0.5);
} else if (pSolo.cellFilter.type[0]=="TopCells") {
- nUMImin = nUMIperCBsorted[max(nCB-1,pSolo.cellFilter.topCells)];
+ nUMImin = nUMIperCBsorted[min(nCB-1,pSolo.cellFilter.topCells)];
};
nUMImin=max(nUMImin,(uint32) 1);//cannot be zero
=====================================
source/Transcript.h
=====================================
@@ -68,7 +68,7 @@ public:
intScore alignScore(char **Read1, char *G, Parameters &P);
int variationAdjust(const Genome &mapGen, char *R);
string generateCigarP(); //generates CIGAR
- void peOverlapSEtoPE(uint* mSta, Transcript &t);
+ void peOverlapSEtoPE(uint* mSta, const Transcript &t);
void extractSpliceJunctions(vector<array<uint64,2>> &sjOut, bool &annotYes);
uint64 chrStartLengthExtended();
=====================================
source/Transcript_alignScore.cpp
=====================================
@@ -5,6 +5,10 @@ intScore Transcript::alignScore(char **Read1, char *G, Parameters &P) {//re-calc
maxScore=0;
nMM=0;
nMatch=0;
+
+ if (nExons==0)
+ return maxScore;
+
char* R=Read1[roStr==0 ? 0:2];
for (uint iex=0;iex<nExons;iex++) {
for (uint ii=0;ii<exons[iex][EX_L];ii++) {
=====================================
source/VERSION
=====================================
@@ -1 +1 @@
-#define STAR_VERSION "2.7.5c"
+#define STAR_VERSION "2.7.6a"
View it on GitLab: https://salsa.debian.org/med-team/rna-star/-/commit/ababb9e70de3862cdd3b1353d0c26b48130eaf70
--
View it on GitLab: https://salsa.debian.org/med-team/rna-star/-/commit/ababb9e70de3862cdd3b1353d0c26b48130eaf70
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20201002/dfb8108d/attachment-0001.html>
More information about the debian-med-commit
mailing list