Set default vals for primaryOnly and canonOnly; Fixed bug leading to different lengths for SEQ and CIGAR when correcting junctions, adding extensive test cases for this.

Author: Dana Elizabeth Wyman

TranscriptClean.py

@@ -9,6 +9,7 @@
 
 # TC Classes
 from transcript import Transcript
+from transcript import check_seq_and_cigar_length
 from spliceJunction import *
 from intronBound import IntronBound
 from optparse import OptionParser
@@ -120,12 +121,13 @@ def getOptions():
     parser.add_option("--primaryOnly", dest ="primaryOnly", action='store_true',
                       help = "If this option is set, TranscriptClean will only \
                       output primary mappings of transcripts (ie it will filter \
-                      out unmapped and multimapped lines from the SAM input.")
+                      out unmapped and multimapped lines from the SAM input.",
+                      default = False)
     parser.add_option("--canonOnly", dest ="canonOnly", action='store_true',
                       help = ("If this option is set, TranscriptClean will "
                       "output only canonical transcripts and transcripts "
                       "containing annotated noncanonical junctions to the "
-                      "clean SAM file at the end of the run."))
+                      "clean SAM file at the end of the run."), default = False)
     parser.add_option("--tmpDir", dest ="tmp_path", 
                       help = ("If you would like the tmp files to be written "
                               "somewhere different than the final output, "
@@ -411,7 +413,7 @@ def correct_transcript(transcript_line, options, refs):
                               options.maxLenIndel, upd_logInfo)
             if ins_TE != "":
                 TE_entries += ins_TE
-            
+ 
             # Deletion correction
             del_TE = correctDeletions(upd_transcript, refs.genome, refs.deletions,
                              options.maxLenIndel, upd_logInfo)
@@ -1212,6 +1214,7 @@ def update_post_ncsj_correction(transcript, splice_jn_num, genome, sjAnnot):
     transcript.jM, transcript.jI = transcript.get_jM_jI_tags_from_sjs() 
     transcript.isCanonical = transcript.recheckCanonical()
     transcript.allJnsAnnotated = transcript.recheckJnsAnnotated()
+
     return
 
 def attempt_jn_correction(transcript, splice_jn_num, genome, ref_donors, 
@@ -1247,17 +1250,19 @@ def attempt_jn_correction(transcript, splice_jn_num, genome, ref_donors,
                                            transcript.POS, splice_jn_num, 
                                            donor, ref_donor.dist, genome, 
                                            transcript.SEQ, transcript.CIGAR)
-
+       
         # Attempt to fix the splice acceptor side
         acceptor = junction.get_splice_acceptor()
         transcript.SEQ, transcript.CIGAR = fix_one_side_of_junction(transcript.CHROM,
                                            transcript.POS, splice_jn_num,
                                            acceptor, ref_acceptor.dist, genome,
                                            transcript.SEQ, transcript.CIGAR)
+        
         # Now, perform updates:
         update_post_ncsj_correction(transcript, splice_jn_num, genome, sjAnnot)
 
-    except:
+    except Exception as e:
+        print(e)
         return False, "Other", combined_dist
 
     return True, "NA", combined_dist 
@@ -1384,6 +1389,9 @@ def fix_one_side_of_junction(chrom, transcript_start, jn_number, intronBound, d,
         newCIGAR = newCIGAR + exonCIGARs[i] + str(intronCIGARs[i]) + "N"
     newCIGAR = newCIGAR + exonCIGARs[-1]
 
+    if not check_seq_and_cigar_length(newSeq, newCIGAR):
+        raise RuntimeError("CIGAR string and sequence are not the same length")
+
     return newSeq, newCIGAR
 
 
@@ -1580,6 +1588,7 @@ def dryRun(sam, options, outfiles):
         write_to_transcript_log(logInfo, tL)
     return 
 
+
 if __name__ == '__main__':
     #pr = cProfile.Profile()
     #pr.enable()

testing_suite/input_files/drosophila_example/bad_correction.sam

@@ -0,0 +1,9 @@
+@SQ	SN:chr2L	LN:23513712
+@SQ	SN:chr2R	LN:25286936
+@SQ	SN:chr3L	LN:28110227
+@SQ	SN:chr3R	LN:32079331
+@SQ	SN:chr4	LN:1348131
+@SQ	SN:chrX	LN:23542271
+@SQ	SN:chrY	LN:3667352
+@PG	ID:minimap2	PN:minimap2	VN:2.15-r905	CL:minimap2 -t 16 -ax splice -uf --secondary=no -C5 /pub/dwyman/TALON-review/talon_on_other_data/refs/drosophila/dmel_r6.fa /pub/dwyman/TALON-review/talon_on_other_data/krizanovic_et_al/dataset7_real_subreads.fastq
+m160713_133433_42182_c101000162550000001823232709161620_s1_p0/121139/11291_13013	0	chr3R	14890311	53	67S116M269N202M10D232M11I151M62N371M61N92M74N284M13D101M6D30M198013N68M36S	*	0	0	AACAGGGTATCAACGCAGCGAGAGTACAGGGATGTGTAACACACCGCCCGAGAGCGACACAACAAAATCTTAAATATTTTTTAGTAATAGGCTATTGAACGAATTAATAAAACGATGGCGTTGCGAGTTCTTGGCGTGACGCAGAATATGCGCTCCTTCATGCGGGGCGTTGACATACAAGTCCGATGGCTGCCACCCTACAGAAGACTCTGCTGAACATACCCGCCACCCAGGTGACCAAGCTGGACAACGGTCTCCGCGTGGCTAGCGAGGATTCCGGTGCCTCTACCGCCACCGTGGGACTGTGGATCGATGCTGGGTCGCGCTCCGAGAACGAGAAGAACAACGGAGTGGCCCACTTCCTGGAGCACATGGCCTCTCACTAGCTTGATGCTTTAAATGCGTATAATATGCGTTTCTGCCACTTTTTAGGGTACCGCCAAGCGCTCGCAAACGGATTTGGAGCTGGAGGTTGAGAACCTGGGCGCCCACTTGAACGCCTACACTTCCAGGGAGCAGACCGTCTTCTACGCCAAGTGTTTGTCCAAGGATGTGCCCAAAGCCGTCGAGATCCTGGCCGACATCATCCAGAACTCCAAGCTGGGTGAGGCGCCACTATAACTGAGATCGCCCGCGAGCGTTCGGTGATTCTGCGCGAGATGCAGGAGGTGGAGAGCAACCTGCAGGAGGTGGTCTTCGATCACCTTCACGCCACCGCCTATCAGGGCACTCCTCTGGGCCAGACCATCCTGGGACCCACCAAGAACATTCAGTCCATTGGCAAGGCCGATCTGACCGACTACATCCAGACCCACTACAAGGCTTCGCGCATCGTTCTGGCCGCTGCTGGCGGCGTCAAGCACGATGATCTGGTCAAGCTTGCCTGCAGCAGCTTGGGTGGCCTGGAGGCCAGTGTGCTGCCCGCGGAGGTGACTCCCTGCCGCTTCACCGGCTCCGAGGTGCGTGTGCGCGACGATTCCCTGCCCCTGGCCCATGTCGCCATCGCCGTCGAAGGCTGCGGCTGGACCGACCAGGACAACATCCCCCTGATGGTGGCCAACACCTTGGTGGGTGCTTGGGATCGTTCCCAAGGCGGTGGCGCCAACAACGCCTCCAACCTGGCCCGTGCCAGCGCTGAGGACAACCTTTGCCACAGCTTCCAATCGTTCAACACATGCTACAAGGACACCGGACTCTGGGGCATTTACTTCGTGTGCGACCCTCTGCAGTGCGAGGATATGCTCTTCAACGTGCAGACCGAGTGGATGCGTCTGTGCACCATGGTTACCGAGGCTGAGGTCGAGCGCGCCAAGAACCTTCTGAAGACCAACATGCTGCTGCAGCTCGACGGCACCACACCCATCTGCGAGGACATTGGCCGCCAGATCCTGTGCTACAACCGCCGCATCCCGCTGCACGAGCTGGAGCAGCGCATCGATGCCGTGAGTGTGGGCAATGTGCGCGACGTCGCAATGAAGTACATCTACGATCGGTGCCCAGCCGTCGCTGCCGTGGGTCCCCCGACTACAACAGAATCCGCTCCTCCATGTACTGGTTGAGGGTTTAAGAGGTTCCGCCCCGCTGTTGTTGGACATTGTAGTTTAATTCAAAAGAAAATTTCAGAGGATCAGTAAACCGATAATACTAAAAATTCAAAGCTCCAACTCCGCGCAGCCCCAAAAAAAATATCGAATAAAAGCTAAAAAAAAAAAAAAAGAAAAAAAAAATAAAGTACCCTGCGTTGATACCAGCTT	*	ms:i:368	AS:i:172	nn:i:0	ts:A:+	tp:A:P	cm:i:15	s1:i:81	s2:i:0	de:f:0.1985	NM:i:1184	MD:Z:109G0A0T0C0A0A0A4A0C0A0A1T1G0G1T0G0C2C0C0C0T2A0G0A0A0G0A0C0T0C0T0G0C0T0G0A0A0C0A0T4G0C0C0A0C0C1A0G0G0T0G0A0C0C1A0G0C0T2A1A0A2G0T0C0T0C1G0C1T1G0C0T1G0C0G0A2A0T0T1C0G0G0T0G1C0T1T0A1C0G0C0C0A0C1G0T0G3C0T1T0G1A1C1A1G0C0T1G0G0T2C1C0T1C0G0A0G1A0C0G0A0G1A1A1C0A0A0C0G0G0A0G0T0G0G2C0A0C0T0T1C0T0G0G0A2A0C0A1G0G1C0T0^TCAAGGTTGG0C0T0C0A0C0T0A4G1T0G0C0T0T0T0A1A0T0G0C0G1A0T0A0A1A0T0G0C0G0T0T0T0C1G0C0C1C0T0T1T0T0A1G0G0T1C3C0A0A0G0C0G0C0T0C1C0A0A0A0C1G0A0T1T1G0A1C2G0A0G0G0T0T0G0A1A0A0C1T0G0G0G0C0G0C0C0C0A0C0T0T0G1A0C0G0C0C0T0A0C0A0C0T0T0C3G0G0A1C0A0G0A1C0G0T0C0T0T0C0T0A0C1C0C0A0A0G0T0G0T0T0T1T0C0C0A0A0G0G0A0T0G0T3C0A0A0A0G0C0C0G0T1G0A1A0T0C0C0T0G0G1C0G0A2T0C1T0C1A0G1A2C0C0A0A1C0T1G1T0G1G0G1T0G0A0G0A0T0C0G0C1C0G0C0G0A1C1T1C0G0G0T0G0A0T0T0C0T0G0C2G0A2T1C0A1G0A0G1T0G0G0A0G0A3A0C0C0T1C0A1G0A0G0G0T1G2T0T1G0A0T0C0A0C1T0T1A0C0G1C1C1G0C0C0T0A0T1A0G0G0G1A0C0T0C1T0C0T0G0G0G0C0C0A0G0A0C0C0A0T2T0G0G0G1C0C0C2C0A0A0G1A0C0A0T0T0C0A1T1C1T0T0G0G0C1A0G0G0C0C0G0A0T0C0T0G2C0G0A1T0A0C0A0T1C0A0G2C0C0A0C0T0A0C0A0A0G2T3C0G1A0T0C0G0T0T3G0C0C2T2T0G0G0C0G2G0T0C1A1C0A1G0A0T0G0A0T0C0T0G0G1C0A0A0G0C0T0T2C0T0G0C0A1C0A0G0C0T0T0G0G1T0G1C1T0G1A1G0C1A1T0G0T2T1C0C0C0G0C0G0G0A0G0G0T0G0A1T0C1C1G0C1G0C0T0T0C0A1C0G1C1C2A1G0T0G0C0G0T0G0T1C0G0C0G0A1G0A0T0T2C0T4C2G1C1A0T0G0T0C0G1C0A0T0C0G0C0C1T0C1A0A1G0C0T1C0G0G0C0T0G0G1C0C0G0A0C0C0A0G0G0A1A0A1A0T0C0C0C0C0C0T1A0T0G0G0T0G0G1C0A0A0C0A0C0C0T2G1G2T0G0C0T0T0G0G0G0A0T0C0G0T0T1C0C0A0A1G0C0G0G0T0G0G0C0G1C0A0A0C0A0A0C0G3C0C0A0A1C2G1C0C0G0T0G0C0C2C0G1T0G0A0G0G0A0C0A0A1C0T0T0T0G0C0C0A2G0C0T0T0C0C0A1T1G0T0T0C0A0A0C3T1C0T0A0C0A0A4A1C2A0C0T0C1G0G0G0G0C0A0T0T1A0C0T0T0C0G0T0G0T0G1G0A1C0C1C0T1C0A1T0G0C0G0A0G0G0A1A0T0G0C0T0C1T2A0C0G0T1C0A0G0A0C0C0G0A0G0T0G0G0A0T0G0C0G0T1T0G2C0A0C0C0A0T1G0T0T0A0C0C0G0A1G1T0G0A0G1T0C0G0A0G0C1C0G2A0A0G0A0A0C0C0T0T1T0G1A0G0A0C1A0A1A0T0G0C0T0G0C0T1C0A1C0T0C0G0A0C0G0G1A1C0A1A0C1C1T0C0T0G0C0G0A2A1A0T0T0G1C0C0G1C0A0G0A0T0C0C0T0G0T0G1T0A1A0A0C0C0G4A0T1C2C0T0G0C0A0C0G0A0G0C0T0G1A0G0C0A1C0G0C0A0T0C0G0A0T1C0C0G0T1A0G0T0G2G2A0A0T0G0T1C0G0C0G1C0G0T0C0G0C0A1T0G0A0A0G0T0A0C3T0A0C1A0T1G1T0G0C0C1A0G0C0C0G0T0C0G0C0T1C1G0T0G0G1T1C0^CGTTGAGAACCTG0C1C0G0A0C0T1C0A0A1A0G0A0A0T0C1G0C1C0C0T1C0A0T0G1A0C0T4G1G1G0T0T1A0A0G0A0G0G0T0T0C1G0C0C0C0C0G0C0T1T0T0G3G1C0A2G0T0A0G0T0T0T2T0T0C2A0A0G0A1A1T0T0T0^CACCCC0C0A0G0A0G0G1T1A1T0A2C0C0G0A0T2T1C1A3A0T2A2G0C0T1C0A0A1T0C1G1G0C1G0C0C0C0C1A1A0A0A4T0C0G0A0A5G0C0T8A6G1	jM:B:c,21,21,21,21,1	jI:B:i,14890427,14890695,14891291,14891352,14891724,14891784,14891877,14891950,14892385,15090397