# HG changeset patch # User drosofff # Date 1481021188 18000 # Node ID 8b3daa745d9bcc00299a4be1ba270081f73b6974 planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/lumpy commit c0bfc4b2215705e1b5fd1d4e60b1d72e5da13c92 diff -r 000000000000 -r 8b3daa745d9b extractSplitReads_BwaMem.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extractSplitReads_BwaMem.py Tue Dec 06 05:46:28 2016 -0500 @@ -0,0 +1,202 @@ +#!/usr/bin/env python + +import sys +import getopt +import string +from optparse import OptionParser +import re + +def extractSplitsFromBwaMem(inFile,numSplits,includeDups,minNonOverlap): + if inFile == "stdin": + data = sys.stdin + else: + data = open(inFile, 'r') + for line in data: + split = 0 + if line[0] == '@': + print line.strip() + continue + samList = line.strip().split('\t') + sam = SAM(samList) + if includeDups==0 and (1024 & sam.flag)==1024: + continue + for el in sam.tags: + if "SA:" in el: + if(len(el.split(";")))<=numSplits: + split = 1 + mate = el.split(",") + mateCigar = mate[3] + mateFlag = int(0) + if mate[2]=="-": mateFlag = int(16) + if split: + read1 = sam.flag & 64 + if read1 == 64: tag = "_1" + else: tag="_2" + samList[0] = sam.query + tag + readCigar = sam.cigar + readCigarOps = extractCigarOps(readCigar,sam.flag) + readQueryPos = calcQueryPosFromCigar(readCigarOps) + mateCigarOps = extractCigarOps(mateCigar,mateFlag) + mateQueryPos = calcQueryPosFromCigar(mateCigarOps) + overlap = calcQueryOverlap(readQueryPos.qsPos,readQueryPos.qePos,mateQueryPos.qsPos,mateQueryPos.qePos) + nonOverlap1 = 1 + readQueryPos.qePos - readQueryPos.qsPos - overlap + nonOverlap2 = 1 + mateQueryPos.qePos - mateQueryPos.qsPos - overlap + mno = min(nonOverlap1, nonOverlap2) + if mno >= minNonOverlap: + print "\t".join(samList) + +#-------------------------------------------------------------------------------------------------- +# functions +#-------------------------------------------------------------------------------------------------- + +class SAM (object): + """ + __very__ basic class for SAM input. + """ + def __init__(self, samList = []): + if len(samList) > 0: + self.query = samList[0] + self.flag = int(samList[1]) + self.ref = samList[2] + self.pos = int(samList[3]) + self.mapq = int(samList[4]) + self.cigar = samList[5] + self.matRef = samList[6] + self.matePos = int(samList[7]) + self.iSize = int(samList[8]) + self.seq = samList[9] + self.qual = samList[10] + self.tags = samList[11:]#tags is a list of each tag:vtype:value sets + self.valid = 1 + else: + self.valid = 0 + self.query = 'null' + + def extractTagValue (self, tagID): + for tag in self.tags: + tagParts = tag.split(':', 2); + if (tagParts[0] == tagID): + if (tagParts[1] == 'i'): + return int(tagParts[2]); + elif (tagParts[1] == 'H'): + return int(tagParts[2],16); + return tagParts[2]; + return None; + +#----------------------------------------------- +cigarPattern = '([0-9]+[MIDNSHP])' +cigarSearch = re.compile(cigarPattern) +atomicCigarPattern = '([0-9]+)([MIDNSHP])' +atomicCigarSearch = re.compile(atomicCigarPattern) + +def extractCigarOps(cigar,flag): + if (cigar == "*"): + cigarOps = [] + elif (flag & 0x0010): + cigarOpStrings = cigarSearch.findall(cigar) + cigarOps = [] + for opString in cigarOpStrings: + cigarOpList = atomicCigarSearch.findall(opString) +# print cigarOpList + # "struct" for the op and it's length + cigar = cigarOp(cigarOpList[0][0], cigarOpList[0][1]) + # add to the list of cigarOps + cigarOps.append(cigar) + cigarOps = cigarOps + cigarOps.reverse() + ##do in reverse order because negative strand## + else: + cigarOpStrings = cigarSearch.findall(cigar) + cigarOps = [] + for opString in cigarOpStrings: + cigarOpList = atomicCigarSearch.findall(opString) + # "struct" for the op and it's length + cigar = cigarOp(cigarOpList[0][0], cigarOpList[0][1]) + # add to the list of cigarOps + cigarOps.append(cigar) +# cigarOps = cigarOps + return(cigarOps) + +def calcQueryPosFromCigar(cigarOps): + qsPos = 0 + qePos = 0 + qLen = 0 + # if first op is a H, need to shift start position + # the opPosition counter sees if the for loop is looking at the first index of the cigar object + opPosition = 0 + for cigar in cigarOps: + if opPosition == 0 and (cigar.op == 'H' or cigar.op == 'S'): + qsPos += cigar.length + qePos += cigar.length + qLen += cigar.length + elif opPosition > 0 and (cigar.op == 'H' or cigar.op == 'S'): + qLen += cigar.length + elif cigar.op == 'M' or cigar.op == 'I': + qePos += cigar.length + qLen += cigar.length + opPosition += 1 + d = queryPos(qsPos, qePos, qLen); + return d + +class cigarOp (object): + """ + sturct to store a discrete CIGAR operations + """ + def __init__(self, opLength, op): + self.length = int(opLength) + self.op = op + +class queryPos (object): + """ + struct to store the start and end positions of query CIGAR operations + """ + def __init__(self, qsPos, qePos, qLen): + self.qsPos = int(qsPos) + self.qePos = int(qePos) + self.qLen = int(qLen) + + +def calcQueryOverlap(s1,e1,s2,e2): + o = 1 + min(e1, e2) - max(s1, s2) + return max(0, o) + +############################################### + +class Usage(Exception): + def __init__(self, msg): + self.msg = msg + +def main(): + + usage = """%prog -i + +extractSplitReads_BwaMem v0.1.0 +Author: Ira Hall +Description: Get split-read alignments from bwa-mem in lumpy compatible format. Ignores reads marked as duplicates. +Works on read or position sorted SAM input. Tested on bwa mem v0.7.5a-r405. + """ + parser = OptionParser(usage) + + parser.add_option("-i", "--inFile", dest="inFile", + help="A SAM file or standard input (-i stdin).", + metavar="FILE") + parser.add_option("-n", "--numSplits", dest="numSplits", default=2, type = "int", + help="The maximum number of split-read mappings to allow per read. Reads with more are excluded. Default=2", + metavar="INT") + parser.add_option("-d", "--includeDups", dest="includeDups", action="store_true",default=0, + help="Include alignments marked as duplicates. Default=False") + parser.add_option("-m", "--minNonOverlap", dest="minNonOverlap", default=20, type = "int", + help="minimum non-overlap between split alignments on the query (default=20)", + metavar="INT") + (opts, args) = parser.parse_args() + if opts.inFile is None: + parser.print_help() + print + else: + try: + extractSplitsFromBwaMem(opts.inFile, opts.numSplits, opts.includeDups, opts.minNonOverlap) + except IOError as err: + sys.stderr.write("IOError " + str(err) + "\n"); + return +if __name__ == "__main__": + sys.exit(main()) diff -r 000000000000 -r 8b3daa745d9b lumpy.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lumpy.xml Tue Dec 06 05:46:28 2016 -0500 @@ -0,0 +1,102 @@ + + find structural variants + + lumpy-sv + samtools + numpy + + + + + lumpy --version + "input.discordants.unsorted.bam" && + samtools view -h input.bam | python $__tool_directory__/extractSplitReads_BwaMem.py -i stdin | samtools view -Sb - > "input.splitters.unsorted.bam" && + samtools sort input.discordants.unsorted.bam > input.discordants.bam && + samtools sort input.splitters.unsorted.bam > input.splitters.bam && + samtools view -r readgroup input.bam + |tail -n +$seq_method.additional_params.samplingValue + |python $__tool_directory__/pairend_distro.py -r 101 -X 4 -N $seq_method.additional_params.samplingValue -o input.lib.histo > meandev.txt && + mean=\$(cat meandev.txt | sed s/mean:// | sed -r s/stdev:.+//) && + stdev=\$(cat meandev.txt | sed -r s/mean:.+stdev://) && + lumpy -mw 4 -tt 0 + -pe id:input.bam,bam_file:input.discordants.bam,histo_file:input.lib.histo,mean:"\$mean",stdev:"\$stdev",read_length:$seq_method.readLength,min_non_overlap:$seq_method.additional_params.min_non_overlap,discordant_z:$seq_method.additional_params.discordant_z,back_distance:$seq_method.additional_params.back_distance,weight:$seq_method.additional_params.weight,min_mapping_threshold:$seq_method.additional_params.min_mapping_threshold + -sr id:input.bam,bam_file:input.splitters.bam,back_distance:$seq_method.additional_params.back_distance,weight:$seq_method.additional_params.weight,min_mapping_threshold:$seq_method.additional_params.min_mapping_threshold > output.vcf && + mv input.discordants.bam $discordants && + mv input.splitters.bam $splits && + mv input.lib.histo $histogram && + mv output.vcf $vcf_call && + rm input.discordants.unsorted.bam input.splitters.unsorted.bam meandev.txt + #end if + #if $seq_method.seq_method_list == "single-read": + samtools view -h input.bam | python $__tool_directory__/extractSplitReads_BwaMem.py -i stdin | samtools view -Sb - > "input.splitters.unsorted.bam" && + lumpy -mw 4 -tt 0 + -sr id:input.bam,bam_file:input.splitters.unsorted.bam,back_distance:$seq_method.additional_params.back_distance,weight:$seq_method.additional_params.weight,min_mapping_threshold:$seq_method.additional_params.min_mapping_threshold > output.vcf && + mv input.splitters.unsorted.bam $splits && + mv output.vcf $vcf_call + #end if + + ]]> + + + + + + + + + + +
+ + + + + + +
+ + +
+ + + +
+ + + + + + + + + seq_method['seq_method_list'] == "paired-end" + + + + seq_method['seq_method_list'] == "paired-end" + + + + + + + + + + + + + + + + + Some help required + + + + 10.1186/gb-2014-15-6-r84 + + diff -r 000000000000 -r 8b3daa745d9b pairend_distro.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pairend_distro.py Tue Dec 06 05:46:28 2016 -0500 @@ -0,0 +1,140 @@ +#!/usr/bin/env python +# (c) 2012 - Ryan M. Layer +# Hall Laboratory +# Quinlan Laboratory +# Department of Computer Science +# Department of Biochemistry and Molecular Genetics +# Department of Public Health Sciences and Center for Public Health Genomics, +# University of Virginia +# rl6sf@virginia.edu + +import sys +import numpy as np +from operator import itemgetter +from optparse import OptionParser + +# some constants for sam/bam field ids +SAM_FLAG = 1 +SAM_REFNAME = 2 +SAM_MATE_REFNAME = 6 +SAM_ISIZE = 8 + +parser = OptionParser() + +parser.add_option("-r", + "--read_length", + type="int", + dest="read_length", + help="Read length") + +parser.add_option("-X", + dest="X", + type="int", + help="Number of stdevs from mean to extend") + +parser.add_option("-N", + dest="N", + type="int", + help="Number to sample") + +parser.add_option("-o", + dest="output_file", + help="Output file") + +parser.add_option("-m", + dest="mads", + type="int", + default=10, + help="Outlier cutoff in # of median absolute deviations (unscaled, upper only)") + +def unscaled_upper_mad(xs): + """Return a tuple consisting of the median of xs followed by the + unscaled median absolute deviation of the values in xs that lie + above the median. + """ + med = np.median(xs) + return med, np.median(xs[xs > med] - med) + + +(options, args) = parser.parse_args() + +if not options.read_length: + parser.error('Read length not given') + +if not options.X: + parser.error('X not given') + +if not options.N: + parser.error('N not given') + +if not options.output_file: + parser.error('Output file not given') + + +required = 97 +restricted = 3484 +flag_mask = required | restricted + +L = [] +c = 0 + +for l in sys.stdin: + if c >= options.N: + break + + A = l.rstrip().split('\t') + flag = int(A[SAM_FLAG]) + refname = A[SAM_REFNAME] + mate_refname = A[SAM_MATE_REFNAME] + isize = int(A[SAM_ISIZE]) + + want = mate_refname == "=" and flag & flag_mask == required and isize >= 0 + if want: + c += 1 + L.append(isize) + +# warn if very few elements in distribution +min_elements = 1000 +if len(L) < min_elements: + sys.stderr.write("Warning: only %s elements in distribution (min: %s)\n" % (len(L), min_elements)) + mean = "NA" + stdev = "NA" + +else: + # Remove outliers + L = np.array(L) + L.sort() + med, umad = unscaled_upper_mad(L) + upper_cutoff = med + options.mads * umad + L = L[L < upper_cutoff] + new_len = len(L) + removed = c - new_len + sys.stderr.write("Removed %d outliers with isize >= %d\n" % + (removed, upper_cutoff)) + c = new_len + + mean = np.mean(L) + stdev = np.std(L) + + start = options.read_length + end = int(mean + options.X*stdev) + + H = [0] * (end - start + 1) + s = 0 + + for x in L: + if (x >= start) and (x <= end): + j = int(x - start) + H[j] = H[ int(x - start) ] + 1 + s += 1 + + f = open(options.output_file, 'w') + + for i in range(end - start): + o = str(i) + "\t" + str(float(H[i])/float(s)) + "\n" + f.write(o) + + + f.close() + +print('mean:' + str(mean) + '\tstdev:' + str(stdev)) diff -r 000000000000 -r 8b3daa745d9b test-data/output.vcf --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.vcf Tue Dec 06 05:46:28 2016 -0500 @@ -0,0 +1,37 @@ +##fileformat=VCFv4.2 +##source=LUMPY +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##INFO= +##ALT= +##ALT= +##ALT= +##ALT= +##ALT= +##ALT= +##FORMAT= +##FORMAT= +##FORMAT= +##FORMAT= +##FORMAT= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT input.bam +hg38_gold_U07000.1 14 1_1 N [hg38_gold_U07000.1:1876[N . . SVTYPE=BND;STRANDS=--:19;EVENT=1;MATEID=1_2;CIPOS=0,0;CIEND=0,2;CIPOS95=0,0;CIEND95=0,0;SU=19;SR=19 GT:SU:SR ./.:19:19 +hg38_gold_U07000.1 1876 1_2 N [hg38_gold_U07000.1:14[N . . SVTYPE=BND;STRANDS=--:19;SECONDARY;EVENT=1;MATEID=1_1;CIPOS=0,2;CIEND=0,0;CIPOS95=0,0;CIEND95=0,0;SU=19;SR=19 GT:SU:SR ./.:19:19 +hg38_gold_U07000.1 10 2_1 N [hg38_gold_U07000.1:1897[N . . SVTYPE=BND;STRANDS=--:19;EVENT=2;MATEID=2_2;CIPOS=-1,0;CIEND=-7,5;CIPOS95=0,1;CIEND95=-2,1;IMPRECISE;SU=19;SR=19 GT:SU:SR ./.:19:19 +hg38_gold_U07000.1 1897 2_2 N [hg38_gold_U07000.1:10[N . . SVTYPE=BND;STRANDS=--:19;SECONDARY;EVENT=2;MATEID=2_1;CIPOS=-7,5;CIEND=-1,0;CIPOS95=-2,1;CIEND95=0,1;IMPRECISE;SU=19;SR=19 GT:SU:SR ./.:19:19 diff -r 000000000000 -r 8b3daa745d9b test-data/sr.input.bam Binary file test-data/sr.input.bam has changed