convert Fasta files into TIGR .contig file

During the process of converting Fasta files with gaps to AFG format, one step required generation of TIGR .contig format. Given below is a simple script to generate a dummy .contig file which has just one read for each contig which is exactly the same as the contig. This is done due to the lack of read tracking in SOAP, SSPACE etc. and also to avoid the huge file size if all the reads are retained.

 #!/usr/bin/perl  
 open FASTA, "< $ARGV[0]" or die "Can't open $ARGV[0] ($!)\n";  
 while($z=<FASTA>){  
  if($z=~ />/)  
  {  
  chomp $z;  
  $header=(split(/>/,$z))[1];  
  $readid=(split(/_/,$header))[1];  
  $seq=;  
  chomp $seq;  
  $seqlen=length $seq;  
  print "##$header 1 $seqlen bases, 00000000 checksum.\n";  
  print "$seq\n";  
  print "#$readid(0) [] $seqlen bases, 00000000 checksum. {$seqlen 0} <1>\n";  
  print "$seq\n";  
  }  
 }  

save it as fasta2contig.pl and use as "perl fasta2contig.pl contig.fa"

Minimus2 larger datasets with Nucmer

While Blat is much more faster than Nucmer, its results can be very different from the results from Nucmer. A version of Minimus2 which splits the data before feeding it to Nucmer and merging it back before processing further is able to give a result which is almost always identical to the results from Nucmer. The changes required to be done are given below:

NUCMER = nucmer
DELTAFILTER = delta-filter
SHOWCOORDS = show-coords
RUNNUCMER = perl runnucmer.pl
#------------------------------------------------------------------------------#

## Building AMOS bank & Dumping reads
10: rm -fr $(BANK)
11: $(BINDIR)/bank-transact -c -z -b $(BANK) -m $(TGT)
12: $(BINDIR)/dumpreads $(BANK) -M $(REFCOUNT) > $(REFSEQ)
13: $(BINDIR)/dumpreads $(BANK) -m $(REFCOUNT) > $(QRYSEQ)

## Getting overlaps
20: $(RUNNUCMER) $(OVERLAP) $(REFSEQ) $(QRYSEQ) $(PREFIX)
21: $(SHOWCOORDS) -H -c -l -o -r -I $(MINID) $(ALIGN) | $(BINDIR)/nucmerAnnotate | egrep 'BEGIN|END|CONTAIN|IDENTITY' > $(COORDS)
22: $(BINDIR)/nucmer2ovl -ignore $(MAXTRIM) -tab $(COORDS) | $(BINDIR)/sort2 > $(OVLTAB)


The code for runnucmer.pl is given below:

#!/usr/bin/perl

$allowedsize=100000000;

$overlap=$ARGV[0];
$mainref=$ARGV[1];
$mainqry=$ARGV[2];
$mainpre=$ARGV[3];
$reflen=0;
$refcount=1;

$qrylen=0;
$qrycount=1;

#open first ref file
open FILE1, ">>",$mainpre."ref".$refcount or die $!;

#reading the refrence sequences
open FILE2, "<",$mainref or die $!;


$seqst_temp="";
$header = ;#reading the first header

while($z = ){
if($z=~ />/)
{#next sequence
$reflen=$reflen+length $seqst_temp;
if($reflen < $allowedsize)
{
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
else
{
$refcount++;
$reflen=0;
close FILE1;
open FILE1, ">>",$mainpre."ref".$refcount or die $!;
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
$header=$z;
$seqst_temp="";
}
else
{
$z =~ s/[\n\t\f\r_0-9\s]//g;
$seqst_temp .= $z;
}
}#end of File while loop

if($reflen < $allowedsize)
{
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
else
{
$refcount++;
$reflen=0;
close FILE1;
open FILE1, ">>",$mainpre."ref".$refcount or die $!;
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
$header=$z;
$seqst_temp="";



close FILE1;
close FILE2;

#open first qry file
open FILE1, ">>",$mainpre."qry".$qrycount or die $!;

#reading the refrence sequences
open FILE2, "<",$mainqry or die $!;


$seqst_temp="";
$header = ;#reading the first header

while($z = ){
if($z=~ />/)
{#next sequence
$qrylen=$qrylen+length $seqst_temp;
if($qrylen < $allowedsize)
{
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
else
{
$qrycount++;
close FILE1;
$qrylen=0;
open FILE1, ">>",$mainpre."qry".$qrycount or die $!;
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
$header=$z;
$seqst_temp="";
}
else
{
$z =~ s/[\n\t\f\r_0-9\s]//g;
$seqst_temp .= $z;
}
}#end of File while loop

if($qrylen < $allowedsize)
{
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
else
{
$qrycount++;
close FILE1;
$qrylen=0;
open FILE1, ">>",$mainpre."qry".$qrycount or die $!;
print FILE1 $header;
print FILE1 $seqst_temp;
print FILE1 "\n";
}
$header=$z;
$seqst_temp="";

close FILE1;
close FILE2;

for($i=1;$i<=$qrycount;$i++){
for($j=1;$j<=$refcount;$j++){
$refseq=$mainpre."ref".$j;
$qryseq=$mainpre."qry".$i;
$prefix=$mainpre."nuc".$j."".$i;
system("/bubo/home/h7/nagarjun/glob/tools/MUMmer3.22/nucmer -maxmatch -c $overlap $refseq $qryseq -p $prefix");
}
}


#merging delta files
open FILE1, ">>",$mainpre.".delta" or die $!;
for($i=1;$i<=$qrycount;$i++){
for($j=1;$j<=$refcount;$j++){
open FILE2, "<",$mainpre."nuc".$j."".$i.".delta" or die $!;
if(($i==1)&&($j==1)){while($a=){print FILE1 $a;}}
else{$a=;$a=;while($a=){print FILE1 $a;}}
}
}

Blat instead of Nucmer in Minimus2

Minimus2 is a good utility to merge assemblies based on the AMOS architecture. However, it has a limitation on the size of the input datasets. This is due to the Nucmer program is its pipline. Apart from Nucmer other steps are able to handle much bigger chunks of data. So instead of using Nucmer i tried using Blat. Although this gave me assemblies, they were not always correct and complete.

May be someday somebody will use it.In the meantime i try to use Nucmer iteratively on smaller chunks and then merge the output before the next step. This is possible as Nucmer is run with the option -maxmatch in Minimus, making it independent of the other sequences in the dataset.

1)The changes to Minimu2.acf in amos/src/Pipeline/:

BINDIR = /usr/local/bin
NUCMER = /usr/local/nucmer
DELTAFILTER = /usr/local/delta-filter
SHOWCOORDS = /usr/local/show-coords
BLAT = blat
BTON = perl psltoNucmer.pl

#-----------------------------

-------------------------------------------------#

## Building AMOS bank & Dumping reads
10: rm -fr $(BANK)
11: $(BINDIR)/bank-transact -c -z -b $(BANK) -m $(TGT)
12: $(BINDIR)/dumpreads $(BANK) -M $(REFCOUNT) > $(REFSEQ)
13: $(BINDIR)/dumpreads $(BANK) -m $(REFCOUNT) > $(QRYSEQ)

## Getting overlaps
20: $(BLAT) $(REFSEQ) $(QRYSEQ) $(BOUT) -fastMap -noHead -minIdentity=98
21: $(BTON) $(BOUT) | $(BINDIR)/nucmerAnnotate | egrep 'BEGIN|END|CONTAIN|IDENTITY' > $(COORDS)
22: $(BINDIR)/nucmer2ovl -ignore $(MAXTRIM) -tab $(COORDS) | $(BINDIR)/sort2 > $(OVLTAB)


2) The code for psltoNucmer.pl is given below:

#!/usr/bin/perl

open FILE2, "<",$ARGV[0] or die $!;

while($z = ){
@values=split(/\t/,$z);
print $values[15]."\t".$values[16]."\t|\t".$values[11]."\t".$values[12]."\t|\t".abs($values[16]-$values[15])."\t".abs($values[12]-$values[11])."\t|\t".(sprintf "%.2f",(($values[0]+$values[2]+$values[3])*100)/abs($values[12]-$values[11]))."\t|\t".$values[14]."\t".$values[10]."\t|\t".(sprintf "%.2f",(abs($values[16]-$values[15])*100/$values[14]))."\t".(sprintf "%.2f",(abs($values[12]-$values[11])*100/$values[10]))."\t|\t".$values[13]."\t".$values[9]."\n";
}#end of File while loop

3) blat output format details can be found here : http://genome.ucsc.edu/goldenPath/help/blatSpec.html

Program to shred contigs into reads with overlap

Most assembly programs seem to have some limit on the length of the reads that they take as input for assembling them. Assembling contigs can be a difficult problem in such cases. Shredding the contigs into smaller reads with a 50% or more of overlap should retain the benefit of having assembled the contigs and use it for further assembly.

Here is a small program to shred contigs into 1kb reads with 500 base pair overlap:

#!/usr/bin/perl

open FILE1, $ARGV[0] or die $!;

while($line1 = ){
my $flag=0,$overlap=500,$length=1000;

chomp($line1);

$line2 = ;
$seqlen=length $line2;
chomp($line2);

while(($seqlen-$flag)>$length){
if(($seqlen-($flag+$overlap))<$length){
$length=$seqlen-$flag;
}
$nextseq=substr $line2,$flag,$length;
print $line1.":".$flag."-".($flag+$length)."\n";
print $nextseq."\n";
$flag=$flag+$overlap;

}#end of seqlen while loop


}#end of file while loop

This script works on a multifasta file with header in the first line and the sequence in the second line.

Handling duplicate accession numbers created by flowsim with sfffile command

Having simulated loads of 454 data with Flowsim, i wanted to assemble everything. Unfortunately Flowsim gives one output .sff file for one input fasta file. In my case there were some 1,00,000 fasta sequences and hence 1,00,000 .sff files.

Newbler (assembler for 454 data) just crashed when asked to assemble the data in these many different files. It worked upto 500 .sff files though.

When all the files where given as input to sfffile to combine into a single file it gave a cryptic "segmentation fault" error. It seemed as though all the loads of data generated by Flowsim would be wasted. As a lost resort i wrote a script to combine the .sff files two at a time.

The script crashed after running for 100 files. The culprit was "Flowsim". It had generated sequences with the same accession number in the different runs. sfffile was failing with "Error: Duplicate accession chr1:56036789..56039961_1148- found. Unable to write SFF file.". After further investigation it became clear that duplciate accession numbers had occured rather frequently. Later versions of Flowsim may have fixed this issue (I used version 0.1).

Since, the data had already been generated it had to "handled". sfffile does not have a flag to ignore duplicate accession numbers. Hence, a simple perl script that checks if sfffile failed in the previous step take care of the duplicate accession numbers. Perl script below:

 #!/usr/bin/perl  
 my $startfile="combined1.sff";  
 for($count=2;$count<8000;$count++){  
 $pcount=$count-1;  
 print "sfffile -o combined".$count.".sff combined".$pcount.".sff myreads.fasta.".$count.".sff \n";  
 $pfilename="combined".$pcount.".sff";  
 if(-e $pfilename){  
 system("sfffile -o combined".$count.".sff combined".$pcount.".sff myreads.fasta.".$count.".sff");  
 $filename="combined".$count.".sff";  
 if(-e $filename){$startfile=$filename;}  
 }  
 else {system("sfffile -o combined".$count.".sff ".$startfile." myreads.fasta.".$count.".sff");}  
 }  

Youtube bug or hidden feature?

Certain Youtube videos which are considered "inappropriate for some users, as flagged by YouTube's user community" are supposed to be hidden from users who are not logged in or who are not above 18 years old. But alas, google seems to have left a bug(or is it a hidden feature--not decipherbale by the underaged?) which allows you to watch any video you want without having to sign in or sign up.

Just give the direct path to the video like "http://www.youtube.com/v/video_id" and the video just opens up in the browser. This seems to have been left open to allow external websites to embed youtube videos. However, one would expect YT with the "brilliance" of google to have blocked embed access to the "marked" videos.

Apart from the flaw in youtube that has been hastily fixed by google, there seems to be minor bugs or rather design decisions that are confusing. Will they just keep patching things or consider a redesign of the model is interesting to watch for.

Stars and stripes in perl

Perl can be made to print in color using TERM::ANSIColor. This perl program tries to generate the stars and stripes with 50 stars and 13 stripes. Although its not perfect in having larger gaps between every second row of stars, it seems to do a pretty ok job. Even the 2/5 ratio between the full and partial stripes is valid. This is the longest surviving version which has managed to complete 50 long years.

Perl does not seem to have the ability to print subscripts or superscripts while printing to the terminal atleast. The uneven spacing can probably be avoided by using subscripts or superscripts.

 #!/usr/bin/perl  
 use Term::ANSIColor qw(:constants);  
 $sixstars=" * * * * * * ";  
 $fivestars=" * * * * * ";  
 $emptystars=" ";  
 $partstripes=" \n";  
 $fullstripes=" \n";  
 print RESET;  
 print BOLD,WHITE,ON_BLUE,$sixstars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$fivestars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$emptystars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$sixstars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$fivestars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$emptystars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$sixstars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$fivestars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$emptystars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$sixstars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$fivestars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$emptystars;  
 print ON_WHITE,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$sixstars;  
 print ON_RED,$partstripes,RESET;  
 print BOLD,WHITE,ON_BLUE,$emptystars;  
 print ON_RED,$partstripes,RESET;  
 for($i=0;$i<3;$i++){  
 print ON_WHITE,$fullstripes,RESET;  
 print ON_WHITE,$fullstripes,RESET;  
 print ON_RED,$fullstripes,RESET;  
 print ON_RED,$fullstripes,RESET;  
 }