view matching.pl @ 47:c75593f79aa9 draft

Uploaded

#!/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";
if(not -d $mapdir){
	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 --max genome_mapped_Mlimit.fa > genome_mapped.bwt  2> run.log`;

#`convert_bowtie_to_blast.pl  genome_mapped.bwt genome_mapped.fa $genome > genome_mapped.bst`;

sub Time{
        my $time=time();
        my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
        $month++;
        $year+=1900;
        if (length($sec) == 1) {$sec = "0"."$sec";}
        if (length($min) == 1) {$min = "0"."$min";}
        if (length($hour) == 1) {$hour = "0"."$hour";}
        if (length($day) == 1) {$day = "0"."$day";}
        if (length($month) == 1) {$month = "0"."$month";}
        #print "$year-$month-$day $hour:$min:$sec\n";
        return("$year-$month-$day-$hour-$min-$sec");
}

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);
}