# HG changeset patch # User Marius van den Beek # Date 1432481381 -7200 # Node ID 942464ea4211777675148e8a5886f20524bb8e90 # Parent 8d604c41010dd7f212560869fa615ec6b38e2e72 remove heads diff -r 8d604c41010d -r 942464ea4211 mismatch_frequencies.py --- 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() - diff -r 8d604c41010d -r 942464ea4211 mismatch_frequencies.xml --- 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 @@ - diff -r 8d604c41010d -r 942464ea4211 test-data/3mismatches_ago2ip_ovary.bam Binary file test-data/3mismatches_ago2ip_ovary.bam has changed diff -r 8d604c41010d -r 942464ea4211 test-data/3mismatches_ago2ip_s2.bam Binary file test-data/3mismatches_ago2ip_s2.bam has changed diff -r 8d604c41010d -r 942464ea4211 test-data/mismatch.pdf Binary file test-data/mismatch.pdf has changed diff -r 8d604c41010d -r 942464ea4211 test-data/mismatch.tab --- 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 diff -r 8d604c41010d -r 942464ea4211 tool_dependencies.xml --- 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 @@ - - - - - - - - - - - -