<commit_before>#include "Statistics.h"
#include "Exceptions.h"
#include "ChromosomalIndex.h"
#include "Log.h"
#include "cmath"
#include "numeric"
#include "algorithm"
#include "BasicStatistics.h"
#include <QVector>
#include "Pileup.h"
#include "FastqFileStream.h"
#include "LinePlot.h"
#include "ScatterPlot.h"
#include "BarPlot.h"
#include "Helper.h"
#include "SampleSimilarity.h"
#include <QFileInfo>
#include <QPair>
#include "Histogram.h"
#include "FilterCascade.h"

QCCollection Statistics::variantList(VariantList variants, bool filter)
{
    QCCollection output;

	//filter variants
	if (filter)
	{
		FilterResult filter_result(variants.count());
		FilterFilterColumnEmpty filter;
		filter.apply(variants, filter_result);
		filter_result.removeFlagged(variants);
	}

    //var_total
	output.insert(QCValue("variant count", variants.count(), "Total number of variants in the target region.", "QC:2000013"));

	//var_perc_dbsnp and high-impact variants
	if (variants.count()==0)
    {
		output.insert(QCValue("known variants percentage", "n/a (no variants)", "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
		output.insert(QCValue("high-impact variants percentage", "n/a (no variants)", "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
	}
	else
	{
		int i_csq = variants.annotationIndexByName("CSQ", true, false);
		if (i_csq==-1)
		{
			output.insert(QCValue("known variants percentage", "n/a (CSQ info field missing)", "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
			output.insert(QCValue("high-impact variants percentage", "n/a (CSQ info field missing)", "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
		}
		else
		{
			double dbsnp_count = 0;
			double high_impact_count = 0;
			for(int i=0; i<variants.count(); ++i)
			{
				if (variants[i].annotations().at(i_csq).contains("|rs")) //works without splitting by transcript
				{
					++dbsnp_count;
				}
				if (variants[i].annotations().at(i_csq).contains("|HIGH|")) //works without splitting by transcript
				{
					++high_impact_count;
				}
			}
			output.insert(QCValue("known variants percentage", 100.0*dbsnp_count/variants.count(), "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
			output.insert(QCValue("high-impact variants percentage", 100.0*high_impact_count/variants.count(), "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
		}
	}

	//homozygous variants
    const int i_gt = variants.annotationIndexByName("GT", true, false);
    if (variants.count()!=0 && i_gt!=-1)
    {
        double hom_count = 0;
        for(int i=0; i<variants.count(); ++i)
        {
			QByteArray geno = variants[i].annotations().at(i_gt);
			if (geno=="1/1" || geno=="1|1")
            {
                ++hom_count;
            }
        }
		output.insert(QCValue("homozygous variants percentage", 100.0*hom_count/variants.count(), "Percentage of variants that are called as homozygous.", "QC:2000016"));
    }
    else
    {
		output.insert(QCValue("homozygous variants percentage", "n/a (GT annotation not found, or no variants)", "Percentage of variants that are called as homozygous.", "QC:2000016"));
    }

    //var_perc_indel / var_ti_tv_ratio
    double indel_count = 0;
    double ti_count = 0;
    double tv_count = 0;
    for(int i=0; i<variants.count(); ++i)
    {
        const Variant& var = variants[i];
		if (var.ref().length()>1 || var.obs().length()>1)
        {
            ++indel_count;
        }
        else if ((var.obs()=="A" && var.ref()=="G") || (var.obs()=="G" && var.ref()=="A") || (var.obs()=="T" && var.ref()=="C") || (var.obs()=="C" && var.ref()=="T"))
        {
            ++ti_count;
        }
        else
        {
            ++tv_count;
        }
    }

    if (variants.count()!=0)
    {
		output.insert(QCValue("indel variants percentage", 100.0*indel_count/variants.count(), "Percentage of variants that are insertions/deletions.", "QC:2000017"));
    }
    else
    {
		output.insert(QCValue("indel variants percentage", "n/a (no variants)", "Percentage of variants that are insertions/deletions.", "QC:2000017"));
    }

    if (tv_count!=0)
    {
		output.insert(QCValue("transition/transversion ratio", ti_count/tv_count , "Transition/transversion ratio of SNV variants.", "QC:2000018"));
    }
    else
    {
		output.insert(QCValue("transition/transversion ratio", "n/a (no variants or tansversions)", "Transition/transversion ratio of SNV variants.", "QC:2000018"));
    }

    return output;
}

QCCollection Statistics::mapping(const BedFile& bed_file, const QString& bam_file, int min_mapq)
{
    //check target region is merged/sorted and create index
    if (!bed_file.isMergedAndSorted())
    {
        THROW(ArgumentException, "Merged and sorted BED file required for coverage details statistics!");
    }
    ChromosomalIndex<BedFile> roi_index(bed_file);

    //open BAM file
	BamReader reader(bam_file);

    //create coverage statistics data structure
    long long roi_bases = 0;
	QHash<int, QMap<int, int> > roi_cov;
    for (int i=0; i<bed_file.count(); ++i)
    {
        const BedLine& line = bed_file[i];

		if (!roi_cov.contains(line.chr().num()))
        {
			roi_cov.insert(line.chr().num(), QMap<int, int>());
        }

        for(int p=line.start(); p<=line.end(); ++p)
        {
			roi_cov[line.chr().num()].insert(p, 0);
        }
        roi_bases += line.length();
    }

    //init counts
    int al_total = 0;
    int al_mapped = 0;
    int al_ontarget = 0;
    int al_dup = 0;
    int al_proper_paired = 0;
    double bases_trimmed = 0;
    double bases_mapped = 0;
    double bases_clipped = 0;
    double insert_size_sum = 0;
	Histogram insert_dist(0, 999, 5);
    long long bases_usable = 0;
    int max_length = 0;
    bool paired_end = false;

    //iterate through all alignments
    BamAlignment al;
	while (reader.getNextAlignment(al))
    {
		//skip secondary alignments
		if (al.isSecondaryAlignment()) continue;

        ++al_total;
		max_length = std::max(max_length, al.length());

        //insert size
		if (al.isPaired())
        {
            paired_end = true;

			if (al.isProperPair())
            {
                ++al_proper_paired;
				int insert_size = std::min(abs(al.insertSize()), 999); //cap insert size at 1000
				insert_size_sum += insert_size;
				insert_dist.inc(insert_size, true);
            }
        }

		if (!al.isUnmapped())
        {
            ++al_mapped;

            //calculate soft/hard-clipped bases
			const int start_pos = al.start();
			const int end_pos = al.end();
			bases_mapped += al.length();
			const QList<CigarOp> cigar_data = al.cigarData();
			foreach(const CigarOp& op, cigar_data)
            {
				if (op.Type==BAM_CSOFT_CLIP || op.Type==BAM_CHARD_CLIP)
                {
					bases_clipped += op.Length;
                }
            }

            //calculate usable bases and base-resolution coverage
			const Chromosome& chr = reader.chromosome(al.chromosomeID());
            QVector<int> indices = roi_index.matchingIndices(chr, start_pos, end_pos);
            if (indices.count()!=0)
            {
                ++al_ontarget;

				if (!al.isDuplicate() && al.mappingQuality()>=min_mapq)
                {
                    foreach(int index, indices)
                    {
                        const int ol_start = std::max(bed_file[index].start(), start_pos);
                        const int ol_end = std::min(bed_file[index].end(), end_pos);
						bases_usable += ol_end - ol_start + 1;
						auto it = roi_cov[chr.num()].lowerBound(ol_start);
						auto end = roi_cov[chr.num()].upperBound(ol_end);
						while (it!=end)
						{
							(*it)++;
							++it;
                        }
                    }
                }
            }
        }

        //trimmed bases (this is not entirely correct if the first alignments are all trimmed, but saves the second pass through the data)
		if (al.length()<max_length)
        {
			bases_trimmed += (max_length - al.length());
        }

		if (al.isDuplicate())
        {
            ++al_dup;
        }
	}

	//calculate coverage depth statistics
	double avg_depth = (double) bases_usable / roi_bases;
	int hist_max = 999;
	int hist_step = 5;
	if (avg_depth>500)
	{
		hist_max += 1000;
		hist_step += 5;
	}
	if (avg_depth>1000)
	{
		hist_max += 1000;
		hist_step += 5;
	}
	Histogram depth_dist(0, hist_max, hist_step);
	QHashIterator<int, QMap<int, int> > it(roi_cov);<commit_msg>Statistics.cpp: do not increase hist_step above 10 to prevent wrong 10x coverage % values for ultradeep sequencing<commit_after>#include "Statistics.h"
#include "Exceptions.h"
#include "ChromosomalIndex.h"
#include "Log.h"
#include "cmath"
#include "numeric"
#include "algorithm"
#include "BasicStatistics.h"
#include <QVector>
#include "Pileup.h"
#include "FastqFileStream.h"
#include "LinePlot.h"
#include "ScatterPlot.h"
#include "BarPlot.h"
#include "Helper.h"
#include "SampleSimilarity.h"
#include <QFileInfo>
#include <QPair>
#include "Histogram.h"
#include "FilterCascade.h"

QCCollection Statistics::variantList(VariantList variants, bool filter)
{
    QCCollection output;

	//filter variants
	if (filter)
	{
		FilterResult filter_result(variants.count());
		FilterFilterColumnEmpty filter;
		filter.apply(variants, filter_result);
		filter_result.removeFlagged(variants);
	}

    //var_total
	output.insert(QCValue("variant count", variants.count(), "Total number of variants in the target region.", "QC:2000013"));

	//var_perc_dbsnp and high-impact variants
	if (variants.count()==0)
    {
		output.insert(QCValue("known variants percentage", "n/a (no variants)", "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
		output.insert(QCValue("high-impact variants percentage", "n/a (no variants)", "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
	}
	else
	{
		int i_csq = variants.annotationIndexByName("CSQ", true, false);
		if (i_csq==-1)
		{
			output.insert(QCValue("known variants percentage", "n/a (CSQ info field missing)", "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
			output.insert(QCValue("high-impact variants percentage", "n/a (CSQ info field missing)", "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
		}
		else
		{
			double dbsnp_count = 0;
			double high_impact_count = 0;
			for(int i=0; i<variants.count(); ++i)
			{
				if (variants[i].annotations().at(i_csq).contains("|rs")) //works without splitting by transcript
				{
					++dbsnp_count;
				}
				if (variants[i].annotations().at(i_csq).contains("|HIGH|")) //works without splitting by transcript
				{
					++high_impact_count;
				}
			}
			output.insert(QCValue("known variants percentage", 100.0*dbsnp_count/variants.count(), "Percentage of variants that are known polymorphisms in the dbSNP database.", "QC:2000014"));
			output.insert(QCValue("high-impact variants percentage", 100.0*high_impact_count/variants.count(), "Percentage of variants with high impact on the protein, i.e. stop-gain, stop-loss, frameshift, splice-acceptor or splice-donor variants.", "QC:2000015"));
		}
	}

	//homozygous variants
    const int i_gt = variants.annotationIndexByName("GT", true, false);
    if (variants.count()!=0 && i_gt!=-1)
    {
        double hom_count = 0;
        for(int i=0; i<variants.count(); ++i)
        {
			QByteArray geno = variants[i].annotations().at(i_gt);
			if (geno=="1/1" || geno=="1|1")
            {
                ++hom_count;
            }
        }
		output.insert(QCValue("homozygous variants percentage", 100.0*hom_count/variants.count(), "Percentage of variants that are called as homozygous.", "QC:2000016"));
    }
    else
    {
		output.insert(QCValue("homozygous variants percentage", "n/a (GT annotation not found, or no variants)", "Percentage of variants that are called as homozygous.", "QC:2000016"));
    }

    //var_perc_indel / var_ti_tv_ratio
    double indel_count = 0;
    double ti_count = 0;
    double tv_count = 0;
    for(int i=0; i<variants.count(); ++i)
    {
        const Variant& var = variants[i];
		if (var.ref().length()>1 || var.obs().length()>1)
        {
            ++indel_count;
        }
        else if ((var.obs()=="A" && var.ref()=="G") || (var.obs()=="G" && var.ref()=="A") || (var.obs()=="T" && var.ref()=="C") || (var.obs()=="C" && var.ref()=="T"))
        {
            ++ti_count;
        }
        else
        {
            ++tv_count;
        }
    }

    if (variants.count()!=0)
    {
		output.insert(QCValue("indel variants percentage", 100.0*indel_count/variants.count(), "Percentage of variants that are insertions/deletions.", "QC:2000017"));
    }
    else
    {
		output.insert(QCValue("indel variants percentage", "n/a (no variants)", "Percentage of variants that are insertions/deletions.", "QC:2000017"));
    }

    if (tv_count!=0)
    {
		output.insert(QCValue("transition/transversion ratio", ti_count/tv_count , "Transition/transversion ratio of SNV variants.", "QC:2000018"));
    }
    else
    {
		output.insert(QCValue("transition/transversion ratio", "n/a (no variants or tansversions)", "Transition/transversion ratio of SNV variants.", "QC:2000018"));
    }

    return output;
}

QCCollection Statistics::mapping(const BedFile& bed_file, const QString& bam_file, int min_mapq)
{
    //check target region is merged/sorted and create index
    if (!bed_file.isMergedAndSorted())
    {
        THROW(ArgumentException, "Merged and sorted BED file required for coverage details statistics!");
    }
    ChromosomalIndex<BedFile> roi_index(bed_file);

    //open BAM file
	BamReader reader(bam_file);

    //create coverage statistics data structure
    long long roi_bases = 0;
	QHash<int, QMap<int, int> > roi_cov;
    for (int i=0; i<bed_file.count(); ++i)
    {
        const BedLine& line = bed_file[i];

		if (!roi_cov.contains(line.chr().num()))
        {
			roi_cov.insert(line.chr().num(), QMap<int, int>());
        }

        for(int p=line.start(); p<=line.end(); ++p)
        {
			roi_cov[line.chr().num()].insert(p, 0);
        }
        roi_bases += line.length();
    }

    //init counts
    int al_total = 0;
    int al_mapped = 0;
    int al_ontarget = 0;
    int al_dup = 0;
    int al_proper_paired = 0;
    double bases_trimmed = 0;
    double bases_mapped = 0;
    double bases_clipped = 0;
    double insert_size_sum = 0;
	Histogram insert_dist(0, 999, 5);
    long long bases_usable = 0;
    int max_length = 0;
    bool paired_end = false;

    //iterate through all alignments
    BamAlignment al;
	while (reader.getNextAlignment(al))
    {
		//skip secondary alignments
		if (al.isSecondaryAlignment()) continue;

        ++al_total;
		max_length = std::max(max_length, al.length());

        //insert size
		if (al.isPaired())
        {
            paired_end = true;

			if (al.isProperPair())
            {
                ++al_proper_paired;
				int insert_size = std::min(abs(al.insertSize()), 999); //cap insert size at 1000
				insert_size_sum += insert_size;
				insert_dist.inc(insert_size, true);
            }
        }

		if (!al.isUnmapped())
        {
            ++al_mapped;

            //calculate soft/hard-clipped bases
			const int start_pos = al.start();
			const int end_pos = al.end();
			bases_mapped += al.length();
			const QList<CigarOp> cigar_data = al.cigarData();
			foreach(const CigarOp& op, cigar_data)
            {
				if (op.Type==BAM_CSOFT_CLIP || op.Type==BAM_CHARD_CLIP)
                {
					bases_clipped += op.Length;
                }
            }

            //calculate usable bases and base-resolution coverage
			const Chromosome& chr = reader.chromosome(al.chromosomeID());
            QVector<int> indices = roi_index.matchingIndices(chr, start_pos, end_pos);
            if (indices.count()!=0)
            {
                ++al_ontarget;

				if (!al.isDuplicate() && al.mappingQuality()>=min_mapq)
                {
                    foreach(int index, indices)
                    {
                        const int ol_start = std::max(bed_file[index].start(), start_pos);
                        const int ol_end = std::min(bed_file[index].end(), end_pos);
						bases_usable += ol_end - ol_start + 1;
						auto it = roi_cov[chr.num()].lowerBound(ol_start);
						auto end = roi_cov[chr.num()].upperBound(ol_end);
						while (it!=end)
						{
							(*it)++;
							++it;
                        }
                    }
                }
            }
        }

        //trimmed bases (this is not entirely correct if the first alignments are all trimmed, but saves the second pass through the data)
		if (al.length()<max_length)
        {
			bases_trimmed += (max_length - al.length());
        }

		if (al.isDuplicate())
        {
            ++al_dup;
        }
	}

	//calculate coverage depth statistics
	double avg_depth = (double) bases_usable / roi_bases;
	int hist_max = 999;
	int hist_step = 5;
	if (avg_depth>500)
	{
		hist_max += 1000;
		hist_step += 5;
	}
	if (avg_depth>1000)
	{
		hist_max += 1000;
	}
	Histogram depth_dist(0, hist_max, hist_step);
	QHashIterator<int, QMap<int, int> > it(roi_cov);
{'fraction_non_alphanumeric': 0.1414620378210556, 'fraction_numerical': 0.02026531188258538, 'mean_word_length': 2.88593989910068, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 35, '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': '6783158', 'n_tokens_mistral': 6695, 'n_tokens_neox': 5802, 'n_words': 1380}