EM.cpp

@@ -612,8 +612,31 @@ int main(int argc, char* argv[]) {
 	fin>>N0>>N1>>N2>>N_tot;
 	fin.close();
 
-	general_assert(N1 > 0, "There are no alignable reads!");
-
+	if (N1 == 0) {
+		printf("Warning: There are no alignable reads!\n");
+		theta.resize(M + 1, 0.0);
+		FILE *fo = NULL;
+		sprintf(thetaF, "%s.theta", statName);
+		fo = fopen(thetaF, "w");
+		fclose(fo);
+		sprintf(modelF, "%s.model", statName);
+		fo = fopen(modelF, "w");
+		fclose(fo);
+		eel.resize(M + 1, 0.0);
+		for (int i = 1; i <= M; ++i) eel[i] = transcripts.getTranscriptAt(i).getLength();
+		double *countv = new double[M + 1];
+		memset(countv, 0, sizeof(double) * (M + 1));
+		writeResultsEM(M, refName, imdName, transcripts, theta, eel, countv, appendNames);
+		if (genBamF) {
+			sprintf(outBamF, "%s.transcript.bam", outName);
+			char command[1005];
+			sprintf(command, "cp %s %s", inpSamF, outBamF);
+			printf("%s\n", command);
+			system(command);
+		}
+		delete[] countv;
+	}
+	else {
 		if ((READ_INT_TYPE)nThreads > N1) nThreads = N1;
 
 		//set model parameters
@@ -642,6 +665,7 @@ int main(int argc, char* argv[]) {
 		case 3 : EM<PairedEndReadQ, PairedEndHit, PairedEndQModel>(); break;
 		default : fprintf(stderr, "Unknown Read Type!\n"); exit(-1);
 		}		
+	}
 
 	time_t b = time(NULL);
 

WriteResults.h

@@ -86,7 +86,8 @@ void calcExpressionValues(int M, const std::vector<double>& theta, const std::ve
 	    frac[i] = theta[i];
 	    denom += frac[i];
 	  }
-	general_assert(denom >= EPSILON, "No alignable reads?!");
+	// general_assert(denom >= EPSILON, "No alignable reads?!");
+	if (denom < EPSILON) denom = 1.0;
 	for (int i = 1; i <= M; i++) frac[i] /= denom;
   
 	//calculate FPKM
@@ -98,6 +99,7 @@ void calcExpressionValues(int M, const std::vector<double>& theta, const std::ve
 	tpm.assign(M + 1, 0.0);
 	denom = 0.0;
 	for (int i = 1; i <= M; i++) denom += fpkm[i];
+	if (denom < EPSILON) denom = 1.0;
 	for (int i = 1; i <= M; i++) tpm[i] = fpkm[i] / denom * 1e6;  
 }
 
@@ -173,7 +175,7 @@ void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts
 		}
 		else {
 			for (int j = b; j < e; j++) {
-				isopct[j] = tpm[j] / gene_tpm[i];
+				isopct[j] = gene_tpm[i] > EPSILON ? tpm[j] / gene_tpm[i] : 0.0;
 				glens[i] += tlens[j] * isopct[j];
 				gene_eels[i] += eel[j] * isopct[j];
 			}
@@ -203,7 +205,7 @@ void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts
 		}
 		else {
 			for (int j = b; j < e; j++) {
-				ta_pct[j] = tpm[j] / trans_tpm[i];
+				ta_pct[j] = trans_tpm[i] > EPSILON ? tpm[j] / trans_tpm[i] : 0.0;
 				trans_lens[i] += tlens[j] * ta_pct[j];
 				trans_eels[i] += eel[j] * ta_pct[j];
 			}
@@ -214,7 +216,7 @@ void writeResultsEM(int M, const char* refName, const char* imdName, Transcripts
 	  for (int i = 0; i < m; i++) 
 	    if (gene_tpm[i] >= EPSILON) {
 	      int b = gt.spAt(i), e = gt.spAt(i + 1);
-	      for (int j = b; j < e; j++) gt_pct[j] = trans_tpm[j] / gene_tpm[i];
+	      for (int j = b; j < e; j++) gt_pct[j] = gene_tpm[i] > EPSILON ? trans_tpm[j] / gene_tpm[i] : 0.0;
 	    }
 	}
 

debian/changelog

+rsem (1.3.2+dfsg-1) unstable; urgency=medium
+
+  * New upstream version
+  * debhelper 12
+  * Standards-Version: 4.4.0
+  * Remove trailing whitespace in debian/changelog
+
+ -- Andreas Tille <tille@debian.org>  Fri, 12 Jul 2019 21:52:15 +0200
+
 rsem (1.3.1+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/compat

-11
+12

debian/control

@@ -4,7 +4,7 @@ Uploaders: Michael R. Crusoe <michael.crusoe@gmail.com>,
            Andreas Tille <tille@debian.org>
 Section: science
 Priority: optional
-Build-Depends: debhelper (>= 11~),
+Build-Depends: debhelper (>= 12~),
                zlib1g-dev,
                libncurses5-dev,
                libboost-dev,
@@ -15,7 +15,7 @@ Build-Depends: debhelper (>= 11~),
 # bowtie is a run-time dependency available on only a few systems, not actually
 # needed for building but it prevents the creation of uninstallable packages
                bowtie | bowtie2
-Standards-Version: 4.2.1
+Standards-Version: 4.4.0
 Vcs-Browser: https://salsa.debian.org/med-team/rsem
 Vcs-Git: https://salsa.debian.org/med-team/rsem.git
 Homepage: http://deweylab.biostat.wisc.edu/rsem/

rsem-calculate-expression

@@ -417,25 +417,11 @@ if (!$is_alignment) {
 	               " --outFileNamePrefix $imdName ";
 	               ##
 	
-#<<<<<<< HEAD
-	 #if ( $gzipped_read_file ) {
-	 #  $command .= ' --readFilesCommand zcat ';
-	 #} elsif ( $bzipped_read_file ) {
-	 #  $command .= ' --readFilesCommand bzcat ';
-	 #}
-	
-	 #if ( $read_type == 0 || $read_type == 1 ) {
-	 #  $command .= " --readFilesIn $mate1_list ";
-	 #} else {
-	 #  $command .= " --readFilesIn $mate1_list $mate2_list";
-	 #}
-#=======
 	    if ( $star_gzipped_read_file ) {
 		$command .= ' --readFilesCommand zcat ';
 	    } elsif ( $star_bzipped_read_file ) {
 		$command .= ' --readFilesCommand bzip2 -c ';
 	    }
-#>>>>>>> master
 	
 	    if ( $read_type == 0 || $read_type == 1 ) {
 		$command .= " --readFilesIn $mate1_list ";
@@ -552,6 +538,16 @@ if ($quiet) { $command .= " -q"; }
 
 &runCommand($command);
 
+my $inpCntF = "$statName.cnt";
+my $local_status = open(INPUT, $inpCntF);
+if ($local_status == 0) { print "Fail to open file $inpF!\n"; exit(-1); }
+my $line = <INPUT>;
+chomp($line);
+my @Ns = split(/ /, $line);
+close(INPUT);
+my $no_aligned = ($Ns[1] == 0);
+
+if (!$no_aligned) {
   $command = "rsem-build-read-index $gap"; 
   if ($read_type == 0) { $command .= " 0 $quiet $imdName\_alignable.fa"; }
   elsif ($read_type == 1) { $command .= " 1 $quiet $imdName\_alignable.fq"; }
@@ -559,6 +555,7 @@ elsif ($read_type == 2) { $command .= " 0 $quiet $imdName\_alignable_1.fa $imdNa
   elsif ($read_type == 3) { $command .= " 1 $quiet $imdName\_alignable_1.fq $imdName\_alignable_2.fq"; }
   else { print "Impossible! read_type is not in [1,2,3,4]!\n"; exit(-1); }
   &runCommand($command);  
+}
 
 my $doesOpen = open(OUTPUT, ">$imdName.mparams");
 if ($doesOpen == 0) { print "Cannot generate $imdName.mparams!\n"; exit(-1); }
@@ -628,6 +625,20 @@ if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
 
 if ($mTime) { $time_start = time(); }
 
+if ($no_aligned) { 
+  print "Since no aligned reads, further steps will not be performed!\n";
+  if (!$keep_intermediate_files) {
+      &runCommand("rm -rf $temp_dir", "Fail to delete the temporary folder!");
+  }
+  if ($mTime) { 
+      open(OUTPUT, ">$sampleName.time");
+      print OUTPUT "Aligning reads: $time_alignment s.\n";
+      print OUTPUT "Estimating expression levels: $time_rsem s.\n";
+      close(OUTPUT);
+  }
+  exit(0); 
+}
+
 if ($calcPME || $calcCI ) {
     $command = "rsem-run-gibbs $refName $imdName $statName $BURNIN $NCV $SAMPLEGAP";
     $command .= " -p $nThreads";