<commit_before>//
// Created by <NAME> on 1/3/17.
//

#include <iostream>
#include <fstream>
#include <random>
#include <algorithm>

#include <OpenMS/CHEMISTRY/Element.h>
#include <OpenMS/CHEMISTRY/AASequence.h>
#include <OpenMS/CHEMISTRY/ResidueDB.h>
#include <OpenMS/CHEMISTRY/ElementDB.h>

#include "FASTAParser.h"

static const ElementDB* elementDB = ElementDB::getInstance();
static const OpenMS::ResidueDB* residueDB = OpenMS::ResidueDB::getInstance();

static std::string AMINO_ACIDS = "ADEFGHIKLNPQRSTVWYCM";
static std::string AMINO_ACIDS_NO_SULFUR = "ADEFGHIKLNPQRSTVWY";
static std::string AMINO_ACIDS_SULFUR = "CM";

std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis_AA(0,1);


std::ofstream* openOutputFiles(std::string base_path, int max_depth)
{
    // create all output files and write header to each
    std::ofstream* outfiles = new std::ofstream[max_depth];
    for (int precursor_isotope = 0; precursor_isotope < max_depth; ++precursor_isotope)
    {
        std::string filename = "Precursor" + std::to_string(precursor_isotope) + ".tab";


        outfiles[precursor_isotope].open(base_path + filename);
        outfiles[precursor_isotope] << "probability" << "\tprecursor.mass" << std::endl;
    }

    return outfiles;
}

void closeOutputFiles(std::ofstream* outfiles, int max_depth)
{
    // close all output files
    for (int precursor_isotope = 0; precursor_isotope < max_depth; ++precursor_isotope)
    {
        outfiles[precursor_isotope].close();
    }
    delete[] outfiles;
}

void write_distribution(const OpenMS::AASequence &p, std::ofstream* outfiles, int max_depth, bool mono)
{
    OpenMS::EmpiricalFormula precursor_ef = p.getFormula();
    OpenMS::IsotopeDistribution precursor_id = precursor_ef.getIsotopeDistribution(0);

    for (int precursor_isotope = 0; precursor_isotope < max_depth && precursor_isotope < precursor_id.size(); ++precursor_isotope)
    {
        double mass = mono ? precursor_ef.getMonoWeight() : precursor_ef.getAverageWeight();

        outfiles[precursor_isotope] << precursor_id.getContainer()[precursor_isotope].second << "\t" << mass << std::endl;
    }
}

void proteome_isotopic_distributions(std::string base_path, std::string fasta_path, float max_mass, int num_sulfurs, int max_depth, bool mono)
{
    std::ofstream* outfiles = openOutputFiles(base_path, max_depth);

    FASTAParser parser(fasta_path, max_mass, 1, 100);
    for (auto itr = parser.begin(); itr != parser.end(); ++itr)
    {
        if (num_sulfurs < 0 || itr->getFormula().getNumberOf(elementDB->getElement("Sulfur")) == num_sulfurs)
        {
            write_distribution(*itr, outfiles, max_depth, mono);
        }
    }

    closeOutputFiles(outfiles, max_depth);
}

OpenMS::AASequence create_random_peptide_sequence(int peptide_length, std::vector<double> aa2prob, int num_sulfurs)
{
    OpenMS::AASequence random_peptide;

    if (num_sulfurs < 0)
    {
        for (int aa_index = 0; aa_index < peptide_length; ++aa_index)
        {
            double rand = dis_AA(gen);
            int index = std::lower_bound(aa2prob.begin(), aa2prob.end(), rand) - aa2prob.begin() - 1;
            random_peptide += residueDB->getResidue(AMINO_ACIDS[index]);
        }
    }
    else
    {
        // for insertion of sulfur containing amino acids
        for (int i = 0; i < num_sulfurs; ++i)
        {
            random_peptide += residueDB->getResidue(AMINO_ACIDS_SULFUR[0]);
        }

        // random amino acid insertion (non Sulfur and Selenium amino acids)
        for (int aa_index = 0; aa_index < peptide_length; ++aa_index)
        {
            double rand = dis_AA(gen);
            int index = std::lower_bound(aa2prob.begin(), aa2prob.end(), rand) - aa2prob.begin() - 1;
            random_peptide += residueDB->getResidue(AMINO_ACIDS_NO_SULFUR[index]);
        }
    }

    return random_peptide;
}


std::map<char, double> getAAProbabilities(std::string fasta_path, bool sulfur)
{
    std::map<char, double> aa2prob;
    for (char aa : AMINO_ACIDS)
    {
        aa2prob[aa] = 0.0;
    }

    std::vector<FASTAFile::FASTAEntry> proteins;
    FASTAFile().load(fasta_path, proteins);

    int count = 0;
    for (Size i = 0; i < proteins.size(); ++i)
    {
        for (int j = 0; j < proteins[i].sequence.size(); ++j)<commit_msg>added stochastic peptide sampling<commit_after>//
// Created by <NAME> on 1/3/17.
//

#include <iostream>
#include <fstream>
#include <random>
#include <algorithm>

#include <OpenMS/CHEMISTRY/Element.h>
#include <OpenMS/CHEMISTRY/AASequence.h>
#include <OpenMS/CHEMISTRY/ResidueDB.h>
#include <OpenMS/CHEMISTRY/ElementDB.h>

#include "FASTAParser.h"

static const ElementDB* elementDB = ElementDB::getInstance();
static const OpenMS::ResidueDB* residueDB = OpenMS::ResidueDB::getInstance();

static std::string AMINO_ACIDS = "ADEFGHIKLNPQRSTVWYCM";
static std::string AMINO_ACIDS_NO_SULFUR = "ADEFGHIKLNPQRSTVWY";
static std::string AMINO_ACIDS_SULFUR = "CM";

std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis_AA(0,1);
std::uniform_int_distribution<> dis_SULFUR(0,1);

std::ofstream* openOutputFiles(std::string base_path, int max_depth)
{
    // create all output files and write header to each
    std::ofstream* outfiles = new std::ofstream[max_depth];
    for (int precursor_isotope = 0; precursor_isotope < max_depth; ++precursor_isotope)
    {
        std::string filename = "Precursor" + std::to_string(precursor_isotope) + ".tab";


        outfiles[precursor_isotope].open(base_path + filename);
        outfiles[precursor_isotope] << "probability" << "\tprecursor.mass" << std::endl;
    }

    return outfiles;
}

void closeOutputFiles(std::ofstream* outfiles, int max_depth)
{
    // close all output files
    for (int precursor_isotope = 0; precursor_isotope < max_depth; ++precursor_isotope)
    {
        outfiles[precursor_isotope].close();
    }
    delete[] outfiles;
}

void write_distribution(const OpenMS::AASequence &p, std::ofstream* outfiles, int max_depth, bool mono)
{
    OpenMS::EmpiricalFormula precursor_ef = p.getFormula();
    OpenMS::IsotopeDistribution precursor_id = precursor_ef.getIsotopeDistribution(0);

    for (int precursor_isotope = 0; precursor_isotope < max_depth && precursor_isotope < precursor_id.size(); ++precursor_isotope)
    {
        double mass = mono ? precursor_ef.getMonoWeight() : precursor_ef.getAverageWeight();

        outfiles[precursor_isotope] << precursor_id.getContainer()[precursor_isotope].second << "\t" << mass << std::endl;
    }
}

void proteome_isotopic_distributions(std::string base_path, std::string fasta_path, float max_mass, int num_sulfurs, int max_depth, bool mono)
{
    std::ofstream* outfiles = openOutputFiles(base_path, max_depth);

    FASTAParser parser(fasta_path, max_mass, 1, 100);
    for (auto itr = parser.begin(); itr != parser.end(); ++itr)
    {
        if (num_sulfurs < 0 || itr->getFormula().getNumberOf(elementDB->getElement("Sulfur")) == num_sulfurs)
        {
            write_distribution(*itr, outfiles, max_depth, mono);
        }
    }

    closeOutputFiles(outfiles, max_depth);
}

OpenMS::AASequence create_random_peptide_sequence(int peptide_length, std::vector<double> aa2prob, int num_sulfurs)
{
    OpenMS::AASequence random_peptide;

    if (num_sulfurs < 0)
    {
        for (int aa_index = 0; aa_index < peptide_length; ++aa_index)
        {
            double rand = dis_AA(gen);
            int index = std::lower_bound(aa2prob.begin(), aa2prob.end(), rand) - aa2prob.begin() - 1;
            random_peptide += residueDB->getResidue(AMINO_ACIDS[index]);
        }
    }
    else
    {
        // for insertion of sulfur containing amino acids
        for (int i = 0; i < num_sulfurs; ++i)
        {
            random_peptide += residueDB->getResidue(AMINO_ACIDS_SULFUR[dis_SULFUR(gen)]);
        }

        // random amino acid insertion (non Sulfur and Selenium amino acids)
        for (int aa_index = 0; aa_index < peptide_length; ++aa_index)
        {
            double rand = dis_AA(gen);
            int index = std::lower_bound(aa2prob.begin(), aa2prob.end(), rand) - aa2prob.begin() - 1;
            random_peptide += residueDB->getResidue(AMINO_ACIDS_NO_SULFUR[index]);
        }
    }

    return random_peptide;
}


std::map<char, double> getAAProbabilities(std::string fasta_path, bool sulfur)
{
    std::map<char, double> aa2prob;
    for (char aa : AMINO_ACIDS)
    {
        aa2prob[aa] = 0.0;
    }

    std::vector<FASTAFile::FASTAEntry> proteins;
    FASTAFile().load(fasta_path, proteins);

    int count = 0;
    for (Size i = 0; i < proteins.size(); ++i)
    {
        for (int j = 0; j < proteins[i].sequence.size(); ++j)
{'fraction_non_alphanumeric': 0.12975699220541037, 'fraction_numerical': 0.009513984410820724, 'mean_word_length': 3.73412913727618, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1516801', 'n_tokens_mistral': 3252, 'n_tokens_neox': 2946, 'n_words': 602}