| Next changeset 1:a006d42dd759 (2016-12-06) |
|
Commit message:
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/lumpy commit c0bfc4b2215705e1b5fd1d4e60b1d72e5da13c92 |
|
added:
extractSplitReads_BwaMem.py lumpy.xml pairend_distro.py test-data/output.vcf test-data/sr.input.bam |
| b |
| 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 <file> + +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() + 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()) |
| b |
| 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 @@ +<tool id="lumpy" name="lumpy-sv" version="0.1"> + <description>find structural variants</description> + <requirements> + <requirement type="package" version="0.2.12">lumpy-sv</requirement> + <requirement type="package" version="1.3.1">samtools</requirement> + <requirement type="package" version="1.11.2">numpy</requirement> + </requirements> + <stdio> + <exit_code range="1:" level="fatal" description="Tool exception" /> + </stdio> + <version_command>lumpy --version</version_command> + <command><![CDATA[ + ln -f -s "$input_file" input.bam && + #if $seq_method.seq_method_list == "paired-end": + samtools view -b -F 1294 input.bam > "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 + + ]]></command> + <!-- basic error handling --> + <inputs> + <param format="bam" name="input_file" type="data" label="One BAM alignment file produced by BWA-mem"/> + <conditional name="seq_method"> + <param help="Paired-end or single-read sequencing" label="Sequencing method" name="seq_method_list" type="select"> + <option selected="True" value="paired-end">Paired-end sequencing</option> + <option value="single-read">Single-read sequencing</option> + </param> + <when value="paired-end"> + <param name="readLength" value="151" type="integer" label="read length" help="e.g. 151 nt" /> + <section name="additional_params" title="Additional Options" expanded="False"> + <param name="samplingValue" value="100000" type="integer" label="number of reads to compute mean and stdev of read length" help="e.g. 10000" /> + <param name="min_non_overlap" value="101" type="integer" label="min_non_overlap" help="e.g. 101" /> + <param name="discordant_z" value="5" type="integer" label="discordant_z" help="e.g. 5" /> + <param name="back_distance" value="10" type="integer" label="back_distance" help="e.g. 10" /> + <param name="weight" value="1" type="integer" label="weight" help="e.g. 1" /> + <param name="min_mapping_threshold" value="20" type="integer" label="min_mapping_threshold" help="e.g. 20" /> + </section> + </when> + <when value="single-read"> + <section name="additional_params" title="Additional Options" expanded="False"> + <param name="back_distance" value="10" type="integer" label="back_distance" help="e.g. 10" /> + <param name="weight" value="1" type="integer" label="weight" help="e.g. 1" /> + <param name="min_mapping_threshold" value="20" type="integer" label="min_mapping_threshold" help="e.g. 20" /> + </section> + </when> + + </conditional> + + </inputs> + + <outputs> + <data format="tabular" name="histogram" type="data" label="${input_file.element_identifier} Fragment size distribution"> + <filter>seq_method['seq_method_list'] == "paired-end"</filter> + </data> + <data format="bam" name="splits" type="data" label="${input_file.element_identifier} Split Reads (Bam format)"/> + <data format="bam" name="discordants" type="data" label="${input_file.element_identifier} Discordant Pairs (Bam format)"> + <filter>seq_method['seq_method_list'] == "paired-end"</filter> + </data> + <data format="vcf" name="vcf_call" type="data" label="${input_file.element_identifier} Variant Calling (vcf format)"/> + </outputs> + + <tests> + <test> + <param name="input_file" value="sr.input.bam" ftype="bam"/> + <param name="seq_method_list" value="single-read" /> + <param name="back_distance" value="10"/> + <param name="weight" value="1" /> + <param name="min_mapping_threshold" value="20" /> + <output name="vcf_call" file="output.vcf" ftype="vcf"/> + </test> + </tests> + + <help> + Some help required + </help> + + <citations> + <citation type="doi">10.1186/gb-2014-15-6-r84</citation> + </citations> +</tool> |
| b |
| 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)) |
| b |
| 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=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant"> +##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles"> +##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> +##INFO=<ID=STRANDS,Number=.,Type=String,Description="Strand orientation of the adjacency in BEDPE format (DEL:+-, DUP:-+, INV:++/--)"> +##INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description="Imprecise structural variation"> +##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants"> +##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants"> +##INFO=<ID=CIPOS95,Number=2,Type=Integer,Description="Confidence interval (95%) around POS for imprecise variants"> +##INFO=<ID=CIEND95,Number=2,Type=Integer,Description="Confidence interval (95%) around END for imprecise variants"> +##INFO=<ID=MATEID,Number=.,Type=String,Description="ID of mate breakends"> +##INFO=<ID=EVENT,Number=1,Type=String,Description="ID of event associated to breakend"> +##INFO=<ID=SECONDARY,Number=0,Type=Flag,Description="Secondary breakend in a multi-line variants"> +##INFO=<ID=SU,Number=.,Type=Integer,Description="Number of pieces of evidence supporting the variant across all samples"> +##INFO=<ID=PE,Number=.,Type=Integer,Description="Number of paired-end reads supporting the variant across all samples"> +##INFO=<ID=SR,Number=.,Type=Integer,Description="Number of split reads supporting the variant across all samples"> +##INFO=<ID=BD,Number=.,Type=Integer,Description="Amount of BED evidence supporting the variant across all samples"> +##INFO=<ID=EV,Number=.,Type=String,Description="Type of LUMPY evidence contributing to the variant call"> +##INFO=<ID=PRPOS,Number=.,Type=String,Description="LUMPY probability curve of the POS breakend"> +##INFO=<ID=PREND,Number=.,Type=String,Description="LUMPY probability curve of the END breakend"> +##ALT=<ID=DEL,Description="Deletion"> +##ALT=<ID=DUP,Description="Duplication"> +##ALT=<ID=INV,Description="Inversion"> +##ALT=<ID=DUP:TANDEM,Description="Tandem duplication"> +##ALT=<ID=INS,Description="Insertion of novel sequence"> +##ALT=<ID=CNV,Description="Copy number variable region"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=SU,Number=1,Type=Integer,Description="Number of pieces of evidence supporting the variant"> +##FORMAT=<ID=PE,Number=1,Type=Integer,Description="Number of paired-end reads supporting the variant"> +##FORMAT=<ID=SR,Number=1,Type=Integer,Description="Number of split reads supporting the variant"> +##FORMAT=<ID=BD,Number=1,Type=Integer,Description="Amount of BED evidence supporting the variant"> +#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 |
| b |
| diff -r 000000000000 -r 8b3daa745d9b test-data/sr.input.bam |
| b |
| Binary file test-data/sr.input.bam has changed |