|
21
|
1 #!/usr/bin/perl -w
|
|
|
2 #Filename:
|
|
|
3 #Author: Tian Dongmei
|
|
|
4 #Email: tiandm@big.ac.cn
|
|
|
5 #Date: 2013/7/19
|
|
|
6 #Modified:
|
|
|
7 #Description:
|
|
|
8 my $version=1.00;
|
|
|
9
|
|
|
10 use strict;
|
|
|
11 use Getopt::Long;
|
|
|
12
|
|
|
13 my %opts;
|
|
|
14 GetOptions(\%opts,"i=s","g=s","index:s","v:i","p:i","r:s","o=s","h");
|
|
|
15 if (!(defined $opts{i} and defined $opts{o} ) || defined $opts{h}) { #necessary arguments
|
|
|
16 &usage;
|
|
|
17 }
|
|
|
18
|
|
|
19 my $filein=$opts{'i'};
|
|
|
20 my $fileout=$opts{'o'};
|
|
|
21 unless ($fileout=~/\/$/) {$fileout.="/";}
|
|
|
22 my $genome=$opts{'g'};
|
|
|
23 my $mis=defined $opts{'v'}? $opts{'v'} : 0;
|
|
|
24 my $hits=defined $opts{'r'}? $opts{'r'} : 25;
|
|
|
25 my $index=defined $opts{'index'} ? $opts{'index'} : "";
|
|
|
26 my $threads=defined $opts{'p'} ? $opts{'p'} : 1;
|
|
|
27
|
|
|
28
|
|
|
29 #my $time=time();
|
|
|
30 #my $mapdir=$fileout."/genome_match_".$time;
|
|
|
31 my $mapdir=$fileout."/genome_match";
|
|
|
32 mkdir $mapdir;
|
|
|
33 chdir $mapdir;
|
|
|
34 ###check genome index
|
|
|
35 if (-s $index.".1.ebwt") {
|
|
|
36 }else{
|
|
|
37 `bowtie-build $genome genome`;
|
|
|
38 $index="genome";
|
|
|
39 }
|
|
|
40
|
|
|
41 ### genome mapping
|
|
|
42 `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`;
|
|
|
43
|
|
|
44 #`convert_bowtie_to_blast.pl genome_mapped.bwt genome_mapped.fa $genome > genome_mapped.bst`;
|
|
|
45
|
|
|
46 sub usage{
|
|
|
47 print <<"USAGE";
|
|
|
48 Version $version
|
|
|
49 Usage:
|
|
|
50 $0 -i -o
|
|
|
51 options:
|
|
|
52 -i input file# input reads fasta/fastq file
|
|
|
53 -g input file# genome file
|
|
|
54 -index file-prefix #(must be indexed by bowtie-build) The parameter
|
|
|
55 string must be the prefix of the bowtie index. For instance, if
|
|
|
56 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
|
|
|
57 the prefix is 'h_sapiens_37_asm'.##can be null
|
|
|
58 -v <int> report end-to-end hits w/ <=v mismatches; ignore qualities,default 0;
|
|
|
59
|
|
|
60 -p/--threads <int> number of alignment threads to launch (default: 1)
|
|
|
61
|
|
|
62 -r int a read is allowed to map up to this number of positions in the genome
|
|
|
63 default is 25
|
|
|
64
|
|
|
65 -o output directory
|
|
|
66
|
|
|
67 -h help
|
|
|
68 USAGE
|
|
|
69 exit(1);
|
|
|
70 }
|
|
|
71
|
