# HG changeset patch
# User mvdbeek
# Date 1432481791 14400
# Node ID 6590be3f8e3f77a0b406a38a2bb3156050b1130c
# Parent ca7b7890ed20b19583a90f9d1ca49a0434e41d5d
planemo upload for repository https://bitbucket.org/drosofff/gedtools/
diff -r ca7b7890ed20 -r 6590be3f8e3f mismatch_frequencies.py
--- a/mismatch_frequencies.py Sun May 24 17:33:33 2015 +0200
+++ b/mismatch_frequencies.py Sun May 24 11:36:31 2015 -0400
@@ -1,16 +1,27 @@
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):
+ 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
@@ -19,27 +30,43 @@
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(
+ self.pysam_alignment,
+ minimal_readlength,
+ maximal_readlength,
+ possible_mismatches
+ )
- def get_mismatches(self, pysam_alignment, minimal_readlength, maximal_readlength):
+ def get_mismatches(self, pysam_alignment, minimal_readlength,
+ maximal_readlength, possible_mismatches):
mismatch_dict = defaultdict(int)
- len_dict={}
- for i in range(minimal_readlength, maximal_readlength+1):
- len_dict[i]=mismatch_dict.copy()
+ 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):
- len_dict[int(alignedread.rlen)]['total_mapped'] += 1
- MD=alignedread.opt('MD')
+ 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):
- 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):
@@ -133,6 +160,7 @@
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):
@@ -140,7 +168,6 @@
else:
return (None, None)
-
def reverseComplement(sequence):
'''do a reverse complement of DNA base.
>>> reverseComplement('ATGC')=='GCAT'
@@ -154,39 +181,66 @@
def barplot(df, library, axes):
df.plot(kind='bar', ax=axes, subplots=False,\
stacked=False, legend='test',\
- 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
-
+ title='Mismatch frequencies for {0}'.format(library))
+
def df_to_tab(df, output):
df.to_csv(output, sep='\t')
-def plot_result(result_dict, args):
- names=args.name
+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=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']
+ library_dict=reduced_dict[library]
df=pd.DataFrame(library_dict)
+ df.drop(['total aligned reads'], inplace=True)
barplot(df, library, axes),
- axes.set_ylabel('Percent of mapped reads with mismatches')
+ 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,
@@ -196,20 +250,32 @@
'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]
- return resultDict
+ references = [MismatchFrequencies(**kw_dict).references for kw_dict in kw_list]
+ return (resultDict, references[0])
def main():
- result_dict=run_MismatchFrequencies(args)
- df=result_dict_to_df(result_dict)
- plot_result(result_dict, args)
- df_to_tab(df, args.output_tab)
+ 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__":
@@ -218,12 +284,17 @@
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()
+
diff -r ca7b7890ed20 -r 6590be3f8e3f mismatch_frequencies.xml
--- a/mismatch_frequencies.xml Sun May 24 17:33:33 2015 +0200
+++ b/mismatch_frequencies.xml Sun May 24 11:36:31 2015 -0400
@@ -1,11 +1,11 @@
-
- Analyze mismatch frequencies in BAM/SAM alignments
-
- pysam
- pandas
- matplotlib
-
- mismatch_frequencies.py --input
+
+ Analyze mismatch frequencies in BAM/SAM alignments
+
+ pysam
+ pandas
+ matplotlib
+
+ mismatch_frequencies.py --input
#for i in $rep
"$i.input_file"
#end for
@@ -17,20 +17,73 @@
--n_mm $number_of_mismatches
--five_p $five_p
--three_p $three_p
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+ --expanded_output_tab $expanded_tab
+ --possible_mismatches $possible_mismatches
+
+
+
+
+
+
+
+ len([False for char in value if not char in " AGCTN"]) == 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+ expanded == "expanded"
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+.. 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
+
+
+
diff -r ca7b7890ed20 -r 6590be3f8e3f test-data/3mismatches_ago2ip_ovary.bam
Binary file test-data/3mismatches_ago2ip_ovary.bam has changed
diff -r ca7b7890ed20 -r 6590be3f8e3f test-data/3mismatches_ago2ip_s2.bam
Binary file test-data/3mismatches_ago2ip_s2.bam has changed
diff -r ca7b7890ed20 -r 6590be3f8e3f test-data/mismatch.pdf
Binary file test-data/mismatch.pdf has changed
diff -r ca7b7890ed20 -r 6590be3f8e3f test-data/mismatch.tab
--- a/test-data/mismatch.tab Sun May 24 17:33:33 2015 +0200
+++ b/test-data/mismatch.tab Sun May 24 11:36:31 2015 -0400
@@ -1,3 +1,3 @@
-library readsize A to C A to G A to T C to A C to G C to T G to A G to C G to T T to A T to C T to G total_mapped
-3mismatches_ago2ip_s2.bam 21 31 5484 69 25 40 137 156 109 188 51 196 29 43881
-3mismatches_ago2ip_ovary.bam 21 293 879 411 452 231 872 845 191 473 384 818 324 138649
+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 ca7b7890ed20 -r 6590be3f8e3f tool_dependencies.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tool_dependencies.xml Sun May 24 11:36:31 2015 -0400
@@ -0,0 +1,12 @@
+
+
+
+
+
+
+
+
+
+
+
+