In silico data sets from BACs for GAII Illumina

[CG&R]

Hi there! Currently I am planning a next-gen seq project with the GAII
of Illumina. But before I dive in, I have one question.

Basically I want to do targetted sequencing of a 300kb region.
The region is notorious for it's copynumber variation and could easily vary in size by 100kb.

Up till now one conventional BAC clone of about 190kb has been paired-end Sanger sequenced. Basically I want to know if Alignment Software can
handle this region filled with repeats.

My question is the following: Is there a software program: where I can give the 190kb sequence as input. And the software then does an in silico cutting/shearing similar as would happen during the actual experiment.
(It would be helpfull if variable fragment sizes can be analysed.)

Then from these fragments, an in silico mate-paired data set is generated, and this is fed back into an allignment program. Basically a new contig is generated, and this is then compared with the actual input sequence.

Any help/suggestions is much appreciated! Thanks

[krobison]

There are some programs out there; see previous thread

http://seqanswers.com/forums/showthread.php?t=806

This my personal one; it needs some more work (doesn't generate quality information & the errors are evenly distributed in the read) but is useful as a baseline.

#!/usr/bin/perl
use strict;
use Bio::SeqIO;
use Statistics:istrib::Normal;

# quick & dirty Illumina read simulator
# Keith Robison. Infinity Pharmaceuticals Inc

# heavily modified from http://wiki.bioinformatics.ucdavis.e...asta_reference

## several flaws in original modified
## 1) proper FASTA reading via Bioperl
## 2) handles multiple input sequences correctly
## 3) reads can come from either strand
## 4) depth of reads set with constant instead of fixed number of reads

my $insertMean=200;
my $insertSd=20;
my $dist = Statistics:istrib::Normal->new(mu => $insertMean,
sigma => $insertSd);

my $coverage=40;
my $readLen=50;
my $bases="ATCG";
my $stem="reads";
if ($ARGV[0] eq '-s') # -s stem where stem will be the beginning of output file names
{
shift(@ARGV); $stem=shift(@ARGV);
}
open(FWD,">$stem.1.fasta");
open(REV,">$stem.2.fasta");
my $sourceSeq=0;
foreach my $arg(@ARGV)
{
my $seqFile=new Bio::SeqIO(-file=>$arg,-format=>"Fasta");
while (my $seq=$seqFile->next_seq())
{
$sourceSeq++;
my $seqLen=$seq->length;
my $nReadsToGenerate=int($seqLen/$readLen * $coverage);
my @fragmentSizes=$dist->rand($nReadsToGenerate);

for (my $i=1; $i<scalar(@fragmentSizes); $i++)
{
my $fragmentSize=$fragmentSizes[$i];
my $pos=int(rand($seqLen-$fragmentSizes[$i]-1));
my $fragment=substr($seq->seq,$pos,$fragmentSize);
if (rand()>=0.5) # reverse complement
{
$fragment=reverse($fragment);
$fragment=~tr/ATCG/TAGC/;
}
for (my $j=0; $j<$readLen; $j++)
{
if (rand()<0.001)
{
substr($fragment,$j,1)=substr($bases,int(rand(4)),1);
}
if (rand()<0.001)
{
substr($fragment,$fragmentSize-$j,1)=substr($bases,int(rand(4)),1);
}
}

print FWD ">$sourceSeq-$i\n",substr($fragment,0,$readLen),"\n";
my $revRead=reverse(substr($fragment,$fragmentSize-$readLen,$readLen));
$revRead=~tr/ATCG/TAGC/;
print REV ">$sourceSeq-$i\n",$revRead,"\n";
}
}
}