diff Plink.pl @ 0:fe39a4677281 draft

Uploaded

--- /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");
+
+}
+
+
+
+
+