| Next changeset 1:3613460e891e (2016-03-23) |
|
Commit message:
planemo upload for repository https://bitbucket.org/drosofff/gedtools/ |
|
added:
mismatch_frequencies.py mismatch_frequencies.xml test-data/3mismatches_ago2ip_ovary.bam test-data/3mismatches_ago2ip_s2.bam test-data/mismatch.pdf test-data/mismatch.tab tool_dependencies.xml |
| b |
| diff -r 000000000000 -r 77de5fc623f9 mismatch_frequencies.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies.py Wed May 27 13:40:23 2015 -0400 |
| [ |
| b'@@ -0,0 +1,300 @@\n+import pysam, re, string\n+import matplotlib.pyplot as plt\n+import pandas as pd\n+import json\n+from collections import defaultdict\n+from collections import OrderedDict\n+import argparse\n+import itertools\n+\n+class MismatchFrequencies:\n+ \'\'\'Iterate over a SAM/BAM alignment file, collecting reads with mismatches. One\n+ class instance per alignment file. The result_dict attribute will contain a\n+ nested dictionary with name, readlength and mismatch count.\'\'\'\n+ def __init__(self, result_dict={}, alignment_file=None, name="name", minimal_readlength=21, \n+ maximal_readlength=21,\n+ number_of_allowed_mismatches=1, \n+ ignore_5p_nucleotides=0, \n+ ignore_3p_nucleotides=0,\n+ possible_mismatches = [\n+ \'AC\', \'AG\', \'AT\',\n+ \'CA\', \'CG\', \'CT\',\n+ \'GA\', \'GC\', \'GT\',\n+ \'TA\', \'TC\', \'TG\'\n+ ]):\n+ \n+ self.result_dict = result_dict\n+ self.name = name\n+ self.minimal_readlength = minimal_readlength\n+ self.maximal_readlength = maximal_readlength\n+ self.number_of_allowed_mismatches = number_of_allowed_mismatches\n+ self.ignore_5p_nucleotides = ignore_5p_nucleotides\n+ self.ignore_3p_nucleotides = ignore_3p_nucleotides\n+ self.possible_mismatches = possible_mismatches\n+ \n+ if alignment_file:\n+ self.pysam_alignment = pysam.Samfile(alignment_file)\n+ self.references = self.pysam_alignment.references #names of fasta reference sequences\n+ result_dict[name]=self.get_mismatches(\n+ self.pysam_alignment, \n+ minimal_readlength, \n+ maximal_readlength,\n+ possible_mismatches\n+ )\n+ \n+ def get_mismatches(self, pysam_alignment, minimal_readlength, \n+ maximal_readlength, possible_mismatches):\n+ mismatch_dict = defaultdict(int)\n+ rec_dd = lambda: defaultdict(rec_dd)\n+ len_dict = rec_dd()\n+ for alignedread in pysam_alignment:\n+ if self.read_is_valid(alignedread, minimal_readlength, maximal_readlength):\n+ chromosome = pysam_alignment.getrname(alignedread.rname)\n+ try:\n+ len_dict[int(alignedread.rlen)][chromosome][\'total valid reads\'] += 1\n+ except TypeError:\n+ len_dict[int(alignedread.rlen)][chromosome][\'total valid reads\'] = 1\n+ MD = alignedread.opt(\'MD\')\n+ if self.read_has_mismatch(alignedread, self.number_of_allowed_mismatches):\n+ (ref_base, mismatch_base)=self.read_to_reference_mismatch(MD, alignedread.seq, alignedread.is_reverse)\n+ if ref_base == None:\n+ continue\n+ else:\n+ for i, base in enumerate(ref_base):\n+ if not ref_base[i]+mismatch_base[i] in possible_mismatches:\n+ continue\n+ try:\n+ len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] += 1\n+ except TypeError:\n+ len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] = 1\n+ return len_dict\n+ \n+ def read_is_valid(self, read, min_readlength, max_readlength):\n+ \'\'\'Filter out reads that are unmatched, too short or\n+ too long or that contian insertions\'\'\'\n+ if read.is_unmapped:\n+ return False\n+ if read.rlen < min_readlength:\n+ return False\n+ if read.rlen > max_readlength:\n+ return False\n+ else:\n+ return True\n+ \n+ def read_has_mismatch(self, read, number_of_allowed_mismatches=1):\n+ \'\'\'keep only reads with one mismatch. Could be simplified\'\'\'\n+'..b'e_mismatches])\n+ line.extend([0])\n+ result.append(line)\n+ df = pd.DataFrame(result, columns=header.split(\'\\t\'))\n+ last_column=3+len(possible_mismatches)\n+ df[\'mismatches/per aligned nucleotides\'] = df.iloc[:,3:last_column].sum(1)/(df.iloc[:,last_column]*df[\'Readlength\'])\n+ return df\n+ \n+def setup_MismatchFrequencies(args):\n+ resultDict=OrderedDict()\n+ kw_list=[{\'result_dict\' : resultDict, \n+ \'alignment_file\' :alignment_file, \n+ \'name\' : name, \n+ \'minimal_readlength\' : args.min, \n+ \'maximal_readlength\' : args.max,\n+ \'number_of_allowed_mismatches\' : args.n_mm,\n+ \'ignore_5p_nucleotides\' : args.five_p, \n+ \'ignore_3p_nucleotides\' : args.three_p,\n+ \'possible_mismatches\' : args.possible_mismatches }\n+ for alignment_file, name in zip(args.input, args.name)]\n+ return (kw_list, resultDict)\n+\n+def nested_dict_to_df(dictionary):\n+ dictionary = {(outerKey, innerKey): values for outerKey, innerDict in dictionary.iteritems() for innerKey, values in innerDict.iteritems()}\n+ df=pd.DataFrame.from_dict(dictionary).transpose()\n+ df.index.names = [\'Library\', \'Readlength\']\n+ return df\n+\n+def run_MismatchFrequencies(args):\n+ kw_list, resultDict=setup_MismatchFrequencies(args)\n+ references = [MismatchFrequencies(**kw_dict).references for kw_dict in kw_list]\n+ return (resultDict, references[0])\n+\n+def main():\n+ result_dict, references = run_MismatchFrequencies(args)\n+ df = format_result_dict(result_dict, references, args.possible_mismatches)\n+ reduced_dict = reduce_result(df, args.possible_mismatches)\n+ plot_result(reduced_dict, args)\n+ reduced_df = nested_dict_to_df(reduced_dict)\n+ df_to_tab(reduced_df, args.output_tab)\n+ if not args.expanded_output_tab == None:\n+ df_to_tab(df, args.expanded_output_tab)\n+ return reduced_dict\n+\n+if __name__ == "__main__":\n+ \n+ parser = argparse.ArgumentParser(description=\'Produce mismatch statistics for BAM/SAM alignment files.\')\n+ parser.add_argument(\'--input\', nargs=\'*\', help=\'Input files in SAM/BAM format\')\n+ parser.add_argument(\'--name\', nargs=\'*\', help=\'Name for input file to display in output file. Should have same length as the number of inputs\')\n+ parser.add_argument(\'--output_pdf\', help=\'Output filename for graph\')\n+ parser.add_argument(\'--output_tab\', help=\'Output filename for table\')\n+ parser.add_argument(\'--expanded_output_tab\', default=None, help=\'Output filename for table\')\n+ parser.add_argument(\'--possible_mismatches\', default=[\n+ \'AC\', \'AG\', \'AT\',\'CA\', \'CG\', \'CT\', \'GA\', \'GC\', \'GT\', \'TA\', \'TC\', \'TG\'\n+ ], nargs=\'+\', help=\'specify mismatches that should be counted for the mismatch frequency. The format is Reference base -> observed base, eg AG for A to G mismatches.\')\n+ parser.add_argument(\'--min\', \'--minimal_readlength\', type=int, help=\'minimum readlength\')\n+ parser.add_argument(\'--max\', \'--maximal_readlength\', type=int, help=\'maximum readlength\')\n+ parser.add_argument(\'--n_mm\', \'--number_allowed_mismatches\', type=int, default=1, help=\'discard reads with more than n mismatches\')\n+ parser.add_argument(\'--five_p\', \'--ignore_5p_nucleotides\', type=int, default=0, help=\'when calculating nucleotide mismatch frequencies ignore the first N nucleotides of the read\')\n+ parser.add_argument(\'--three_p\', \'--ignore_3p_nucleotides\', type=int, default=1, help=\'when calculating nucleotide mismatch frequencies ignore the last N nucleotides of the read\')\n+ #args = parser.parse_args([\'--input\', \'3mismatches_ago2ip_s2.bam\', \'3mismatches_ago2ip_ovary.bam\',\'--possible_mismatches\',\'AC\',\'AG\', \'CG\', \'TG\', \'CT\',\'--name\', \'Siomi1\', \'Siomi2\' , \'--five_p\', \'3\',\'--three_p\',\'3\',\'--output_pdf\', \'out.pdf\', \'--output_tab\', \'out.tab\', \'--expanded_output_tab\', \'expanded.tab\', \'--min\', \'20\', \'--max\', \'22\'])\n+ args = parser.parse_args()\n+ reduced_dict = main()\n+\n+\n' |
| b |
| diff -r 000000000000 -r 77de5fc623f9 mismatch_frequencies.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies.xml Wed May 27 13:40:23 2015 -0400 |
| [ |
| @@ -0,0 +1,89 @@ +<tool id="mismatch_frequencies" name="Mismatch Frequencies" version="0.0.9" hidden="false" > + <description>Analyze mismatch frequencies in BAM/SAM alignments</description> + <requirements> + <requirement type="package" version="0.7.7">pysam</requirement> + <requirement type="package" version="0.14.1">pandas</requirement> + <requirement type="package" version="1.2.1">matplotlib</requirement> + </requirements> + <command interpreter="python">mismatch_frequencies.py --input + #for i in $rep + "$i.input_file" + #end for + --name + #for i in $rep + "$i.input_file.name" + #end for + --output_pdf $output_pdf --output_tab $output_tab --min $min_length --max $max_length + --n_mm $number_of_mismatches + --five_p $five_p + --three_p $three_p + --expanded_output_tab $expanded_tab + --possible_mismatches $possible_mismatches + </command> + <inputs> + <repeat name="rep" title="alignment files"> + <param name="input_file" type="data" format="bam,sam" label="Alignment file" help="The input alignment file(s) for which to analyze the mismatches."/> + </repeat> + <param name="number_of_mismatches" label="Maximum number of allowed mismatches per read" help="Discard reads with more than the chosen number of mismatches from the frequency calculation" type="integer" value="3"/> + <param name="possible_mismatches" label="Specify mismatches that should be counted" help="Ignores mismatches that are not listed" type="text" value="AC AG AT CA CG CT GA GC GT TA TC TG"> + <validator type="expression" message="Allowed values are AGCTN, seperated by space.">len([False for char in value if not char in " AGCTN"]) == 0</validator> + </param> + <param name="min_length" label="Minumum read length to analyse" type="integer" value="21"/> + <param name="max_length" label="Maximum read length to analyse" type="integer" value="21"/> + <param name="five_p" label="Ignore mismatches in the first N nucleotides of a read" type="integer" value="0"/> + <param name="three_p" label="Ignore mismatches in the last N nucleotides of a read" help="useful to discriminate between tailing events and editing events" type="integer" value="3"/> + <param help="Output expanded tabular format" label="Nucleotide mismatches per reference sequence" name="expanded" type="select"> + <option select="true" value="false">No</option> + <option value="expanded">Yes</option> + </param> + </inputs> + <outputs> + <data format="tabular" name="output_tab" /> + <data format="tabular" name="expanded_tab"> + <filter> expanded == "expanded"</filter> + </data> + <data format="pdf" name="output_pdf" /> + </outputs> + <tests> + <test> + <param name="rep_0|input_file" value="3mismatches_ago2ip_s2.bam" ftype="bam" /> + <param name="rep_1|input_file" value="3mismatches_ago2ip_ovary.bam" ftype="bam" /> + <param name="number_of_mismatches" value="1" /> + <param name="min_length" value="21" /> + <param name="max_length" value="21" /> + <param name="three_p" value="0" /> + <param name="five_p" value="0" /> + <output name="tabular" file="mismatch.tab" ftype="tabular"/> + <!-- + <output name="pdf" file="mismatch.pdf" ftype="pdf"/> + --> + </test> + </tests> + <help> + +.. class:: infomark + + +***What it does*** + +This tool reconstitues for each aligned read of an alignment file in SAM/BAM format whether +a mismatch is annotated in the MD tag, and if that is the case counts the identity of the +mismatch relative to the reference sequence. The output is a PDF document with the calculated +frequency for each mismatch that occured relative to the total number of valid reads and a table +with the corresponding values. Read length can be limited to a specific read length, and 5 prime and +3 prime-most nucleotides of a read can be ignored. + +---- + +.. class:: warningmark + +***Warning*** + +This tool skips all read that have insertions and has been tested only with bowtie and bowtie2 +generated alignment files. + +Written by Marius van den Beek, m.vandenbeek at gmail . com + </help> + <citations> + </citations> +</tool> |
| b |
| diff -r 000000000000 -r 77de5fc623f9 test-data/3mismatches_ago2ip_ovary.bam |
| b |
| Binary file test-data/3mismatches_ago2ip_ovary.bam has changed |
| b |
| diff -r 000000000000 -r 77de5fc623f9 test-data/3mismatches_ago2ip_s2.bam |
| b |
| Binary file test-data/3mismatches_ago2ip_s2.bam has changed |
| b |
| diff -r 000000000000 -r 77de5fc623f9 test-data/mismatch.pdf |
| b |
| Binary file test-data/mismatch.pdf has changed |
| b |
| diff -r 000000000000 -r 77de5fc623f9 test-data/mismatch.tab --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/mismatch.tab Wed May 27 13:40:23 2015 -0400 |
| b |
| @@ -0,0 +1,3 @@ +Library Readlength AC AG AT CA CG CT GA GC GT TA TC TG total aligned reads +3mismatches_ago2ip_ovary.bam 21 380 1214 524 581 278 1127 1032 239 595 483 973 394 138649 +3mismatches_ago2ip_s2.bam 21 48 6503 106 68 46 173 222 144 220 90 232 40 43881 |
| b |
| diff -r 000000000000 -r 77de5fc623f9 tool_dependencies.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Wed May 27 13:40:23 2015 -0400 |
| b |
| @@ -0,0 +1,12 @@ +<?xml version="1.0"?> +<tool_dependency> + <package name="pysam" version="0.7.7"> + <repository changeset_revision="b62538c8c664" name="package_pysam_0_7_7" owner="iuc" toolshed="https://toolshed.g2.bx.psu.edu" /> + </package> + <package name="pandas" version="0.14.1"> + <repository changeset_revision="18e65cea168d" name="package_pandas_0_14" owner="iuc" toolshed="https://toolshed.g2.bx.psu.edu" /> + </package> + <package name="matplotlib" version="1.2.1"> + <repository changeset_revision="a03ee94316b5" name="package_matplotlib_1_2" owner="iuc" toolshed="https://toolshed.g2.bx.psu.edu" /> + </package> +</tool_dependency> |