<div dir="ltr">Hello,<div><ol><li>I run MACS2 to find out were peaks are located (see attached).<br></li><li>I wrote a script (see attached) which:<br></li><ol><li>read the peaks information </li><li>extract the sequence from <i>genome.fasta</i> with help of the peak information<br></li><li>the extracted sequences are saved in a new FASTA file.</li></ol></ol><div>The <i>MACS_peaks.csv</i> file contains a column <i>fold_enrichmen</i>t and the script safe this information in <i>peak.fold_enrichment.</i></div><div><br></div><div>How is it possible to sort the output sequences with highest peak.fold_enrichment?</div></div><div><br></div><div><div><i>#!/usr/bin/env python</i></div><div><i>from Bio import SeqIO</i></div><div><i>from Bio.SeqRecord import SeqRecord</i></div><div><i>from collections import defaultdict</i></div><div><i>from collections import OrderedDict</i></div><div><i>from Bio.Alphabet import generic_dna</i></div><div><i>from Bio.Seq import Seq # BioPython Seq object</i></div><div><i>from Bio import SeqIO # BioPython SeqIO for manipulating fasta</i></div><div><i>import itertools</i></div><div><i><br></i></div><div><i><br></i></div><div><i>class PeakInfo:</i></div><div><i><br></i></div><div><i>    start = 0</i></div><div><i>    end = 0</i></div><div><i>    len = 0</i></div><div><i>    pileup = 0.0</i></div><div><i>    pvalue = 0.0</i></div><div><i>    fold_enrichment = 0.0</i></div><div><i>    qvalue = 0.0</i></div><div><i>    peak_name = ""</i></div><div><i><br></i></div><div><i><br></i></div><div><i>def extract_peaks():</i></div><div><i><br></i></div><div><i>    with open('MACS2_peaks.csv') as f:</i></div><div><i>        peaks = defaultdict(PeakInfo)</i></div><div><i>        for line in itertools.islice(f, 26, None):  # start=17, stop=None</i></div><div><i><br></i></div><div><i>            try:</i></div><div><i>                parts = [part.strip() for part in line.split(',')]</i></div><div><i>            except IndexError:</i></div><div><i>                continue</i></div><div><i><br></i></div><div><i>            peak = PeakInfo()</i></div><div><i>            peak.start = int(parts[1])</i></div><div><i>            peak.end = int(parts[2])</i></div><div><i>            peak.len = int(parts[3])</i></div><div><i>            peak.pvalue = float(parts[4])</i></div><div><i>            peak.pileup = float(parts[5])</i></div><div><i>            peak.fold_enrichment = float(parts[6])</i></div><div><i>            peak.qvalue = float(parts[7])</i></div><div><i>            peak.peak_name = parts[8]</i></div><div><i>            peaks[parts[0]] = peak</i></div><div><i><br></i></div><div><i>    return peaks</i></div><div><i><br></i></div><div><i><br></i></div><div><i>def extract_peak_seqs(fasta_file_name, out_file_name, peaks):</i></div><div><i><br></i></div><div><i>    out_file = open(out_file_name, "w")</i></div><div><i>    records = [] # list of SeqRecord to be written</i></div><div><i><br></i></div><div><i>    for record in SeqIO.parse(open(fasta_file_name, "rU"), "fasta"):</i></div><div><i>        # an additional check to see if there are contigs with no uniref90 reference genes to begin with</i></div><div><i>        if not <a href="http://record.id">record.id</a> in peaks:</i></div><div><i>            continue</i></div><div><i><br></i></div><div><i>        peak = peaks[<a href="http://record.id">record.id</a>]</i></div><div><i><br></i></div><div><i>        peak.fold_enrichment</i></div><div><i><br></i></div><div><i>        records.append(SeqRecord(Seq(str(record.seq[peak.start-1:peak.end]), generic_dna),</i></div><div><i>                        id=<a href="http://record.id">record.id</a> + "_" + peak.peak_name, description=""))</i></div><div><i><br></i></div><div><i>    SeqIO.write(records, out_file, "fasta")</i></div><div><i><br></i></div><div><i>    out_file.close()</i></div><div><i><br></i></div><div><i><br></i></div><div><i>if __name__ == "__main__":</i></div><div><i><br></i></div><div><i>    peaks = extract_peaks()</i></div><div><i><br></i></div><div><i>    extract_peak_seqs("genome.fasta", "peaks.fasta", peaks)</i></div><div><br></div></div><div><br></div><div>Thank you in advance.</div><div><br></div><div>Mic</div></div>