Fix flair-quantify error when -output_bam is not specified (fixes #433)

Author: diekhans

src/flair/count_sam_transcripts.py

@@ -124,7 +124,7 @@ def check_splicesites(coveredpos, exonpos, tstart, tend):
 def get_matchvals(args, md):
 	matchvals = []
 	if args.stringent or args.check_splice:
-		mdblocks = re.findall('\d+|\D+', md)  # ['531', '^CCAGGTGAGCCGCCCGCG', '50', 'G', '2031']
+		mdblocks = re.findall(r'\d+|\D+', md)  # ['531', '^CCAGGTGAGCCGCCCGCG', '50', 'G', '2031']
 		for b in mdblocks:
 			if b[0] != '^':
 				if b.isnumeric():
@@ -234,17 +234,20 @@ def parsesam(args, transcripttoexons):
 				##for transcriptome alignment, always take rightmost side on transcript
 				if args.remove_internal_priming:
 					notinternalpriming = remove_internal_priming.removeinternalpriming(read.reference_name,
-																				   read.reference_start,
-																				   read.reference_end, False,
-																				   genome, None, transcripttoexons,
-																				   args.intprimingthreshold,
-																				   args.intprimingfracAs)
+													   read.reference_start,
+													   read.reference_end, False,
+													   genome, None, transcripttoexons,
+													   args.intprimingthreshold,
+													   args.intprimingfracAs)
 				else: notinternalpriming = True
 				if notinternalpriming:
 					pos = read.reference_start
-					alignscore = read.get_tag('AS')
+					try:
+						alignscore = read.get_tag('AS')
+						mdtag = read.get_tag('MD')
+					except KeyError as ex:
+						raise Exception(f"Missing AS or MD tag in alignment of '{read.query_name}' in '{args.sam.name}'") from ex
 					cigar = read.cigartuples
-					mdtag = read.get_tag('MD')
 					tlen = samfile.get_reference_length(transcript)
 					if lastread and readname != lastread:
 						assignedt = getbesttranscript(curr_transcripts, args, transcripttoexons)

src/flair/flair_quantify.py

@@ -1,5 +1,6 @@
 #! /usr/bin/env python3
 
+import os
 import sys
 import re
 import argparse
@@ -88,7 +89,7 @@ def quantify(isoform_sequences=''):
 			readFileRoot = tempfile.NamedTemporaryFile().name
 			if args.temp_dir != '':
 				if not os.path.isdir(args.temp_dir):
-					if subprocess.call(['mkdir', args.temp_dir]):
+					if subprocess.call(['mkdir', '-p', args.temp_dir]):
 						sys.stderr.write('Could not make temporary directory {}\n'.format(args.temp_dir))
 						return 1
 				readFileRoot = args.temp_dir + '/' + readFileRoot[readFileRoot.rfind('/')+1:]
@@ -101,9 +102,7 @@ def quantify(isoform_sequences=''):
 
 	for num, sample in enumerate(samData, 0):
 		sys.stderr.write('Step 1/3. Aligning sample %s_%s, %s/%s \n' % (sample[0], sample[2], num+1, len(samData)))
-		if args.output_bam: mm2_command = ['minimap2', '--MD', '-a', '-N', '4', '-t', str(args.t), args.i, sample[-2]]
-		else: mm2_command = ['minimap2', '-a', '-N', '4', '-t', str(args.t), args.i, sample[-2]]
-		# mm2_command = ['minimap2', '-a', '-N', '4', '-t', str(args.t), args.i, sample[-2]]
+		mm2_command = ['minimap2', '--MD', '-a', '-N', '4', '-t', str(args.t), args.i, sample[-2]]
 
 		# TODO: Replace this with proper try/except Exception as ex
 		try:
@@ -140,8 +139,7 @@ def quantify(isoform_sequences=''):
 		if args.generate_map or args.output_bam:
 			count_cmd += ['--generate_map', args.o+'.'+sample+'.'+group+'.isoform.read.map.txt']
 
-		if subprocess.call(count_cmd):
-			return 1
+		subprocess.check_call(count_cmd)
 		for line in open(samOut+'.counts.txt'):
 			line = line.rstrip().split('\t')
 			iso, numreads = line[0], line[1]
@@ -206,5 +204,6 @@ def quantify(isoform_sequences=''):
 	return args.o+'.counts.tsv'
 
 if __name__ == '__main__':
-	quantify()
-
+	# FIXME: need proper error handling
+	status = quantify()
+	os.exit(status)

test/Makefile

@@ -146,7 +146,7 @@ test-collapse-support : mkdirs ${GENOME_FA}
 ###
 # FLAIR QUANTIFY
 ###
-quantify-tests: test-quantify test-quantify-help
+quantify-tests: test-quantify test-quantify-nobam test-quantify-help
 
 test-quantify : mkdirs
 	rm -rf $O/$@.tmp
@@ -155,6 +155,12 @@ test-quantify : mkdirs
 	diff $E/$@.tpm.tsv $O/$@.tpm.tsv
 	diff $E/$@.A1.A.flair.aligned.readtotranscript.txt $O/$@.A1.A.flair.aligned.readtotranscript.txt
 
+# without map or bam output
+test-quantify-nobam : mkdirs
+	rm -rf $O/$@.tmp
+	flair quantify -r $(READS_MANIFEST) -i $(ISOFORMS_FA) --isoform_bed $(ISOFORMS_BED) --temp_dir $O/$@.tmp --tpm --sample_id_only -o $O/$@
+	diff $E/$@.tpm.tsv $O/$@.tpm.tsv
+
 test-quantify-help: mkdirs
 	flair quantify --help >&$O/$@.out
 	diff $E/$@.out $O/$@.out

test/expected/test-quantify-nobam.tpm.tsv

@@ -0,0 +1,11 @@
+ID	A1	A2	A3	B1	B2	B3
+ENST00000225792.10_ENSG00000108654.15	395061.7283950617	364161.84971098264	345679.012345679	263157.8947368421	315217.39130434784	376470.5882352941
+ENST00000256078.9_ENSG00000133703.12	0.0	5780.346820809248	6172.839506172839	47846.889952153106	27173.91304347826	5882.35294117647
+ENST00000311936.8_ENSG00000133703.12	111111.11111111111	92485.54913294797	98765.43209876542	114832.53588516745	103260.86956521739	94117.64705882352
+ENST00000348547.6_ENSG00000125991.19	111111.11111111111	104046.24277456648	123456.79012345678	172248.8038277512	135869.5652173913	117647.0588235294
+ENST00000357394.8_ENSG00000125991.19	37037.03703703704	23121.387283236993	55555.555555555555	33492.82296650718	43478.260869565216	52941.17647058823
+ENST00000451605.1_ENSG00000125991.19	74074.07407407407	167630.0578034682	98765.43209876542	138755.98086124402	146739.13043478262	94117.64705882352
+ENST00000557334.5_ENSG00000133703.12	18518.51851851852	0.0	24691.358024691355	43062.2009569378	86956.52173913043	23529.41176470588
+HISEQ:1287:HKCG7BCX3:1:1101:12283:18224-0_ENSG00000125991.19	6172.839506172839	0.0	6172.839506172839	0.0	0.0	5882.35294117647
+HISEQ:1287:HKCG7BCX3:1:1103:5649:96402-0_ENSG00000108654.15	135802.46913580247	86705.20231213873	141975.3086419753	76555.02392344497	48913.04347826087	135294.11764705883
+HISEQ:1287:HKCG7BCX3:1:1106:15192:70510-0_ENSG00000108654.15	111111.11111111111	156069.36416184972	98765.43209876542	110047.84688995215	92391.30434782608	94117.64705882352