<commit_before>#include "Variation.h"
#include "streamFuns.h"
#include "SequenceFuns.h"
#include "TimeFunctions.h"
#include "serviceFuns.h"

Variation::Variation (Parameters* Pin)
{
    P=Pin;
    if (!P->var.yes)
    {
        yes=false;
        return;
    };
    
    yes=true;
    
    varOutFileName=P->outFileNamePrefix+"Variation.out";
    varOutStream.open(varOutFileName);
    
    vcfFile=P->var.vcfFile;
    loadVCF(vcfFile);

};

void scanVCF(ifstream& vcf, bool load, Parameters* P, SNP& snp)
{
    snp.N=0;
    uint nlines=0;
    while (true)
    {
        string chr,id, ref, alt, dummy, sample;
        uint pos;
        nlines++;
        
        vcf >> chr;
        if (!vcf.good())
        {
            break;
        };
        
        if (chr.at(0)!='#')
        {
            vcf >> pos >> id >> ref >> alt >> dummy >> dummy >> dummy >> dummy >> sample;
            
            vector <string> altV(3);
            
            if (ref.size()==1 && splitString(alt,',',altV)==1)
            {
                altV.insert(altV.begin(),ref);//add ref to the beginning
                
                if (P->chrNameIndex.count(chr)==0) {//chr not in Genome
                    if (!load)
                    {
                        P->inOut->logMain << "WARNING: while processing varVCFfile file=" << P->var.vcfFile <<": chromosome '"<<chr<<"' not found in Genome fasta file\n";
                    };
                } else if (sample.size()<3)
                {
                    //undefined genotype
                } else if (sample.size()>3 && sample.at(3)!=':')
                {
                    if (!load)
                    {
                        P->inOut->logMain << "WARNING: while processing varVCFfile file=" << P->var.vcfFile <<": more than 2 alleles per sample for line number "<<nlines<<"\n";
                    };
                } else if (sample.at(0)=='0' && sample.at(2)=='0')
                {    
                    //both alleles are reference, no need to do anything                    
                } else if (altV.at( atoi(&sample.at(0)) ).at(0)==ref.at(0) && altV.at( atoi(&sample.at(2)) ).at(0)==ref.at(0))
                {
                    //both alleles are reference, no need to do anything
                } else
                {
                    if (load)
                    {
                        snp.loci[snp.N]=pos-1+P->chrStart[P->chrNameIndex[chr]];
                        snp.nt[snp.N][0]=convertNt01234( ref.at(0) );
                        snp.nt[snp.N][1]=convertNt01234( altV.at( atoi(&sample.at(0)) ).at(0) );
                        snp.nt[snp.N][2]=convertNt01234( altV.at( atoi(&sample.at(2)) ).at(0) );
                        
                        //cout << chr <<"\t"<< pos <<"\t"<< ref.at(0) <<"\t"<< altV.at( atoi(&sample.at(0)) ).at(0) <<"\t"<< altV.at( atoi(&sample.at(2)) ).at(0) <<"\n";

                    };
                    snp.N++;
                };
            };
        };
        getline(vcf,dummy);
    };
};


void Variation::loadVCF(string fileIn)
{
    time_t rawTime;
    time(&rawTime);
    P->inOut->logMain     << timeMonthDayTime(rawTime) <<" ..... Loading variations VCF GTF\n" <<flush;
    *P->inOut->logStdOut  << timeMonthDayTime(rawTime) <<" ..... Loading variations VCF GTF\n" <<flush;
        
    ifstream & vcf = ifstrOpen(fileIn, ERROR_OUT, "SOLUTION: check the path and permissions of the VCF file: "+fileIn, P);
    scanVCF(vcf, false, P, snp);
    
    snp.loci=new uint[snp.N];
    snp.nt.resize(snp.N);
//     snp.nt1=new char[snp.N][3];
    
//     for (int ii=0,ii<3)
//     {
//         snp.nt[
    vcf.clear();
    vcf.seekg(0,ios::beg);
    scanVCF(vcf, true, P, snp);
    vcf.close();

    time(&rawTime);
    P->inOut->logMain << timeMonthDayTime(rawTime) <<" ..... Loaded VCF data, found "<<snp.N<< " SNPs"<<endl;
    
    uint *s1=new uint[2*snp.N];
    for (uint ii=0;ii<snp.N; ii++) 
    {
        s1[2*ii]=snp.loci[ii];
        s1[2*ii+1]=ii;
    };
    
    qsort((void*)s1, snp.N, 2*sizeof(uint), funCompareUint1);
    
    vector< array<char,3> > nt1=snp.nt;
    for (uint ii=0;ii<snp.N; ii++) 
    {
        snp.loci[ii]=s1[2*ii];
        snp.nt[ii]=nt1.at(s1[2*ii+1]);
    };
    //sort SNPs by coordinate
    time(&rawTime);
    P->inOut->logMain << timeMonthDayTime(rawTime) <<" ..... Finished sorting VCF data"<<endl;
};

void SNP::snpOnBlocks(uint blockStart, uint blockL, int blockShift, vector<vector<array<int,2>>> &snpV)
{
    int32 isnp=binarySearch1b <uint> (blockStart, loci, N);
    while (isnp<N && loci[isnp]<(blockStart+blockL))
    {
        for (int ii=0;ii<2;ii++)
        {
            if (nt[isnp][ii+1]!=nt[isnp][0])
            {//allele different from reference
                array<int,2> snp1;
                snp1[0]=(int) (loci[isnp]-blockStart)+blockShift;
                snp1[1]=(int) nt[isnp][ii+1];
                snpV[ii].push_back(snp1);
            };
        };
        ++isnp;
    };
};

vector<vector<array<int,2>>> Variation::sjdbSnp(uint sjStart, uint sjEnd, uint sjdbOverhang)
{          
    vector<vector<array<int,2>>> snpV(2);
    
    if (!yes)
    {//no variation
        return snpV;
    };

    snp.snpOnBlocks(sjStart-sjdbOverhang, sjdbOverhang, 0,            snpV);
    snp.snpOnBlocks(sjEnd+1,              sjdbOverhang, sjdbOverhang, snpV);

    if (snpV.at(0).empty() && snpV.at(1).empty())
    {
        snpV.pop_back();
    } else if (snpV.at(0) == snpV.at(1))
    {
        snpV.pop_back();
    };
     
    return snpV;
};<commit_msg>Fixed a bug that prevent correct operation when inserting junctions without Variation.<commit_after>#include "Variation.h"
#include "streamFuns.h"
#include "SequenceFuns.h"
#include "TimeFunctions.h"
#include "serviceFuns.h"

Variation::Variation (Parameters* Pin)
{
    P=Pin;
    if (!P->var.yes)
    {
        yes=false;
        return;
    };
    
    yes=true;
    
    varOutFileName=P->outFileNamePrefix+"Variation.out";
    varOutStream.open(varOutFileName);
    
    vcfFile=P->var.vcfFile;
    loadVCF(vcfFile);

};

void scanVCF(ifstream& vcf, bool load, Parameters* P, SNP& snp)
{
    snp.N=0;
    uint nlines=0;
    while (true)
    {
        string chr,id, ref, alt, dummy, sample;
        uint pos;
        nlines++;
        
        vcf >> chr;
        if (!vcf.good())
        {
            break;
        };
        
        if (chr.at(0)!='#')
        {
            vcf >> pos >> id >> ref >> alt >> dummy >> dummy >> dummy >> dummy >> sample;
            
            vector <string> altV(3);
            
            if (ref.size()==1 && splitString(alt,',',altV)==1)
            {
                altV.insert(altV.begin(),ref);//add ref to the beginning
                
                if (P->chrNameIndex.count(chr)==0) {//chr not in Genome
                    if (!load)
                    {
                        P->inOut->logMain << "WARNING: while processing varVCFfile file=" << P->var.vcfFile <<": chromosome '"<<chr<<"' not found in Genome fasta file\n";
                    };
                } else if (sample.size()<3)
                {
                    //undefined genotype
                } else if (sample.size()>3 && sample.at(3)!=':')
                {
                    if (!load)
                    {
                        P->inOut->logMain << "WARNING: while processing varVCFfile file=" << P->var.vcfFile <<": more than 2 alleles per sample for line number "<<nlines<<"\n";
                    };
                } else if (sample.at(0)=='0' && sample.at(2)=='0')
                {    
                    //both alleles are reference, no need to do anything                    
                } else if (altV.at( atoi(&sample.at(0)) ).at(0)==ref.at(0) && altV.at( atoi(&sample.at(2)) ).at(0)==ref.at(0))
                {
                    //both alleles are reference, no need to do anything
                } else
                {
                    if (load)
                    {
                        snp.loci[snp.N]=pos-1+P->chrStart[P->chrNameIndex[chr]];
                        snp.nt[snp.N][0]=convertNt01234( ref.at(0) );
                        snp.nt[snp.N][1]=convertNt01234( altV.at( atoi(&sample.at(0)) ).at(0) );
                        snp.nt[snp.N][2]=convertNt01234( altV.at( atoi(&sample.at(2)) ).at(0) );
                        
                        //cout << chr <<"\t"<< pos <<"\t"<< ref.at(0) <<"\t"<< altV.at( atoi(&sample.at(0)) ).at(0) <<"\t"<< altV.at( atoi(&sample.at(2)) ).at(0) <<"\n";

                    };
                    snp.N++;
                };
            };
        };
        getline(vcf,dummy);
    };
};


void Variation::loadVCF(string fileIn)
{
    time_t rawTime;
    time(&rawTime);
    P->inOut->logMain     << timeMonthDayTime(rawTime) <<" ..... Loading variations VCF GTF\n" <<flush;
    *P->inOut->logStdOut  << timeMonthDayTime(rawTime) <<" ..... Loading variations VCF GTF\n" <<flush;
        
    ifstream & vcf = ifstrOpen(fileIn, ERROR_OUT, "SOLUTION: check the path and permissions of the VCF file: "+fileIn, P);
    scanVCF(vcf, false, P, snp);
    
    snp.loci=new uint[snp.N];
    snp.nt.resize(snp.N);
//     snp.nt1=new char[snp.N][3];
    
//     for (int ii=0,ii<3)
//     {
//         snp.nt[
    vcf.clear();
    vcf.seekg(0,ios::beg);
    scanVCF(vcf, true, P, snp);
    vcf.close();

    time(&rawTime);
    P->inOut->logMain << timeMonthDayTime(rawTime) <<" ..... Loaded VCF data, found "<<snp.N<< " SNPs"<<endl;
    
    uint *s1=new uint[2*snp.N];
    for (uint ii=0;ii<snp.N; ii++) 
    {
        s1[2*ii]=snp.loci[ii];
        s1[2*ii+1]=ii;
    };
    
    qsort((void*)s1, snp.N, 2*sizeof(uint), funCompareUint1);
    
    vector< array<char,3> > nt1=snp.nt;
    for (uint ii=0;ii<snp.N; ii++) 
    {
        snp.loci[ii]=s1[2*ii];
        snp.nt[ii]=nt1.at(s1[2*ii+1]);
    };
    //sort SNPs by coordinate
    time(&rawTime);
    P->inOut->logMain << timeMonthDayTime(rawTime) <<" ..... Finished sorting VCF data"<<endl;
};

void SNP::snpOnBlocks(uint blockStart, uint blockL, int blockShift, vector<vector<array<int,2>>> &snpV)
{
    int32 isnp=binarySearch1b <uint> (blockStart, loci, N);
    while (isnp<N && loci[isnp]<(blockStart+blockL))
    {
        for (int ii=0;ii<2;ii++)
        {
            if (nt[isnp][ii+1]!=nt[isnp][0])
            {//allele different from reference
                array<int,2> snp1;
                snp1[0]=(int) (loci[isnp]-blockStart)+blockShift;
                snp1[1]=(int) nt[isnp][ii+1];
                snpV[ii].push_back(snp1);
            };
        };
        ++isnp;
    };
};

vector<vector<array<int,2>>> Variation::sjdbSnp(uint sjStart, uint sjEnd, uint sjdbOverhang)
{          
    vector<vector<array<int,2>>> snpV(2);
    
    if (!yes)
    {//no variation, return 1 empty element 
        vector<vector<array<int,2>>> snpV1(1);
        return snpV1;
    };

    snp.snpOnBlocks(sjStart-sjdbOverhang, sjdbOverhang, 0,            snpV);
    snp.snpOnBlocks(sjEnd+1,              sjdbOverhang, sjdbOverhang, snpV);

    if (snpV.at(0).empty() && snpV.at(1).empty())
    {
        snpV.pop_back();
    } else if (snpV.at(0) == snpV.at(1))
    {
        snpV.pop_back();
    };
     
    return snpV;
};
{'fraction_non_alphanumeric': 0.17521746848925973, 'fraction_numerical': 0.01819634297887449, 'mean_word_length': 1.6890214797136038, 'pattern_counts': {'":': 4, '<': 184, '<?xml version=': 0, '>': 116, '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': '2931076', 'n_tokens_mistral': 4156, 'n_tokens_neox': 3979, 'n_words': 723}