mismatch_frequencies: mismatch_frequencies.py annotate

annotate mismatch_frequencies.py @ 26:18772b1eb16a draft

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

author	mvdbeek
date	Mon, 29 Jun 2015 05:52:19 -0400
parents	6590be3f8e3f
children	bb3d962e8cb1

rev	line source
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	1 import pysam, re, string
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	2 import matplotlib.pyplot as plt
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	3 import pandas as pd
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	4 import json
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	5 from collections import defaultdict
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	6 from collections import OrderedDict
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	7 import argparse
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	8 import itertools
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	9
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	10 class MismatchFrequencies:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	11 '''Iterate over a SAM/BAM alignment file, collecting reads with mismatches. One
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	12 class instance per alignment file. The result_dict attribute will contain a
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	13 nested dictionary with name, readlength and mismatch count.'''
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	14 def __init__(self, result_dict={}, alignment_file=None, name="name", minimal_readlength=21,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	15 maximal_readlength=21,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	16 number_of_allowed_mismatches=1,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	17 ignore_5p_nucleotides=0,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	18 ignore_3p_nucleotides=0,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	19 possible_mismatches = [
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	20 'AC', 'AG', 'AT',
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	21 'CA', 'CG', 'CT',
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	22 'GA', 'GC', 'GT',
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	23 'TA', 'TC', 'TG'
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	24 ]):
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	25
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	26 self.result_dict = result_dict
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	27 self.name = name
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	28 self.minimal_readlength = minimal_readlength
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	29 self.maximal_readlength = maximal_readlength
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	30 self.number_of_allowed_mismatches = number_of_allowed_mismatches
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	31 self.ignore_5p_nucleotides = ignore_5p_nucleotides
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	32 self.ignore_3p_nucleotides = ignore_3p_nucleotides
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	33 self.possible_mismatches = possible_mismatches
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	34
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	35 if alignment_file:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	36 self.pysam_alignment = pysam.Samfile(alignment_file)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	37 self.references = self.pysam_alignment.references #names of fasta reference sequences
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	38 result_dict[name]=self.get_mismatches(
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	39 self.pysam_alignment,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	40 minimal_readlength,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	41 maximal_readlength,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	42 possible_mismatches
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	43 )
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	44
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	45 def get_mismatches(self, pysam_alignment, minimal_readlength,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	46 maximal_readlength, possible_mismatches):
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	47 mismatch_dict = defaultdict(int)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	48 rec_dd = lambda: defaultdict(rec_dd)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	49 len_dict = rec_dd()
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	50 for alignedread in pysam_alignment:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	51 if self.read_is_valid(alignedread, minimal_readlength, maximal_readlength):
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	52 chromosome = pysam_alignment.getrname(alignedread.rname)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	53 try:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	54 len_dict[int(alignedread.rlen)][chromosome]['total valid reads'] += 1
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	55 except TypeError:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	56 len_dict[int(alignedread.rlen)][chromosome]['total valid reads'] = 1
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	57 MD = alignedread.opt('MD')
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	58 if self.read_has_mismatch(alignedread, self.number_of_allowed_mismatches):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	59 (ref_base, mismatch_base)=self.read_to_reference_mismatch(MD, alignedread.seq, alignedread.is_reverse)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	60 if ref_base == None:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	61 continue
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	62 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	63 for i, base in enumerate(ref_base):
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	64 if not ref_base[i]+mismatch_base[i] in possible_mismatches:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	65 continue
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	66 try:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	67 len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] += 1
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	68 except TypeError:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	69 len_dict[int(alignedread.rlen)][chromosome][ref_base[i]+mismatch_base[i]] = 1
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	70 return len_dict
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	71
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	72 def read_is_valid(self, read, min_readlength, max_readlength):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	73 '''Filter out reads that are unmatched, too short or
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	74 too long or that contian insertions'''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	75 if read.is_unmapped:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	76 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	77 if read.rlen < min_readlength:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	78 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	79 if read.rlen > max_readlength:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	80 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	81 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	82 return True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	83
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	84 def read_has_mismatch(self, read, number_of_allowed_mismatches=1):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	85 '''keep only reads with one mismatch. Could be simplified'''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	86 NM=read.opt('NM')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	87 if NM <1: #filter out reads with no mismatch
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	88 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	89 if NM >number_of_allowed_mismatches: #filter out reads with more than 1 mismtach
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	90 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	91 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	92 return True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	93
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	94 def mismatch_in_allowed_region(self, readseq, mismatch_position):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	95 '''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	96 >>> M = MismatchFrequencies()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	97 >>> readseq = 'AAAAAA'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	98 >>> mismatch_position = 2
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	99 >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	100 True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	101 >>> M = MismatchFrequencies(ignore_3p_nucleotides=2, ignore_5p_nucleotides=2)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	102 >>> readseq = 'AAAAAA'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	103 >>> mismatch_position = 1
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	104 >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	105 False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	106 >>> readseq = 'AAAAAA'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	107 >>> mismatch_position = 4
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	108 >>> M.mismatch_in_allowed_region(readseq, mismatch_position)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	109 False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	110 '''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	111 mismatch_position+=1 # To compensate for starting the count at 0
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	112 five_p = self.ignore_5p_nucleotides
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	113 three_p = self.ignore_3p_nucleotides
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	114 if any([five_p > 0, three_p > 0]):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	115 if any([mismatch_position <= five_p,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	116 mismatch_position >= (len(readseq)+1-three_p)]): #Again compensate for starting the count at 0
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	117 return False
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	118 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	119 return True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	120 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	121 return True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	122
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	123 def read_to_reference_mismatch(self, MD, readseq, is_reverse):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	124 '''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	125 This is where the magic happens. The MD tag contains SNP and indel information,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	126 without looking to the genome sequence. This is a typical MD tag: 3C0G2A6.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	127 3 bases of the read align to the reference, followed by a mismatch, where the
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	128 reference base is C, followed by 10 bases aligned to the reference.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	129 suppose a reference 'CTTCGATAATCCTT'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	130 \|\|\| \|\| \|\|\|\|\|\|
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	131 and a read 'CTTATATTATCCTT'.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	132 This situation is represented by the above MD tag.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	133 Given MD tag and read sequence this function returns the reference base C, G and A,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	134 and the mismatched base A, T, T.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	135 >>> M = MismatchFrequencies()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	136 >>> MD='3C0G2A7'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	137 >>> seq='CTTATATTATCCTT'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	138 >>> result=M.read_to_reference_mismatch(MD, seq, is_reverse=False)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	139 >>> result[0]=="CGA"
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	140 True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	141 >>> result[1]=="ATT"
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	142 True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	143 >>>
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	144 '''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	145 search=re.finditer('[ATGC]',MD)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	146 if '^' in MD:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	147 print 'WARNING insertion detected, mismatch calling skipped for this read!!!'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	148 return (None, None)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	149 start_index=0 # refers to the leading integer of the MD string before an edited base
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	150 current_position=0 # position of the mismatched nucleotide in the MD tag string
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	151 mismatch_position=0 # position of edited base in current read
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	152 reference_base=""
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	153 mismatched_base=""
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	154 for result in search:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	155 current_position=result.start()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	156 mismatch_position=mismatch_position+1+int(MD[start_index:current_position]) #converts the leading characters before an edited base into integers
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	157 start_index=result.end()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	158 reference_base+=MD[result.end()-1]
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	159 mismatched_base+=readseq[mismatch_position-1]
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	160 if is_reverse:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	161 reference_base=reverseComplement(reference_base)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	162 mismatched_base=reverseComplement(mismatched_base)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	163 mismatch_position=len(readseq)-mismatch_position-1
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	164 if mismatched_base=='N':
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	165 return (None, None)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	166 if self.mismatch_in_allowed_region(readseq, mismatch_position):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	167 return (reference_base, mismatched_base)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	168 else:
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	169 return (None, None)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	170
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	171 def reverseComplement(sequence):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	172 '''do a reverse complement of DNA base.
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	173 >>> reverseComplement('ATGC')=='GCAT'
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	174 True
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	175 >>>
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	176 '''
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	177 sequence=sequence.upper()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	178 complement = string.maketrans('ATCGN', 'TAGCN')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	179 return sequence.upper().translate(complement)[::-1]
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	180
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	181 def barplot(df, library, axes):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	182 df.plot(kind='bar', ax=axes, subplots=False,\
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	183 stacked=False, legend='test',\
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	184 title='Mismatch frequencies for {0}'.format(library))
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	185
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	186 def df_to_tab(df, output):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	187 df.to_csv(output, sep='\t')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	188
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	189 def reduce_result(df, possible_mismatches):
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	190 '''takes a pandas dataframe with full mismatch details and
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	191 summarises the results for plotting.'''
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	192 alignments = df['Alignment_file'].unique()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	193 readlengths = df['Readlength'].unique()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	194 combinations = itertools.product(*[alignments, readlengths]) #generate all possible combinations of readlength and alignment files
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	195 reduced_dict = {}
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	196 frames = []
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	197 last_column = 3+len(possible_mismatches)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	198 for combination in combinations:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	199 library_subset = df[df['Alignment_file'] == combination[0]]
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	200 library_readlength_subset = library_subset[library_subset['Readlength'] == combination[1]]
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	201 sum_of_library_and_readlength = library_readlength_subset.iloc[:,3:last_column+1].sum()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	202 if not reduced_dict.has_key(combination[0]):
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	203 reduced_dict[combination[0]] = {}
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	204 reduced_dict[combination[0]][combination[1]] = sum_of_library_and_readlength.to_dict()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	205 return reduced_dict
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	206
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	207 def plot_result(reduced_dict, args):
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	208 names=reduced_dict.keys()
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	209 nrows=len(names)/2+1
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	210 fig = plt.figure(figsize=(16,32))
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	211 for i,library in enumerate (names):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	212 axes=fig.add_subplot(nrows,2,i+1)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	213 library_dict=reduced_dict[library]
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	214 df=pd.DataFrame(library_dict)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	215 df.drop(['total aligned reads'], inplace=True)
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	216 barplot(df, library, axes),
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	217 axes.set_ylabel('Mismatch count / all valid reads * readlength')
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	218 fig.savefig(args.output_pdf, format='pdf')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	219
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	220 def format_result_dict(result_dict, chromosomes, possible_mismatches):
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	221 '''Turn nested dictionary into preformatted tab seperated lines'''
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	222 header = "Reference sequence\tAlignment_file\tReadlength\t" + "\t".join(
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	223 possible_mismatches) + "\ttotal aligned reads"
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	224 libraries = result_dict.keys()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	225 readlengths = result_dict[libraries[0]].keys()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	226 result = []
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	227 for chromosome in chromosomes:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	228 for library in libraries:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	229 for readlength in readlengths:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	230 line = []
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	231 line.extend([chromosome, library, readlength])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	232 try:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	233 line.extend([result_dict[library][readlength][chromosome].get(mismatch, 0) for mismatch in possible_mismatches])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	234 line.extend([result_dict[library][readlength][chromosome].get(u'total valid reads', 0)])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	235 except KeyError:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	236 line.extend([0 for mismatch in possible_mismatches])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	237 line.extend([0])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	238 result.append(line)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	239 df = pd.DataFrame(result, columns=header.split('\t'))
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	240 last_column=3+len(possible_mismatches)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	241 df['mismatches/per aligned nucleotides'] = df.iloc[:,3:last_column].sum(1)/(df.iloc[:,last_column]*df['Readlength'])
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	242 return df
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	243
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	244 def setup_MismatchFrequencies(args):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	245 resultDict=OrderedDict()
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	246 kw_list=[{'result_dict' : resultDict,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	247 'alignment_file' :alignment_file,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	248 'name' : name,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	249 'minimal_readlength' : args.min,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	250 'maximal_readlength' : args.max,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	251 'number_of_allowed_mismatches' : args.n_mm,
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	252 'ignore_5p_nucleotides' : args.five_p,
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	253 'ignore_3p_nucleotides' : args.three_p,
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	254 'possible_mismatches' : args.possible_mismatches }
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	255 for alignment_file, name in zip(args.input, args.name)]
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	256 return (kw_list, resultDict)
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	257
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	258 def nested_dict_to_df(dictionary):
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	259 dictionary = {(outerKey, innerKey): values for outerKey, innerDict in dictionary.iteritems() for innerKey, values in innerDict.iteritems()}
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	260 df=pd.DataFrame.from_dict(dictionary).transpose()
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	261 df.index.names = ['Library', 'Readlength']
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	262 return df
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	263
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	264 def run_MismatchFrequencies(args):
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	265 kw_list, resultDict=setup_MismatchFrequencies(args)
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	266 references = [MismatchFrequencies(**kw_dict).references for kw_dict in kw_list]
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	267 return (resultDict, references[0])
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	268
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	269 def main():
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	270 result_dict, references = run_MismatchFrequencies(args)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	271 df = format_result_dict(result_dict, references, args.possible_mismatches)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	272 reduced_dict = reduce_result(df, args.possible_mismatches)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	273 plot_result(reduced_dict, args)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	274 reduced_df = nested_dict_to_df(reduced_dict)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	275 df_to_tab(reduced_df, args.output_tab)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	276 if not args.expanded_output_tab == None:
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	277 df_to_tab(df, args.expanded_output_tab)
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	278 return reduced_dict
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	279
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	280 if __name__ == "__main__":
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	281
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	282 parser = argparse.ArgumentParser(description='Produce mismatch statistics for BAM/SAM alignment files.')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	283 parser.add_argument('--input', nargs='*', help='Input files in SAM/BAM format')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	284 parser.add_argument('--name', nargs='*', help='Name for input file to display in output file. Should have same length as the number of inputs')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	285 parser.add_argument('--output_pdf', help='Output filename for graph')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	286 parser.add_argument('--output_tab', help='Output filename for table')
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	287 parser.add_argument('--expanded_output_tab', default=None, help='Output filename for table')
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	288 parser.add_argument('--possible_mismatches', default=[
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	289 'AC', 'AG', 'AT','CA', 'CG', 'CT', 'GA', 'GC', 'GT', 'TA', 'TC', 'TG'
6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	290 ], 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.')
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	291 parser.add_argument('--min', '--minimal_readlength', type=int, help='minimum readlength')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	292 parser.add_argument('--max', '--maximal_readlength', type=int, help='maximum readlength')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	293 parser.add_argument('--n_mm', '--number_allowed_mismatches', type=int, default=1, help='discard reads with more than n mismatches')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	294 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')
848d799e6fe8 Uploaded mvdbeek parents: diff changeset	295 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')
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	296 #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'])
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	297 args = parser.parse_args()
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	298 reduced_dict = main()
12 848d799e6fe8 Uploaded mvdbeek parents: diff changeset	299
24 6590be3f8e3f planemo upload for repository https://bitbucket.org/drosofff/gedtools/ mvdbeek parents: 22 diff changeset	300

Mercurial > repos > mvdbeek > mismatch_frequencies

annotate mismatch_frequencies.py @ 26:18772b1eb16a draft