Mercurial > repos > mvdbeek > mismatch_frequencies
changeset 8:2e041a1564ad
Uploaded
| author | mvdbeek |
|---|---|
| date | Mon, 02 Feb 2015 10:33:21 -0500 |
| parents | 270681625775 |
| children | 9e4311a3270e |
| files | mismatch_frequencies.py mismatch_frequencies.xml mismatch_frequencies/.hg_archival.txt mismatch_frequencies/mismatch_frequencies.py mismatch_frequencies/mismatch_frequencies.xml mismatch_frequencies/test-data/3mismatches_ago2ip_ovary.bam mismatch_frequencies/test-data/3mismatches_ago2ip_s2.bam mismatch_frequencies/test-data/mismatch.pdf mismatch_frequencies/test-data/mismatch.tab mismatch_frequencies/tool_dependencies.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 | 15 files changed, 319 insertions(+), 314 deletions(-) [+] |
line wrap: on
line diff
--- a/mismatch_frequencies.py Wed Jan 28 13:12:56 2015 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,229 +0,0 @@ -import pysam, re, string -import matplotlib.pyplot as plt -import pandas as pd -from collections import defaultdict -from collections import OrderedDict -import argparse - -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): - - 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 - - if alignment_file: - self.pysam_alignment = pysam.Samfile(alignment_file) - result_dict[name]=self.get_mismatches(self.pysam_alignment, minimal_readlength, maximal_readlength) - - def get_mismatches(self, pysam_alignment, minimal_readlength, maximal_readlength): - mismatch_dict = defaultdict(int) - 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): - len_dict[int(alignedread.rlen)]['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 - 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''' - if read.is_unmapped: - return False - if read.rlen < min_readlength: - return False - if read.rlen > max_readlength: - return False - else: - return True - - def read_has_mismatch(self, read, number_of_allowed_mismatches=1): - '''keep only reads with one mismatch. Could be simplified''' - NM=read.opt('NM') - if NM <1: #filter out reads with no mismatch - return False - if NM >number_of_allowed_mismatches: #filter out reads with more than 1 mismtach - return False - else: - return True - - def mismatch_in_allowed_region(self, readseq, mismatch_position): - ''' - >>> M = MismatchFrequencies() - >>> readseq = 'AAAAAA' - >>> mismatch_position = 2 - >>> M.mismatch_in_allowed_region(readseq, mismatch_position) - True - >>> M = MismatchFrequencies(ignore_3p_nucleotides=2, ignore_5p_nucleotides=2) - >>> readseq = 'AAAAAA' - >>> mismatch_position = 1 - >>> M.mismatch_in_allowed_region(readseq, mismatch_position) - False - >>> readseq = 'AAAAAA' - >>> mismatch_position = 4 - >>> M.mismatch_in_allowed_region(readseq, mismatch_position) - False - ''' - mismatch_position+=1 # To compensate for starting the count at 0 - five_p = self.ignore_5p_nucleotides - three_p = self.ignore_3p_nucleotides - if any([five_p > 0, three_p > 0]): - if any([mismatch_position <= five_p, - mismatch_position >= (len(readseq)+1-three_p)]): #Again compensate for starting the count at 0 - return False - else: - return True - else: - return True - - def read_to_reference_mismatch(self, MD, readseq, is_reverse): - ''' - This is where the magic happens. The MD tag contains SNP and indel information, - without looking to the genome sequence. This is a typical MD tag: 3C0G2A6. - 3 bases of the read align to the reference, followed by a mismatch, where the - reference base is C, followed by 10 bases aligned to the reference. - suppose a reference 'CTTCGATAATCCTT' - ||| || |||||| - and a read 'CTTATATTATCCTT'. - This situation is represented by the above MD tag. - Given MD tag and read sequence this function returns the reference base C, G and A, - and the mismatched base A, T, T. - >>> M = MismatchFrequencies() - >>> MD='3C0G2A7' - >>> seq='CTTATATTATCCTT' - >>> result=M.read_to_reference_mismatch(MD, seq, is_reverse=False) - >>> result[0]=="CGA" - True - >>> result[1]=="ATT" - True - >>> - ''' - search=re.finditer('[ATGC]',MD) - if '^' in MD: - print 'WARNING insertion detected, mismatch calling skipped for this read!!!' - return (None, None) - start_index=0 # refers to the leading integer of the MD string before an edited base - current_position=0 # position of the mismatched nucleotide in the MD tag string - mismatch_position=0 # position of edited base in current read - reference_base="" - mismatched_base="" - for result in search: - current_position=result.start() - mismatch_position=mismatch_position+1+int(MD[start_index:current_position]) #converts the leading characters before an edited base into integers - start_index=result.end() - reference_base+=MD[result.end()-1] - mismatched_base+=readseq[mismatch_position-1] - 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): - return (reference_base, mismatched_base) - else: - return (None, None) - -def reverseComplement(sequence): - '''do a reverse complement of DNA base. - >>> reverseComplement('ATGC')=='GCAT' - True - >>> - ''' - sequence=sequence.upper() - complement = string.maketrans('ATCGN', 'TAGCN') - return sequence.upper().translate(complement)[::-1] - -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): - 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=result_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) - barplot(df, library, axes), - axes.set_ylabel('Mismatch count / all valid reads * 100') - fig.savefig(args.output_pdf, format='pdf') - -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} - for alignment_file, name in zip(args.input, args.name)] - return (kw_list, resultDict) - -def run_MismatchFrequencies(args): - kw_list, resultDict=setup_MismatchFrequencies(args) - [MismatchFrequencies(**kw_dict) for kw_dict in kw_list] - return resultDict - -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) - -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('--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() - main() -
--- a/mismatch_frequencies.xml Wed Jan 28 13:12:56 2015 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -<tool id="mismatch_frequencies" name="Mismatch Frequencies" version="0.0.5" 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 - --name - #for i in $rep - "$i.input_file.name" - #end for - --output_pdf $output_pdf --output_tab $output_tab --min $min_length --max $max_length - --n_mm $number_of_mismatches - --five_p $five_p - --three_p $three_p - </command> - <inputs> - <repeat name="rep" title="alignment files" min="1"> - <param name="input_file" type="data" format="bam,sam" label="Alignment file" help="The input alignment file(s) for which you want to calculate mismatch frequencies."/> - </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="tabular" name="output_tab" /> - <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" /> - <output name="tabular" file="mismatch.tab" ftype="tabular"/> - </test> - </tests> - <help> - -.. 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> -</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies/.hg_archival.txt Mon Feb 02 10:33:21 2015 -0500 @@ -0,0 +1,5 @@ +repo: e8ebe5132737d76a94ec92eda7d1fc1059bffb88 +node: 27068162577556368c6e40a43e87bd71be98a713 +branch: default +latesttag: null +latesttagdistance: 6
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies/mismatch_frequencies.py Mon Feb 02 10:33:21 2015 -0500 @@ -0,0 +1,229 @@ +import pysam, re, string +import matplotlib.pyplot as plt +import pandas as pd +from collections import defaultdict +from collections import OrderedDict +import argparse + +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): + + 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 + + if alignment_file: + self.pysam_alignment = pysam.Samfile(alignment_file) + result_dict[name]=self.get_mismatches(self.pysam_alignment, minimal_readlength, maximal_readlength) + + def get_mismatches(self, pysam_alignment, minimal_readlength, maximal_readlength): + mismatch_dict = defaultdict(int) + 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): + len_dict[int(alignedread.rlen)]['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 + 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''' + if read.is_unmapped: + return False + if read.rlen < min_readlength: + return False + if read.rlen > max_readlength: + return False + else: + return True + + def read_has_mismatch(self, read, number_of_allowed_mismatches=1): + '''keep only reads with one mismatch. Could be simplified''' + NM=read.opt('NM') + if NM <1: #filter out reads with no mismatch + return False + if NM >number_of_allowed_mismatches: #filter out reads with more than 1 mismtach + return False + else: + return True + + def mismatch_in_allowed_region(self, readseq, mismatch_position): + ''' + >>> M = MismatchFrequencies() + >>> readseq = 'AAAAAA' + >>> mismatch_position = 2 + >>> M.mismatch_in_allowed_region(readseq, mismatch_position) + True + >>> M = MismatchFrequencies(ignore_3p_nucleotides=2, ignore_5p_nucleotides=2) + >>> readseq = 'AAAAAA' + >>> mismatch_position = 1 + >>> M.mismatch_in_allowed_region(readseq, mismatch_position) + False + >>> readseq = 'AAAAAA' + >>> mismatch_position = 4 + >>> M.mismatch_in_allowed_region(readseq, mismatch_position) + False + ''' + mismatch_position+=1 # To compensate for starting the count at 0 + five_p = self.ignore_5p_nucleotides + three_p = self.ignore_3p_nucleotides + if any([five_p > 0, three_p > 0]): + if any([mismatch_position <= five_p, + mismatch_position >= (len(readseq)+1-three_p)]): #Again compensate for starting the count at 0 + return False + else: + return True + else: + return True + + def read_to_reference_mismatch(self, MD, readseq, is_reverse): + ''' + This is where the magic happens. The MD tag contains SNP and indel information, + without looking to the genome sequence. This is a typical MD tag: 3C0G2A6. + 3 bases of the read align to the reference, followed by a mismatch, where the + reference base is C, followed by 10 bases aligned to the reference. + suppose a reference 'CTTCGATAATCCTT' + ||| || |||||| + and a read 'CTTATATTATCCTT'. + This situation is represented by the above MD tag. + Given MD tag and read sequence this function returns the reference base C, G and A, + and the mismatched base A, T, T. + >>> M = MismatchFrequencies() + >>> MD='3C0G2A7' + >>> seq='CTTATATTATCCTT' + >>> result=M.read_to_reference_mismatch(MD, seq, is_reverse=False) + >>> result[0]=="CGA" + True + >>> result[1]=="ATT" + True + >>> + ''' + search=re.finditer('[ATGC]',MD) + if '^' in MD: + print 'WARNING insertion detected, mismatch calling skipped for this read!!!' + return (None, None) + start_index=0 # refers to the leading integer of the MD string before an edited base + current_position=0 # position of the mismatched nucleotide in the MD tag string + mismatch_position=0 # position of edited base in current read + reference_base="" + mismatched_base="" + for result in search: + current_position=result.start() + mismatch_position=mismatch_position+1+int(MD[start_index:current_position]) #converts the leading characters before an edited base into integers + start_index=result.end() + reference_base+=MD[result.end()-1] + mismatched_base+=readseq[mismatch_position-1] + 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): + return (reference_base, mismatched_base) + else: + return (None, None) + +def reverseComplement(sequence): + '''do a reverse complement of DNA base. + >>> reverseComplement('ATGC')=='GCAT' + True + >>> + ''' + sequence=sequence.upper() + complement = string.maketrans('ATCGN', 'TAGCN') + return sequence.upper().translate(complement)[::-1] + +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): + 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=result_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) + barplot(df, library, axes), + axes.set_ylabel('Mismatch count / all valid reads * 100') + fig.savefig(args.output_pdf, format='pdf') + +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} + for alignment_file, name in zip(args.input, args.name)] + return (kw_list, resultDict) + +def run_MismatchFrequencies(args): + kw_list, resultDict=setup_MismatchFrequencies(args) + [MismatchFrequencies(**kw_dict) for kw_dict in kw_list] + return resultDict + +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) + +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('--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() + main() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies/mismatch_frequencies.xml Mon Feb 02 10:33:21 2015 -0500 @@ -0,0 +1,70 @@ +<tool id="mismatch_frequencies" name="Mismatch Frequencies" version="0.0.5" 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 + --name + #for i in $rep + "$i.input_file.name" + #end for + --output_pdf $output_pdf --output_tab $output_tab --min $min_length --max $max_length + --n_mm $number_of_mismatches + --five_p $five_p + --three_p $three_p + </command> + <inputs> + <repeat name="rep" title="alignment files" min="1"> + <param name="input_file" type="data" format="bam,sam" label="Alignment file" help="The input alignment file(s) for which you want to calculate mismatch frequencies."/> + </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="tabular" name="output_tab" /> + <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" /> + <output name="tabular" file="mismatch.tab" ftype="tabular"/> + </test> + </tests> + <help> + +.. 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> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies/test-data/mismatch.tab Mon Feb 02 10:33:21 2015 -0500 @@ -0,0 +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 valid reads +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mismatch_frequencies/tool_dependencies.xml Mon Feb 02 10:33:21 2015 -0500 @@ -0,0 +1,12 @@ +<?xml version="1.0"?> +<tool_dependency> + <package name="pysam" version="0.8.1"> + <repository changeset_revision="6b6843e15541" name="package_pysam_0_8_1" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> + </package> + <package name="pandas" version="0.14"> + <repository changeset_revision="21afd61aae1e" name="package_pandas_0_14" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> + </package> + <package name="matplotlib" version="1.4"> + <repository changeset_revision="38e91928f905" name="package_matplotlib_1_4" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> + </package> +</tool_dependency>
--- a/test-data/mismatch.tab Wed Jan 28 13:12:56 2015 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -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 valid reads -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
--- a/tool_dependencies.xml Wed Jan 28 13:12:56 2015 +0100 +++ /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="a7f103854ad5" name="package_pysam_0_7_7" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> - <package name="pandas" version="0.14"> - <repository changeset_revision="21afd61aae1e" name="package_pandas_0_14" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> - <package name="matplotlib" version="1.4"> - <repository changeset_revision="6424ce261dab" name="package_matplotlib_1_4" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> -</tool_dependency>