Perl & Genomics

#!/usr/bin/perl
 
print STDERR "
  [Write instructions here]
";
 
($inputsequence, $coverage, $readlength) = @ARGV;
 
if ($readlength==0) {
	die " FATAL ERROR:\n You have to specify three parameters.\n";
}
 
open(I, "$inputsequence") || die " FATAL ERROR:\n Unable to load sequence from $inputsequence.\n";
 
while (<I>) {
	chomp($_);
 
	if ($_=~/>/) {
		$countheaders++;
	} else {
		$sequence.=$_;
	}
}
 
if ($countheaders>1) {
	print STDERR "Warning: you gave me a multifasta with $countheaders sequences.\n I merged them :)\n";
}
 
$genomesize = length($sequence);
 
$readsnumber = int($genomesize*$coverage/$readlength);
 
print STDERR "
Genome size    $genomesize
Reads number   $readsnumber
";
for ($i=1; $i<=$readsnumber; $i++) {
	$position = int(rand($genomesize-$readlength));
	$readsequence = substr($sequence, $position, $readlength);
	print ">Sequence$i-pos:$position\n$readsequence\n";
}