Mercurial > repos > dereeper > plink
diff Plink.pl @ 0:fe39a4677281 draft
Uploaded
author | dereeper |
---|---|
date | Fri, 05 Aug 2016 09:46:55 -0400 |
parents | |
children | d6a7be1b5adb |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Plink.pl Fri Aug 05 09:46:55 2016 -0400 @@ -0,0 +1,242 @@ + +#!/usr/bin/perl + +use strict; +use Getopt::Long; +use Bio::SeqIO; + +my $usage = qq~Usage:$0 <args> [<opts>] + +where <args> are: + + -i, --input <VCF input> + -o, --out <Output basename> + + <opts> are: + + -s, --samples <Samples to be analyzed. Comma separated list> + -c, --chromosomes <List of chromosomes to be analyzed.> + -e, --export <Output format (VCF/freq/plink. Default: VCF> + -f, --frequency <Minimum MAF. Default: 0.001> + -m, --max_freq <Maximum MAF. Default: 0.5> + -a, --allow_missing <Allowed missing data proportion per site. Must be comprised between 0 and 1. Default: 1> + -t, --type <Type of polymorphisms to keep (ALL/SNP). Default: ALL> + -b, --bounds <Lower bound and upper bound for a range of sites to be processed (start,end). Default: 1, 100000000> + -r, --remove_filt <Remove all sites with a FILTER flag other than PASS (true/false). Default: false> + -d, --distance <Thin sites so that no two sites are within the specified distance from one another. Default: 0> +~; +$usage .= "\n"; + +my ($input,$out); + +my $PLINK_EXE = "plink"; + + +#my $indel_size_max = 500; +#my $indel_size_min = 1; +my $frequency_max = 0.5; +my $frequency_min = 0.001; +my $pos_max = 100000000000; +my $pos_min = 0; +my $filter_snp_type = "all"; +my $remove_filt = "False"; + +my $missing_data = 1; +my $export = "VCF"; +my $type = "ALL"; +my $bounds; +my $samples; +my $chromosomes; +my $thin; + +GetOptions( + "input=s" => \$input, + "out=s" => \$out, + "samples=s" => \$samples, + "chromosomes=s" => \$chromosomes, + "frequency=s" => \$frequency_min, + "max_freq=s" => \$frequency_max, + "allow_missing=s"=> \$missing_data, + "export=s" => \$export, + "type=s" => \$type, + "bounds=s" => \$bounds, + "remove_filt=s" => \$remove_filt, + "distance=s" => \$thin +); + + +die $usage + if ( !$input || !$out); + +if ($samples && $samples =~/^([\w\,\-\.]+)\s*$/){ + $samples = $1; +} +elsif ($samples){ + die "Error: Samples must be a comma separated list of string\n"; +} +if ($bounds && $bounds =~/^([\d\,]+)\s*$/){ + $bounds = $1; +} +elsif($bounds){ + die "Error: Bounds must be a comma separated list of integers\n"; +} + +my $minfreq_cmd = ""; +if ($frequency_min && $frequency_min > 0 && $frequency_min =~/^([\d\.]+)\s*$/){ + $frequency_min = $1; + $minfreq_cmd = "--maf $frequency_min"; +} +elsif ($frequency_min == 0){ + $minfreq_cmd = ""; +} +elsif ($frequency_min){ + die "Error: frequency must be an integer\n"; +} +if ($thin && $thin =~/^([\d\.]+)\s*$/){ + $thin = $1; +} +elsif ($thin){ + die "Error: frequency must be an integer\n"; +} +my $maxfreq_cmd = ""; +if ($frequency_max && $frequency_max =~/^([\d\.]+)\s*$/){ + $frequency_max = $1; + if ($frequency_max < 0.5){ + $maxfreq_cmd = "--max-maf $frequency_max"; + } +} +elsif($frequency_max){ + die "Error: frequency must be an integer\n"; +} +if ($missing_data =~/^([\d\.]+)\s*$/){ + $missing_data = $1; + #$missing_data = 1 - $missing_data; +} +elsif ($missing_data){ + die "Error: Missing data must be an integer\n"; +} +if ($export && $export =~/^([\w]+)\s*$/){ + $export = $1; +} +elsif($export){ + die "Error: Export must be a string\n"; +} +if ($type && $type =~/^([\w]+)\s*$/){ + $type = $1; +} +elsif($type){ + die "Error: Type must be a string\n"; +} + + +my @dnasamples; +if ($samples) +{ + @dnasamples = split(",",$samples); +} +my @boundaries; +if ($bounds) +{ + @boundaries = split(",",$bounds); +} + + +my $experiment = "chromosomes"; +my $table = ""; +my %genes; +my @snp_ids; +my @snp_ids_and_positions; +my @snp_ids_and_positions_all; +my $gene; +my $snp_num = 0; +my %ref_sequences; +my %snps_of_gene; + +my $indiv_cmd = ""; +if (@dnasamples) +{ + if (scalar @dnasamples > 1) + { + open(my $S,">$out.samples"); + foreach my $samp(@dnasamples){ + print $S "$samp $samp\n"; + } + close($S); + $indiv_cmd = "--keep $out.samples "; + } + else + { + $indiv_cmd = "--indv " . join(" --indv ",@dnasamples); + } +} + +my $chrom_cmd = ""; +if ($chromosomes) +{ + $chrom_cmd = "--chr ".$chromosomes +} + +my $export_cmd = "--recode vcf-iid"; +if ($export eq "bcf"){ + $export_cmd = "--recode bcf"; +} +if ($export eq "freq"){ + $export_cmd = "--freq"; +} +if ($export eq "plink"){ + $export_cmd = "--make-bed"; +} +if ($export eq "bed"){ + $export_cmd = "--make-bed"; +} + + +my $bounds_cmd = ""; +if (@boundaries && $chrom_cmd=~/\w/ && $chrom_cmd !~/,/) +{ + $bounds_cmd = "--from-bp $boundaries[0] --to-bp $boundaries[1]"; +} + + + +my $type_cmd = ""; +if ($type eq "SNP") +{ + $type_cmd = "--snps-only"; +} + +my $filt_cmd = ""; +if ($remove_filt eq "true") +{ + $filt_cmd = "--remove-filtered-all"; +} + +my $thin_cmd = ""; +if ($thin){ + $thin_cmd = "--bp-space $thin"; +} + +#my $bcf_input = $input; +#$bcf_input =~s/vcf/bcf/g; +my $bcf_input; +my $bed_input = $input; +$bed_input =~s/\.bed//g; + +if (-e "$bed_input.bed"){ + system("$PLINK_EXE --bfile $bed_input --out $out $type_cmd $export_cmd $chrom_cmd $indiv_cmd $minfreq_cmd $maxfreq_cmd --geno $missing_data $thin_cmd $bounds_cmd --allow-extra-chr 1>$out.plink.stdout 2>$out.plink.stderr"); + # for first 1000 SNPs + system("$PLINK_EXE --bfile $bed_input --out $out.recode $type_cmd --recode vcf-fid $chrom_cmd $indiv_cmd $minfreq_cmd $maxfreq_cmd --geno $missing_data $thin_cmd $bounds_cmd --allow-extra-chr --thin-count 800 1>$out.2.plink.stdout 2>$out.2.plink.stderr"); +} +elsif (-e $bcf_input){ + system("$PLINK_EXE --bcf $bcf_input --out $out $type_cmd $export_cmd $chrom_cmd $indiv_cmd $minfreq_cmd $maxfreq_cmd --geno $missing_data $thin_cmd $bounds_cmd --allow-extra-chr 1>$out.plink.stdout 2>$out.plink.stderr"); +} +else +{ + system("$PLINK_EXE --vcf $input --out $out $type_cmd $export_cmd $chrom_cmd $indiv_cmd $minfreq_cmd $maxfreq_cmd --geno $missing_data $thin_cmd $bounds_cmd --allow-extra-chr 1>$out.3.plink.stdout 2>$out.3.plink.stderr"); + +} + + + + +