Hi everyone,
I have a file with many predicted miRNAs. I need to perform a randomization test to identify which of these miRNAs are highly probable. This test is to randomize each predicted miRNA 1000 times and calculate each randomized sequence's MFE value (this can be easily done by RNAfold). My current problem is how to generate 1000 radomized sequences for each miRNA? I make an example below.
Predicted miRNAs file
>miR1
augcgugaccguaugcuac
>miR2
uuuggugcguagucguacg
............
>miR100
auaugagucguacguacgu
Radomized sequences file
>miR1_1
ugcggaccguaugcuacau
>miR1_2
ugcggaccuugcuacauga
...........
>miR1_1000
ugcggaccguaugcuacua
.............
............
>miR100_1000
auaugagucgacguacu
Could anyone being familar with perl help to solve this problem? For the next steps including calculating MFE etc, I can do them myself. But I believe someone who know well RNAfold and miRNA prediction can produce a pipline for this work. I attached a Nucleic Acid Research reference link here http://nar.oxfordjournals.org/content/37/suppl_1/D111.full . In the method section, when searing RNAFOLD, you can find the authors' method. THANK YOU in advance!
2 answers
Just a quick advice. You should use the pre-miRNA sequence (the typical hairpin-like structure) and not the mature miRNA to fold the secondary structure. You can use Randfold ( http://bioinformatics.oxfordjournals.org/content/20/17/2911 ) to do that.
download Randfold : http://bioinformatics.psb.ugent.be/supplementary_data/erbon/nov2003/
#!/usr/bin/perl
use strict;
use warnings;
use List::Util 'shuffle';
my $rep = 1000; # permutation per mirna
$/ = "\n>";
while (<>) {
s/>//g;
my ($id, $seq) = split (/\n/, $_);
my @seq = split (//, $seq);
my %seen = ();
for (my $n=1; $n<=$rep; $n++) {
my @rand_seq = shuffle(@seq);
my $new_seq = join "", @rand_seq;
next if (defined $seen{$new_seq}); # skip sequences already generated
next if ($new_seq eq $seq); # skip if both are the same miRNA
print ">$id\_$n\n$new_seq\n";
$seen{$new_seq} = 1;
}
}
Save as any perl script and run as: perl randMiRNA.pl < mirna.fasta > random_mirnas.fasta
THANKS a lot for your script. brentp's suggestion is good, but I don't think to test a random sequence from intergenic and intron regions, because I am testing the second structure of already predicted miRNAs. From previous publications, they didn't test random genomic sequences. Anyway, your discussions are very good. JC, I have one more question: in your code you use
List::Util 'shuffle';
Would you please explain briefly about shuffle? Do i need to install biosquid? THANK YOU VERY MUCH!
List::Util is a core module in Perl, you don't need to install it. The "shuffle" function returns a list (array) in random order, I'm guessing it uses Fisher-Yates permutation but you can check the main algorithm in the source code: http://search.cpan.org/~pevans/Scalar-List-Utils-1.35/lib/List/Util.pm#shuffle_LIST
Log in to answer this question.
Instead of looking for a script to generate random sequences, may be you can use "off the shelf" tools to generate random sequences from your input sequence. Take a look at biosquid package and especially shuffle. (these are Debian/Ubuntu packages and I am not sure if you can find an alternative for other distro or windows and I've not tried it myself)
Yes; I'd use EMBOSS shuffleseq, which can easily be run from a (Bio)Perl (or other) script if required.