<commit_before>// Package depth parallelizes calls to samtools depths and outputs:
// 1) $prefix.callable.bed that contains collapsed per-base regions of NO/LOW/or CALLABLE coverage.
// where low is < MinCov.
// 2) $prefix.depth.bed that contains the average depth for each window interval specified by WindowSize.
// TODO: output gc-content in depth windows.
package depth

import (
	"bufio"
	"bytes"
	"fmt"
	"io"
	"log"
	"os"
	"regexp"
	"runtime"
	"strconv"
	"strings"

	arg "github.com/alexflint/go-arg"
	"github.com/brentp/faidx"
	"github.com/brentp/gargs/process"
	"github.com/brentp/xopen"
	"github.com/fatih/color"
)

type dargs struct {
	WindowSize   int    `arg:"-w,help:window size in which to calculate high-depth regions"`
	MaxMeanDepth int    `arg:"-m,help:windows with depth > than this are high-depth. The default reports the depth of all regions."`
	Ordered      bool   `args:"-o,help:force output to be in same order as input even with -p."`
	Q            int    `arg:"-Q,help:mapping quality cutoff"`
	Chrom        string `arg:"-c,help:optional chromosome to limit analysis"`
	MinCov       int    `arg:"help:minimum depth considered callable"`
	Stats        bool   `arg:"-s,help:report sequence stats [GC CpG masked] for each window"`
	Reference    string `arg:"-r,required,help:path to reference fasta"`
	Processes    int    `arg:"-p,help:number of processors to parallelize."`
	Bed          string `arg:"-b,help:optional file of positions or regions to restrict depth calculations."`
	Prefix       string `arg:"required,help:prefix for output files depth.bed and callable.bed"`
	Bam          string `arg:"positional,required,help:bam for which to calculate depth"`
}

// we echo the region first so the callback knows the full extents even if there is NOTE
// coverage for part of it.
const command = "echo %s; samtools depth -Q %d -d %d -r %s --reference %s %s"

// this is the size in basepairs of the genomic chunks for parallelization.
const step = 5000000

var exitCode = 0

func pcheck(e error) {
	if e != nil {
		log.Fatal(e)
	}
}

func min(a, b int) int {
	if a < b {
		return a
	}
	return b
}

func max(a, b int) int {
	if a > b {
		return a
	}
	return b
}

// match chrom:start-end and chrom\tstart\tend
var re = regexp.MustCompile("(.+?)[:\t](\\d+)([\\-\t])(\\d+).*?")

func chromStartEndFromLine(line []byte) (string, int, int) {
	ret := re.FindSubmatch(line)
	if len(ret) != 5 {
		log.Fatal("couldn't get region from line", string(line))
	}
	chrom, start, isep, end := ret[1], ret[2], ret[3], ret[4]
	// convert from bed to chrom:start-end region so add 1 to start
	istart, err := strconv.Atoi(string(start))
	if err != nil {
		log.Fatal(err)
	}
	if bytes.Equal(isep, []byte{'-'}) {
		istart--
	}
	iend, err := strconv.Atoi(string(end))
	if err != nil {
		log.Fatal(err)
	}
	return string(chrom), max(istart, 0), iend
}

func regionFromLine(line []byte) string {
	chrom, start, end := chromStartEndFromLine(line)
	// convert from bed to chrom:start-end region so add 1 to start
	return fmt.Sprintf("%s:%d-%d", chrom, start+1, end)
}

// when the user specified a Bed file of regions for coverage, this is used.
func genFromBed(ch chan string, args dargs) {
	rdr, err := xopen.Ropen(args.Bed)
	pcheck(err)
	for {
		line, err := rdr.ReadBytes('\n')
		if err == io.EOF {
			break
		}
		pcheck(err)
		if len(line) == 0 {
			continue
		}
		region := regionFromLine(line)
		ch <- fmt.Sprintf(command, region, args.Q, args.MaxMeanDepth+1000,
			region, args.Reference, args.Bam)
	}
	close(ch)
}

func genCommands(args dargs) chan string {
	ch := make(chan string)
	if args.Bed != "" {
		go genFromBed(ch, args)
		return ch
	}
	go func() {

		rdr, err := xopen.Ropen(args.Reference + ".fai")
		pcheck(err)
		for {
			line, err := rdr.ReadString('\n')
			if err == io.EOF {
				break
			}
			pcheck(err)
			toks := strings.Split(line, "\t")

			chrom := toks[0]
			if args.Chrom != "" && chrom != args.Chrom {
				continue
			}
			length, err := strconv.Atoi(toks[1])
			pcheck(err)
			for i := 0; i < length; i += step {
				region := fmt.Sprintf("%s:%d-%d", chrom, i, min(i+step, length))
				ch <- fmt.Sprintf(command, region, args.Q, args.MaxMeanDepth+1000,
					region, args.Reference, args.Bam)
			}
		}
		close(ch)
	}()
	return ch
}

// Main is run from the dispatcher
func Main() {

	args := dargs{WindowSize: 250,
		MaxMeanDepth: 0,
		MinCov:       4,
		Q:            1}
	p := arg.MustParse(&args)
	if args.Prefix == "" {<commit_msg>remove fasta from depth call as it's not used.<commit_after>// Package depth parallelizes calls to samtools depths and outputs:
// 1) $prefix.callable.bed that contains collapsed per-base regions of NO/LOW/or CALLABLE coverage.
// where low is < MinCov.
// 2) $prefix.depth.bed that contains the average depth for each window interval specified by WindowSize.
// TODO: output gc-content in depth windows.
package depth

import (
	"bufio"
	"bytes"
	"fmt"
	"io"
	"log"
	"os"
	"regexp"
	"runtime"
	"strconv"
	"strings"

	arg "github.com/alexflint/go-arg"
	"github.com/brentp/faidx"
	"github.com/brentp/gargs/process"
	"github.com/brentp/xopen"
	"github.com/fatih/color"
)

type dargs struct {
	WindowSize   int    `arg:"-w,help:window size in which to calculate high-depth regions"`
	MaxMeanDepth int    `arg:"-m,help:windows with depth > than this are high-depth. The default reports the depth of all regions."`
	Ordered      bool   `args:"-o,help:force output to be in same order as input even with -p."`
	Q            int    `arg:"-Q,help:mapping quality cutoff"`
	Chrom        string `arg:"-c,help:optional chromosome to limit analysis"`
	MinCov       int    `arg:"help:minimum depth considered callable"`
	Stats        bool   `arg:"-s,help:report sequence stats [GC CpG masked] for each window"`
	Reference    string `arg:"-r,required,help:path to reference fasta"`
	Processes    int    `arg:"-p,help:number of processors to parallelize."`
	Bed          string `arg:"-b,help:optional file of positions or regions to restrict depth calculations."`
	Prefix       string `arg:"required,help:prefix for output files depth.bed and callable.bed"`
	Bam          string `arg:"positional,required,help:bam for which to calculate depth"`
}

// we echo the region first so the callback knows the full extents even if there is NOTE
// coverage for part of it.
const command = "echo %s; samtools depth -Q %d -d %d -r %s %s"

// this is the size in basepairs of the genomic chunks for parallelization.
const step = 5000000

var exitCode = 0

func pcheck(e error) {
	if e != nil {
		log.Fatal(e)
	}
}

func min(a, b int) int {
	if a < b {
		return a
	}
	return b
}

func max(a, b int) int {
	if a > b {
		return a
	}
	return b
}

// match chrom:start-end and chrom\tstart\tend
var re = regexp.MustCompile("(.+?)[:\t](\\d+)([\\-\t])(\\d+).*?")

func chromStartEndFromLine(line []byte) (string, int, int) {
	ret := re.FindSubmatch(line)
	if len(ret) != 5 {
		log.Fatal("couldn't get region from line", string(line))
	}
	chrom, start, isep, end := ret[1], ret[2], ret[3], ret[4]
	// convert from bed to chrom:start-end region so add 1 to start
	istart, err := strconv.Atoi(string(start))
	if err != nil {
		log.Fatal(err)
	}
	if bytes.Equal(isep, []byte{'-'}) {
		istart--
	}
	iend, err := strconv.Atoi(string(end))
	if err != nil {
		log.Fatal(err)
	}
	return string(chrom), max(istart, 0), iend
}

func regionFromLine(line []byte) string {
	chrom, start, end := chromStartEndFromLine(line)
	// convert from bed to chrom:start-end region so add 1 to start
	return fmt.Sprintf("%s:%d-%d", chrom, start+1, end)
}

// when the user specified a Bed file of regions for coverage, this is used.
func genFromBed(ch chan string, args dargs) {
	rdr, err := xopen.Ropen(args.Bed)
	pcheck(err)
	for {
		line, err := rdr.ReadBytes('\n')
		if err == io.EOF {
			break
		}
		pcheck(err)
		if len(line) == 0 {
			continue
		}
		region := regionFromLine(line)
		ch <- fmt.Sprintf(command, region, args.Q, args.MaxMeanDepth+1000,
			region, args.Bam)
	}
	close(ch)
}

func genCommands(args dargs) chan string {
	ch := make(chan string)
	if args.Bed != "" {
		go genFromBed(ch, args)
		return ch
	}
	go func() {

		rdr, err := xopen.Ropen(args.Reference + ".fai")
		pcheck(err)
		for {
			line, err := rdr.ReadString('\n')
			if err == io.EOF {
				break
			}
			pcheck(err)
			toks := strings.Split(line, "\t")

			chrom := toks[0]
			if args.Chrom != "" && chrom != args.Chrom {
				continue
			}
			length, err := strconv.Atoi(toks[1])
			pcheck(err)
			for i := 0; i < length; i += step {
				region := fmt.Sprintf("%s:%d-%d", chrom, i, min(i+step, length))
				ch <- fmt.Sprintf(command, region, args.Q, args.MaxMeanDepth+1000,
					region, args.Bam)
			}
		}
		close(ch)
	}()
	return ch
}

// Main is run from the dispatcher
func Main() {

	args := dargs{WindowSize: 250,
		MaxMeanDepth: 0,
		MinCov:       4,
		Q:            1}
	p := arg.MustParse(&args)
	if args.Prefix == "" {
{'fraction_non_alphanumeric': 0.14386633774388877, 'fraction_numerical': 0.008297824624355236, 'mean_word_length': 3.2390684410646386, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 7, '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': '601267', 'n_tokens_mistral': 3478, 'n_tokens_neox': 3204, 'n_words': 990}