<commit_before>#!/usr/bin/perl
## <NAME>, Illinois Tech - 2016
## Version 1.3

use strict;
use warnings;
use Getopt::Long qw(GetOptions);

## User-defined environment variables.
## Change these to reflect your settings or leave blank (e.g. $bwa = '';) if set in $PATH
my $samtools = '/opt/samtools-1.3.1/bin/';	## Path to samtools 1.3.1+ - http://www.htslib.org/
my $bcftools = '/opt/bcftools-1.3.1/';		## Path to bcftools 1.3.1+ - http://www.htslib.org/
my $bwa = '/usr/bin/';				## Path to BWA - http://bio-bwa.sourceforge.net/
my $bowtie2 = '/opt/bowtie2-2.2.9/';		## Path to Bowtie2 - http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
my $hisat2 = '/opt/hisat2-2.0.4/';		## Path to HISAT2 - https://ccb.jhu.edu/software/hisat2/index.shtml
my $freebayes = '/opt/freebayes/bin/';		## Path to FreeBayes -  https://github.com/ekg/freebayes
my $varscan = '/opt/VarScan/';			## Path to VarScan2 jar file to use -  https://github.com/dkoboldt/varscan
my $mash = '/opt/Mash/bin/';			## Path to Mash - https://github.com/marbl/Mash

## Define which VarScan2 jar file to use
my $varjar = 'VarScan.v2.4.2.jar';

## Defining options
my $usage = "
USAGE = perl get_SNPs.pl [options]

EXAMPLE (simple): get_SNPs.pl --fasta *.fasta --fastq *.fastq
EXAMPLE (advanced): get_SNPs.pl --fasta *.fasta --fastq *.fastq --mapper bowtie2 --caller freebayes --threads 16

OPTIONS:

## Genetic distances
--mh	Evaluate genetic distances using Mash (Ondov et al. DOI: 10.1186/s13059-016-0997-x)
--out	Output file name [default: Mash.txt]
--sort	Sort Mash output by decreasing order of similarity

## Mapping options
--fasta		Reference genome(s) in fasta file
--fastq		Fastq reads to be mapped against reference(s)
--mapper	Read mapping tool: bwa, bowtie2 or hisat2 [default: bowtie2]
--caller	Variant caller: varscan2, bcftools or freebayes [default: varscan2]
--algo		BWA mapping algorithm:  bwasw, mem, samse [default: bwasw]
--threads	Number of processing threads [default: 16]
--bam		Keeps BAM files generated
--sam		Keeps SAM files generated; SAM files can be quite large

## VarScan2 parameters (see http://dkoboldt.github.io/varscan/using-varscan.html)
--indel				Calculates indels 	## Runs mpileup2indel
--mc (--min-coverage)		[default: 15]		## Minimum read depth at a position to make a call
--mr (--min-reads2)		[default: 5]		## Minimum supporting reads at a position to call variants
--maq (--min-avg-qual)		[default: 28]		## Minimum base quality at a position to count a read
--mvf (--min-var-freq)		[default: 0.2]		## Minimum variant allele frequency threshold
--mhom (--min-freq-for-hom)	[default: 0.75]		## Minimum frequency to call homozygote
--pv (--p-value)		[default: 1e-02]	## P-value threshold for calling variants 
--sf (--strand-filter)		[default: 0]		## 0 or 1; 1 ignores variants with >90% support on one strand

## FreeBayes/BCFtools (see https://github.com/ekg/freebayes/; https://samtools.github.io/bcftools/bcftools.html)
--ploidy			[default: 1]		## Change ploidy (if needed)
";

die "$usage\n" unless@ARGV;

## Genetic distances
my $mh = '';
my $out = 'Mash.txt';
my $sort = '';
## Mapping
my $mapper = 'bowtie2';
my $caller = 'varscan2';
my $algo = 'bwasw';
my $threads = 16;
my $bam = '';
my $sam = '';
my @fasta;
my @fastq;
## VarScan
my $indel = '';
my $mc = 15;
my $mr = 5;
my $maq = 28;
my $mvf = '0.2';
my $mhom = '0.75';
my $pv = '1e-02';
my $sf = 0;
## FreeBayes/BCFtools
my $ploidy = 1;

GetOptions(
	'mh' => \$mh,
	'out=s' => \$out,
	'sort' => \$sort,
	'mapper=s' => \$mapper,
	'caller=s' => \$caller,
	'algo=s' => \$algo,
	'threads=i' => \$threads,
	'bam' => \$bam,
	'sam' => \$sam,
	'fasta=s@{1,}' => \@fasta,
	'fastq=s@{1,}' => \@fastq,
	'indel' => \$indel,
	'mc|min-coverage=i' => \$mc,
	'mr|min-reads2=i' => \$mr,
	'maq|min-avg-qual=i' => \$maq,
	'mvf|min-var-freq=s' => \$mvf,
	'mhom|min-freq-for-hom=s' => \$mhom,
	'pv|p-value=s' => \$pv,
	'sf|strand-filter=i' => \$sf,
	'ploidy=i' => \$ploidy,
);

## Timestamps and logs
my $start = localtime();
my $tstart = time;
my $todo = scalar(@fastq)*scalar(@fasta);
open LOG, ">time.$mapper.$caller.log"; ## Keep track of running time
print LOG "Mapping/SNP calling started on: $start\n";
print LOG "A total of $todo pairwise comparisons will be performed\n";
open MAP, ">>mapping.$mapper.log"; ## Keep track of read mapper STDERR messages
my $comparison = 0;

## Genetic distances with MASH

if ($mh){
	system "echo Running Mash genetic distance analysis...";
	system "$mash"."mash sketch -o reference @fasta";
	system "$mash"."mash dist reference.msh @fasta > $out";
	if ($sort){
		system "echo Sorting out Mash results - See $out.sorted";
		system "sort -gk3 $out > $out.sorted";
	}<commit_msg>Mash support<commit_after>#!/usr/bin/perl
## <NAME>, Illinois Tech - 2016
## Version 1.3

use strict;
use warnings;
use Getopt::Long qw(GetOptions);

## User-defined environment variables.
## Change these to reflect your settings or leave blank (e.g. $bwa = '';) if set in $PATH
my $samtools = '/opt/samtools-1.3.1/bin/';	## Path to samtools 1.3.1+ - http://www.htslib.org/
my $bcftools = '/opt/bcftools-1.3.1/';		## Path to bcftools 1.3.1+ - http://www.htslib.org/
my $bwa = '/usr/bin/';				## Path to BWA - http://bio-bwa.sourceforge.net/
my $bowtie2 = '/opt/bowtie2-2.2.9/';		## Path to Bowtie2 - http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
my $hisat2 = '/opt/hisat2-2.0.4/';		## Path to HISAT2 - https://ccb.jhu.edu/software/hisat2/index.shtml
my $freebayes = '/opt/freebayes/bin/';		## Path to FreeBayes -  https://github.com/ekg/freebayes
my $varscan = '/opt/VarScan/';			## Path to VarScan2 jar file to use -  https://github.com/dkoboldt/varscan
my $mash = '/opt/Mash/bin/';			## Path to Mash - https://github.com/marbl/Mash

## Define which VarScan2 jar file to use
my $varjar = 'VarScan.v2.4.2.jar';

## Defining options
my $usage = "
USAGE = perl get_SNPs.pl [options]

EXAMPLE (simple): get_SNPs.pl --fasta *.fasta --fastq *.fastq
EXAMPLE (advanced): get_SNPs.pl --fasta *.fasta --fastq *.fastq --mapper bowtie2 --caller freebayes --threads 16

OPTIONS:

## Genetic distances
--mh	Evaluate genetic distances using Mash (Ondov et al. DOI: 10.1186/s13059-016-0997-x)
--out	Output file name [default: Mash.txt]
--sort	Sort Mash output by decreasing order of similarity

## Mapping options
--fasta		Reference genome(s) in fasta file
--fastq		Fastq reads to be mapped against reference(s)
--mapper	Read mapping tool: bwa, bowtie2 or hisat2 [default: bowtie2]
--caller	Variant caller: varscan2, bcftools or freebayes [default: varscan2]
--algo		BWA mapping algorithm:  bwasw, mem, samse [default: bwasw]
--threads	Number of processing threads [default: 16]
--bam		Keeps BAM files generated
--sam		Keeps SAM files generated; SAM files can be quite large

## VarScan2 parameters (see http://dkoboldt.github.io/varscan/using-varscan.html)
--indel				Calculates indels 	## Runs mpileup2indel
--mc (--min-coverage)		[default: 15]		## Minimum read depth at a position to make a call
--mr (--min-reads2)		[default: 5]		## Minimum supporting reads at a position to call variants
--maq (--min-avg-qual)		[default: 28]		## Minimum base quality at a position to count a read
--mvf (--min-var-freq)		[default: 0.2]		## Minimum variant allele frequency threshold
--mhom (--min-freq-for-hom)	[default: 0.75]		## Minimum frequency to call homozygote
--pv (--p-value)		[default: 1e-02]	## P-value threshold for calling variants 
--sf (--strand-filter)		[default: 0]		## 0 or 1; 1 ignores variants with >90% support on one strand

## FreeBayes/BCFtools (see https://github.com/ekg/freebayes/; https://samtools.github.io/bcftools/bcftools.html)
--ploidy			[default: 1]		## Change ploidy (if needed)
";

die "$usage\n" unless@ARGV;

## Genetic distances
my $mh = '';
my $out = 'Mash.txt';
my $sort = '';
## Mapping
my $mapper = 'bowtie2';
my $caller = 'varscan2';
my $algo = 'bwasw';
my $threads = 16;
my $bam = '';
my $sam = '';
my @fasta;
my @fastq;
## VarScan
my $indel = '';
my $mc = 15;
my $mr = 5;
my $maq = 28;
my $mvf = '0.2';
my $mhom = '0.75';
my $pv = '1e-02';
my $sf = 0;
## FreeBayes/BCFtools
my $ploidy = 1;

GetOptions(
	'mh' => \$mh,
	'out=s' => \$out,
	'sort' => \$sort,
	'mapper=s' => \$mapper,
	'caller=s' => \$caller,
	'algo=s' => \$algo,
	'threads=i' => \$threads,
	'bam' => \$bam,
	'sam' => \$sam,
	'fasta=s@{1,}' => \@fasta,
	'fastq=s@{1,}' => \@fastq,
	'indel' => \$indel,
	'mc|min-coverage=i' => \$mc,
	'mr|min-reads2=i' => \$mr,
	'maq|min-avg-qual=i' => \$maq,
	'mvf|min-var-freq=s' => \$mvf,
	'mhom|min-freq-for-hom=s' => \$mhom,
	'pv|p-value=s' => \$pv,
	'sf|strand-filter=i' => \$sf,
	'ploidy=i' => \$ploidy,
);

## Timestamps and logs
my $start = localtime();
my $tstart = time;
my $todo = scalar(@fastq)*scalar(@fasta);
open LOG, ">time.$mapper.$caller.log"; ## Keep track of running time
print LOG "Mapping/SNP calling started on: $start\n";
print LOG "A total of $todo pairwise comparisons will be performed\n";
open MAP, ">>mapping.$mapper.log"; ## Keep track of read mapper STDERR messages
my $comparison = 0;

## Genetic distances with MASH
if ($mh){
	system "echo Running Mash genetic distance analysis...";
	system "$mash"."mash sketch -o reference @fasta";
	system "$mash"."mash dist reference.msh @fasta > $out";
	if ($sort){
		system "echo Sorting out Mash results - See $out.sorted";
		system "sort -gk3 $out > $out.sorted";
	}
{'fraction_non_alphanumeric': 0.18345477925417916, 'fraction_numerical': 0.02400342906129447, 'mean_word_length': 3.6409746394828444, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 57, 'https://': 12, 'lorem ipsum': 0, 'www.': 4, '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': '6662753', 'n_tokens_mistral': 3968, 'n_tokens_neox': 3662, 'n_words': 1051}