|
37
|
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 use File::Basename;
|
|
|
13
|
|
|
14 my %opts;
|
|
47
|
15 GetOptions(\%opts,"i=s","ref=s","index:s","v:i","p:i","o=s","h");
|
|
37
|
16 if (!(defined $opts{i} and defined $opts{o} ) || defined $opts{h}) { #necessary arguments
|
|
|
17 &usage;
|
|
|
18 }
|
|
|
19
|
|
|
20 my $filein=$opts{'i'};
|
|
|
21 my $dir=$opts{'o'};
|
|
|
22 unless ($dir=~/\/$/) {$dir.="/";}
|
|
|
23 my $rfam=$opts{'ref'};
|
|
|
24 my $mis=defined $opts{'v'}? $opts{'v'} : 0;
|
|
|
25 my $index=defined $opts{'index'} ? $opts{'index'} : "";
|
|
|
26 my $threads=defined $opts{'p'} ? $opts{'p'} : 1;
|
|
|
27
|
|
|
28 if (not -d $dir) {
|
|
|
29 mkdir $dir;
|
|
|
30 }
|
|
|
31
|
|
|
32
|
|
|
33 my $time=Time();
|
|
47
|
34
|
|
|
35 my $mapdir=$dir."/rfam_match";
|
|
37
|
36 if(not -d $mapdir){
|
|
|
37 mkdir $mapdir;
|
|
|
38 }
|
|
|
39 chdir $mapdir;
|
|
|
40 ###check genome index
|
|
|
41 if (-s $index.".1.ebwt") {
|
|
|
42 }else{
|
|
|
43 &checkACGT($rfam);
|
|
|
44 `bowtie-build $rfam $rfam`;
|
|
|
45 $index="$rfam";
|
|
|
46 }
|
|
|
47 ### genome mapping
|
|
44
|
48 `bowtie -v $mis -f -p $threads -k 1 $index $filein --al rfam_mapped.fa --un rfam_not_mapped.fa > rfam_mapped.bwt 2> run.log`;
|
|
37
|
49
|
|
|
50 sub checkACGT{
|
|
|
51 my $string;
|
|
|
52 open IN,"<$rfam";
|
|
|
53 while (my $aline=<IN>) {
|
|
|
54 if ($aline!~/^>/) {
|
|
|
55 $aline=~s/U/T/gi;
|
|
|
56 }
|
|
|
57 $string .=$aline;
|
|
|
58 }
|
|
|
59 close IN;
|
|
|
60 $rfam=basename($rfam);
|
|
|
61 open OUT, ">$rfam";
|
|
|
62 print OUT $string;
|
|
|
63 close OUT;
|
|
|
64 }
|
|
|
65
|
|
|
66 sub Time{
|
|
|
67 my $time=time();
|
|
|
68 my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
|
|
|
69 $month++;
|
|
|
70 $year+=1900;
|
|
|
71 if (length($sec) == 1) {$sec = "0"."$sec";}
|
|
|
72 if (length($min) == 1) {$min = "0"."$min";}
|
|
|
73 if (length($hour) == 1) {$hour = "0"."$hour";}
|
|
|
74 if (length($day) == 1) {$day = "0"."$day";}
|
|
|
75 if (length($month) == 1) {$month = "0"."$month";}
|
|
|
76 #print "$year-$month-$day $hour:$min:$sec\n";
|
|
|
77 return("$year-$month-$day-$hour-$min-$sec");
|
|
|
78 }
|
|
|
79 sub usage{
|
|
|
80 print <<"USAGE";
|
|
|
81 Version $version
|
|
|
82 Usage:
|
|
|
83 $0 -i -o
|
|
|
84 options:
|
|
|
85 -i input file# input reads fasta/fastq file
|
|
|
86 -ref input file# rfam file, which do not contain miRNAs
|
|
|
87 -index file-prefix #(must be indexed by bowtie-build) The parameter
|
|
|
88 string must be the prefix of the bowtie index. For instance, if
|
|
|
89 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
|
|
|
90 the prefix is 'h_sapiens_37_asm'.##can be null
|
|
|
91 -v <int> report end-to-end hits w/ <=v mismatches; ignore qualities,default 0;
|
|
|
92
|
|
|
93 -p/--threads <int> number of alignment threads to launch (default: 1)
|
|
|
94
|
|
|
95 -o output directory
|
|
|
96 -h help
|
|
|
97 USAGE
|
|
|
98 exit(1);
|
|
|
99 }
|
|
|
100
|
