diff matching.pl @ 21:9dcffd531c76 draft

Uploaded
author big-tiandm
date Wed, 05 Nov 2014 21:09:35 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matching.pl	Wed Nov 05 21:09:35 2014 -0500
@@ -0,0 +1,71 @@
+#!/usr/bin/perl -w
+#Filename:
+#Author: Tian Dongmei
+#Email: tiandm@big.ac.cn
+#Date: 2013/7/19
+#Modified:
+#Description: 
+my $version=1.00;
+
+use strict;
+use Getopt::Long;
+
+my %opts;
+GetOptions(\%opts,"i=s","g=s","index:s","v:i","p:i","r:s","o=s","h");
+if (!(defined $opts{i} and defined $opts{o} ) || defined $opts{h}) { #necessary arguments
+&usage;
+}
+
+my $filein=$opts{'i'};
+my $fileout=$opts{'o'};
+unless ($fileout=~/\/$/) {$fileout.="/";}
+my $genome=$opts{'g'};
+my $mis=defined $opts{'v'}? $opts{'v'} : 0;
+my $hits=defined $opts{'r'}? $opts{'r'} : 25;
+my $index=defined $opts{'index'} ? $opts{'index'} : "";
+my $threads=defined $opts{'p'} ? $opts{'p'} : 1;
+
+
+#my $time=time();
+#my $mapdir=$fileout."/genome_match_".$time;
+my $mapdir=$fileout."/genome_match";
+mkdir $mapdir;
+chdir $mapdir;
+###check genome index
+if (-s $index.".1.ebwt") {
+}else{
+	`bowtie-build $genome genome`;
+	$index="genome";
+}
+
+### genome mapping
+`bowtie -v $mis -f -p $threads -m $hits -a --best --strata $index $filein --al genome_mapped.fa --un genome_not_mapped.fa > genome_mapped.bwt 2> run.log`;
+
+#`convert_bowtie_to_blast.pl  genome_mapped.bwt genome_mapped.fa $genome > genome_mapped.bst`;
+
+sub usage{
+print <<"USAGE";
+Version $version
+Usage:
+$0 -i -o
+options:
+-i input file# input reads fasta/fastq file
+-g input file# genome file
+-index file-prefix #(must be indexed by bowtie-build) The parameter
+                string must be the prefix of the bowtie index. For instance, if
+                the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
+                the prefix is 'h_sapiens_37_asm'.##can be null
+-v <int>           report end-to-end hits w/ <=v mismatches; ignore qualities,default 0;
+
+-p/--threads <int> number of alignment threads to launch (default: 1)
+
+-r int          a read is allowed to map up to this number of positions in the genome
+                default is 25 
+
+-o output directory 
+
+-h help
+USAGE
+exit(1);
+}
+