comparison matching.pl @ 21:9dcffd531c76 draft

Uploaded

21:9dcffd531c76

#!/usr/bin/perl -w
#Filename:
#Author: Tian Dongmei
#Email: [e-mail]
#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);
}