Linxe's eye

I needed to test some sequence for known TFBS, I was looking for: (1) a public database of TFBS and (2) a program to search the motifs in a sequence(s). First I found Jaspar as an alternative to the commercial TransFac DB , I downloaded the matrices for vertebrate. The second problem was to find a program to search the motifs in any sequence, I did not find any updated program (many are published, but few websites are active nowadays). So, I decided to write my own program, I share the code if anyone can use it: #!/usr/bin/perl -w use strict; use Getopt::Long; =head1 NAME jasparScan.pl =head1 USAGE jasparScan.pl -f FASTA -m MATRIX [-s MODELS] [-c CUT_OFF] OPTIONS: -h --help Help screen -f --fasta Fasta file to scan -m --matrix Jaspar motif matrix -s --models List of models, separated by ':' [Def: all] -c --cutoff Similarity cut-off [Def: 0.90] =head1 DESCRIPTION Script to map a TF binding site matrix (from Jaspar DB) into a fasta sequence. The simila...

Sometimes I need to sample large data sets, so I randomly select some lines in the file. My files are generally text-records, one record by line, then I wrote this small script to do the task: #!/usr/bin/perl -w use strict; =head1 NAME selectRandomLines.pl =head1 DESCRIPTION Select random lines in a file. =cut $ARGV[2] or die "Usage: selectRandomLines.pl TOTAL_LINES_IN_FILE NUM_LINES_WANTED FILE_NAME\n"; my $total = shift @ARGV; # Total lines in the file my $want = shift @ARGV; # Total lines to select my $file = shift @ARGV; # The file my $line = 0; # Line counter my %select = randomSelect($total, $want); # Hash with selected lines open FILE, "$file" or die "cannot open $file\n"; while (<FILE>) { print "$_" if (defined $select{$line++}); } close FILE; =head1 SUBROUTINES randomSelect() CALL: randomSelect(TOTAL_ELEMENTS, TOTAL_WANTED) [NUM, NUM] RETURN: %s [HASH] =cut sub randomSelect { my $t = shift @_; my $n = shift...

I was coding to query a SQLite database using Perl::DBI, I have a small code like this: #!/usr/bin/perl -w use strict; use DBI; my $dbh = DBI->connect("dbi:SQLite:dbname=db_file", "", "") or die "error in connection\n"; my $sth = $dbh->prepare($sql) or die "cannot prepare SQL\n"; $sth->execute(); while (my @data = $sth->fetchrow_array()) { #process data } $sth->finish; $dbh->disconnect; But, every time I execute this script I obtain this warning: closing dbh with active statement handles .... blah, blah, blah I double checked the code, the documentation and ran some debug examples, I obtained the expected result, so why this message? Finally, I found the solution in the PerlMonks website, the problem is a bad status return in the module when you use close() method, so the simple solution is to " undef $sth " at the end: $sth->finish; undef $sth; $dbh->disconnect;

Last day, M asks me how to compact the microRNA mature fasta-file from miRbase , the problem is many miRNA are conservated across species, so the fasta file with all the mature sequences has a lot of redundancy and he wants unique sequences to screen some sequences. This is a simple Perl solution, the first step is to read and collect the sequences, the trick is to save in a hash using the nucleic sequences as "key" and append the names as "value". After this, simply extract the "keys" and parse a little the "values" to have a common identifier. Code: #!/usr/bin/perl -w use strict; =head1 NAME uniq_mir.pl FASTA =head1 DESCRIPTION Filter the multifasta file of miRbase (mature.fa) into unique miR by sequence. =cut $ARGV[0] or die "Usage: uniq_mir.pl FASTA\n"; my $name = ''; my %seqs = (); open F, "$ARGV[0]" or die "cannot open $ARGV[0]\n"; while ( <F>) { chomp; if (/^>/) { s/>//; my @line = ...

I know it! Original in XKCD: http://xkcd.com/224/ What's your favorite programmer computer cartoon?: http://stackoverflow.com/questions/84556/whats-your-favorite-programmer-cartoon

#!/usr/bin/perl -w use strict; =head1 NAME oligoGenerator.pl =head1 DESCRIPTION Perl script to generate all possible combinations of size k using an alphabet @a, we use function recursion. My intention is to create all possible oligonucleotides (DNA alphabet or ACGT) but can be extended to any other field using a different alphabet. Output also can be printed in other forms, you can put other delimiters in the push function or in the final array printed. =cut my ($k) = 4; # definition of the word size my @a = qw/A C G T/; # definition of the alphabet my @words = createWords($k, @a);# main function print join("\n", @words), "\n"; # print output sub createWords { my $k = shift @_; $k--; my @old = @_; my @new = (); if ($k < 1) { return @old; } else { foreach my $e (@old) { foreach my $n (@a) { push @new, "$e$n"; # add new element } } createWords($k, @new); # recursion c...

Last week my good friend M asks me to create a script to draw some sequences, the main problem is to visually see differences in transcript orientation (sense and antisense), some time ago I created similar tools for mapping short sequences like 454 pyro or siRNAs. Now the big problem is to have a good assembly, many regions can extend a "contig" to the left or the right, so internal coordinates change every time you add a new element, with short sequences you have a target well defined and the extension is relative to it. After thinking a little, I decide not extend my code to the assembly (yes, I'm lazy) and use Phrap to perform the job. A parser read the output and extract the alignment information. Second problem was the fasta naming convention, the original fasta file is a mixed names of other sequences, so I decided to create and tag an unique ID, "seq_##". Extending the problem, M asks me to be possible to add more sequences, based in the sequence homolo...

Do you know what is a Perl Golf problem? It's a general problem formulated and you try to resolve with a minimal number of characters in a perl script, who writes less win. Some times is a good habit to see, admire and think in this beautiful pearls. Commonly there are a lot in the Perl Monks website. Today I was looking for a more simple and effective subroutine to translate a DNA/RNA sequence into the corresponding peptide version using the typical genetic code, I used the typical solution with a hash storing the code and call the sequence in block with substr or pop/shift. I found this solutions in a Perl Golf challenge : # Typical solution hashing the codes: sub f0 { #by tadman my %g = ( # . - Stop 'UAA'=>'.','UAG'=>'.','UGA'=>'.', # A - Alanine 'GCU'=>'A','GCC'=>'A','GCA'=>'A','GCG'=>'A', # C - Cysteine '...

I have used Perl for my projects for more than 6 years, when I used the Arabidopsis genome it was easy to load a full chromosome into a string variable, but now I work mining bigger genomes, like human or mouse, and I had the habit to pass direct variable or references to subroutines, bad idea. The point is when you are using many subs to perform calculations inside the sequences, the memory increase exponentially, one of my scripts takes the sequences, divide into blocks and perform calculation by each block, and store in a hash ( O-O programming), first time I made a mistake and create an infinite loop and sucks the entire memory, including the SWAP, our server has 32 GB in RAM and 32 GB in SWAP, fortunately the server survive to my error. Next I fixed it, and now the script require about 3 GB , so much ... After think a little I decide to not pass nothing to the subs, declare globally access for principal variables and this reduces the memory needed to 300 MB for the same proc...

Yesterday I was disturbed of my coding-state to answer a basic problem, which's the best font to display a sequences alignment? First introduce the concept of sequence alignment, if you want to compare 2 or more sequences in their composition of nucleic bases (for DNA or RNA) or amino acid (for peptids ), you must search similar blocks between all the context and try to fix the rest as possible. Many algorithms has been created to perform this task, with different approaches, if you want to know more, please go to Wikipedia or bioinformatics books. Ok , let's consider we have our alignment, really this is a text file with a special structure to easy visualize the blocks found. Some people make the mistake to cut-and-paste in a text document without change the font style, a font is just an image to represent letters or symbols, some are very artistic but generally the font does not conserve a proportional aspect in size, so the beautiful alignment changes its proportion and vi...

These days I had working on web interfaces for some Perl scripts, the server is an 8-cores Linux box with CentOS with Apache as web server, but some tests didn't work well and some probably fall into an infinite loops. The problem was I can not kill these process, because I must not know the admin password and I can not have access for the web server, the scripts are created in my home and when activate the process is submit by the apache user. So the solution was a script which ask the PID and send a termination signal. Here is the code, maybe someone can use it: #!/tools64/bin/perl -w use strict; =head1 NAME reaperl.pl =head1 DESCRIPTION Web service to stop some uncontrolled process created by the Perl::CGI. Simply presents a web form to indicate the PID and call a system termination. =cut use CGI::Pretty qw/:standard/; print header; print start_html(-title=>'reaperl.pl'); print start_form; print p('Kill this process: ', textfield(-name=>'proc')); pr...