galaxy-dist commit e2ba0e9c6852: Updating sam_indel_filter tool; adding NGS: Indel Analysis section; adding indel_analysis, indel_table, and indel_sam2interval tools
# HG changeset patch -- Bitbucket.org # Project galaxy-dist # URL http://bitbucket.org/galaxy/galaxy-dist/overview # User Kelly Vincent <kpvincent@bx.psu.edu> # Date 1279084192 14400 # Node ID e2ba0e9c6852f2acba9d7119c09b703d3bc954be # Parent f31c41763836603f42da474c3723d27143085cf8 Updating sam_indel_filter tool; adding NGS: Indel Analysis section; adding indel_analysis, indel_table, and indel_sam2interval tools --- /dev/null +++ b/test-data/indel_table_in1.interval @@ -0,0 +1,6 @@ +chrM 300 301 D - 6 +chrM 303 304 D - 19 +chrM 359 360 D - 1 +chrM 410 411 D - 1 +chrM 435 436 D - 1 +chrM 753 754 I A 1 --- a/test-data/sam_indel_filter_in2.sam +++ b/test-data/sam_indel_filter_in2.sam @@ -1,12 +1,12 @@ -081017-and-081020:1:6:774:1836 0 PHIX174 4973 37 26M1I9M * 0 0 GCTTAAAGCTACCAGTTATATGGCTGTTTGGTTTTT IIIIIIIIIIIIIIIIIIIIII@III/IE;%II;I= XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:33C1 -081017-and-081020:1:6:1193:793 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIII7IIIDIIIIIIII,=(>%II? XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 -081017-and-081020:1:8:753:970 0 PHIX174 2974 37 29M1I6M * 0 0 GTGGCGCCATGTCTAAATTGTTTGGAGGCGGGTCAA IIIIIIIIIIIIIIII4IIIIII3I&IIIII*%%&' XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:35 -081017-and-081020:1:17:361:871 0 PHIX174 4739 37 30M1I5M * 0 0 TGAGTATGGTACAGCTAATGGCCGTCTTTATTTTCC IIIIIIIIIIIIIIIIIIIGIIIII%II&'III%/# XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:28C6 -081017-and-081020:1:18:1164:1678 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIIIIIIIII;IIII1I0I)II.I- XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 -081017-and-081020:1:20:754:1256 16 PHIX174 4772 37 4M1I31M * 0 0 CCCTGGCGGTGCATTTTATGCGGACACTTCCTACAG &II(IIIII3IIIII7II,IIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:2A32 -081017-and-081020:1:24:1326:917 0 PHIX174 188 37 25M1D11M * 0 0 TTCGCCATCAACTAACGATTCTGTCAAACCTGACGC IIIIIIIIIIIIIIIIIIIIIIIIII2/&II>'IEI XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:25^A3A7 -081017-and-081020:1:25:1466:511 16 PHIX174 1179 37 21M1D15M * 0 0 CTATTGACTCTACTGTAGACATTTTACTTTTTATGT :I<=IIIIII5IIGI5IIIIIIIIIIIIIIIIIIII XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:21^T15 -081017-and-081020:1:27:1267:1275 0 PHIX174 3716 37 28M1D8M * 0 0 TCATCAGCAAACGCAGAATCAGCGGTATGCTCTTCT IIIIIIIIIIIIIIIII;IIIIIII87III%I(@I. XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:28^G8 -081017-and-081020:1:94:1649:147 16 PHIX174 2755 37 15M2D21M * 0 0 TATACCGTCAAGGACTGTGACTATTGACGTCCTTCC 4IIIIII@I7IIIIIIIIIIIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:15^TG21 -081017-and-081020:1:95:74:43 0 PHIX174 4038 37 29M1I6M * 0 0 ATCGAGGCTCTTAAACCTGCTATTGAGGCTTTTTGG IIIIIIIIIIIIIIIIICI;8I,I>IIIIII1I%5& XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 -081017-and-081020:1:95:58:307 16 PHIX174 3859 37 2M1I33M * 0 0 ACAAATGTCTGGAAAGACGGTAAAGCTGATGGTATT I&*IIIIIIII;IIIBIII>IICIIIIIIFII:III XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:5T29 +081017-and-081020:1:6:774:1836 0 PHIX174 4973 37 26M1I9M * 0 0 GCTTAAAGCTACCAGTTATATGGCTGTTTGGTTTTT IIIIIIIIIIIIIIIIIIIIII@III/IE;%II;I= XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:33C1 +081017-and-081020:1:6:1193:793 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIII7IIIDIIIIIIII,=(>%II? XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:8:753:970 0 PHIX174 2974 37 29M1I6M * 0 0 GTGGCGCCATGTCTAAATTGTTTGGAGGCGGGTCAA IIIIIIIIIIIIIIII4IIIIII3I&IIIII*%%&' XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:35 +081017-and-081020:1:17:361:871 0 PHIX174 4739 37 30M1I5M * 0 0 TGAGTATGGTACAGCTAATGGCCGTCTTTATTTTCC IIIIIIIIIIIIIIIIIIIGIIIII%II&'III%/# XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:28C6 +081017-and-081020:1:18:1164:1678 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIIIIIIIII;IIII1I0I)II.I- XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:20:754:1256 16 PHIX174 4772 37 4M1I31M * 0 0 CCCTGGCGGTGCATTTTATGCGGACACTTCCTACAG &II(IIIII3IIIII7II,IIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:2A32 +081017-and-081020:1:24:1326:917 0 PHIX174 188 37 25M1D11M * 0 0 TTCGCCATCAACTAACGATTCTGTCAAACCTGACGC IIIIIIIIIIIIIIIIIIIIIIIIII2/&II>'IEI XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:25^A3A7 +081017-and-081020:1:25:1466:511 16 PHIX174 1179 37 21M1D15M * 0 0 CTATTGACTCTACTGTAGACATTTTACTTTTTATGT :I<=IIIIII5IIGI5IIIIIIIIIIIIIIIIIIII XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:21^T15 +081017-and-081020:1:27:1267:1275 0 PHIX174 3716 37 28M1D8M * 0 0 TCATCAGCAAACGCAGAATCAGCGGTATGCTCTTCT IIIIIIIIIIIIIIIII;IIIIIII87III%I(@I. XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:28^G8 +081017-and-081020:1:94:1649:147 16 PHIX174 2755 37 15M2D21M * 0 0 TATACCGTCAAGGACTGTGACTATTGACGTCCTTCC 4IIIIII@I7IIIIIIIIIIIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:15^TG21 +081017-and-081020:1:95:74:43 0 PHIX174 4038 37 29M1I6M * 0 0 ATCGAGGCTCTTAAACCTGCTATTGAGGCTTTTTGG IIIIIIIIIIIIIIIIICI;8I,I>IIIIII1I%5& XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:95:58:307 16 PHIX174 3859 37 2M1I33M * 0 0 ACAAATGTCTGGAAAGACGGTAAAGCTGATGGTATT I&*IIIIIIII;IIIBIII>IICIIIIIIFII:III XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:5T29 --- /dev/null +++ b/test-data/indel_table_out1.interval @@ -0,0 +1,15 @@ +chrM 300 301 D - 14 6 0 8 +chrM 303 304 D - 41 19 0 22 +chrM 359 360 D - 1 1 0 0 +chrM 410 411 D - 4 1 1 2 +chrM 435 436 D - 2 1 0 1 +chrM 525 526 D - 1 0 0 1 +chrM 714 715 D - 1 0 1 0 +chrM 753 754 I A 2 1 0 1 +chrM 958 959 D - 1 0 0 1 +chrM 995 996 D - 3 0 0 3 +chrM 995 997 D - 1 0 1 0 +chrM 1168 1169 I A 1 0 1 0 +chrM 1168 1169 I C 1 0 0 1 +chrM 1274 1297 D - 1 0 1 0 +chrM 1296 1297 D - 2 0 1 1 --- /dev/null +++ b/test-data/indel_analysis_out4.interval @@ -0,0 +1,5 @@ +ref 13 14 T 2 13.33 14.29 +ref 20 21 GG 1 8.33 8.33 +ref 22 23 A 1 8.33 11.11 +ref 24 25 G 1 11.11 14.29 +ref 25 26 T 1 12.50 14.29 --- /dev/null +++ b/tools/indels/indel_analysis.py @@ -0,0 +1,240 @@ +#!/usr/bin/env python + +""" +Given an input sam file, provides analysis of the indels.. + +usage: %prog [options] [input3 sum3[ input4 sum4[ input5 sum5[...]]]] + -i, --input=i: The sam file to analyze + -t, --threshold=t: The deletion frequency threshold + -I, --out_ins=I: The interval output file showing insertions + -D, --out_del=D: The interval output file showing deletions +""" + +import re, sys +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse + + +def stop_err( msg ): + sys.stderr.write( '%s\n' % msg ) + sys.exit() + +def add_to_ref_pos( ref_pos, pos, bases ): + """ + Adds the bases and counts to the ref_pos dict + """ + for j, base in enumerate( bases ): + try: + ref_pos[ pos + j ][ base ] += 1 + except KeyError: + try: + ref_pos[ pos + j ][ base ] = 1 + except KeyError: + ref_pos[ pos + j ] = { base: 1 } + +def __main__(): + #Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + # prep output files + out_ins = open( options.out_ins, 'wb' ) + out_del = open( options.out_del, 'wb' ) + # pattern + pat = re.compile( '(^(?P<lmatch>\d+)M(?P<ins_del_width>\d+)(?P<ins_del>[ID])(?P<rmatch>\d+)M$)|(^(?P<match_width>\d+)M$)' ) + pat_multi = re.compile( '(\d+[MIDNSHP])(\d+[MIDNSHP])(\d+[MIDNSHP])+' ) + # for tracking occurences at each pos of ref + ref_pos = {} + indels = {} + num_reads = {} + multi_indel_lines = 0 + # go through all lines in input file + for i,line in enumerate( open( options.input, 'rb' ) ): + if line and not line.startswith( '#' ) and not line.startswith( '@' ) : + split_line = line.split( '\t' ) + chrom = split_line[2].strip() + pos = int( split_line[3].strip() ) + cigar = split_line[5].strip() + bases = split_line[9].strip() + # if not an indel or match, exit + if chrom == '*': + continue + # find matches like 3M2D7M or 7M3I10M + matches = '' + m = pat.search( cigar ) + # unprocessable CIGAR + if not m: + m = pat_multi.search( cigar ) + # skip this line if no match + if not m: + continue + # account for multiple indels or operations we don't process + else: + multi_indel_lines += 1 + # get matching parts for the indel or full match if matching + else: + if not ref_pos.has_key( chrom ): + ref_pos[ chrom ] = {} + indels[ chrom ] = { 'D': {}, 'I': {} } + if not num_reads.has_key( chrom ): + num_reads[ chrom ] = {} + parts = m.groupdict() + if parts[ 'match_width' ] or ( parts[ 'lmatch' ] and parts[ 'ins_del_width' ] and parts[ 'rmatch' ] ): + match = parts + # see if matches meet filter requirements + if match: + # match/mismatch + if parts[ 'match_width' ]: + add_to_ref_pos( ref_pos[ chrom ], pos, bases ) + for i, base in enumerate( bases ): + try: + num_reads[ chrom ][ i + pos ] += 1 + except KeyError: + num_reads[ chrom ][ i + pos ] = 1 + # indel + else: + # pieces of CIGAR string + left = int( match[ 'lmatch' ] ) + middle = int( match[ 'ins_del_width' ] ) + right = int( match[ 'rmatch' ] ) + left_bases = bases[ : left ] + if match[ 'ins_del' ] == 'I': + middle_bases = bases[ left : left + middle ] + else: + middle_bases = '' + right_bases = bases[ -right : ] + start = pos + left + # add data to ref_pos dict for match/mismatch bases on left and on right + add_to_ref_pos( ref_pos[ chrom ], pos, left_bases ) + for i, base in enumerate( left_bases ): + try: + num_reads[ chrom ][ i + pos ] += 1 + except KeyError: + num_reads[ chrom ][ i + pos ] = 1 + if match[ 'ins_del' ] == 'I': + add_to_ref_pos( ref_pos[ chrom ], start, right_bases ) + indel_pos = start + else: + add_to_ref_pos( ref_pos[ chrom ], start + middle, right_bases ) + indel_pos = start + middle + for i, base in enumerate( right_bases ): + try: + num_reads[ chrom ][ i + indel_pos ] += 1 + except KeyError: + num_reads[ chrom ][ i + indel_pos ] = 1 + # for insertions, count instances of particular inserted bases + if match[ 'ins_del' ] == 'I': + if indels[ chrom ][ 'I' ].has_key( start ): + if indels[ chrom ][ 'I' ][ start ].has_key( middle_bases ): + indels[ chrom ][ 'I' ][ start ][ middle_bases ] += 1 + else: + indels[ chrom ][ 'I' ][ start ][ middle_bases ] = 1 + else: + indels[ chrom ][ 'I' ][ start ] = { middle_bases: 1 } + # for deletions, count number of deletions bases + else: + if indels[ chrom ][ 'D' ].has_key( start ): + if indels[ chrom ][ 'D' ][ start ].has_key( middle ): + indels[ chrom ][ 'D' ][ start ][ middle ] += 1 + else: + indels[ chrom ][ 'D' ][ start ][ middle ] = 1 + else: + indels[ chrom ][ 'D' ][ start ] = { middle: 1 } + # compute deletion frequencies and insertion frequencies for checking against threshold + freqs = {} + ins_freqs = {} + chroms = ref_pos.keys() + chroms.sort() + for chrom in chroms: + freqs[ chrom ] = {} + ins_freqs[ chrom ] = {} + poses = num_reads[ chrom ].keys() + poses.sort() + for pos in poses: + # all reads touching this particular position + freqs[ chrom ][ pos ] = {} + sum_counts = 0.0 + sum_counts_end = 0.0 + # get basic counts (match/mismatch) + if num_reads[ chrom ].has_key( pos ): + sum_counts += float( num_reads[ chrom ][ pos ] ) + try: + sum_counts_end += float( num_reads[ chrom ][ pos + 1 ] ) + except KeyError: + pass + # add deletions also touching this position + try: + sum_counts += float( sum( indels[ chrom ][ 'D' ][ pos ].values() ) ) + try: + sum_counts_end += float( sum( indels[ chrom ][ 'D' ][ pos + 1 ].values() ) ) + except KeyError: + pass + for d in indels[ chrom ][ 'D' ][ pos ].keys(): + freqs[ chrom ][ pos ][ '-' * d ] = indels[ chrom ][ 'D' ][ pos ][ d ] / sum_counts + except KeyError: + pass + # calculate actual frequencies + # deletions + freqs[ chrom ][ pos ][ 'total' ] = sum_counts + for base in ref_pos[ chrom ][ pos ].keys(): + try: + prop = float( ref_pos[ chrom ][ pos ][ base ] ) / sum_counts + freqs[ chrom ][ pos ][ base ] = prop + except ZeroDivisionError: + freqs[ chrom ][ pos ][ base ] = 0.0 + try: + for d in indels[ chrom ][ 'D' ][ pos ].keys(): + freqs[ chrom ][ pos ][ '-' * d ] = indels[ chrom ][ 'D' ][ pos ][ d ] / sum_counts + except KeyError: + pass + # insertions + try: + for bases in indels[ chrom ][ 'I' ][ pos ].keys(): + prop_start = indels[ chrom ][ 'I' ][ pos ][ bases ] / ( indels[ chrom ][ 'I' ][ pos ][ bases ] + sum_counts ) + try: + prop_end = indels[ chrom ][ 'I' ][ pos ][ bases ] / sum_counts_end + except ZeroDivisionError: + prop_end = 0.0 + try: + ins_freqs[ chrom ][ pos ][ bases ] = [ prop_start, prop_end ] + except KeyError: + ins_freqs[ chrom ][ pos ] = { bases: [ prop_start, prop_end ] } + except KeyError, e: + pass + # output to files if meet threshold requirement + threshold = float( options.threshold ) + #out_del.write( '#Chrom\tStart\tEnd\t#Del\t#Reads\t%TotReads\n' ) + #out_ins.write( '#Chrom\tStart\tEnd\tInsBases\t#Reads\t%TotReadsAtStart\t%ReadsAtEnd\n' ) + for chrom in ref_pos.keys(): + # deletions file + poses = indels[ chrom ][ 'D' ].keys() + poses.sort() + for pos in poses: + start = pos + dels = indels[ chrom ][ 'D' ][ start ].keys() + dels.sort() + for d in dels: + end = start + d + prop = freqs[ chrom ][ start ][ '-' * d ] + if prop > threshold : + out_del.write( '%s\t%s\t%s\t%s\t%s\t%.2f\n' % ( chrom, start, end, d, indels[ chrom ][ 'D' ][ pos ][ d ], 100.0 * prop ) ) + # insertions file + poses = indels[ chrom ][ 'I' ].keys() + poses.sort() + for pos in poses: + start = pos + end = pos + 1 + ins_bases = indels[ chrom ][ 'I' ][ start ].keys() + ins_bases.sort() + for bases in ins_bases: + prop_start = ins_freqs[ chrom ][ start ][ bases ][0] + prop_end = ins_freqs[ chrom ][ start ][ bases ][1] + if prop_start > threshold or prop_end > threshold: + out_ins.write( '%s\t%s\t%s\t%s\t%s\t%.2f\t%.2f\n' % ( chrom, start, end, bases, indels[ chrom ][ 'I' ][ start ][ bases ], 100.0 * prop_start, 100.0 * prop_end ) ) + # close out files + out_del.close() + out_ins.close() + # if skipped lines because of more than one indel, output message + if multi_indel_lines > 0: + sys.stdout.write( '%s alignments were skipped because they contained more than one indel or had unhandled operations (N/S/H/P).' % multi_indel_lines ) + +if __name__=="__main__": __main__() --- /dev/null +++ b/test-data/indel_sam2interval_out1.interval @@ -0,0 +1,6 @@ +ref 133 134 I C 1 +ref 256 258 D - 1 +ref 48819784 48819785 I A 1 +ref 87824726 87824727 I G 1 +ref 188841437 188841438 I A 1 +ref 190341171 190341172 D - 1 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -3,7 +3,7 @@ <section name="Get Data" id="getext"><tool file="data_source/upload.xml"/><tool file="data_source/ucsc_tablebrowser.xml" /> - <tool file="data_source/ucsc_tablebrowser_test.xml" /> + <tool file="data_source/ucsc_tablebrowser_test.xml" /><tool file="data_source/ucsc_tablebrowser_archaea.xml" /><tool file="data_source/bx_browser.xml" /><tool file="data_source/microbial_import.xml" /> @@ -23,7 +23,7 @@ <tool file="data_source/encode_db.xml" /><tool file="data_source/epigraph_import.xml" /><tool file="data_source/epigraph_import_test.xml" /> - <tool file="data_source/hbvar.xml" /> + <tool file="data_source/hbvar.xml" /><tool file="validation/fix_errors.xml" /></section><section name="Send Data" id="send"> @@ -112,7 +112,7 @@ <tool file="extract/phastOdds/phastOdds_tool.xml" /></section><section name="Operate on Genomic Intervals" id="bxops"> - <tool file="new_operations/intersect.xml" /> + <tool file="new_operations/intersect.xml" /><tool file="new_operations/subtract.xml" /><tool file="new_operations/merge.xml" /><tool file="new_operations/concat.xml" /> @@ -213,7 +213,7 @@ <tool file="fastx_toolkit/fastx_quality_statistics.xml" /><tool file="fastx_toolkit/fastq_quality_boxplot.xml" /><tool file="fastx_toolkit/fastx_nucleotides_distribution.xml" /> - <tool file="metag_tools/split_paired_reads.xml" + <tool file="metag_tools/split_paired_reads.xml" <tool file="fastx_toolkit/fastq_to_fasta.xml" /> --><label text="Roche-454 data" id="454" /> @@ -255,6 +255,12 @@ <tool file="metag_tools/megablast_xml_parser.xml" /><tool file="sr_mapping/PerM.xml" /></section> + <section name="NGS: Indel Analysis" id="indel_analysis"> + <tool file="indels/sam_indel_filter.xml" /> + <tool file="indels/indel_sam2interval.xml" /> + <tool file="indels/indel_table.xml" />" + <tool file="indels/indel_analysis.xml" /> + </section><section name="NGS: Expression Analysis" id="ngs-rna-tools"><label text="RNA-seq" id="rna_seq" /><tool file="ngs_rna/tophat_wrapper.xml" /> @@ -267,7 +273,6 @@ </section><section name="NGS: SAM Tools" id="samtools"><tool file="samtools/sam_bitwise_flag_filter.xml" /> - <tool file="samtools/sam_indel_filter.xml" /><tool file="samtools/sam2interval.xml" /><tool file="samtools/sam_to_bam.xml" /><tool file="samtools/sam_merge.xml" /> @@ -309,12 +314,12 @@ <tool file="rgenetics/rgManQQ.xml"/></section><!-- - TODO: uncomment the following EMBOSS section whenever - moving to test, but comment it in .sample to eliminate - it from buildbot functional tests since these tools + TODO: uncomment the following EMBOSS section whenever + moving to test, but comment it in .sample to eliminate + it from buildbot functional tests since these tools rarely change. --> -<!-- +<!-- <section name="EMBOSS" id="EMBOSSLite"><tool file="emboss_5/emboss_antigenic.xml" /><tool file="emboss_5/emboss_backtranseq.xml" /> @@ -329,13 +334,13 @@ <tool file="emboss_5/emboss_chips.xml" /><tool file="emboss_5/emboss_cirdna.xml" /><tool file="emboss_5/emboss_codcmp.xml" /> - <tool file="emboss_5/emboss_coderet.xml" /> + <tool file="emboss_5/emboss_coderet.xml" /><tool file="emboss_5/emboss_compseq.xml" /> - <tool file="emboss_5/emboss_cpgplot.xml" /> + <tool file="emboss_5/emboss_cpgplot.xml" /><tool file="emboss_5/emboss_cpgreport.xml" /><tool file="emboss_5/emboss_cusp.xml" /><tool file="emboss_5/emboss_cutseq.xml" /> - <tool file="emboss_5/emboss_dan.xml" /> + <tool file="emboss_5/emboss_dan.xml" /><tool file="emboss_5/emboss_degapseq.xml" /><tool file="emboss_5/emboss_descseq.xml" /><tool file="emboss_5/emboss_diffseq.xml" /> @@ -351,7 +356,7 @@ <tool file="emboss_5/emboss_etandem.xml" /><tool file="emboss_5/emboss_extractfeat.xml" /><tool file="emboss_5/emboss_extractseq.xml" /> - <tool file="emboss_5/emboss_freak.xml" /> + <tool file="emboss_5/emboss_freak.xml" /><tool file="emboss_5/emboss_fuzznuc.xml" /><tool file="emboss_5/emboss_fuzzpro.xml" /><tool file="emboss_5/emboss_fuzztran.xml" /> @@ -372,7 +377,7 @@ <tool file="emboss_5/emboss_merger.xml" /><tool file="emboss_5/emboss_msbar.xml" /><tool file="emboss_5/emboss_needle.xml" /> - <tool file="emboss_5/emboss_newcpgreport.xml" /> + <tool file="emboss_5/emboss_newcpgreport.xml" /><tool file="emboss_5/emboss_newcpgseek.xml" /><tool file="emboss_5/emboss_newseq.xml" /><tool file="emboss_5/emboss_noreturn.xml" /> @@ -400,7 +405,7 @@ <tool file="emboss_5/emboss_revseq.xml" /><tool file="emboss_5/emboss_seqmatchall.xml" /><tool file="emboss_5/emboss_seqret.xml" /> - <tool file="emboss_5/emboss_showfeat.xml" /> + <tool file="emboss_5/emboss_showfeat.xml" /><tool file="emboss_5/emboss_shuffleseq.xml" /><tool file="emboss_5/emboss_sigcleave.xml" /><tool file="emboss_5/emboss_sirna.xml" /> @@ -422,7 +427,7 @@ <tool file="emboss_5/emboss_water.xml" /><tool file="emboss_5/emboss_wobble.xml" /><tool file="emboss_5/emboss_wordcount.xml" /> - <tool file="emboss_5/emboss_wordmatch.xml" /> - </section> + <tool file="emboss_5/emboss_wordmatch.xml" /> + </section> --></toolbox> --- a/tools/samtools/sam_indel_filter.xml +++ /dev/null @@ -1,77 +0,0 @@ -<tool id="sam_indel_filter" name="Filter SAM" version="1.0.0"> - <description>for indels</description> - <command interpreter="python"> - sam_indel_filter.py - --input=$input1 - --quality_threshold=$quality_threshold - --adjacent_bases=$adjacent_bases - --output=$out_file1 - </command> - <inputs> - <param format="sam" name="input1" type="data" label="Select dataset to filter" /> - <param name="quality_threshold" type="integer" value="40" label="Quality threshold for adjacent bases" help="Takes Phred value assuming Sanger scale; usually between 0 and 40, but up to 93" /> - <param name="adjacent_bases" type="integer" value="1" label="The number of adjacent bases to match on either side of the indel" help="If one side is shorter than this width, it will still match if the long-enough side matches" /> - </inputs> - <outputs> - <data format="sam" name="out_file1" /> - </outputs> - <tests> - <test> - <param name="input1" value="sam_indel_filter_in1.sam" ftype="sam"/> - <param name="quality_threshold" value="14"/> - <param name="adjacent_bases" value="2"/> - <output name="out_file1" file="sam_indel_filter_out1.sam" ftype="sam"/> - </test> - <test> - <param name="input1" value="sam_indel_filter_in1.sam" ftype="sam"/> - <param name="quality_threshold" value="29"/> - <param name="adjacent_bases" value="5"/> - <output name="out_file1" file="sam_indel_filter_out2.sam" ftype="sam"/> - </test> - <test> - <param name="input1" value="sam_indel_filter_in2.sam" ftype="sam"/> - <param name="quality_threshold" value="7"/> - <param name="adjacent_bases" value="1"/> - <output name="out_file1" file="sam_indel_filter_out3.sam" ftype="sam"/> - </test> - </tests> - <help> - -**What it does** - -Allows extracting indels from SAM. Currently it can handle SAM with alignments with only one insertion or one deletion, and will skip that alignment if it encounters one with more than one indel. It matches CIGAR strings (column 6 in the SAM file) like 5M3I5M or 4M2D10M, so there must be a match or mismatch of sufficient length on either side of the indel. - ------ - -**Example** - -Suppose you have the following:: - - r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - r770 181 ref 116 0 24M = 72131356 0 TTGGTGCGCGCGGTTGAGGGTTGG $$(#%%#$%#%####$%%##$### - r1945 177 ref 41710908 0 23M 190342418 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 - r3671 117 ref 190342418 0 24M = 190342418 0 CTGGCGTTCTCGGCGTGGATGGGT #####$$##$#%#%%###%$#$## - r3671 153 ref 190342418 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - r3824 117 ref 80324999 0 24M = 80324999 0 TCCAGTCGCGTTGTTAGGTTCGGA #$#$$$#####%##%%###**#+/ - r3824 153 ref 80324999 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - r4795 81 ref 26739130 0 23M 57401793 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 - r4795 161 ref 57401793 37 23M 26739130 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 - r4800 16 ref 241 255 15M1D8M * 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII - r5377 170 ref 59090793 37 23M 26739130 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 - r5612 151 ref 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - - -To select only alignments with indels, you need to determine the minimum quality you want the adjacent bases to have, as well as the number of adjacent bases to check. If you set the quality threshold to 47 and the number of bases to check to 2, you will get the following output:: - - r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - r4800 16 ref 241 255 15M1D8M * 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII - r5612 151 ref 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 - - -For more information on SAM, please consult the `SAM format description`__. - -.. __: http://www.ncbi.nlm.nih.gov/pubmed/19505943 - - - </help> -</tool> --- a/tools/samtools/sam_indel_filter.py +++ /dev/null @@ -1,107 +0,0 @@ -#!/usr/bin/env python - -""" -Allows user to filter out non-indels from SAM. - -usage: %prog [options] - -i, --input=i: Input SAM file to be filtered - -q, --quality_threshold=q: Minimum quality value for adjacent bases - -a, --adjacent_bases=a: Number of adjacent bases on each size to check qualities - -o, --output=o: Filtered output SAM file -""" - -import re, sys -from galaxy import eggs -import pkg_resources; pkg_resources.require( "bx-python" ) -from bx.cookbook import doc_optparse - - -def stop_err( msg ): - sys.stderr.write( '%s\n' % msg ) - sys.exit() - -def __main__(): - #Parse Command Line - options, args = doc_optparse.parse( __doc__ ) - # prep output file - output = open( options.output, 'wb' ) - # patterns - pat_indel = re.compile( '(?P<before_match>(\d+[MNSHP])*)(?P<lmatch>\d+)M(?P<ins_del_width>\d+)(?P<ins_del>[ID])(?P<rmatch>\d+)M(?P<after_match>(\d+[MNSHP])*)' ) - pat_matches = re.compile( '(\d+[MIDNSHP])+' ) - try: - qual_thresh = int( options.quality_threshold ) + 33 - if qual_thresh < 33 or qual_thresh > 126: - raise ValueError - except ValueError: - stop_err( 'Your quality threshold should be an integer between 0 and 93, inclusive.' ) - try: - adj_bases = int( options.adjacent_bases ) - if adj_bases < 1: - raise ValueError - except ValueError: - stop_err( 'The number of adjacent bases should be an integer greater than 1.' ) - # record lines skipped because of more than one indel - multi_indel_lines = 0 - # go through all lines in input file - for i,line in enumerate(open( options.input, 'rb' )): - if i > 1000: - break - if line and not line.startswith( '#' ) and not line.startswith( '@' ) : - split_line = line.split( '\t' ) - cigar = split_line[5] - # find all possible matches, like 3M2D7M and 7M3I10M in 3M2D7M3I10M - cigar_copy = cigar[:] - matches = [] - while len( cigar_copy ) >= 6: # nMnInM or nMnDnM - m = pat_indel.search( cigar_copy ) - if not m: - break - else: - parts = m.groupdict() - pre_left = 0 - if m.start() > 0: - pre_left_groups = pat_matches.search( cigar_copy[ : m.start() ] ) - if pre_left_groups: - for pl in pre_left_groups.groups(): - if pl.endswith( 'M' ) or pl.endswith( 'S' ) or pl.endswith( 'P' ): - pre_left += pl[:-1] - parts[ 'pre_left' ] = pre_left - matches.append( parts ) - cigar_copy = cigar_copy[ len( parts[ 'lmatch' ] ) + 1 : ] - # see if matches meet filter requirements - if len( matches ) > 1: - multi_indel_lines += 1 - elif len( matches ) == 1: - pre_left = int( matches[0][ 'pre_left' ] ) - left = int( matches[0][ 'lmatch' ] ) - right = int( matches[0][ 'rmatch' ] ) - if matches[0][ 'ins_del' ] == 'D': - middle = int( matches[0][ 'ins_del_width' ] ) - else: - middle = 0 - # if there are enough adjacent bases to check, then do so - if left >= adj_bases and right >= adj_bases: - quals = split_line[10] - left_quals = quals[ pre_left : pre_left + left ][ -adj_bases : ] - middle_quals = quals[ pre_left + left : pre_left + left + middle ] - right_quals = quals[ pre_left + left + middle : pre_left + left + middle + right ][ : adj_bases ] - qual_thresh_met = True - for l in left_quals: - if ord( l ) < qual_thresh: - qual_thresh_met = False - break - if qual_thresh_met: - for r in right_quals: - if ord( r ) < qual_thresh: - qual_thresh_met = False - break - # if filter reqs met, output line - if qual_thresh_met: - output.write( line ) - # close out file - output.close() - # if skipped lines because of more than one indel, output message - if multi_indel_lines > 0: - sys.stdout.write( '%s alignments were skipped because they contained more than one indel.' % multi_indel_lines ) - -if __name__=="__main__": __main__() --- /dev/null +++ b/test-data/indel_analysis_in1.sam @@ -0,0 +1,22 @@ +r770 89 ref 6 37 7M1I5M = 0 0 TCGATCTTCATAG !0//110AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r1124 113 ref 4 0 23M = 0 0 CATCGTTCTGTTAGATCTACGTA PQRVUMNXYRPUXYXWXSOSZ]M XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r1231 69 * 0 0 * * 0 0 AGACCGGGCGGGGTGGCGTTCGGT %##+'#######%###$#$##$(# +r1563 133 * 0 0 * * 0 0 GTTCGTGGCCGGTGGGTGTTTGGG ###$$#$#$&#####$'$#$###$ +r1789 177 ref 6 0 17M = 0 0 TCGATCGCTTAGTTCTC SQQWZY]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r3671 153 ref 10 37 6M1I6M = 0 0 TCTCTTTAGGTCT /<<!"0/////// XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r3824 153 ref 5 37 8M1I7M = 0 0 ATCGATGTTCTTAGAT 4;6//11!"100110/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r4800 16 ref 7 255 5M2D6M = 0 0 CGATCTTTGAT IIIIIIIIIIC XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r5612 151 ref 5 37 8M1D9M = 0 0 ATCTATCTTTTGATCTC /<<!"0/4/*/7//B0/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r5929 151 ref 11 37 3M1I10M = 0 0 CTCCTTAGCTCTCC /<<!"0/4//7//0 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r6743 69 * 0 0 * * 0 0 TGCCGTGTCTTGCTAACGCCGATT #'#$$#$###%%##$$$$###### +r9145 115 ref 11 0 19M = 0 -1 CTCTTAGCTCTCCGAATTAG 7753:<5#"4!&=9518A> XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:137 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 +r11770 89 ref 10 37 10M2I8M = 0 0 TCTCTTAGATGGCTCCGTAT 00/02!!0/120210AA4/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r13671 153 ref 1 37 12M1I12M = 0 0 TCGCATCGATCTCCTTAGATCTCCG /<""<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r13762 133 * 0 0 * * 0 0 TGGGTGGATGTGTTGTCGTTCATG #$#$###$#$#######$#$#### +r13824 153 ref 13 37 9M1I7M = 0 0 CATAGATCTACCGGATT 4;6//11!"11100110 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r24800 16 ref 3 255 15M2D9M = 0 0 GCATCTATCTGATAGCTCCGAATT IIIIIIIII45"CCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r25612 151 ref 1 37 9M1D5M = 0 0 TCGCATCGACTCTT 0/4/*/7//00/1C XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r25786 151 ref 21 37 4M1I7M = 0 0 TGCGTTATTGGG <!"0/70/BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r27899 69 * 0 0 * * 0 0 CTGCGTGTTGGTGTCTACTGGGGT #%#'##$#$##&%#%$$$%#%#'# +r29192 133 * 0 0 * * 0 0 GTGCGTCGGGGAGGGTGCTGTCGG ######%#$%#$$###($###&&% +r29962 16 ref 20 37 4M1I7M = 0 0 CTCCGGTATGAGG <!"0/70/7BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 --- /dev/null +++ b/test-data/indel_table_in2.interval @@ -0,0 +1,6 @@ +chrM 410 411 D - 1 +chrM 714 715 D - 1 +chrM 995 997 D - 1 +chrM 1168 1169 I A 1 +chrM 1274 1297 D - 1 +chrM 1296 1297 D - 1 --- a/test-data/sam_indel_filter_out1.sam +++ b/test-data/sam_indel_filter_out1.sam @@ -1,3 +1,3 @@ -1378_28_770 89 chr11.nib:1-134452384 72131356 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_28_770 89 chr11.nib:1-134452384 72131356 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 1378_69_800 16 chr11.nib:1-125234658 241 255 15M1D8M * 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII -1378_72_1612 151 chrY.nib:1-124295114 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_72_1612 151 chrY.nib:1-124295114 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 --- /dev/null +++ b/tools/indels/indel_table.xml @@ -0,0 +1,122 @@ +<tool id="indel_table" name="Indel Analysis Table" version="1.0.0"> + <description>for combining indel interval data</description> + <command interpreter="python"> + indel_table.py + --input1=$input1 + --sum1=$sum1 + --input2=$input2 + --sum2=$sum2 + --output=$output1 + #for $i in $inputs + ${i.input} + ${i.sum} + #end for + </command> + <inputs> + <param format="interval" name="input1" type="data" label="Select first file to add" /> + <param name="sum1" type="boolean" checked="true" truevalue="true" falsevalue="false" label="Include first file's totals in overall total" /> + <param format="interval" name="input2" type="data" label="Select second file to add" /> + <param name="sum2" type="boolean" checked="true" truevalue="true" falsevalue="false" label="Include second file's totals in overall total" /> + <repeat name="inputs" title="Input Files"> + <param name="input" label="Add file" type="data" format="interval" /> + <param name="sum" type="boolean" checked="true" truevalue="true" falsevalue="false" label="Include file's totals in overall total" /> + </repeat> + </inputs> + <outputs> + <data format="interval" name="output1" /> + </outputs> + <tests> + <test> + <param name="input1" value="indel_table_in1.interval" ftype="interval" /> + <param name="sum1" value="true"/> + <param name="input2" value="indel_table_in2.interval" ftype="interval" /> + <param name="sum2" value="true" /> + <param name="input" value="indel_table_in3.interval" ftype="interval" /> + <param name="sum" value="true" /> + <output name="output1" file="indel_table_out1.interval" ftype="interval" /> + </test> + </tests> + <help> + +**What it does** + +Creates a table allowing for analysis and comparison of indel data. Combines any number of interval files that have been produced by the tool that converts indel SAM data to interval format. Includes overall total counts for all or some files. The tool has the option to not include a given file's counts in the total column. This could be useful for combined data if the counts for certain indels might be included more than once. + +The exact columns of the output will depend on the columns of the input. Here is the detailed specification of the output columns:: + + Column Description + ------------------------------- ---------------------------------------------------------------------------------- + 1 ... m "Indel" All the "indel" columns, which contain the info that will be checked for equality + m + 1 Total Occurrences Total number of occurrences of this indel across all (included) files + m + 2 Occurrences for File 1 Number of occurrences of this indel for first file + m + 3 Occurrences for File 2 Number of occurrences of this indel for second file + [m + ...] [...] [Number of occurrences of this indel for ... file] + +The most likely columns would be from the output of the Convert SAM to Interval/BED tool, so: Chromosome, Start position, End position, I/D (Insertion/Deletion), -/<base(s)> (Deletion/Inserted base(s)), Total Occurrences (across files), Occurrences for File 1, Occurrences for File 2, etc. See below for an example. + + +----- + +**Example** + +Suppose you have the following 4 files:: + + chrM 300 301 D - 6 + chrM 303 304 D - 19 + chrM 359 360 D - 1 + chrM 410 411 D - 1 + chrM 435 436 D - 1 + + chrM 410 411 D - 1 + chrM 714 715 D - 1 + chrM 995 997 D - 1 + chrM 1168 1169 I A 1 + chrM 1296 1297 D - 1 + + chrM 300 301 D - 8 + chrM 525 526 D - 1 + chrM 958 959 D - 1 + chrM 995 996 D - 3 + chrM 1168 1169 I C 1 + chrM 1296 1297 D - 1 + + chrM 303 304 D - 22 + chrM 410 411 D - 1 + chrM 435 436 D - 1 + chrM 714 715 D - 1 + chrM 753 754 I A 1 + chrM 1168 1169 I A 1 + +and the fifth file:: + + chrM 303 304 D - 22 + chrM 410 411 D - 2 + chrM 435 436 D - 1 + chrM 714 715 D - 2 + chrM 753 754 I A 1 + chrM 995 997 D - 1 + chrM 1168 1169 I A 2 + chrM 1296 1297 D - 1 + +The following will be produced if you include the first four files in the sum, but not the fifth:: + + chrM 300 301 D - 14 6 0 8 0 0 + chrM 303 304 D - 41 19 0 0 22 22 + chrM 359 360 D - 1 1 0 0 0 0 + chrM 410 411 D - 3 1 1 0 1 2 + chrM 435 436 D - 2 1 0 0 1 2 + chrM 525 526 D - 1 0 0 1 0 0 + chrM 714 715 D - 2 0 1 0 1 2 + chrM 753 754 I A 1 0 0 0 1 1 + chrM 958 959 D - 1 0 0 1 0 0 + chrM 995 996 D - 3 0 0 3 0 0 + chrM 995 997 D - 1 0 1 0 0 1 + chrM 1168 1169 I A 2 0 1 0 1 2 + chrM 1168 1169 I C 1 0 0 1 0 0 + chrM 1296 1297 D - 2 0 1 1 0 1 + +The first numeric column includes the total or the next four columns, but not the fifth. + + + </help> +</tool> --- /dev/null +++ b/tools/indels/sam_indel_filter.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python + +""" +Allows user to filter out non-indels from SAM. + +usage: %prog [options] + -i, --input=i: Input SAM file to be filtered + -q, --quality_threshold=q: Minimum quality value for adjacent bases + -a, --adjacent_bases=a: Number of adjacent bases on each size to check qualities + -o, --output=o: Filtered output SAM file +""" + +import re, sys +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse + + +def stop_err( msg ): + sys.stderr.write( '%s\n' % msg ) + sys.exit() + +def __main__(): + #Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + # prep output file + output = open( options.output, 'wb' ) + # patterns + pat_indel = re.compile( '(?P<before_match>(\d+[MNSHP])*)(?P<lmatch>\d+)M(?P<ins_del_width>\d+)(?P<ins_del>[ID])(?P<rmatch>\d+)M(?P<after_match>(\d+[MNSHP])*)' ) + pat_matches = re.compile( '(\d+[MIDNSHP])+' ) + try: + qual_thresh = int( options.quality_threshold ) + 33 + if qual_thresh < 33 or qual_thresh > 126: + raise ValueError + except ValueError: + stop_err( 'Your quality threshold should be an integer between 0 and 93, inclusive.' ) + try: + adj_bases = int( options.adjacent_bases ) + if adj_bases < 1: + raise ValueError + except ValueError: + stop_err( 'The number of adjacent bases should be an integer greater than 1.' ) + # record lines skipped because of more than one indel + multi_indel_lines = 0 + # go through all lines in input file + for i,line in enumerate(open( options.input, 'rb' )): + if line and not line.startswith( '#' ) and not line.startswith( '@' ) : + split_line = line.split( '\t' ) + cigar = split_line[5] + # find all possible matches, like 3M2D7M and 7M3I10M in 3M2D7M3I10M + cigar_copy = cigar[:] + matches = [] + while len( cigar_copy ) >= 6: # nMnInM or nMnDnM + m = pat_indel.search( cigar_copy ) + if not m: + break + else: + parts = m.groupdict() + if parts[ 'lmatch' ] and parts[ 'ins_del_width' ] and parts[ 'rmatch' ]: + pre_left = 0 + if m.start() > 0: + pre_left_groups = pat_matches.search( cigar_copy[ : m.start() ] ) + if pre_left_groups: + for pl in pre_left_groups.groups(): + if pl.endswith( 'M' ) or pl.endswith( 'S' ) or pl.endswith( 'P' ): + pre_left += pl[:-1] + parts[ 'pre_left' ] = pre_left + matches.append( parts ) + cigar_copy = cigar_copy[ len( parts[ 'lmatch' ] ) + 1 : ] + # see if matches meet filter requirements + if len( matches ) > 1: + multi_indel_lines += 1 + elif len( matches ) == 1: + pre_left = int( matches[0][ 'pre_left' ] ) + left = int( matches[0][ 'lmatch' ] ) + right = int( matches[0][ 'rmatch' ] ) + if matches[0][ 'ins_del' ] == 'D': + middle = int( matches[0][ 'ins_del_width' ] ) + else: + middle = 0 + # if there are enough adjacent bases to check, then do so + if left >= adj_bases and right >= adj_bases: + quals = split_line[10] + left_quals = quals[ pre_left : pre_left + left ][ -adj_bases : ] + middle_quals = quals[ pre_left + left : pre_left + left + middle ] + right_quals = quals[ pre_left + left + middle : pre_left + left + middle + right ][ : adj_bases ] + qual_thresh_met = True + for l in left_quals: + if ord( l ) < qual_thresh: + qual_thresh_met = False + break + if qual_thresh_met: + for r in right_quals: + if ord( r ) < qual_thresh: + qual_thresh_met = False + break + # if filter reqs met, output line + if qual_thresh_met: + output.write( line ) + # close out file + output.close() + # if skipped lines because of more than one indel, output message + if multi_indel_lines > 0: + sys.stdout.write( '%s alignments were skipped because they contained more than one indel.' % multi_indel_lines ) + +if __name__=="__main__": __main__() --- /dev/null +++ b/test-data/indel_sam2interval_out2.bed @@ -0,0 +1,4 @@ +ref 133 134 +ref 48819784 48819785 +ref 87824726 87824727 +ref 188841437 188841438 --- /dev/null +++ b/test-data/indel_sam2interval_out3.bed @@ -0,0 +1,2 @@ +ref 256 258 +ref 190341171 190341172 --- a/test-data/sam_indel_filter_out3.sam +++ b/test-data/sam_indel_filter_out3.sam @@ -1,10 +1,10 @@ -081017-and-081020:1:6:774:1836 0 PHIX174 4973 37 26M1I9M * 0 0 GCTTAAAGCTACCAGTTATATGGCTGTTTGGTTTTT IIIIIIIIIIIIIIIIIIIIII@III/IE;%II;I= XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:33C1 -081017-and-081020:1:6:1193:793 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIII7IIIDIIIIIIII,=(>%II? XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 -081017-and-081020:1:8:753:970 0 PHIX174 2974 37 29M1I6M * 0 0 GTGGCGCCATGTCTAAATTGTTTGGAGGCGGGTCAA IIIIIIIIIIIIIIII4IIIIII3I&IIIII*%%&' XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:35 -081017-and-081020:1:18:1164:1678 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIIIIIIIII;IIII1I0I)II.I- XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 -081017-and-081020:1:20:754:1256 16 PHIX174 4772 37 4M1I31M * 0 0 CCCTGGCGGTGCATTTTATGCGGACACTTCCTACAG &II(IIIII3IIIII7II,IIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:2A32 -081017-and-081020:1:24:1326:917 0 PHIX174 188 37 25M1D11M * 0 0 TTCGCCATCAACTAACGATTCTGTCAAACCTGACGC IIIIIIIIIIIIIIIIIIIIIIIIII2/&II>'IEI XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:25^A3A7 -081017-and-081020:1:25:1466:511 16 PHIX174 1179 37 21M1D15M * 0 0 CTATTGACTCTACTGTAGACATTTTACTTTTTATGT :I<=IIIIII5IIGI5IIIIIIIIIIIIIIIIIIII XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:21^T15 -081017-and-081020:1:27:1267:1275 0 PHIX174 3716 37 28M1D8M * 0 0 TCATCAGCAAACGCAGAATCAGCGGTATGCTCTTCT IIIIIIIIIIIIIIIII;IIIIIII87III%I(@I. XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:28^G8 -081017-and-081020:1:94:1649:147 16 PHIX174 2755 37 15M2D21M * 0 0 TATACCGTCAAGGACTGTGACTATTGACGTCCTTCC 4IIIIII@I7IIIIIIIIIIIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:15^TG21 -081017-and-081020:1:95:74:43 0 PHIX174 4038 37 29M1I6M * 0 0 ATCGAGGCTCTTAAACCTGCTATTGAGGCTTTTTGG IIIIIIIIIIIIIIIIICI;8I,I>IIIIII1I%5& XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:6:774:1836 0 PHIX174 4973 37 26M1I9M * 0 0 GCTTAAAGCTACCAGTTATATGGCTGTTTGGTTTTT IIIIIIIIIIIIIIIIIIIIII@III/IE;%II;I= XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:33C1 +081017-and-081020:1:6:1193:793 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIII7IIIDIIIIIIII,=(>%II? XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:8:753:970 0 PHIX174 2974 37 29M1I6M * 0 0 GTGGCGCCATGTCTAAATTGTTTGGAGGCGGGTCAA IIIIIIIIIIIIIIII4IIIIII3I&IIIII*%%&' XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:35 +081017-and-081020:1:18:1164:1678 0 PHIX174 4971 37 27M1I8M * 0 0 CCGCTTAAAGCTACCAGTTATATGGCTGGTTGTTTT IIIIIIIIIIIIIIIIIIIII;IIII1I0I)II.I- XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 +081017-and-081020:1:20:754:1256 16 PHIX174 4772 37 4M1I31M * 0 0 CCCTGGCGGTGCATTTTATGCGGACACTTCCTACAG &II(IIIII3IIIII7II,IIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:2A32 +081017-and-081020:1:24:1326:917 0 PHIX174 188 37 25M1D11M * 0 0 TTCGCCATCAACTAACGATTCTGTCAAACCTGACGC IIIIIIIIIIIIIIIIIIIIIIIIII2/&II>'IEI XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:25^A3A7 +081017-and-081020:1:25:1466:511 16 PHIX174 1179 37 21M1D15M * 0 0 CTATTGACTCTACTGTAGACATTTTACTTTTTATGT :I<=IIIIII5IIGI5IIIIIIIIIIIIIIIIIIII XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:21^T15 +081017-and-081020:1:27:1267:1275 0 PHIX174 3716 37 28M1D8M * 0 0 TCATCAGCAAACGCAGAATCAGCGGTATGCTCTTCT IIIIIIIIIIIIIIIII;IIIIIII87III%I(@I. XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:1 MD:Z:28^G8 +081017-and-081020:1:94:1649:147 16 PHIX174 2755 37 15M2D21M * 0 0 TATACCGTCAAGGACTGTGACTATTGACGTCCTTCC 4IIIIII@I7IIIIIIIIIIIIIIIIIIIIIIIIII XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:15^TG21 +081017-and-081020:1:95:74:43 0 PHIX174 4038 37 29M1I6M * 0 0 ATCGAGGCTCTTAAACCTGCTATTGAGGCTTTTTGG IIIIIIIIIIIIIIIIICI;8I,I>IIIIII1I%5& XT:A:U NM:i:2 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:31G3 --- /dev/null +++ b/tools/indels/indel_analysis.xml @@ -0,0 +1,150 @@ +<tool id="indel_analysis" name="Indel Analysis" version="1.0.0"> + <description></description> + <command interpreter="python"> + indel_analysis.py + --input=$input1 + --threshold=$threshold + --out_ins=$out_ins + --out_del=$out_del + </command> + <inputs> + <param format="sam" name="input1" type="data" label="Select sam file to analyze" /> + <param name="threshold" type="float" value="0.015" size="5" label="Frequency threshold" help="Cutoff" /> + </inputs> + <outputs> + <data format="interval" name="out_del" /> + <data format="interval" name="out_ins" /> + </outputs> + <tests> + <test> + <param name="input1" value="indel_analysis_in1.sam" ftype="sam"/> + <param name="threshold" value="0.017"/> + <output name="out_del" file="indel_analysis_out1.interval" ftype="interval"/> + <output name="out_ins" file="indel_analysis_out2.interval" ftype="interval"/> + </test> + <test> + <param name="input1" value="indel_analysis_in2.sam" ftype="sam"/> + <param name="threshold" value="0.08"/> + <output name="out_del" file="indel_analysis_out3.interval" ftype="interval"/> + <output name="out_ins" file="indel_analysis_out4.interval" ftype="interval"/> + </test> + </tests> + <help> + +**What it does** + +Given an input sam file, this tool provides analysis of the indels. It filters out matches that do not meet the frequency threshold. The way this frequency of occurence is calculated is different for deletions and insertions. The CIGAR string's "M" can indicate an exact match or a mismatch. For SAM containing the following bits of information (assuming the reference "ACTGCTCGAT"):: + + CHROM POS CIGAR SEQ + ref 3 2M1I3M TATCTC + ref 2 3M1D2M ATGTC + ref 4 2M2I3M GTTGAAG + ref 1 2M2D2M ACCT + ref 2 3M1I2M TCCATC + ref 7 4M CTAT + ref 5 5M CGTGA + +The following totals would be calculated:: + + Counts for chromosome "ref", where "-" indicates a deletion and + "+" an insertion + ---------------------------------------------------------------- + POS BASE NUMREADS DELPROPCALC DELPROP INSPROPCALC INSPROP + 1 A 1 1/1 1.00 -- -- + 2 A 1 1/2 0.50 -- -- + C 1 1/2 0.50 -- -- + 3 T 3 3/4 0.75 -- -- + -- 1 1/4 0.25 -- -- + 4 A 1 1/5 0.20 -- -- + G 2 2/5 0.40 -- -- + C 1 1/5 0.20 -- -- + - 1 1/5 0.20 -- -- + 5 C 4 4/5 0.80 -- -- + T 1 1/5 0.20 -- -- + +A 1 --- --- 1/6 0.17 + +T 1 --- --- 1/6 0.17 + 6 A 1 1/6 0.17 -- -- + C 1 1/6 0.17 -- -- + T 4 4/6 0.67 -- -- + +TG 1 --- --- 1/7 0.14 + 7 A 1 1/6 0.17 -- -- + C 3 3/6 0.50 -- -- + G 1 1/6 0.17 -- -- + T 1 1/6 0.17 -- -- + 8 G 2 2/3 0.67 -- -- + T 1 1/3 0.33 -- -- + 9 A 2 2/2 1.00 -- -- + 10 T 1 1/1 1.00 -- -- + +Note that the way these are calculated may not be immediately clear. First, the basic total number of reads at a given position is the number of reads with a particular base plus the number of reads with that a deletion at that given position. Note that deletions of two bases and one base would be counted separately. Insertions are not counted in this total, which is used to calculate the proportion of occurrences of each individual base and deletion. For position 4 above, the reference base is G, and there are 2 occurrences of it along with one each of mismatching bases A and C. Also, there is on 1-base deletion. So there are a total of 5 matches/mismatches/deletions, and the proportions for each base are either 1/5 = 0.20 or 2/5 = 0.40, and for the deletion it is 1/5 = 0.20. Insertions are slightly more complicated. Each insertion is regarded individually, and the total number of occurrences of that insertion is divided by the sum of the number of its occurrences and the b asic total. So for position 5, there are a total of 5 matches/mismatches/deletions, and two insertions that each occur once, so each has a insertion has a proportion of 1/6 = 0.17. + +The DELPROP or INSPROP needs to be greater than the threshold frequency specified by the user. + +The output varies for deletions and insertions, though for both, the first three columns are chromosome, start position, and end position. + +Columns in the deletions file:: + + Column Description + ---------------------------- ------------------------------------------------------------------------------------ + 1 Chrom Chromosome + 2 Start Starting position + 3 End Ending position + 4 Number of Deleted Base(s) The number of bases deleted at Start position + 5 Frequency Percentage Frequency of this exact deletion (2 and 1 are mutually exclusive), as percentage (%) + +Columns in the insertions file:: + + Column Description + -------------------------- ------------------------------------------------------------------------------------------------------------- + 1 Chrom Chromosome + 2 Start Starting position + 3 End Ending position (always Start + 1 for insertions) + 4 Inserted Base(s) The exact base(s) inserted at Start position + 5 Freq. Perc. at Start Frequency of this exact insertion given Start position ("GG" and "G" are considered distinct), as percentage (%) + 6 Freq. Perc. at End Frequency of this exact insertion given End position ("GG" and "G" are considered distinct), as percentage (%) + +Before using this tool, you probably will want to use the Filter SAM for indels tool to filter out indels on bases with insufficient quality scores, but this is not required. + + +----- + +**Example** + +If you set the threshold to 0.0 and have the following SAM file:: + + r770 89 ref 6 37 7M1I5M = 0 0 TGGATCTTCATAG !0//110AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r1124 113 ref 4 0 23M = 0 0 CATCGTTCTGTTAGATCTACGTA PQRVUMNXYRPUXYXWXSOSZ]M XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r1789 177 ref 6 0 17M = 0 0 TCGATCGCTTAGTTCTC SQQWZY]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r3671 153 ref 10 37 6M1I6M = 0 0 TCTCTTTAGGTCT /<<!"0/////// XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r3824 153 ref 5 37 8M1I7M = 0 0 ATTGATGTTCTTAGAT 4;6//11!"100110/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r4800 16 ref 7 255 5M2D6M = 0 0 CGATCTTTGAT IIIIIIIIIIC XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5612 151 ref 5 37 8M1D9M = 0 0 ATCTATCTTTTGATCTC /<<!"0/4/*/7//B0/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5612 151 ref 11 37 3M1I10M = 0 0 CTCCTTAGCTCTCC /<<!"0/4//7//0 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r9145 115 ref 11 0 19M = 0 -1 CTCTTAGCTCTCCGAATTAG 7753:<5#"4!&=9518A> XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:137 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 + r11770 89 ref 10 37 10M2I8M = 0 0 TCTCTTAGATGGCTCCGTAT 00/02!!0/120210AA4/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r13671 153 ref 1 37 12M1I12M = 0 0 TCGCATCGATCTCCGTAGATCTCCG /<""<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r13824 153 ref 13 37 9M1I7M = 0 0 CATAGATCTACCGGATT 4;6//11!"11100110 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r24800 16 ref 3 255 15M2D9M = 0 0 GCATCTATCTGATAGCTCCGAATT IIIIIIIII45"CCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r25612 151 ref 1 37 9M1D5M = 0 0 TCGCATCGACTCTT 0/4/*/7//00/1C XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r25612 151 ref 21 37 4M1I7M = 0 0 TGCGTTATTGGG <!"0/70/BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r29962 16 ref 20 37 4M1I7M = 0 0 CTCCGGTATGAGG <!"0/70/7BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + +The following will be produced (deletions file followed by insertions file):: + + ref 10 11 1 1 9.09 + ref 12 14 2 1 7.69 + ref 13 14 1 1 7.69 + ref 18 20 2 1 8.33 + + ref 13 14 C 1 6.25 6.67 + ref 13 14 T 2 12.50 13.33 + ref 14 15 C 1 6.67 7.14 + ref 16 17 T 1 7.14 7.69 + ref 20 21 GG 1 8.33 8.33 + ref 22 23 A 1 8.33 11.11 + ref 24 25 G 1 11.11 12.50 + ref 25 26 T 1 12.50 14.29 + + + </help> +</tool> --- /dev/null +++ b/test-data/indel_analysis_out2.interval @@ -0,0 +1,8 @@ +ref 13 14 C 1 7.14 7.14 +ref 13 14 T 2 13.33 14.29 +ref 14 15 C 1 6.67 7.14 +ref 16 17 T 1 7.14 7.69 +ref 20 21 GG 1 8.33 8.33 +ref 22 23 A 1 8.33 11.11 +ref 24 25 G 1 11.11 14.29 +ref 25 26 T 1 12.50 14.29 --- /dev/null +++ b/test-data/indel_analysis_in2.sam @@ -0,0 +1,16 @@ +r770 89 ref 6 37 7M1I5M = 0 0 TCGATCTTCATAG !0//110AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r1124 113 ref 4 0 23M = 0 0 CATCGTTCTGTTAGATCTACGTA PQRVUMNXYRPUXYXWXSOSZ]M XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r1789 177 ref 6 0 17M = 0 0 TCGATCGCTTAGTTCTC SQQWZY]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r3671 153 ref 10 37 6M1I6M = 0 0 TCTCTTTAGGTCT /<<!"0/////// XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r3824 153 ref 5 37 8M1I7M = 0 0 ATCGATGTTCTTAGAT 4;6//11!"100110/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r4800 16 ref 7 255 5M2D6M = 0 0 CGATCTTTGAT IIIIIIIIIIC XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r5612 151 ref 5 37 8M1D9M = 0 0 ATCTATCTTTTGATCTC /<<!"0/4/*/7//B0/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r5929 151 ref 11 37 3M1I10M = 0 0 CTCCTTAGCTCTCC /<<!"0/4//7//0 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r9145 115 ref 11 0 19M = 0 -1 CTCTTAGCTCTCCGAATTAG 7753:<5#"4!&=9518A> XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:137 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 +r11770 89 ref 10 37 10M2I8M = 0 0 TCTCTTAGATGGCTCCGTAT 00/02!!0/120210AA4/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r13671 153 ref 1 37 12M1I12M = 0 0 TCGCATCGATCTCCTTAGATCTCCG /<""<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r13824 153 ref 13 37 9M1I7M = 0 0 CATAGATCTACCGGATT 4;6//11!"11100110 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r24800 16 ref 3 255 15M2D9M = 0 0 GCATCTATCTGATAGCTCCGAATT IIIIIIIII45"CCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r25612 151 ref 1 37 9M1D5M = 0 0 TCGCATCGACTCTT 0/4/*/7//00/1C XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r25786 151 ref 21 37 4M1I7M = 0 0 TGCGTTATTGGG <!"0/70/BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r29962 16 ref 20 37 4M1I7M = 0 0 CTCCGGTATGAGG <!"0/70/7BC01 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 --- a/test-data/sam_indel_filter_in1.sam +++ b/test-data/sam_indel_filter_in1.sam @@ -1,15 +1,15 @@ -1378_28_770 89 chr11.nib:1-134452384 72131356 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_28_770 89 chr11.nib:1-134452384 72131356 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 1378_28_770 181 chr11.nib:1-134452384 72131356 0 13M2I6M1D5M = 72131356 0 TTGGTGCGCGCGGTTGAGGGTTGG $$(#%%#$%#%####$%%##$### -1378_33_1945 113 chr2.nib:1-242951149 181247988 0 23M chr12.nib:1-132349534 41710908 0 GAGAGAGAGAGAGAGAGAGAGAG PQRVUMNXYRPUXYXWXSOSZ]M XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 -1378_33_1945 177 chr12.nib:1-132349534 41710908 0 23M chr2.nib:1-242951149 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 -1378_35_263 115 chr16.nib:1-88827254 19671878 0 23M = 19671877 -1 AGAGAGAGAGAGAGAGAGAGTCT 77543:<55#"4!&=964518A> XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:137 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 -1378_35_263 179 chr16.nib:1-88827254 19671877 0 13M1I4M1I5M = 19671878 1 GAGAGAGAGAGAGAGAGAGAGTC LE7402DD34FL:27AKE>;432 XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:265 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +1378_33_1945 113 chr2.nib:1-242951149 181247988 0 23M chr12.nib:1-132349534 41710908 0 GAGAGAGAGAGAGAGAGAGAGAG PQRVUMNXYRPUXYXWXSOSZ]M XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +1378_33_1945 177 chr12.nib:1-132349534 41710908 0 23M chr2.nib:1-242951149 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +1378_35_263 115 chr16.nib:1-88827254 19671878 0 23M = 19671877 -1 AGAGAGAGAGAGAGAGAGAGTCT 77543:<55#"4!&=964518A> XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:137 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 +1378_35_263 179 chr16.nib:1-88827254 19671877 0 13M1I4M1I5M = 19671878 1 GAGAGAGAGAGAGAGAGAGAGTC LE7402DD34FL:27AKE>;432 XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:265 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 1378_51_1671 117 chr2.nib:1-242951149 190342418 0 24M = 190342418 0 CTGGCGTTCTCGGCGTGGATGGGT #####$$##$#%#%%###%$#$## -1378_51_1671 153 chr2.nib:1-242951149 190342418 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_51_1671 153 chr2.nib:1-242951149 190342418 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 1378_56_324 117 chr2.nib:1-242951149 80324999 0 24M = 80324999 0 TCCAGTCGCGTTGTTAGGTTCGGA #$#$$$#####%##%%###**#+/ -1378_56_324 153 chr2.nib:1-242951149 80324999 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 -1378_67_1795 81 chr16.nib:1-88827254 26739130 0 23M chrY.nib:1-57772954 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 -1378_67_1795 161 chrY.nib:1-57772954 57401793 37 23M chr16.nib:1-88827254 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +1378_56_324 153 chr2.nib:1-242951149 80324999 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_67_1795 81 chr16.nib:1-88827254 26739130 0 23M chrY.nib:1-57772954 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 +1378_67_1795 161 chrY.nib:1-57772954 57401793 37 23M chr16.nib:1-88827254 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 1378_69_800 16 chr11.nib:1-125234658 241 255 15M1D8M * 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII -1378_69_1777 170 chrX.nib:1-59090954 59090793 37 23M chr16.nib:1-88827254 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 -1378_72_1612 151 chrY.nib:1-124295114 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +1378_69_1777 170 chrX.nib:1-59090954 59090793 37 23M chr16.nib:1-88827254 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +1378_72_1612 151 chrY.nib:1-124295114 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 --- /dev/null +++ b/test-data/indel_analysis_out1.interval @@ -0,0 +1,4 @@ +ref 10 11 1 1 9.09 +ref 12 14 2 1 7.69 +ref 13 14 1 1 7.69 +ref 18 20 2 1 8.33 --- /dev/null +++ b/tools/indels/indel_sam2interval.xml @@ -0,0 +1,121 @@ +<tool id="indel_sam2interval" name="Convert SAM to interval/BED" version="1.0.0"> + <description>for indels</description> + <command interpreter="python"> + indel_sam2interval.py + --input=$input1 + --include_base=$include_base + --collapse=$collapse + --int_out=$output1 + #if $ins_out.include_ins_out == "true" + --bed_ins_out=$output2 + #else + --bed_ins_out="None" + #end if + #if $del_out.include_del_out == "true" + --bed_del_out=$output3 + #else + --bed_del_out="None" + #end if + </command> + <inputs> + <param format="sam" name="input1" type="data" label="Select dataset to convert" /> + <param name="include_base" type="boolean" checked="true" truevalue="true" falsevalue="false" label="Include the relevant base(s) for each insertion (and a dash (-) for deletions)" /> + <param name="collapse" type="boolean" checked="true" truevalue="true" falsevalue="false" label="Collapse repeated locations onto single line with counts" /> + <conditional name="ins_out"> + <param name="include_ins_out" type="select" label="Include insertions output bed file?"> + <option value="true">Yes</option> + <option value="false">No</option> + </param> + <when value="true" /> + <when value="false" /> + </conditional> + <conditional name="del_out"> + <param name="include_del_out" type="select" label="Include insertions output bed file?"> + <option value="true">Yes</option> + <option value="false">No</option> + </param> + <when value="true" /> + <when value="false" /> + </conditional> + </inputs> + <outputs> + <data format="interval" name="output1" /> + <data format="bed" name="output2"> + <filter>ins_out[ "include_ins_out" ] = "true"</filter> + </data> + <data format="bed" name="output3"> + <filter>del_out[ "include_del_out" ] = "true"</filter> + </data> + </outputs> + <tests> + <test> + <param name="input1" value="indel_sam2interval_in1.sam" ftype="sam"/> + <param name="include_base" value="true"/> + <param name="collapse" value="true"/> + <param name="include_ins_out" value="true" /> + <param name="include_del_out" value="true" /> + <output name="output1" file="indel_sam2interval_out1.interval" ftype="interval"/> + <output name="output2" file="indel_sam2interval_out2.bed" ftype="bed"/> + <output name="output3" file="indel_sam2interval_out3.bed" ftype="bed"/> + </test> + </tests> + <help> + +**What it does** + +Given a SAM file containing indels, converts these to an interval file with a column indicating whether it is an insertion or a deletion, and then also can create a BED file for each type (one for insertions, one for deletions). The interval file can be combined with other like files to create a table useful for analysis with the Indel Analysis Table tool. The BED files can be useful for visualizing the reads. + +----- + +**Example** + +Suppose you have the following:: + + r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r780 181 ref 4567 0 24M = 72131356 0 TTGGTGCGCGCGGTTGAGGGTTGG $$(#%%#$%#%####$%%##$### XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r1231 69 * 0 0 * * 0 0 AGACCGGGCGGGGTGGCGTTCGGT %##+'#######%###$#$##$(# + r1563 133 * 0 0 * * 0 0 GTTCGTGGCCGGTGGGTGTTTGGG ###$$#$#$&#####$'$#$###$ + r1945 177 ref 71908 0 23M 190342418 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r3671 117 ref 31903418 0 24M = 190342418 0 CTGGCGTTCTCGGCGTGGATGGGT #####$$##$#%#%%###%$#$## XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r3673 153 ref 48819768 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r3824 117 ref 80729921 0 24M = 80324999 0 TCCAGTCGCGTTGTTAGGTTCGGA #$#$$$#####%##%%###**#+/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r3911 153 ref 87824718 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r4795 81 ref 126739130 0 23M 57401793 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 + r4797 161 ref 57401793 37 23M 26739130 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r4800 16 ref 241 255 15M2D8M = 1550 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5377 170 ref 52090793 37 23M 26739130 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r5612 151 ref 190341152 37 19M1D4M = 190342418 0 TCTAACTTAGCCTCATAATGCTAA /<<!"0/4/*/7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5623 151 ref 188841418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r7899 69 * 0 0 * * 0 0 CTGCGTGTTGGTGTCTACTGGGGT #%#'##$#$##&%#%$$$%#%#'# + r9192 133 * 0 0 * * 0 0 GTGCGTCGGGGAGGGTGCTGTCGG ######%#$%#$$###($###&&% + + +The following three files will be produced (Interval, Insertions BED and Deletions BED):: + + ref 133 134 I C 1 + ref 256 258 D - 1 + ref 48819784 48819785 I A 1 + ref 87824726 87824727 I G 1 + ref 188841437 188841419 I T 1 + ref 190341171 190341172 D - 1 + + ref 133 134 + ref 48819784 48819785 + ref 87824726 87824727 + ref 188841437 188841419 + + ref 256 258 + ref 190341171 190341172 + + + + + + +For more information on SAM, please consult the `SAM format description`__. + +.. __: http://www.ncbi.nlm.nih.gov/pubmed/19505943 + + + </help> +</tool> --- /dev/null +++ b/tools/indels/indel_sam2interval.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python + +""" +Allows user to filter out non-indels from SAM. + +usage: %prog [options] + -i, --input=i: The input SAM file + -u, --include_base=u: Whether or not to include the base for insertions + -c, --collapse=c: Wheter to collapse multiple occurrences of a location with counts shown + -o, --int_out=o: The interval output file for the converted SAM file + -b, --bed_ins_out=b: The bed output file with insertions only for the converted SAM file + -d, --bed_del_out=d: The bed output file with deletions only for the converted SAM file +""" + +import re, sys +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse + + +def stop_err( msg ): + sys.stderr.write( '%s\n' % msg ) + sys.exit() + +def numeric_sort( text1, text2 ): + """ + For two items containing space-separated text, compares equivalent pieces + numerically if both numeric or as text otherwise + """ + pieces1 = text1.split() + pieces2 = text2.split() + if len( pieces1 ) == 0: + return 1 + if len( pieces2 ) == 0: + return -1 + for i, pc1 in enumerate( pieces1 ): + if i == len( pieces2 ): + return 1 + if not pieces2[i].isdigit(): + if pc1.isdigit(): + return -1 + else: + if pc1 > pieces2[i]: + return 1 + elif pc1 < pieces2[i]: + return -1 + else: + if not pc1.isdigit(): + return 1 + else: + if int( pc1 ) > int( pieces2[i] ): + return 1 + elif int( pc1 ) < int( pieces2[i] ): + return -1 + if i < len( pieces2 ) - 1: + return -1 + return 0 + +def __main__(): + #Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + + # open up output files + output = open( options.int_out, 'wb' ) + if options.bed_ins_out != 'None': + output_bed_ins = open( options.bed_ins_out, 'wb' ) + else: + output_bed_ins = None + if options.bed_del_out != 'None': + output_bed_del = open( options.bed_del_out, 'wb' ) + else: + output_bed_del = None + + # the pattern to match, assuming just one indel per cigar string + pat_indel = re.compile( '(?P<lmatch>\d+)M(?P<ins_del_width>\d+)(?P<ins_del>[ID])(?P<rmatch>\d+)M' ) + pat_multi = re.compile( '(\d+[MIDNSHP])(\d+[MIDNSHP])(\d+[MIDNSHP])+' ) + + # go through all lines in input file + out_data = [] + multi_indel_lines = 0 + for line in open( options.input, 'rb' ): + if line and not line.startswith( '#' ) and not line.startswith( '@' ) : + split_line = line.split( '\t' ) + if split_line < 12: + continue + # grab relevant pieces + cigar = split_line[5] + pos = int( split_line[3] ) + chr = split_line[2] + base_string = split_line[9] + # parse cigar string + m = pat_indel.search( cigar ) + if not m: + m = pat_multi.search( cigar ) + # skip this line if no match + if not m: + continue + # account for multiple indels or operations we don't process + else: + multi_indel_lines += 1 + continue + else: + match = m.groupdict() + left = int( match[ 'lmatch' ] ) + middle = int( match[ 'ins_del_width' ] ) + middle_type = match[ 'ins_del' ] + bases = base_string[ left : left + middle ] + # calculate start and end positions, and output to insertion or deletion file + start = left + pos + if middle_type == 'D': + end = start + middle + d = [ chr, start, end, 'D' ] + if options.include_base == "true": + d.append( '-' ) + out_data.append( tuple( d ) ) + if output_bed_del: + output_bed_del.write( '%s\t%s\t%s\n' % ( chr, start, end ) ) + else: + end = start + 1#+ middle + d = [ chr, start, end, 'I' ] + if options.include_base == "true": + d.append( bases ) + out_data.append( tuple( d ) ) + if output_bed_ins: + output_bed_ins.write( '%s\t%s\t%s\n' % ( chr, start, end ) ) + # output to interval file + if options.collapse == 'true': + out_dict = {} + # first collapse and get counts + for data in out_data: + location = ' '.join( [ '%s' % d for d in data ] ) + try: + out_dict[ location ].append( data ) + except KeyError: + out_dict[ location ] = [ data ] + locations = out_dict.keys() + locations.sort( numeric_sort ) + for loc in locations: + output.write( '%s\t%s\n' % ( '\t'.join( [ '%s' % d for d in out_dict[ loc ][0] ] ), len( out_dict[ loc ] ) ) ) + else: + for data in out_data: + output.write( '%s\n' % '\t'.join( [ '%s' % d for d in data ] ) ) + + # cleanup, close files + if output_bed_ins: + output_bed_ins.close() + if output_bed_del: + output_bed_del.close() + output.close() + + # if skipped lines because of more than one indel, output message + if multi_indel_lines > 0: + sys.stdout.write( '%s alignments were skipped because they contained more than one indel or had unhandled operations (N/S/H/P).' % multi_indel_lines ) + +if __name__=="__main__": __main__() --- /dev/null +++ b/tools/indels/indel_table.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python + +""" +Combines several interval files containing indels with counts. All input files need to have the same number of columns. + +usage: %prog [options] [input3 sum3[ input4 sum4[ input5 sum5[...]]]] + -1, --input1=1: The first input file + -s, --sum1=s: Whether or not to include the totals from first file in overall total + -2, --input2=2: The second input file + -S, --sum2=S: Whether or not to include the totals from second file in overall total + -o, --output=o: The interval output file for the combined files +""" + +import re, sys +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.cookbook import doc_optparse + + +def stop_err( msg ): + sys.stderr.write( '%s\n' % msg ) + sys.exit() + +def numeric_sort( text1, text2 ): + """ + For two items containing space-separated text, compares equivalent pieces + numerically if both numeric or as text otherwise + """ + pieces1 = text1.split() + pieces2 = text2.split() + if len( pieces1 ) == 0: + return 1 + if len( pieces2 ) == 0: + return -1 + for i, pc1 in enumerate( pieces1 ): + if i == len( pieces2 ): + return 1 + if not pieces2[i].isdigit(): + if pc1.isdigit(): + return -1 + else: + if pc1 > pieces2[i]: + return 1 + elif pc1 < pieces2[i]: + return -1 + else: + if not pc1.isdigit(): + return 1 + else: + if int( pc1 ) > int( pieces2[i] ): + return 1 + elif int( pc1 ) < int( pieces2[i] ): + return -1 + if i < len( pieces2 ) - 1: + return -1 + return 0 + +def __main__(): + # Parse Command Line + options, args = doc_optparse.parse( __doc__ ) + inputs = [ options.input1, options.input2 ] + includes = [ options.sum1, options.sum2 ] + inputs.extend( [ a for i, a in enumerate( args ) if i % 2 == 0 ] ) + includes.extend( [ a for i, a in enumerate( args ) if i % 2 == 1 ] ) + num_cols = 0 + counts = {} + # read in data from all files and get total counts + try: + for i, input in enumerate( inputs ): + for line in open( input, 'rb' ): + sp_line = line.strip().split( '\t' ) + # set num_cols on first pass + if num_cols == 0: + if len( sp_line ) < 4: + raise Exception, 'There need to be at least 4 columns in the file: Chrom, Start, End, and Count' + num_cols = len( sp_line ) + # deal with differing number of columns + elif len( sp_line ) != num_cols: + raise Exception, 'All of the files need to have the same number of columns (current %s != %s of first line)' % ( len( sp_line ), num_cols ) + # get actual counts for each indel + indel = '\t'.join( sp_line[:-1] ) + try: + count = int( sp_line[-1] ) + except ValueError, e: + raise Exception, 'The last column of each file must be numeric, with the count of the number of instances of that indel: %s' % str( e ) + # total across all included files + if includes[i] == "true": + try: + counts[ indel ]['tot'] += count + except ( IndexError, KeyError ): + counts[ indel ] = { 'tot': count } + # counts for ith file + counts[ indel ][i] = count + except Exception, e: + stop_err( 'Failed to read all input files:\n%s' % str( e ) ) + # output combined results to table file + try: + output = open( options.output, 'wb' ) + count_keys = counts.keys() + count_keys.sort( numeric_sort ) + for indel in count_keys: + count_out = [ str( counts[ indel ][ 'tot' ] ) ] + for i in range( len( inputs ) ): + try: + count_out.append( str( counts[ indel ][i] ) ) + except KeyError: + count_out.append( '0' ) + output.write( '%s\t%s\n' % ( indel, '\t'.join( count_out ) ) ) + output.close() + except Exception, e: + stop_err( 'Failed to output data: %s' % str( e ) ) + +if __name__=="__main__": __main__() --- /dev/null +++ b/test-data/indel_analysis_out3.interval @@ -0,0 +1,2 @@ +ref 10 11 1 1 9.09 +ref 18 20 2 1 8.33 --- /dev/null +++ b/test-data/indel_table_in3.interval @@ -0,0 +1,10 @@ +chrM 300 301 D - 8 +chrM 303 304 D - 22 +chrM 410 411 D - 2 +chrM 435 436 D - 1 +chrM 525 526 D - 1 +chrM 753 754 I A 1 +chrM 958 959 D - 1 +chrM 995 996 D - 3 +chrM 1168 1169 I C 1 +chrM 1296 1297 D - 1 --- /dev/null +++ b/test-data/indel_sam2interval_in1.sam @@ -0,0 +1,17 @@ +r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r780 181 ref 4567 0 24M = 72131356 0 TTGGTGCGCGCGGTTGAGGGTTGG $$(#%%#$%#%####$%%##$### XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r1231 69 * 0 0 * * 0 0 AGACCGGGCGGGGTGGCGTTCGGT %##+'#######%###$#$##$(# +r1563 133 * 0 0 * * 0 0 GTTCGTGGCCGGTGGGTGTTTGGG ###$$#$#$&#####$'$#$###$ +r1945 177 ref 71908 0 23M 190342418 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r3671 117 ref 31903418 0 24M = 190342418 0 CTGGCGTTCTCGGCGTGGATGGGT #####$$##$#%#%%###%$#$## +r3673 153 ref 48819768 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r3824 117 ref 80729921 0 24M = 80324999 0 TCCAGTCGCGTTGTTAGGTTCGGA #$#$$$#####%##%%###**#+/ +r3911 153 ref 87824718 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r4795 81 ref 126739130 0 23M 57401793 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 +r4797 161 ref 57401793 37 23M 26739130 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r4800 16 ref 241 255 15M2D8M = 1550 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII +r5377 170 ref 52090793 37 23M 26739130 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 +r5612 151 ref 190341152 37 19M1D4M = 190342418 0 TCTAACTTAGCCTCATAATGCTAA /<<!"0/4/*/7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r5623 151 ref 188841418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 +r7899 69 * 0 0 * * 0 0 CTGCGTGTTGGTGTCTACTGGGGT #%#'##$#$##&%#%$$$%#%#'# +r9192 133 * 0 0 * * 0 0 GTGCGTCGGGGAGGGTGCTGTCGG ######%#$%#$$###($###&&% --- /dev/null +++ b/tools/indels/sam_indel_filter.xml @@ -0,0 +1,77 @@ +<tool id="sam_indel_filter" name="Filter SAM" version="1.0.0"> + <description>for indels</description> + <command interpreter="python"> + sam_indel_filter.py + --input=$input1 + --quality_threshold=$quality_threshold + --adjacent_bases=$adjacent_bases + --output=$out_file1 + </command> + <inputs> + <param format="sam" name="input1" type="data" label="Select dataset to filter" /> + <param name="quality_threshold" type="integer" value="40" label="Quality threshold for adjacent bases" help="Takes Phred value assuming Sanger scale; usually between 0 and 40, but up to 93" /> + <param name="adjacent_bases" type="integer" value="1" label="The number of adjacent bases to match on either side of the indel" help="If one side is shorter than this width, the read will be excluded" /> + </inputs> + <outputs> + <data format="sam" name="out_file1" /> + </outputs> + <tests> + <test> + <param name="input1" value="sam_indel_filter_in1.sam" ftype="sam"/> + <param name="quality_threshold" value="14"/> + <param name="adjacent_bases" value="2"/> + <output name="out_file1" file="sam_indel_filter_out1.sam" ftype="sam"/> + </test> + <test> + <param name="input1" value="sam_indel_filter_in1.sam" ftype="sam"/> + <param name="quality_threshold" value="29"/> + <param name="adjacent_bases" value="5"/> + <output name="out_file1" file="sam_indel_filter_out2.sam" ftype="sam"/> + </test> + <test> + <param name="input1" value="sam_indel_filter_in2.sam" ftype="sam"/> + <param name="quality_threshold" value="7"/> + <param name="adjacent_bases" value="1"/> + <output name="out_file1" file="sam_indel_filter_out3.sam" ftype="sam"/> + </test> + </tests> + <help> + +**What it does** + +Allows extracting indels from SAM produced by BWA. Currently it can handle SAM with alignments that have only one insertion or one deletion, and will skip that alignment if it encounters one with more than one indel. It matches CIGAR strings (column 6 in the SAM file) like 5M3I5M or 4M2D10M, so there must be a match or mismatch of sufficient length on either side of the indel. + +----- + +**Example** + +Suppose you have the following:: + + r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r770 181 ref 116 0 24M = 72131356 0 TTGGTGCGCGCGGTTGAGGGTTGG $$(#%%#$%#%####$%%##$### + r1945 177 ref 41710908 0 23M 190342418 181247988 0 AGAGAGAGAGAGAGAGAGAGAGA SQQWZYURVYWX]]YXTSY]]ZM XT:A:R CM:i:0 SM:i:0 AM:i:0 X0:i:163148 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r3671 117 ref 190342418 0 24M = 190342418 0 CTGGCGTTCTCGGCGTGGATGGGT #####$$##$#%#%%###%$#$## + r3671 153 ref 190342418 37 16M1I6M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0///////00/!!0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r3824 117 ref 80324999 0 24M = 80324999 0 TCCAGTCGCGTTGTTAGGTTCGGA #$#$$$#####%##%%###**#+/ + r3824 153 ref 80324999 37 8M1I14M = 80324999 0 TTTAGCCCGAAATGCCTAGAGCA 4;6//11!"11100110////00 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r4795 81 ref 26739130 0 23M 57401793 57401793 0 TGGCATTCCTGTAGGCAGAGAGG AZWWZS]!"QNXZ]VQ]]]/2]] XT:A:R CM:i:2 SM:i:0 AM:i:0 X0:i:3 X1:i:0 XM:i:2 XO:i:0 XG:i:0 MD:Z:23 + r4795 161 ref 57401793 37 23M 26739130 26739130 0 GATCACCCAGGTGATGTAACTCC ]WV]]]]WW]]]]]]]]]]PU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r4800 16 ref 241 255 15M1D8M = 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5377 170 ref 59090793 37 23M 26739130 26739130 0 TATCAATAAGGTGATGTAACTCG ]WV]ABAWW]]]]]P]P//GU]] XT:A:U CM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23 + r5612 151 ref 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + + +To select only alignments with indels, you need to determine the minimum quality you want the adjacent bases to have, as well as the number of adjacent bases to check. If you set the quality threshold to 47 and the number of bases to check to 2, you will get the following output:: + + r770 89 ref 116 37 17M1I5M = 72131356 0 CACACTGTGACAGACAGCGCAGC 00/02!!0//1200210AA44/1 XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r4800 16 ref 241 255 15M1D8M = 0 0 CGTGGCCGGCGGGCCGAAGGCAT IIIIIIIIIICCCCIII?IIIII XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + r5612 151 ref 190342418 37 19M1I3M = 190342418 0 TCTAACTTAGCCTCATAATAGCT /<<!"0/4//7//00/BC0121/ XT:A:U CM:i:2 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:1 XO:i:1 XG:i:1 MD:Z:22 + + +For more information on SAM, please consult the `SAM format description`__. + +.. __: http://www.ncbi.nlm.nih.gov/pubmed/19505943 + + + </help> +</tool>
participants (1)
-
commits-noreply＠bitbucket.org