← Back to branch summary

~ubuntu-branches/ubuntu/vivid/cctools/vivid

« back to all changes in this revision

Viewing changes to sand/test/filter_verification/gen_random_reads.pl

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2011-05-07 09:05:00 UTC
  • Revision ID: james.westby@ubuntu.com-20110507090500-lqpmdtwndor6e7os
Tags: upstream-3.3.2
Import upstream version 3.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
sand/test/filter_verification/gen_random_reads.pl
 
1
#!/usr/bin/perl
 
2
# Assume we already have two 10,000-base DNA sequences who are reverse
 
3
# complement to each other, generate random 100-base reads from these two
 
4
# sequences. The two sequences are expected to be in file "rando.seq" and
 
5
# "random_revcom.seq", as generated by the "gen_randome_sequence.pl" script.
 
6
# The created reads will be saved in file "random.fa" in the FASTA format.
 
7
# "random_revcom.seq" respectively.
 
8
use strict;
 
9
use warnings;
 
10
 
 
11
my $sequence_length = 100000;
 
12
my $read_length = 100;
 
13
my $seq_file = "random.seq";
 
14
my $seq_revcom_file = "random_revcom.seq";
 
15
my $fasta_file = "random.fa";
 
16
 
 
17
my $random_number;
 
18
my $char;
 
19
my $i;
 
20
my $seq;
 
21
my $seq_revcom;
 
22
my $kmer;
 
23
my $num_of_reads;
 
24
 
 
25
# Check sequence_length argument
 
26
my $numArgs;
 
27
$numArgs = $#ARGV+1;
 
28
if ($numArgs == 1) {
 
29
        $sequence_length = $ARGV[0];
 
30
}
 
31
 
 
32
open SEQFILE, "$seq_file" or die $!;
 
33
$seq = <SEQFILE>;
 
34
print length($seq) . " bases are read from $seq_file.\n";
 
35
close(SEQFILE);
 
36
 
 
37
open SEQFILE, "$seq_revcom_file" or die $!;
 
38
$seq_revcom = <SEQFILE>;
 
39
print length($seq_revcom) . " bases are read from $seq_revcom_file.\n";
 
40
close(SEQFILE);
 
41
 
 
42
 
 
43
open FASTAFILE, ">$fasta_file" or die $!;
 
44
$num_of_reads = ($sequence_length / $read_length) * 10;
 
45
for($i = 0; $i < $num_of_reads; $i++) {
 
46
        $random_number = int(rand($sequence_length - $read_length + 1));
 
47
        $kmer = substr($seq, $random_number, $read_length);
 
48
        print FASTAFILE ">$random_number\n$kmer\n";
 
49
}
 
50
 
 
51
for($i = 0; $i < $num_of_reads; $i++) {
 
52
        $random_number = int(rand($sequence_length - $read_length + 1));
 
53
        $kmer = substr($seq_revcom, $random_number, $read_length);
 
54
        print FASTAFILE ">rc$random_number\n$kmer\n";
 
55
}
 
56
close(FASTAFILE);
 
57
 
 
58
$num_of_reads = $num_of_reads * 2;
 
59
print "$num_of_reads random reads are generated successfully!\n";
 
60
exit(0);