Mercurial > repos > big-tiandm > mirplant2
changeset 11:861eb4fbfbd0 draft
Uploaded
author | big-tiandm |
---|---|
date | Fri, 25 Jul 2014 05:20:44 -0400 |
parents | 5ac42e277c76 |
children | dc5a29826c7d |
files | matching.pl |
diffstat | 1 files changed, 90 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/matching.pl Fri Jul 25 05:20:44 2014 -0400 @@ -0,0 +1,90 @@ +#!/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","time: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(); +if (defined $opts{'time'}) { + $time=$opts{'time'}; +} + +my $mapdir=$fileout."/genome_match_".$time; +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`; + +#`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); +} +