mismatch_frequencies: mismatch_frequencies.py comparison

comparison mismatch_frequencies.py @ 19:f7da7f3e2c98

planemo upload for repository https://bitbucket.org/drosofff/gedtools/

author	mvdbeek
date	Sun, 24 May 2015 11:12:24 -0400
parents	848d799e6fe8
children	942464ea4211

comparison

equal deleted inserted replaced

-:e3d950e52e38
+:f7da7f3e2c98
 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,
+def __init__(self, result_dict={}, alignment_file=None, name="name", minimal_readlength=21,
-number_of_allowed_mismatches=1, ignore_5p_nucleotides=0, ignore_3p_nucleotides=0):
+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)
-result_dict[name]=self.get_mismatches(self.pysam_alignment, minimal_readlength, maximal_readlength)
+self.references = self.pysam_alignment.references #names of fasta reference sequences
+result_dict[name]=self.get_mismatches(
-def get_mismatches(self, pysam_alignment, minimal_readlength, maximal_readlength):
+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)
-len_dict={}
+rec_dd = lambda: defaultdict(rec_dd)
-for i in range(minimal_readlength, maximal_readlength+1):
+len_dict = rec_dd()
-len_dict[i]=mismatch_dict.copy()
 for alignedread in pysam_alignment:
 if self.read_is_valid(alignedread, minimal_readlength, maximal_readlength):
-len_dict[int(alignedread.rlen)]['total valid reads'] += 1
+chromosome = pysam_alignment.getrname(alignedread.rname)
-MD=alignedread.opt('MD')
+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):
-len_dict[int(alignedread.rlen)][ref_base[i]+' to '+mismatch_base[i]] += 1
+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'''
 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 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 plot_result(result_dict, args):
+def reduce_result(df, possible_mismatches):
-names=args.name
+'''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=result_dict[library]
+library_dict=reduced_dict[library]
-for length in library_dict.keys():
-for mismatch in library_dict[length]:
-if mismatch == 'total valid reads':
-continue
-library_dict[length][mismatch]=library_dict[length][mismatch]/float(library_dict[length]['total valid reads'])*100
-del library_dict[length]['total valid reads']
 df=pd.DataFrame(library_dict)
+df.drop(['total aligned reads'], inplace=True)
 barplot(df, library, axes),
-axes.set_ylabel('Mismatch count / all valid reads * 100')
+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}
+'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)
-[MismatchFrequencies(**kw_dict) for kw_dict in kw_list]
+references = [MismatchFrequencies(**kw_dict).references for kw_dict in kw_list]
-return resultDict
+return (resultDict, references[0])
 def main():
-result_dict=run_MismatchFrequencies(args)
+result_dict, references = run_MismatchFrequencies(args)
-df=result_dict_to_df(result_dict)
+df = format_result_dict(result_dict, references, args.possible_mismatches)
-plot_result(result_dict, args)
+reduced_dict = reduce_result(df, args.possible_mismatches)
-df_to_tab(df, args.output_tab)
+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.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(['--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()
-main()
+reduced_dict = main()

Mercurial > repos > mvdbeek > mismatch_frequencies

comparison mismatch_frequencies.py @ 19:f7da7f3e2c98