--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mismatch_frequencies.py	Wed May 27 13:40:23 2015 -0400
@@ -0,0 +1,300 @@
+import pysam, re, string
+import matplotlib.pyplot as plt
+import pandas as pd
+import json
+from collections import defaultdict
+from collections import OrderedDict
+import argparse
+import itertools
+
+class MismatchFrequencies:
+    '''Iterate over a SAM/BAM alignment file, collecting reads with mismatches. One
+    class instance per alignment file. The result_dict attribute will contain a
+    nested dictionary with name, readlength and mismatch count.'''
+    def __init__(self, result_dict={}, alignment_file=None, name="name", minimal_readlength=21, 
+                 maximal_readlength=21,
+                 number_of_allowed_mismatches=1, 
+                 ignore_5p_nucleotides=0, 
+                 ignore_3p_nucleotides=0,
+                 possible_mismatches = [
+                        'AC', 'AG', 'AT',
+                        'CA', 'CG', 'CT',
+                        'GA', 'GC', 'GT',
+                        'TA', 'TC', 'TG'
+                ]):
+    
+        self.result_dict = result_dict
+        self.name = name
+        self.minimal_readlength = minimal_readlength
+        self.maximal_readlength = maximal_readlength
+        self.number_of_allowed_mismatches = number_of_allowed_mismatches
+        self.ignore_5p_nucleotides = ignore_5p_nucleotides
+        self.ignore_3p_nucleotides = ignore_3p_nucleotides
+        self.possible_mismatches = possible_mismatches
+        
+        if alignment_file:
+            self.pysam_alignment = pysam.Samfile(alignment_file)
+            self.references = self.pysam_alignment.references #names of fasta reference sequences
+            result_dict[name]=self.get_mismatches(
+                self.pysam_alignment, 
+                minimal_readlength, 
+                maximal_readlength,
+                possible_mismatches
+            )
+    
+    def get_mismatches(self, pysam_alignment, minimal_readlength, 
+                       maximal_readlength, possible_mismatches):
+        mismatch_dict = defaultdict(int)
+        rec_dd = lambda: defaultdict(rec_dd)
+        len_dict = rec_dd()
+        for alignedread in pysam_alignment:
+            if self.read_is_valid(alignedread, minimal_readlength, maximal_readlength):
+                chromosome = pysam_alignment.getrname(alignedread.rname)
+                try:
+                    len_dict[int(alignedread.rlen)][chromosome]['total valid reads'] += 1
+                except TypeError:
+                    len_dict[int(alignedread.rlen)][chromosome]['total valid reads'] = 1
+                MD = alignedread.opt('MD')
+                if self.read_has_mismatch(alignedread, self.number_of_allowed_mismatches):
+                    (ref_base, mismatch_base)=self.read_to_reference_mismatch(MD, alignedread.seq, alignedread.is_reverse)
+                    if ref_base == None:
+                            continue
+                    else:
+                        for i, base in enumerate(ref_base):
+                            if not ref_base[i]+mismatch_base[i] in possible_mismatches:
+                                continue
+                            try:
+                                len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] += 1
+                            except TypeError:
+                                len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] = 1
+        return len_dict
+    
+    def read_is_valid(self, read, min_readlength, max_readlength):
+        '''Filter out reads that are unmatched, too short or
+        too long or that contian insertions'''
+        if read.is_unmapped:
+            return False
+        if read.rlen < min_readlength:
+            return False
+        if read.rlen > max_readlength:
+            return False
+        else:
+            return True
+    
+    def read_has_mismatch(self, read, number_of_allowed_mismatches=1):
+        '''keep only reads with one mismatch. Could be simplified'''
+        NM=read.opt('NM')
+        if NM <1: #filter out reads with no mismatch
+            return False
+        if NM >number_of_allowed_mismatches: #filter out reads with more than 1 mismtach
+            return False
+        else:
+            return True
+        
+    def mismatch_in_allowed_region(self, readseq, mismatch_position):
+        '''
+        >>> M = MismatchFrequencies()
+        >>> readseq = 'AAAAAA'
+        >>> mismatch_position = 2
+        >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
+        True
+        >>> M = MismatchFrequencies(ignore_3p_nucleotides=2, ignore_5p_nucleotides=2)
+        >>> readseq = 'AAAAAA'
+        >>> mismatch_position = 1
+        >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
+        False
+        >>> readseq = 'AAAAAA'
+        >>> mismatch_position = 4
+        >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
+        False
+        '''
+        mismatch_position+=1 # To compensate for starting the count at 0
+        five_p = self.ignore_5p_nucleotides
+        three_p = self.ignore_3p_nucleotides
+        if any([five_p > 0, three_p > 0]):
+            if any([mismatch_position <= five_p, 
+                    mismatch_position >= (len(readseq)+1-three_p)]): #Again compensate for starting the count at 0
+                return False
+            else:
+                return True
+        else:
+            return True
+            
+    def read_to_reference_mismatch(self, MD, readseq, is_reverse):
+        '''
+        This is where the magic happens. The MD tag contains SNP and indel information,
+        without looking to the genome sequence. This is a typical MD tag: 3C0G2A6.
+        3 bases of the read align to the reference, followed by a mismatch, where the
+        reference base is C, followed by 10 bases aligned to the reference. 
+        suppose a reference 'CTTCGATAATCCTT'
+                             |||  || ||||||
+                 and a read 'CTTATATTATCCTT'. 
+        This situation is represented by the above MD tag. 
+        Given MD tag and read sequence this function returns the reference base C, G and A, 
+        and the mismatched base A, T, T.
+        >>> M = MismatchFrequencies()
+        >>> MD='3C0G2A7'
+        >>> seq='CTTATATTATCCTT'
+        >>> result=M.read_to_reference_mismatch(MD, seq, is_reverse=False)
+        >>> result[0]=="CGA"
+        True
+        >>> result[1]=="ATT"
+        True
+        >>> 
+        '''
+        search=re.finditer('[ATGC]',MD)
+        if '^' in MD:
+            print 'WARNING insertion detected, mismatch calling skipped for this read!!!'
+            return (None, None)
+        start_index=0 # refers to the leading integer of the MD string before an edited base
+        current_position=0 # position of the mismatched nucleotide in the MD tag string
+        mismatch_position=0 # position of edited base in current read 
+        reference_base=""
+        mismatched_base=""
+        for result in search:
+            current_position=result.start()
+            mismatch_position=mismatch_position+1+int(MD[start_index:current_position]) #converts the leading characters before an edited base into integers
+            start_index=result.end()
+            reference_base+=MD[result.end()-1]
+            mismatched_base+=readseq[mismatch_position-1]
+        if is_reverse:
+            reference_base=reverseComplement(reference_base)
+            mismatched_base=reverseComplement(mismatched_base)
+            mismatch_position=len(readseq)-mismatch_position-1
+        if mismatched_base=='N':
+            return (None, None)
+        if self.mismatch_in_allowed_region(readseq, mismatch_position):
+            return (reference_base, mismatched_base)
+        else:
+            return (None, None)
+
+def reverseComplement(sequence):
+    '''do a reverse complement of DNA base.
+    >>> reverseComplement('ATGC')=='GCAT'
+    True
+    >>> 
+    '''
+    sequence=sequence.upper()
+    complement = string.maketrans('ATCGN', 'TAGCN')
+    return sequence.upper().translate(complement)[::-1]
+
+def barplot(df, library, axes):
+    df.plot(kind='bar', ax=axes, subplots=False,\
+            stacked=False, legend='test',\
+            title='Mismatch frequencies for {0}'.format(library))
+    
+def df_to_tab(df, output):
+    df.to_csv(output, sep='\t')
+
+def reduce_result(df, possible_mismatches):
+    '''takes a pandas dataframe with full mismatch details and
+    summarises the results for plotting.'''
+    alignments = df['Alignment_file'].unique()
+    readlengths = df['Readlength'].unique()
+    combinations = itertools.product(*[alignments, readlengths]) #generate all possible combinations of readlength and alignment files
+    reduced_dict = {}
+    frames = []
+    last_column = 3+len(possible_mismatches)
+    for combination in combinations:
+        library_subset = df[df['Alignment_file'] == combination[0]]
+        library_readlength_subset = library_subset[library_subset['Readlength'] == combination[1]]
+        sum_of_library_and_readlength = library_readlength_subset.iloc[:,3:last_column+1].sum()
+        if not reduced_dict.has_key(combination[0]):
+            reduced_dict[combination[0]] = {}
+        reduced_dict[combination[0]][combination[1]] = sum_of_library_and_readlength.to_dict()
+    return reduced_dict
+
+def plot_result(reduced_dict, args):
+    names=reduced_dict.keys()
+    nrows=len(names)/2+1
+    fig = plt.figure(figsize=(16,32))
+    for i,library in enumerate (names):
+        axes=fig.add_subplot(nrows,2,i+1)
+        library_dict=reduced_dict[library]
+        df=pd.DataFrame(library_dict)
+        df.drop(['total aligned reads'], inplace=True)
+        barplot(df, library, axes),
+        axes.set_ylabel('Mismatch count / all valid reads * readlength')
+    fig.savefig(args.output_pdf, format='pdf')    
+
+def format_result_dict(result_dict, chromosomes, possible_mismatches):
+    '''Turn nested dictionary into preformatted tab seperated lines'''
+    header = "Reference sequence\tAlignment_file\tReadlength\t" + "\t".join(
+        possible_mismatches) + "\ttotal aligned reads"
+    libraries = result_dict.keys()
+    readlengths = result_dict[libraries[0]].keys()
+    result = []
+    for chromosome in chromosomes:
+        for library in libraries:
+            for readlength in readlengths:
+                line = []
+                line.extend([chromosome, library, readlength])
+                try:
+                    line.extend([result_dict[library][readlength][chromosome].get(mismatch, 0) for mismatch in possible_mismatches])
+                    line.extend([result_dict[library][readlength][chromosome].get(u'total valid reads', 0)])
+                except KeyError:
+                    line.extend([0 for mismatch in possible_mismatches])
+                    line.extend([0])
+                result.append(line)
+    df = pd.DataFrame(result, columns=header.split('\t'))
+    last_column=3+len(possible_mismatches)
+    df['mismatches/per aligned nucleotides'] = df.iloc[:,3:last_column].sum(1)/(df.iloc[:,last_column]*df['Readlength'])
+    return df
+  
+def setup_MismatchFrequencies(args):
+    resultDict=OrderedDict()
+    kw_list=[{'result_dict' : resultDict, 
+             'alignment_file' :alignment_file, 
+             'name' : name, 
+             'minimal_readlength' : args.min, 
+             'maximal_readlength' : args.max,
+             'number_of_allowed_mismatches' : args.n_mm,
+             'ignore_5p_nucleotides' : args.five_p, 
+             'ignore_3p_nucleotides' : args.three_p,
+             'possible_mismatches' : args.possible_mismatches }
+             for alignment_file, name in zip(args.input, args.name)]
+    return (kw_list, resultDict)
+
+def nested_dict_to_df(dictionary):
+    dictionary = {(outerKey, innerKey): values for outerKey, innerDict in dictionary.iteritems() for innerKey, values in innerDict.iteritems()}
+    df=pd.DataFrame.from_dict(dictionary).transpose()
+    df.index.names = ['Library', 'Readlength']
+    return df
+
+def run_MismatchFrequencies(args):
+    kw_list, resultDict=setup_MismatchFrequencies(args)
+    references = [MismatchFrequencies(**kw_dict).references for kw_dict in kw_list]
+    return (resultDict, references[0])
+
+def main():
+    result_dict, references = run_MismatchFrequencies(args)
+    df = format_result_dict(result_dict, references, args.possible_mismatches)
+    reduced_dict = reduce_result(df, args.possible_mismatches)
+    plot_result(reduced_dict, args)
+    reduced_df = nested_dict_to_df(reduced_dict)
+    df_to_tab(reduced_df, args.output_tab)
+    if not args.expanded_output_tab == None:
+        df_to_tab(df, args.expanded_output_tab)
+    return reduced_dict
+
+if __name__ == "__main__":
+    
+    parser = argparse.ArgumentParser(description='Produce mismatch statistics for BAM/SAM alignment files.')
+    parser.add_argument('--input', nargs='*', help='Input files in SAM/BAM format')
+    parser.add_argument('--name', nargs='*', help='Name for input file to display in output file. Should have same length as the number of inputs')
+    parser.add_argument('--output_pdf', help='Output filename for graph')
+    parser.add_argument('--output_tab', help='Output filename for table')
+    parser.add_argument('--expanded_output_tab', default=None, help='Output filename for table')
+    parser.add_argument('--possible_mismatches', default=[
+            'AC', 'AG', 'AT','CA', 'CG', 'CT', 'GA', 'GC', 'GT', 'TA', 'TC', 'TG'
+        ], 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.')
+    parser.add_argument('--min', '--minimal_readlength', type=int, help='minimum readlength')
+    parser.add_argument('--max', '--maximal_readlength', type=int, help='maximum readlength')
+    parser.add_argument('--n_mm', '--number_allowed_mismatches', type=int, default=1, help='discard reads with more than n mismatches')
+    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')
+    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')
+    #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'])
+    args = parser.parse_args()
+    reduced_dict = main()
+
+

--- /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>

Binary file test-data/3mismatches_ago2ip_ovary.bam has changed

Binary file test-data/3mismatches_ago2ip_s2.bam has changed

Binary file test-data/mismatch.pdf has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/mismatch.tab	Wed May 27 13:40:23 2015 -0400
@@ -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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tool_dependencies.xml	Wed May 27 13:40:23 2015 -0400
@@ -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>