Mercurial > repos > mingchen0919 > split_multifasta
changeset 0:efd5c022b54d draft
planemo upload
| author | mingchen0919 |
|---|---|
| date | Mon, 09 Apr 2018 12:27:49 -0400 |
| parents | |
| children | 8ef62ca3938b |
| files | getopt_specification.csv helper.R split.pl split_fasta.Rmd split_fasta.sh split_fasta.xml split_fasta_render.R split_multifasta.pl |
| diffstat | 8 files changed, 658 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/getopt_specification.csv Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,7 @@ +short flag,argument mask,data type,variable name +o,1,character,report +d,1,character,report.files_path +s,1,character,sink_message +A,1,character,fasta_input +B,1,character,number +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/helper.R Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,28 @@ +#' \code{getopt_specification_matrix} returns a getopt specification matrix. +#' +#' @param specification_file a cvs file within the \code{galaxy_tool_directory} which stores getopt specification matrix data. +#' The first column are short flags, the second column are argument masks, the third column +#' is data types. The fourth column are variable names used in the tool XML. These three columns are required. +#' @param gtg_name the name of a running GTG. +getopt_specification_matrix = function(specification_file, gtg_name = 'gtg', tool_dir = Sys.getenv('TOOL_DIR')) { + df = read.csv(paste0(tool_dir, '/', specification_file), + header = TRUE, stringsAsFactors = FALSE) + # check if there are duplicated short flags + short_flags = df[, 1] + if (length(unique(short_flags)) < length(short_flags)) { + cat('----Duplicated short flags found ----\n') + cat('short flags: ', df[, 1][duplicated(df[, 1])], '\n') + stop('Duplicated short flags are not allowed.') + } + + # use short flags to generate long flags + long_flags = paste0('X_', df[, 1]) + + # specification matrix + df2 = data.frame(long_flags = long_flags, + short_flags = df[, 1], + argument_mask = df[, 2], + data_type = df[, 3]) + + as.matrix(df2) +} \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split.pl Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,52 @@ +#!/usr/bin/perl + + + if($ARGV[0] eq "" || $ARGV[1] eq ""){ + die "\n\t Usage : perl <thisScript.pl> <file to be split> <number of partitions> \n\n"; + } + + + $homfile = $ARGV[0]; + $numOfFiles = $ARGV[1]; + + + system("grep -c '^>' $homfile > out"); + open IN, "out" || die "File not found - 2\n"; + $numOfSeqs = <IN>; + close IN; + + print "Number of seqs is $numOfSeqs\n"; + my $numPerFile = $numOfSeqs/$numOfFiles; + print "Num per File is $numPerFile\n"; + + open IN, $homfile || die "File not found - 1\n"; + $lineIn = <IN>; + + for($i = 1; $i <= $numOfFiles; $i++){ + print "$i\n"; + open FILE, ">".$homfile.".".$i || die "Can't open file"; + print FILE $lineIn; + $seqs = 1; + $lineIn = <IN>; + while(defined $lineIn && $seqs < $numPerFile){ + print FILE $lineIn; + if ($lineIn =~ /^>/) { $seqs++; } + $lineIn = <IN>; + } + while(defined $lineIn && $lineIn !~ /^>/){ + print FILE $lineIn; + $lineIn = <IN>; + } + close FILE; + } + $i = $i -1; + open FILE, ">>".$homfile.".".$i; + while ($lineIn = <IN>){ + print FILE $lineIn; + } + close FILE; + + close IN; + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split_fasta.Rmd Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,47 @@ +--- +title: 'FASTA splitter' +output: + html_document: + highlight: pygments +--- + +```{r setup, include=FALSE, warning=FALSE, message=FALSE} +knitr::opts_chunk$set(echo = TRUE, error = TRUE) +``` + + +```{bash} +# build job-script +mkdir -p ${WORKING_DIR}/fasta_files + +# single-end.sh +cat <<EOF >${X_d}/job-script.sh +${X_t}/split_multifasta.pl \\ + --input_file=${X_A} \\ + --seqs_per_file=${X_B} \\ + --output_dir=${WORKING_DIR}/fasta_files > ${X_d}/fasta_splitter-log.txt 2>&1 +EOF +``` + +```{bash, 'run jobs', echo=FALSE} +# run job script, always use absolute path. +# we want to run all jobs within the working path. +sh ${X_d}/job-script.sh +``` + +```{r, 'display output directory contents', results='asis', echo=FALSE} +## after the job is done, we list all files from the output directory. +## full relative path to the output directory needs to be displayed. + +cat('##All output files') +cat('\n\n') +all_files = list.files(path = opt$X_d, + full.names = TRUE, + recursive = TRUE) + +for (f in sub(opt$X_d, '.', all_files) ) { + cat('* [', f, '](', f, ')\n') +} +cat('\n') +``` +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split_fasta.sh Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,9 @@ +export TOOL_DIR='${__tool_directory__}' && + +Rscript '${__tool_directory__}/'split_fasta_render.R + + -o '$report' + -d '$report.files_path' + -s '$sink_message' + -A '$fasta_input' + -B '$number'
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split_fasta.xml Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,41 @@ +<tool id="statonlab_fasta_splitter" name="statonlab_fasta_splitter" version="1.0.0"> + <description>Split a single FASTA file into multiple smaller FASTA files
 + </description> + <requirements> + <requirement type="package" version="1.15.0.6-0">pandoc</requirement><requirement type="package" version="1.20.0">r-getopt</requirement><requirement type="package" version="1.6">r-rmarkdown</requirement><requirement type="package" version="5.22.2.1">perl</requirement></requirements> + <stdio> + <regex match="XXX" source="stderr" level="warning" description="Check the warnings_and_errors.txt file for more details."/></stdio> + <command><![CDATA[export TOOL_DIR='${__tool_directory__}' && + +Rscript '${__tool_directory__}/'split_fasta_render.R + + -o '$report' + -d '$report.files_path' + -s '$sink_message' + -A '$fasta_input' + -B '$number' +]]></command> + <inputs> + <param type="data" name="fasta_input" label="FASTA input" optional="False" multiple="False"/><param type="integer" name="number" label="Number" help="The number of files the FASTA file will be split into" optional="False" value="200" min="1"/></inputs> + <outputs> + <data name="report" format="html" label="${tool.name} report on ${on_string}" hidden="false"/><data name="sink_message" format="txt" label="log on ${on_string}" from_work_dir="warnings_and_errors.txt" hidden="false"/><collection name="list_collection" type="list" label="${tool.name} on ${on_string}"><discover_datasets pattern="__name_and_ext__" directory="fasta_files" visible="true"/></collection></outputs> + <citations> + <citation type="bibtex"><![CDATA[ + @article{allaire2016rmarkdown, + title={rmarkdown: Dynamic Documents for R, 2016}, + author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff + and Wickham, Hadley and Atkins, Aron and Hyndman, Rob}, + journal={R package version 0.9}, + volume={6}, + year={2016} + } + ]]></citation><citation type="bibtex"><![CDATA[ + @book{xie2015dynamic, + title={Dynamic Documents with R and knitr}, + author={Xie, Yihui}, + volume={29}, + year={2015}, + publisher={CRC Press} + } + ]]></citation></citations> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split_fasta_render.R Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,52 @@ +##============ Sink warnings and errors to a file ============== +## use the sink() function to wrap all code within it. +##============================================================== +zz = file('warnings_and_errors.txt') +sink(zz) +sink(zz, type = 'message') + +#------------import libraries-------------------- +options(stringsAsFactors = FALSE) + +library(getopt) +library(rmarkdown) +#------------------------------------------------ + + +#------------get arguments into R-------------------- +# load helper function +source(paste0(Sys.getenv('TOOL_DIR'), '/helper.R')) +# import getopt specification matrix from a csv file +opt = getopt(getopt_specification_matrix('getopt_specification.csv')) +opt$X_t = Sys.getenv('TOOL_DIR') +working_dir = getwd() +Sys.setenv(WORKING_DIR = working_dir) +#---------------------------------------------------- + + +#-----------using passed arguments in R +# to define system environment variables--- +do.call(Sys.setenv, opt[-1]) +#---------------------------------------------------- + +#---------- often used variables ---------------- +# OUTPUT_DIR: path to the output associated directory, which stores all outputs +# TOOL_DIR: path to the tool installation directory +OUTPUT_DIR = opt$X_d +TOOL_DIR = opt$X_t +OUTPUT_REPORT = opt$X_o +RMD_NAME = 'split_fasta.Rmd' + +# create the output associated directory to store all outputs +dir.create(OUTPUT_DIR, recursive = TRUE) + +#-----------------render Rmd-------------- +render(paste0(TOOL_DIR, '/', RMD_NAME), output_file = OUTPUT_REPORT) +#------------------------------------------ + +#==============the end============== + + +##--------end of code rendering .Rmd templates---------------- +sink() +##=========== End of sinking output============================= \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/split_multifasta.pl Mon Apr 09 12:27:49 2018 -0400 @@ -0,0 +1,422 @@ +#!/usr/bin/perl + +#BEGIN{foreach (@INC) {s/\/usr\/local\/packages/\/local\/platform/}}; +#use lib (@INC,$ENV{"PERL_MOD_DIR"}); +#no lib "$ENV{PERL_MOD_DIR}/i686-linux"; +#no lib "."; + +=head1 NAME + +split_multifasta.pl - split a single FASTA file containing multiple sequences into separate files. + +=head1 SYNOPSIS + +USAGE: split_multifasta.pl + --input_file=/path/to/some_file.fsa + --output_dir=/path/to/somedir + [ --output_list=/path/to/somefile.list + --output_subdir_size=1000 + --output_subdir_prefix=fasta + --seqs_per_file=1 + --compress_output=1 + ] + + + split_multifasta.pl --in snapdmel.aa --output_dir=./ --f=snaa --seqs_per_file=1000 + +=head1 OPTIONS + +B<--input_file,-i> + The input multi-fasta file to split. + +B<--output_dir,-o> + The directory to which the output files will be written. + +B<--output_list,-s> + Write a list file containing the paths of each of the regular output files. This may be useful + for later scripts that can accept a list as input. + +B<--output_file_prefix,-f> + If defined, each file created will have this string prepended to its name. This is ignored unless + writing multiple sequences to each output file using the --seqs_per_file option with a value greater + than 1, else each file created will just be a number. + +B<--output_subdir_size,-u> + If defined, this script will create numbered subdirectories in the output directory, each + containing this many sequences files. Once this limit is reached, another subdirectory + is created. + +B<--output_subdir_prefix,-p> + To be used along with --output_subdir_size, this allows more control of the names of the + subdirectories created. Rather than just incrementing numbers (like 10), each subdirectory + will be named with this prefix (like prefix10). + +B<--compress_output,-c> + Output fasta files will be gzipped when written. + +B<--debug,-d> + Debug level. Use a large number to turn on verbose debugging. + +B<--log,-l> + Log file + +B<--help,-h> + This help message + +=head1 DESCRIPTION + +This script is used to split a single FASTA file containing multiple sequences into separate +files containing one sequence each. + +=head1 INPUT + +The input is defined with --input_file and should be a single fasta file. File extensions are +ignored. When creating this multi-entry FASTA file, one should take care to make the first +*word* after the > symbol a unique value, as it will be used as the file name for that sequence. +For example: + + >gi53791237 Tragulus javanicus p97bcnt gene for p97Bcnt + ACAGGAGAAGAGACTGAAGAGACACGTTCAGGAGAAGAGCAAGAGAAGCCTAAAGAAATGCAAGAAGTTA + AACTCACCAAATCACTTGTTGAAGAAGTCAGGTAACATGACATTCACAAACTTCAAAACTAGTTCTTTAA + AAAGGAACATCTCTCTTTTAATATGTATGCATTATTAATTTATTTACTCATTGGCGTGGAGGAGGAAATG + + >gi15387669 Corynebacterium callunae pCC1 plasmid + ATGCATGCTAGTGTGGTGAGTATGAGCACACACATTCATGGGCACCGCCGGGGTGCAGGGGGGCTTGCCC + CTTGTCCATGCGGGGTGTGGGGCTTGCCCCGCCGATAGAGACCGGCCACCACCATGGCACCCGGTCGCGG + GGTGATCGGCCACCACCACCGCCCCCGGCCACTCTCCCCCTGTCTAGGCCATATTTCAGGCCGTCCACTG + +Whitespace is ignored within the input file. See the OUTPUT section for more on creation of +output files. + +=head1 OUTPUT + +The name of each output sequence file is pulled from the FASTA header of that sequence. The +first *word* after the > symbol will be used as the file name, along with the extension .fsa. +The word is defined as all the text after the > symbol up to the first whitespace. + +If the above example were your input file, two files would be created: + + gi53791237.fsa + gi15387669.fsa + +Any characters other than a-z A-Z 0-9 . _ - in the ID will be changed into an +underscore. This only occurs in the file name; the original FASTA header within the file +will be unmodified. + +You can pass a path to the optional --output_list to create a text file containing the full paths +to each of the FASTA files created by this script. + +Two other optional arguments, --output_subdir_size and --output_subdir_prefix, can be used +on input sets that are too large to write out to one directory. This depends on the limitations +of your file system, but you usually don't want 100,000 files written in the same directory. + +If you have an FASTA file containing 95000 sequences, and use the following option: + + --output_dir=/some/path + --output_subdir_size=30000 + +The following will be created: + + directory file count + --------------------------------- + /some/path/1/ 30000 + /some/path/2/ 30000 + /some/path/3/ 30000 + /some/path/4/ 5000 + +If you choose to create a list file (and you probably want to), it will contain these proper paths. + +You may not want the subdirectories to simply be numbers, as above, so you can use the +--output_subdir_prefix option. For example: + + --output_dir=/some/path + --output_subdir_size=30000 + --output_subdir_prefix=fasta + +The following will be created: + + directory file count + --------------------------------- + /some/path/fasta1/ 30000 + /some/path/fasta2/ 30000 + /some/path/fasta3/ 30000 + /some/path/fasta4/ 5000 + +Finally, you can write multiple sequences to each output file using the --seqs_per_file option, which +can be used along with --outupt_subdir_size and --output_subdir_prefix. The main difference to note +is that, if you use --seqs_per_file, the fasta file created will no longer be named using values +taken from the header, since it will contain multiple headers. Instead, the file will simply be +named using sequential numbers starting at 1 (like 1.fsa). For example: + + --output_dir=/some/path + --output_subdir_size=3000 + --output_subdir_prefix=fasta + --seqs_per_file=10 + +The following will be created: + + directory file count + --------------------------------- + /some/path/fasta1/ 3000 + /some/path/fasta2/ 3000 + /some/path/fasta3/ 3000 + /some/path/fasta4/ 500 + +=head1 CONTACT + + Joshua Orvis + jorvis@tigr.org + +=cut + +use strict; +use Getopt::Long; +# qw(:config no_ignore_case no_auto_abbrev pass_through); +use Pod::Usage; +# BEGIN { +# use Ergatis::Logger; +# } + +my %options = (); +my $results = GetOptions (\%options, + 'input_file|i=s', + 'output_dir|o=s', + 'output_file_prefix|f=s', + 'output_list|s=s', + 'output_subdir_size|u=s', + 'output_subdir_prefix|p=s', + 'seqs_per_file|n|e=s', + 'compress_output|c=s', + 'log|l=s', + 'debug=s', + 'help|h') || pod2usage(); + +# my $logfile = $options{'log'} || Ergatis::Logger::get_default_logfilename(); +# my $logger = new Ergatis::Logger('LOG_FILE'=>$logfile, +# 'LOG_LEVEL'=>$options{'debug'}); +# $logger = $logger->get_logger(); + + +my $logfile = $options{'log'} || "log.file"; +my $logger = new logger('LOG_FILE'=>$logfile, + 'LOG_LEVEL'=>$options{'debug'}); + +## display documentation +if( $options{'help'} ){ + pod2usage( {-exitval => 0, -verbose => 2, -output => \*STDERR} ); +} + +## make sure everything passed was peachy +&check_parameters(\%options); + +## open the list file if one was passed +my $listfh; +if (defined $options{output_list}) { + open($listfh, ">$options{output_list}") || $logger->logdie("couldn't create $options{output_list} list file"); +} + +my $first = 1; +my $seq = ''; +my $header; + +my $sfh; + +## load the sequence file +if ($options{'input_file'} =~ /\.(gz|gzip)$/) { + open ($sfh, "<:gzip", $options{'input_file'}) + || $logger->logdie("can't open sequence file:\n$!"); +} else { + open ($sfh, "<$options{'input_file'}") + || $logger->logdie("can't open sequence file:\n$!"); +} + +my $sub_dir = 1; +my $seq_file_count = 0; + +## keep track of how many sequences are in the current output file +my $seqs_in_file = 0; +my $group_filename_prefix = 1; + +## holds the output file handle +my $ofh; + +while (<$sfh>) { + ## if we find a header line ... + if (/^\>(.*)/) { + + ## write the previous sequence before continuing with this one + unless ($first) { + &writeSequence(\$header, \$seq); + + ## reset the sequence + $seq = ''; + } + + $first = 0; + $header = $1; + + ## else we've found a sequence line + } else { + ## skip it if it is just whitespace + next if (/^\s*$/); + + ## record this portion of the sequence + $seq .= $_; + } +} + +## don't forget the last sequence +&writeSequence(\$header, \$seq); + +exit; + +sub check_parameters { + my $options = shift; + + ## make sure input_file and output_dir were passed + unless ( $options{input_file} && $options{output_dir} ) { + $logger->logdie("You must pass both --input_file and --output_dir"); + } + + ## make sure input_file exists + if (! -e $options{input_file} ) { + if ( -e "$options{input_file}.gz" ) { + $options{input_file} .= '.gz'; + } else { + $logger->logdie("the input file passed ($options{input_file}) cannot be read or does not exist"); + } + } + + ## make sure the output_dir exists + if (! -e "$options{output_dir}") { + $logger->logdie("the output directory passed could not be read or does not exist"); + } + + ## seqs_per_file, if passed, must be at least one + if (defined $options{seqs_per_file} && $options{seqs_per_file} < 1) { + $logger->logdie("seq_per_file setting cannot be less than one"); + } + + ## handle some defaults + $options{output_subdir_size} = 0 unless ($options{output_subdir_size}); + $options{output_subdir_prefix} = '' unless ($options{output_subdir_prefix}); + $options{seqs_per_file} = 1 unless ($options{seqs_per_file}); + $options{output_file_prefix} = '' unless ($options{output_file_prefix}); +} + +sub writeSequence { + my ($header, $seq) = @_; + + ## the id used to write the output file will be the first thing + ## in the header up to the first whitespace. get that. + $$header =~ /^(\S+)/ || $logger->logdie( "can't pull out an id on header $$header" ); + my $id = $1; + + ## because it is going to be the filename, we're going to take out the characters that are bad form to use + ## legal characters = a-z A-Z 0-9 - . _ + $id =~ s/[^a-z0-9\-_.]/_/gi; + + my $dirpath; + + ## if we're writing more than one sequence to a file, change the id from + ## fasta header to the current group file name + if ($options{seqs_per_file} > 1) { + $id = $group_filename_prefix; + + ## did the user ask for a file prefix? + if ( $options{output_file_prefix} ) { + $id = $options{output_file_prefix} . $id; + } + } + + + ## the path depends on whether we are using output subdirectories + if ($options{output_subdir_size}) { + $dirpath = "$options{'output_dir'}/$options{output_subdir_prefix}$sub_dir"; + } else { + $dirpath = "$options{'output_dir'}"; + } + + ## did the user ask for a file prefix? + my $filepath = "$dirpath/$id.fsa"; + + ## take any // out of the filepath + $filepath =~ s|/+|/|g; + + ## write the sequence + $logger->debug("Writing sequence to $filepath") if ($logger->is_debug()); + + ## open a new output file if we need to + ## if we're writing multiple sequences per file, we only open a new + ## one when $seqs_in_file = 0 (first sequence) + if ($seqs_in_file == 0) { + + ## if the directory we want to write to doesn't exist yet, create it + mkdir($dirpath) unless (-e $dirpath); + + + if ($options{'compress_output'}) { + open ($ofh, ">:gzip", $filepath.".gz") + || $logger->logdie("can't create '$filepath.gz':\n$!"); + } else { + open ($ofh, ">$filepath") || $logger->logdie("can't create '$filepath':\n$!"); + + } + $seq_file_count++; + + ## add the file we just wrote to the list, if we were asked to + if (defined $options{output_list}) { + print $listfh "$filepath\n"; + } + } + + ## if we're doing output subdirs and hit our size limit, increment to the next dir + if ($options{output_subdir_size} && $options{output_subdir_size} == $seq_file_count) { + $seq_file_count = 0; + $sub_dir++; + } + + ## write the sequence + print $ofh ">$$header\n$$seq\n"; + $seqs_in_file++; + + ## if we hit the limit of how many we want in each file, set the next file name and + ## reset the count of seqs within the file + if ($options{seqs_per_file} == $seqs_in_file) { + $seqs_in_file = 0; + $group_filename_prefix++; + } +} + + +package logger; + +sub new { + my $packname= shift; + my %args= @_; + my $self= \%args; + bless($self,$packname); + return $self; +} + +sub get_logger { + my $self= shift; + return $self; +} + +sub logdie { + my $self= shift; + die @_; +} + +sub debug { + my $self= shift; + warn @_; +} + +sub is_debug { + shift->{LOG_LEVEL} || 0; +} + + +1;