Mercurial > repos > mvdbeek > mismatch_frequencies

changeset 22:942464ea4211

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
remove heads
author Marius van den Beek <m.vandenbeek@gmail.com>
date Sun, 24 May 2015 17:29:41 +0200
parents 8d604c41010d
children ca7b7890ed20
files 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
diffstat 7 files changed, 65 insertions(+), 204 deletions(-) [+]
line wrap: on
line diff
--- a/mismatch_frequencies.py	Sun May 24 17:26:57 2015 +0200
+++ b/mismatch_frequencies.py	Sun May 24 17:29:41 2015 +0200
@@ -1,27 +1,16 @@
 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'
-                ]):
+    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):
     
         self.result_dict = result_dict
         self.name = name
@@ -30,43 +19,27 @@
         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
-            )
+            result_dict[name]=self.get_mismatches(self.pysam_alignment, minimal_readlength, maximal_readlength)
     
-    def get_mismatches(self, pysam_alignment, minimal_readlength, 
-                       maximal_readlength, possible_mismatches):
+    def get_mismatches(self, pysam_alignment, minimal_readlength, maximal_readlength):
         mismatch_dict = defaultdict(int)
-        rec_dd = lambda: defaultdict(rec_dd)
-        len_dict = rec_dd()
+        len_dict={}
+        for i in range(minimal_readlength, maximal_readlength+1):
+            len_dict[i]=mismatch_dict.copy()
         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')
+                len_dict[int(alignedread.rlen)]['total_mapped'] += 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
+                            len_dict[int(alignedread.rlen)][ref_base[i]+' to '+mismatch_base[i]] += 1
         return len_dict
     
     def read_is_valid(self, read, min_readlength, max_readlength):
@@ -160,7 +133,6 @@
         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):
@@ -168,6 +140,7 @@
         else:
             return (None, None)
 
+
 def reverseComplement(sequence):
     '''do a reverse complement of DNA base.
     >>> reverseComplement('ATGC')=='GCAT'
@@ -181,66 +154,39 @@
 def barplot(df, library, axes):
     df.plot(kind='bar', ax=axes, subplots=False,\
             stacked=False, legend='test',\
-            title='Mismatch frequencies for {0}'.format(library))
-    
+            title='Mismatches in TE small RNAs from {0}'.format(library))
+  
+def result_dict_to_df(result_dict):
+    mismatches = []
+    libraries = []
+    for mismatch, library in result_dict.iteritems():
+        mismatches.append(mismatch)
+        libraries.append(pd.DataFrame.from_dict(library, orient='index'))
+    df=pd.concat(libraries, keys=mismatches)
+    df.index.names = ['library', 'readsize']
+    return df
+
 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()
+def plot_result(result_dict, args):
+    names=args.name
     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]
+        library_dict=result_dict[library]
+        for length in library_dict.keys():
+            for mismatch in library_dict[length]:
+                if mismatch == 'total_mapped':
+                    continue
+                library_dict[length][mismatch]=library_dict[length][mismatch]/float(library_dict[length]['total_mapped'])*100
+            del library_dict[length]['total_mapped']
         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')
+        axes.set_ylabel('Percent of mapped reads with mismatches')
     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, 
@@ -250,32 +196,20 @@
              '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 }
+             'ignore_3p_nucleotides' : args.three_p} 
              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])
+    [MismatchFrequencies(**kw_dict) for kw_dict in kw_list]
+    return resultDict
 
 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
+    result_dict=run_MismatchFrequencies(args)
+    df=result_dict_to_df(result_dict)
+    plot_result(result_dict, args)
+    df_to_tab(df, args.output_tab)
 
 if __name__ == "__main__":
     
@@ -284,17 +218,12 @@
     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(['--input', '3mismatches_ago2ip.bam', '2mismatch.bam', '--name', 'Siomi1', 'Siomi2' , '--five_p', '3','--three_p','3','--output_pdf', 'out.pdf', '--output_tab', 'out.tab', '--min', '21', '--max', '21'])
     args = parser.parse_args()
-    reduced_dict = main()
+    main()
 
-
--- a/mismatch_frequencies.xml	Sun May 24 17:26:57 2015 +0200
+++ b/mismatch_frequencies.xml	Sun May 24 17:29:41 2015 +0200
@@ -1,11 +1,11 @@
-<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 
+<tool id="mismatch_frequencies" name="Mismatch Frequencies" version="0.0.3" 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">pandas</requirement>
+    	    <requirement type="package" version="1.4">matplotlib</requirement>
+        </requirements>
+	<command interpreter="python">mismatch_frequencies.py --input 
 		#for i in $rep
 			"$i.input_file" 
 		#end for
@@ -17,73 +17,20 @@
                  --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="fasta" 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>
+        </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="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"/>
+	</inputs>
+        <outputs>
+                <data format="pdf" name="output_pdf" />
+                <data format="tabular" name="output_tab" />
+        </outputs>
 
-.. 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
--- a/test-data/mismatch.tab	Sun May 24 17:26:57 2015 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,3 +0,0 @@
-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
--- a/tool_dependencies.xml	Sun May 24 17:26:57 2015 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-<?xml version="1.0"?>
-<tool_dependency>
-    <package name="pysam" version="0.7.7">
-        <repository changeset_revision="ca10c522f37e" name="package_pysam_0_7_7" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
-    </package>
-    <package name="pandas" version="0.14.1">
-        <repository changeset_revision="e27b7cf19fef" name="package_pandas_0_14" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
-    </package>
-    <package name="matplotlib" version="1.2.1">
-        <repository changeset_revision="dc29f1f60887" name="package_matplotlib_1_2" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
-    </package>
-</tool_dependency>