genomics

Tool and Paper Watch

2019-03-14T00:13:00.003-07:00

Welcome back! Here I have decided to update the tools available for each of the sections:
Genomics; Transcriptomics; Databases; Statistics; Data analysis:

Genomics:
Assembly:

Annotation: Updating and improving genome annotation: https://academic.oup.com/nar/article/47/2/559/5213907
Comparative Genomics: https://nph.onlinelibrary.wiley.com/doi/full/10.1111/nph.15687
Genome compartmentalization: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4824034/

Transcriptomics:
Analysis:
1. DEE2: Digital Expression Explorer 2: Is a database containing analyzed RNAseq data from multiple platforms: http://dee2.io/about.html
2. RNAseqAnalysis pipeline and resources: https://github.com/crazyhottommy/RNA-seq-analysis#databases
3. Best Practises for detecting germline variations in human genomes; https://www.nature.com/articles/s41587-019-0054-x
4. Deep learning with empirical RNAseq evidence (DART) : https://github.com/Xinglab/DARTS
https://www.nature.com/articles/s41592-019-0351-9

Metagenome assemblies
https://academic.oup.com/bib/advance-article/doi/10.1093/bib/bbz025/5371422

Databases:

Statistics & R:
Cola: A General Framework for Consensus and Hierarchical partitioning: https://github.com/jokergoo/cola

Fastq format converter in Fastx

2013-07-06T06:02:00.001-07:00

Fastx is a useful tool for quality clipping, getting quality statistics etc. Just running with -i and -o parameters, it will complain if you have file in sanger fastq format. there is an undocumented parameter that takes the quality offset. For instance if you want to run a sanger fastq file, then do the following:
./fastx_quality_stats -i tmp -Q 33 -o tmp.out
But if the file is in Solexa format, then you dont have to specify the -Q option.

However, there is always a need to have a script in your tool box that can convert one fastq format to another. I found a nice script online probably written by Heng Li

It takes about 30 mins to run in a single processor mode running on the head node of a cluster with 49 GB memory took approximately 1 hour for a 25 GB input file.

#!/usr/bin/perl -w

# Author: lh3
#
# fq_all2std.pl was originally distributed as part of the Maq software package (http://maq.sourceforge.net/)
# and is distributed under the GNU Public License (GPL) http://www.gnu.org/copyleft/gpl.html.
#
# Modified: Lance Parsons
# Added illumina2std function
# Fixed CR/LF conversion in sol2std
# Version: 0.1.7_lp

use strict;
use warnings;
use Getopt::Std;

my $usage = qq(
Usage:   fq_all2std.pl <command> <in.txt>

Command: scarf2std      Convert SCARF format to the standard/Sanger FASTQ
         fqint2std      Convert FASTQ-int format to the standard/Sanger FASTQ
         sol2std        Convert Solexa/Illumina FASTQ (Pipline 1.0 and below) to the standard FASTQ
         illumina2std Convert Illumina FASTQ (Pipeline v1.3) to the standard FASTQ
         std2sol        Convert standard FASTQ to Solexa/Illumina FASTQ (simplified)
         std2illumina Convert standard FASTQ to Illumina FASTQ (Pipeline v1.3)
         fa2std         Convert FASTA to the standard FASTQ
         seqprb2std     Convert .seq and .prb files to the standard FASTQ
         fq2fa          Convert various FASTQ-like format to FASTA
         export2sol     Convert Solexa export format to Solexa FASTQ
         export2std     Convert Solexa export format to Sanger FASTQ
         csfa2std       Convert AB SOLiD read format to Sanger FASTQ
         std2qual       Convert standard FASTQ to .seq+.qual
         instruction    Explanation to different format
         example        Show examples of various formats

Note:    Read/quality sequences MUST be presented in one line.
\n);

die($usage) if ( @ARGV < 1 );

# Solexa->Sanger quality conversion table
my @conv_table;
my @rev_conv_table;
for ( -64 .. 64 ) {
 $conv_table[ $_ + 64 ] = chr( int( 33 + 10 * log( 1 + 10**( $_ / 10.0 ) ) / log(10) + .499 ) );
 $rev_conv_table[ int( 33 + 10 * log( 1 + 10**( $_ / 10.0 ) ) / log(10) + .499 ) ] = chr( $_ + 64 );
}

# parsing command line
my $cmd      = shift;
my %cmd_hash = (
 scarf2std    => \&scarf2std,
 fqint2std    => \&fqint2std,
 sol2std      => \&sol2std,
 fa2std       => \&fa2std,
 fq2fa        => \&fq2fa,
 example      => \&example,
 instruction  => \&instruction,
 export2sol   => \&export2sol,
 export2std   => \&export2std,
 csfa2std     => \&csfa2std,
 seqprb2std   => \&seqprb2std,
 std2sol      => \&std2sol,
 illumina2std => \&illumina2std,
 std2illumina => \&std2illumina,
 std2qual     => \&std2qual
);
if ( defined( $cmd_hash{$cmd} ) ) {
 &{ $cmd_hash{$cmd} };
}
else {
 die("** Unrecognized command $cmd");
}

sub fa2std {
 my %opts = ( q => 25 );
 getopts( 'q:', \%opts );
 die("Usage: fq_all2std.pl fa2std [-q $opts{q}] <in.fa>\n") if ( -t STDIN && @ARGV == 0 );
 my $q = chr( $opts{q} + 33 );
 while (<>) {
  if (/^>(\S+)/) {
   print "\@$1\n";
   $_ = <>;
   print "$_+\n", $q x ( length($_) - 1 ), "\n";
  }
 }
}

sub csfa2std {
 my %opts = ( q => 25, Q => '', l => 0 );
 getopts( 'q:Q:l:', \%opts );
 die( "
Usage:   fq_all2std.pl csfa2std [options] <in.csfa>\n
Options: -q INT     default base quality [$opts{q}]
         -Q FILE    quality file [null]
         -l INT     output read length, 0 for auto [$opts{l}]

Note:    For paired-end alignment, Maq requires two sequence files as the
         input. The n-th read in the first file should forms a read pair with
         the n-th read in the second file. However, SOLiD reads may be
         singletons and therefore further prepocessing is needed.
\n" ) if ( -t STDIN && @ARGV == 0 );
 my ( $fh, $name, $seq );
 my $len = $opts{l};
 my $q   = chr( $opts{q} + 33 );
 if ( $opts{Q} ) {
  open( $fh, $opts{Q} ) || die("** fail to open quality file '$opts{Q}'");
 }
 while (1) {
  while (<>) { last if (/^>/); }
  last unless ($_);
  /^>(\S+)/;
  $name = $1;
  $_ = substr( <>, 2 );
  chomp;
  tr/0123./ACGTN/;
  $seq = $_;

  if ($fh) {    # .qual file is available
   while (<$fh>) { last if (/^>(\S+)/); }
   /^>(\S+)/;
   die("** unmatched seq-qual name: '$name' ne '$1'") unless ( $1 eq $name );
   $_ = <$fh>;
   s/(\s*\d+\s*)/chr(int($1) + 33)/eg;
   $_ = substr( $_, 1 );
  }
  else {
   $_ = $q x length($seq);
  }
  if ( $name =~ /^(\S+)_F\d$/ ) {    # change read name for maq
   $name = "$1/1";
  }
  elsif ( $name =~ /^(\S+)_R\d$/ ) {
   $name = "$1/2";
  }
  if ($len) {                        # chop the sequence if required
   $seq = substr( $seq, 1, $len );
   $_   = substr( $_,   1, $len );
  }
  print "\@$name\n$seq\n+\n$_\n";
 }
 close($fh) if ($fh);
}

sub fq2fa {
 while (<>) {
  if (/^@(\S+)/) {
   print ">$1\n";
   $_ = <>;
   print;
   <>;
   <>;
  }
 }
}

sub scarf2std {
 while (<>) {
  my @t = split( ':', $_ );
  my $name = join( '_', @t[ 0 .. 4 ] );
  print "\@$name\n$t[5]\n+\n";
  my $qual = '';
  @t = split( /\s/, $t[6] );
  $qual .= $conv_table[ $_ + 64 ] for (@t);
  print "$qual\n";
 }
}

sub seqprb2std {
 die("Usage: fq_all2std.pl seqprb2std <in.seq.txt> <in.prb.txt>\n") if ( @ARGV != 2 );
 my ( $fhs, $fhq );
 open( $fhs, $ARGV[0] ) || die;
 open( $fhq, $ARGV[1] ) || die;
 while (<$fhs>) {
  my @t = split;
  my $name = join( ":", @t[ 0 .. 3 ] );
  $t[4] =~ tr/./N/;
  print "\@$name\n$t[4]\n+\n";
  $_ = <$fhq>;
  @t = split;
  my $q   = '';
  my $max = -100;

  for ( 0 .. $#t ) {
   $max = $t[$_] if ( $t[$_] > $max );
   if ( ( $_ & 0x3 ) == 3 ) {
    $q .= $conv_table[ $max + 64 ];
    $max = -100;
   }
  }
  print "$q\n";
 }
 close($fhs);
 close($fhq);
}

sub export2sol {
 while (<>) {
  chomp;
  my @t = split( "\t", $_ );
  my $output = *STDOUT;
  if ( $t[21] ne 'Y' ) {
   $output = *STDERR;
  }

  my $x = ( defined( $t[7] ) && ( $t[7] == 1 || $t[7] == 2 ) ) ? "/$t[7]" : '';
  $t[0] =~ s/^(SLXA|HWI)-//;
  $t[0] =~ s/_Human//i;
  $t[0] =~ s/_PhiX//i;
  $t[0] =~ s/_R1//;
  my $rn_head = ( $t[0] =~ /(^[A-Z]+\d+_\d+)/ ) ? $1 : ( $t[1] ? "$t[0]_$t[1]" : $t[0] );
  print {$output} "\@$rn_head:$t[2]:$t[3]:$t[4]:$t[5]$x\n$t[8]\n+\n$t[9]\n";
 }
}

sub export2std {
 while (<>) {
  chomp;
  my @t = split( "\t", $_ );
  my $output = *STDOUT;
  if ( $t[21] ne 'Y' ) {
   $output = *STDERR;
  }
  my $x = ( defined( $t[7] ) && ( $t[7] eq 1 || $t[7] eq 2 ) ) ? "/$t[7]" : '';
  $t[0] =~ s/^SLXA-//;
  $t[0] =~ s/_Human//i;
  $t[0] =~ s/_PhiX//i;
  $t[0] =~ s/_R1//;
  my $rn_head = ( $t[0] =~ /(^[A-Z]+\d+_\d+)/ ) ? $1 : ( $t[1] ? "$t[0]_$t[1]" : $t[0] );
  print {$output} "\@$rn_head:$t[2]:$t[3]:$t[4]:$t[5]$x\n$t[8]\n";
  my @s = split( '', $t[9] );
  my $qual = '';
  $qual .= $conv_table[ ord($_) ] for (@s);
  print {$output} "+\n$qual\n";
 }
}

sub fqint2std {
 while (<>) {
  if (/^@/) {
   print;
   $_ = <>;
   print;
   $_ = <>;
   $_ = <>;
   my @t    = split;
   my $qual = '';
   $qual .= $conv_table[ $_ + 64 ] for (@t);
   print "+\n$qual\n";
  }
 }
}

sub sol2std {
 my $max = 0;
 while (<>) {
  if (/^@/) {
   print;
   $_ = <>;
   print;
   $_ = <>;
   $_ = <>;

   # Added to eliminate carriage return conversion
   chomp;
   my @t = split( '', $_ );
   my $qual = '';
   $qual .= $conv_table[ ord($_) ] for (@t);
   print "+\n$qual\n";
  }
 }
}

sub illumina2std {
 my $max = 0;
 while (<>) {
  if (/^@/) {
   print;
   $_ = <>;
   print;
   $_ = <>;
   $_ = <>;

   # Added to eliminate carriage return conversion
   chomp;
   my @t = split( '', $_ );
   my $qual = '';
   $qual .= chr( ord($_) - 31 ) for (@t);
   print "+\n$qual\n";
  }
 }
}

sub std2illumina {
 my $max = 0;
 while (<>) {
  if (/^@/) {
   print;
   $_ = <>;
   print;
   $_ = <>;
   $_ = <>;
   chomp;
   my @t = split( '', $_ );
   my $qual = '';
   $qual .= chr( ord($_) + 31 ) for (@t);
   print "+\n$qual\n";
  }
 }
}

sub std2sol {
 my $max = 0;
 while (<>) {
  if (/^@/) {
   print;
   $_ = <>;
   print;
   $_ = <>;
   $_ = <>;
   chomp;
   
   #tr/!-]/@-|/; 
   #print "+\n$_\n";
   
   my @t = split( '', $_ );
   my $qual = '';
   $qual .= $rev_conv_table[ ord($_) ] for (@t);
   print "+\n$qual\n";
  }
 }
}

sub std2qual {
 die("Usage fq_all2std.pl std2qual <out.prefix> <in.fastq>\n") if ( @ARGV == 0 );
 my $pre = shift(@ARGV);
 my ( $fhs, $fhq );
 open( $fhs, ">$pre.seq" )  || die;
 open( $fhq, ">$pre.qual" ) || die;
 while (<>) {
  s/^@/>/;
  print $fhs $_;
  print $fhq $_;
  $_ = <>;
  print $fhs $_;
  <>;
  $_ = <>;
  s/([!-~])/" ".(ord($1)-33)/eg;
  $_ = substr( $_, 1 );
  print $fhq $_;
 }
 close($fhs);
 close($fhq);
}

sub instruction {

 print "
FASTQ format is first used in the Sanger Institute, and therefore
we take the Sanger specification as the standard FASTQ. Although
Solexa/Illumina reads file looks pretty much like the standard
FASTQ, they are different in that the qualities are scaled
differently. In the quality string, if you can see a character
with its ASCII code higher than 90, probably your file is in the
Solexa/Illumina format.

Sometimes we also use an integer, instead of a single character,
to explicitly show the qualities. In that case, negative
qualities indicates that Solexa/Illumina qualities are used.

";

}

sub example {
 my $exam_scarf = '
USI-EAS50_1:4:2:710:120:GTCAAAGTAATAATAGGAGATTTGAGCTATTT:23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 19 23 23 23 18 23 23 23
USI-EAS50_1:4:2:690:87:GTTTTTTTTTTTCTTTCCATTAATTTCCCTTT:23 23 23 23 23 23 23 23 23 23 23 23 12 23 23 23 23 23 16 23 23 9 18 23 23 23 12 23 18 23 23 23
USI-EAS50_1:4:2:709:32:GAGAAGTCAAACCTGTGTTAGAAATTTTATAC:23 23 23 23 23 23 23 23 20 23 23 23 23 23 23 23 23 23 23 23 23 12 23 18 23 23 23 23 23 23 23 23
USI-EAS50_1:4:2:886:890:GCTTATTTAAAAATTTACTTGGGGTTGTCTTT:23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
USI-EAS50_1:4:2:682:91:GGGTTTCTAGACTAAAGGGATTTAACAAGTTT:23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 20 23 23 23 23 23 23 23 23 23 23 23 18 23 23 23 23
USI-EAS50_1:4:2:663:928:GAATTTGTTTGAAGAGTGTCATGGTCAGATCT:23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23 23
';

 my $exam_fqint = '
@4_1_912_360
AAGGGGCTAGAGAAACACGTAATGAAGGGAGGACTC
+4_1_912_360
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 21 40 40 40 40 40 40 40 40 40 26 40 40 14 39 40 40
@4_1_54_483
TAATAAATGTGCTTCCTTGATGCATGTGCTATGATT
+4_1_54_483
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 16 40 40 40 28 40 40 40 40 40 40 16 40 40 5 40 40
@4_1_537_334
ATTGATGATGCTGTGCACCTAGCAAGAAGTTGCATA
+4_1_537_334
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 21 29 40 40 33 40 40 33 40 40 33 31 40 40 40 40 18 26 40 -2
@4_1_920_361
AACGGCACAATCCAGGTTGATGCCTACGGCGGGTAC
+4_1_920_361
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 39 40 40 40 40 40 40 40 40 31 40 40 40 40 40 40 15 5 -1 3
@4_1_784_155
AATGCATGCTTCGAATGGCATTCTCTTCAATCACGA
+4_1_784_155
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 31 40 40 40 40 40
@4_1_595_150
AAAGACGTGGCCAGATGGGTGGCCAAGTGCCCGACT
+4_1_595_150
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 30 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 14 40 40
';

 my $exam_sol = '
@SLXA-B3_649_FC8437_R1_1_1_610_79
GATGTGCAATACCTTTGTAGAGGAA
+SLXA-B3_649_FC8437_R1_1_1_610_79
YYYYYYYYYYYYYYYYYYWYWYYSU
@SLXA-B3_649_FC8437_R1_1_1_397_389
GGTTTGAGAAAGAGAAATGAGATAA
+SLXA-B3_649_FC8437_R1_1_1_397_389
YYYYYYYYYWYYYYWWYYYWYWYWW
@SLXA-B3_649_FC8437_R1_1_1_850_123
GAGGGTGTTGATCATGATGATGGCG
+SLXA-B3_649_FC8437_R1_1_1_850_123
YYYYYYYYYYYYYWYYWYYSYYYSY
@SLXA-B3_649_FC8437_R1_1_1_362_549
GGAAACAAAGTTTTTCTCAACATAG
+SLXA-B3_649_FC8437_R1_1_1_362_549
YYYYYYYYYYYYYYYYYYWWWWYWY
@SLXA-B3_649_FC8437_R1_1_1_183_714
GTATTATTTAATGGCATACACTCAA
+SLXA-B3_649_FC8437_R1_1_1_183_714
YYYYYYYYYYWYYYYWYWWUWWWQQ
';

 print qq(
solexa
======
$exam_sol
scarf
=====
$exam_scarf
fqint
=====
$exam_fqint
);
}

Not mine - I found this csfasta to fastq conversion script that worked best...

2013-07-04T06:55:00.002-07:00

#!/usr/bin/perl -w
use strict;
use Getopt::Long;

my $usage = qq(
This script is used for converting SOLiD/ABi FASTQ to Standard(Sanger) FASTQ
Usage : perl $0 [Option]
Option: -seq seq file
-qual qual file
-fastq FASTQ file
-o output file
-help Help Information
Example: perl $0 -fastq solid.fastq -o std.fastq

Author : BENM <BinxiaoFeng\@gmail.com>
Version: 1.2
Date : 2009-10-06
Update : 2012-03-27
\n);
my ($Seq,$Qual,$Fastq,$Output,$Header,$Help);
my %opts;
GetOptions
(
\%opts,
"seq:s"=>\$Seq,
"qual:s"=>\$Qual,
"fastq:s"=>\$Fastq,
"o:s"=>\$Output,
"header:s"=>\$Header,
"help"=>\$Help
);

die($usage) if (((!defined $Seq)||(!defined $Qual))&&(!defined $Fastq)&&(!defined $Output)||($Help));

# Solexa->Sanger quality conversion table
my @conv_table;
for (-64..64) {
$conv_table[$_+64] = chr(int(33 + 10*log(1+10**($_/10.0))/log(10)+.499));
}

# SOLiD color code
my @code = ([0,1,2,3],[1,0,3,2],[2,3,0,1],[3,2,1,0]);
my @bases = qw(A C G T);
my %decode = ();
foreach my $i(0..3)
{
foreach my $j(0..3)
{
$decode{$code[$i]->[$j]} -> {$bases[$i]} = $bases[$j];
}
}

open (OUT,">$Output") || die "Fail to create FASTA file:$Output for writing\n";
if (defined $Fastq)
{
open (IN,$Fastq) || die "Fail to open FASTQ file:$Fastq for reading\n";
while (<IN>)
{
if (/^@/)
{
print OUT $_;
$_ = <IN>;
s/\s+$//;
my $seq = ($_=~/\d+/) ? col2base($_) : $_;
print OUT "$seq\n";
$_ = <IN>; print OUT $_; $_ = <IN>;
s/\s+$//;
my @t = split;
my $qual = '';
if ($_=~/^\d+\s*/)
{
my @t=split;
$qual .= $conv_table[$_+64] for (@t); #intfq: SOLiD -> Standard
}
else
{
$qual .= $_;
}
substr($qual,0,1,'') if (length($qual)>length($seq));
print OUT "$qual\n";
}
}
close IN;
}
elsif ((defined $Seq)&&(defined $Qual))
{
my ($Title,$Name,$seq,$qual)=("","","","");
open (IN1,$Seq) || die "Fail to open Seq file:$Seq for reading\n";
open (IN2,$Qual) || die "Fail to open Quality file:$Qual for reading\n";
while(<IN2>)
{
if ($_=~/^\#\s+Title\:\s+(\S+)/)
{
$Title=$1;
my $in1=<IN1>;
while($in1 !~ /^\#\s+Title\:\s+(\S+)/)
{
$in1=<IN1>;
}
}
else
{
if ($_=~/^#\s\w+/)
{
<IN2>;<IN1>;
}
elsif (/^[\@\>](\S+)/)
{
chomp $seq;
chomp $qual;
substr($qual,0,1,"") if (length($qual)>length($seq));
warn "The length of quality is not match with seqence!\nSeq: $seq\nQual: $qual\n" if (length($qual)!= length($seq));
print OUT "$Name\n$seq\n+\n$qual\n" if ($seq ne "")&&($qual ne "");
($seq,$qual)=("","");
$Name="@".$Title.$1;
$_ = <IN1>;
if (/^[\@\>]/)
{
$_ = <IN1>;
}
my $tmp_seq;
while (($_ ne "")&&($_ !~ /^[\@\>]/)&&(!eof))
{
s/\s+$//;
$tmp_seq.=$_;
$_=<IN1>;
chomp $_;
}
$tmp_seq .= $_ if ((eof)&&($_ !~ /^[\@\>]/));
chomp $tmp_seq;
$seq = (($tmp_seq=~/\d+/)) ? col2base($tmp_seq) : $tmp_seq;
$_ = <IN2>;
s/^\s+//;
s/\s+$//;
if (($_=~/\d+\s+\d+/)||($_=~/^\d+$/))
{
$qual .= " " if ($qual=~/\d+$/);
my @t = split /\s+/,$_;
$qual .= $conv_table[$_+64] for (@t);
}
else
{
$qual .= $_;
}
}
else
{
s/^\s+//;
s/\s+$//;
if (($_=~/\d+\s+\d+/)||($_=~/^\d+$/))
{
$qual .= " " if ($qual=~/\d+$/);
my @t = split /\s+/,$_;
$qual .= $conv_table[$_+64] for (@t);
}
else
{
$qual .= $_;
}
}
}
}

substr($qual,0,1,"") if (length($qual)>length($seq));
print OUT "$Name\n$seq\n+\n$qual\n" if ($seq ne "")&&($qual ne "");
close IN1;
close IN2;
}
close OUT;

sub col2base
{
my $col = shift;
my @colors = split '',$col;
my $string = $colors[0];
if($string !~ /[acgt]/i){
warn "$col has no header base\n";
return 0;
}
my $last_base = $string;
my $current_base = '';
for(my $i=1;$i<@colors;$i++)
{
if (($last_base=~/N/i)&&(($colors[$i] eq ".")||($colors[$i]==5)))
{
$current_base = $bases[int(rand(@bases))];
}
else
{
$current_base = (exists $decode{$colors[$i]}->{$last_base}) ? $decode{$colors[$i]}->{$last_base} : "N";
}
$string .= $current_base;
$last_base = $current_base;
}
substr($string,0,1,"");
return $string;
}

Installing Samtools and BEDTools

2013-06-23T06:51:00.001-07:00

Trouble comes as a package when you start your lab. In my effort to resurrect our old database Eumicrobedb.org, I face new challenges every day. What earlier seemed to a be a cake walk, does not look that easy. I have installed so many open source programs in my last stint that I lost count of it. When running this software, I realize that many more dependencies need to be installed. Not to mention my troubles dealing with a crashed server and subsequent loss of all installations:((((

recently I was stuck with easy installations samtools-0.1.19 and bedtools-2.17.0 for more than 3 days. Earlier versions compiled fine in my SUSE linux machine but the same sources give trouble in Red Hat Linux.

While compiling BEDTools, it complained about "undefined reference to `gzopen64'". Searching the forums I figured out that it is complaining about zlib. I checked, it was there in my server and also in the path. However, finally changing the makefile with the following compiled the program.

Look for the line in Makefile 'export LIBS'

and add the following:
export LIBS = YOUR_PATH/libz.so.1.2.7 -lz

do a make clean and do make
It should compile fine.

Samtools:

With Samtools My problem lasted for a very long time. It always exited with error
samtools error bam_import.c:76: undefined reference to gzopen64

I tried re-installing zlib, added zlib path to LD_LIBRARY_PATH, installed latest version zlib, nothing worked.
Finally I changed the line in Makefile that says :

CFLAGS= -g -Wall -O2 to
CFLAGS= -g -Wall -O2 -L /usr/local/lib (This is where my zlib libraries are located)

viola.... It got installed.

NOTE: Install as regular user. May be later you can copy the binary files to the system
directories.

Installing CEGMA...

2013-05-19T02:24:00.001-07:00

Installing CEGMA can be a bit of a challenge for most of the users including me! I had installed it earlier, but this time it really gave me a very hard time in RH 6.2. Apart from installing a set of pre-requisitites, installation of genewise is bit of a challenge - especially so because now its status is archival. The INSTALL instruction that comes with the package is not very useful..

Genewise has several pre-requisites:

glib -> Provides low level data structure support for C programming. Can be found here:
http://www.linuxfromscratch.org/blfs/view/svn/general/glib2.html

glib has several pre-requisites such as
1. pkg-config: This package provides tools for passing include and library path so that packages can be built using configure and make utilities.

This package can be found here: http://pkgconfig.freedesktop.org/releases/pkg-config-0.28.tar.gz

2. libffi: This library provides high level portable interface for programmers to call any function at run time and can be found at: ftp://sourceware.org/pub/libffi/libffi-3.0.13.tar.gz

This package also requires a patch that can be found at:

http://www.linuxfromscratch.org/patches/blfs/svn/libffi-3.0.13-includedir-1.patch

3. Python2.7.5 if it is not already there: http://www.python.org/ftp/python/2.7.5/Python-2.7.5.tar.xz

4. PCRE library: This contains perl like regular expression libraries:
http://downloads.sourceforge.net/pcre/pcre-8.32.tar.bz2

All these installation will go smoothly, unless you have some problems with your system. Once done, go
Download wise2.4.1 source from http://wise.sourcearchive.com/downloads/2.4.1/

Do the following:

go to wise-2.4.1/src/HMMer2
then replace 'getline' to 'my_getline' in sqio.c
replace 'isnumber' in src/models/phasemodel.c into 'isdigit'

Then go and check the makefile under each directory under src and relace 'glib-config --libs' to 'pkg-config --libs glib-2.0' and also glib-config --cflags' to 'pkg-config --cflags glib-2.0'
using the following command:

find ./ -type f -name "makefile" -exec sed -i.old 's/glib-config --cflags/pkg-config --cflags glib-2.0/g' "{}" +;
and

find ./ -type f -name "makefile" -exec sed -i.old 's/glib-config --libs/pkg-config --libs glib-2.0/g' "{}" +;

Also set path for wiseconfigdir to wisecfg file that is somewhere inside the distribution.

Then do

make all

It works fine for RH6.2

Make sure blast program and Hmmer3 is on your path. Also install geneid in your system as that is also a pre-requisite. Once installed, first set the environment variables for the following:

CEGMA; CEGMATMP; PERL5LIB

PATH=$PATH:$HOME/bin:/usr/adadata/bin:/bin
export WISECONFIGDIR=/usr/adadata/wise2.4.0/wisecfg
export CEGMA=/usr/adadata/cegma_v2.4.010312
export CEGMATMP=/usr/adadata/cegma_v2.4.010312/tmp/
export PERL5LIB=$CEGMA/lib

Now finally Run CEGMA with the default parameters:

bin/cegma --genome sample/sample.dna --protein sample/sample.prot -o sample

References: http://korflab.ucdavis.edu/datasets/cegma/README

Creating Network File System for your servers

2013-04-18T23:50:00.001-07:00

Network File system is a great way to use your resources efficiently especially when they are connected to each other by network. There is a great article here [http://www.tldp.org/HOWTO/NFS-HOWTO/] that guides you about how to achieve this. For the impatient, I am putting a shorter step by step version here, so that you can quickly achieve it.

As a pre-requisite you need to have root access in all your servers you are trying to connect.

1. edit /etc/exports file to have the following entries:
/usr/data nnn.n.n.n(rw,sync)
/usr/data mmm.m.m.m(rw,sync)

where /usr/data in the local machine is the directory you want to share with machines with ip nnn.n.n.n and mmm.m.m.m. Here rw means read write, sync is to synchronise.
You can change this file in all the servers you want to share items with.

2. make entries in your /etc/hosts file which should read something like:
127.0.0.1 xxx.localdomain localhost
nnn.n.n.n yyy.localdomain -> where your machine has the name 'xxx' and the ip address of the machine you are trying to have in network is nnn.n.n.n and the name of that machine is: 'yyy'. This way you can keep adding the number of machines in the network that you wish to add.

3. In newer machines you can make the following entries into your /etc/hosts.deny file:
portmap:ALL
lockd:ALL
mountd:ALL
rquotad:ALL
statd:ALL
All these commands are harmless, and nothing happens if your system does not have a particular protocol. But it is a great way to protect your machine.
4. In /etc/hosts.allow file add the following lines:
service:nnn.n.n.n , mmm.m.m.m
lockd:nnn.n.n.n , mmm.m.m.m
rquotad:nnn.n.n.n , mmm.m.m.m
mountd:nnn.n.n.n , mmm.m.m.m
statd:nnn.n.n.n , mmm.m.m.m
5. Run rpcinfo -p to check if portmapper, rquotad, mountd, nfs, amd, status, nlockmgr is running.

(where nnn.n.n.n and mmm.m.m.m are the other two machines in your network you wish to give nfs access to)

6. Run exportfs -ra to make nfsd to read /etc/exports

7. Check /proc/filesystems and see if there is an entry on nfs. If it is there then you are good to go

8. Check if rpc.statd and rpc.lockd is running by doing ps aux | grep statd or ps aux | grep lockd

9. Finally mount the file system. Make sure you have mount package is installed in your computer.
be a root in machine xxx.x.x.x and you want to have /usr/data drive of machine nnn.n.n.n to be mounted in /share/data in your local file system -> do the following:
mount nnn.n.n.n:/usr/data /share/data
and keep adding mount drives for other machines as well.

Do the same for other machines if you want to have drives from your machine xxx.x.x.x to be mounted else where. Now finally you are done!!

Happy creating network file system....

Training Augustus gene predictor for your organism

2013-03-22T17:55:00.001-07:00

Lately I have been asked by multiple people to solve the training problem of Augustus for their organism data. Although I have done it earlier, this time, I faced unusually long time in solving this issue. Augustus wants you to provide a training dataset either in genbank format or in gff format. It typically wants users to remove duplicate genes; proteins with 70% identity, in order to get over the over fitting problem. Also, remove any kind of overlapping genes from your training files. If you could not find a script that can do it efficiently, please contact me directly I will provide you one (sorry, I am yet to set my source code page up).

I have generated a genbank file from my gff and the annotation files (again script can be sent upon request). Howver, few things I have overlooked and they were leaving empty lines between 2 genomic LOCUS entries. If empty space is left then the randomSplitter.pl that is used to split the gb training file into 2 random parts fails. You either get all your data in .test file or all of them in .train file. Suppose you decide to move on with your .train file with all the data and .test as a empty file. Surprises that you will end up with std::bad_alloc() error from augustus, that has nothing to do with the memory of your system. So, if you have an autogenerated genbank file with empty lines, open the file in vim editor and get rid of empty lines using :g/^$/d.

Once done, try splitting your *.gb file with randomSplit.pl. This will do the magic. Now try to train your program and it should run file.

Here is a check list of commands what you should be doing to train for Augustus.

perl scripts/new_species.pl --species=sojae
perl scripts/randomSplit.pl psv5.gb 30
./bin/etraining --species=sojae psv5.gb.train
# mock prediction
./bin/augustus --species=sojae psv5.gb.test | tee firsttest.out
grep -A 22 Evaluation firsttest.out

This is what you get when you run the last command. This is a pretty good indicator that your training is almost good to go.

******* Evaluation of gene prediction *******

---------------------------------------------\
| sensitivity | specificity |
---------------------------------------------|
nucleotide level | 0.861 | 0.674 |
---------------------------------------------/

--------------------------------------------------------------------------------
--------------------------\
| #pred | #anno | | FP = false pos. | FN = false neg. |
| |
| total/ | total/ | TP |--------------------|--------------------|
sensitivity | specificity |
| unique | unique | | part | ovlp | wrng | part | ovlp | wrng |
| |
--------------------------------------------------------------------------------
--------------------------|
| | | | 24559 | 20733 |
| |
exon level | 34213 | 30387 | 9654 | ------------------ | ------------------ |
0.318 | 0.282 |
| 34213 | 30387 | | 8709 | 4167 | 11683 | 9312 | 5502 | 5919 |
| |
--------------------------------------------------------------------------------
--------------------------/

----------------------------------------------------------------------------\
transcript | #pred | #anno | TP | FP | FN | sensitivity | specificity |
----------------------------------------------------------------------------|
gene level | 10348 | 12700 | 2547 | 7801 | 10153 | 0.201 | 0.246 |
----------------------------------------------------------------------------/

Sorry this output is jumbled but in essence it tells that 2547 genes were predicted exactly. Although exon and gene level results does not look that good, you could actually try and train the software using various types of exons that are accurate.

perl scripts/optimize_augustus.pl psv5.gb.train # This will take a while....

Resetting forgotten password in RH enterprise Linux version

2012-12-04T21:59:00.002-08:00

We have installed the latest version Redhat.Enterprise.Linux.v6.UPDATE.2.X86_64 on one of our servers. Probably the grub prompt time setting is set to 0, so we never see a grub prompt when the computer booted. Numerous posts indicate that one has to set the run level to 1 during boot time to set the kernel settings to single user mode, then just type a password. Howver, in our machine, we could not get that prompt at all.

Yesterday, one of my students, Dhritiman Jana was in the lab and he spent around few hours of his time in solving the issue. Following are his note on how to get to the grub prompt:
1. During boot time tap any key and enter into the grub prompt.
2. Goto kernel entry and press e to edit. After "rhgb quiet" enter "1"
to set the runlevel to single user mode.
3. Press b to boot the edited kernel entry.
4. Reset Password by using #passwd

Hope this solution works for all!

When your web server gets attacked.....

2012-11-20T07:26:00.001-08:00

This post may be a mixture of rants and some quick tips on how to save your web server when you suspect there is some malicious attack on it and you have limited admin access and support.

It becomes imperative to create a world writable directory if you have a web server that plots, maps data on user demand. I would like to hear more from people if they knew how to do this otherwise. I have been looking for solutions into this problem, but did not find any. Now I have settled for something that makes your web program clean its own mess.

Suppose say you have a world writable directory called /tmp/XX and some hacker has been able to put in some malicious code into it and is being able to run it. Now my solution for that is to delete the content of this folder after n hours. Just put in a small command inside your web program that runs frequently, and when that runs, the world writable folder will be cleaned.

# To delete all files older than 6 hours in /tmp directory

$command="find /tmp -type f -name \"*.*\" -mmin +360 -exec rm {} \\;";
system($command);

Will probably partly solve your problem.

Creating login/password less ssh connection with multiple remote machines

2012-11-06T01:48:00.000-08:00

Many times while working in multiple machines, you may like to automate certain processes where the programs can directly access information/data from another server effortlessly - without a password. This can be done using ssh-keygen protocol.

So, what happens here is; you have a local machine, lets call it 'A' and you have a remote machine, lets call it as 'B'. You have an account in 'B' and that is say 'myname'. Everytime you log into that machine using ssh, you have to do something like:

$ ssh myname@B

$Password:

$myname@B:

In order to directly log into a machine without a password, you have to generate a pair of keys; called as a public key and a private key. The public key is the public information and the private key is only known to your local machine i.e; 'A'. You can use ssh-keygen to create a pair of keys in a given time. This is how you should proceed:

$ ssh-keygen -t rsa

Generating public/private rsa key pair.
Enter file in which to save the key (/home/Sucheta/.ssh/id_rsa): /home/Sucheta/.
ssh/iicb_rsa [ Remember here to enter a new file name, else it will over write any other id_rsa file that you may have saved earlier for any other computer ]
Enter passphrase (empty for no passphrase): [Enter a paraphrase that is > 4 characters long. This is essential if your computer has more than one pair of public private keys for more than one remote server]
Enter same passphrase again:
Your identification has been saved in /home/Sucheta/.ssh/iicb_rsa.
Your public key has been saved in /home/Sucheta/.ssh/iicb_rsa.pub.
The key fingerprint is:
95:13:96:1b:66:ef:36:74:25:76:05:23:64:58:bb:94 Sucheta@Sucheta-PC
The key's randomart image is:
+--[ RSA 2048]----+
|          o== o.o|
|         .*+ +o.o|
|         o++E. + |
|         ..oo.. |
|        S o..   |
|            +    |
|           . .   |
|                 |
|                 |
+-----------------+

Then do:

$ssh-copy-id myname@B

[This command will append your public key in the ~/.ssh/authorized_keys file in the remote host. You can also do this manually by logging back to your remote computer and copy pasting your public key in the 'authorized_keys' file. Make sure your public key is copy pasted in one single line.

Another thing to remember is, depending on the OS and version, the file that needs to have the public key in the remote machine may be different. In order to confirm that it is indeed called as "authorized_keys" do the following:
[root@Apala ssh]# cat /etc/ssh/sshd_config | grep Keys
# HostKeys for protocol version 2
#AuthorizedKeysFile .ssh/authorized_keys
#AuthorizedKeysCommand none
#AuthorizedKeysCommandRunAs nobody

This tells you indeed the file that stores public key in your remote computer is named as authorized_keys .

Next time you try to create another passwordless connection to another computer, just repeat the above steps. And always remember to write the public and private key into different files, else it will overwrite contents into id_rsa and id_rsa.pub file. Add a passphrase too.

One more important thing to remember is to check the file permission for "authorized_keys" file in the remote machine. Always set it to 700.

Using this, you can also automate file transfer by an sftp or any other remote ftp protocols

Some of the Query strategies for Oomycetes Genomes in FungiDB

2012-08-08T00:48:00.002-07:00

Querying Oomycetes genomes in FungiDB 2.0

FungiDB 2.0 released with 6 Oomycetes genomes:

Phytophthora sojae V5.0

Phytophthora ramorum V1.0

Phytophthora capsici V11.0

Phytophthora infestans V4.0

Pythium ultimum V 2.0

Hyaloperonospora arabidopsidis V8.3

Naming conventions:

Gene ids:

There are 2 types of standards followed in naming of oomycetes gene ids in FungiDB. In most of the cases, the gene ids are represented by a 5 letter word followed by an underscore, followed by the gene identifier devised by the sequencing center e.g; Physo_517720, Phyra_74442, Hyaar_813319, Phyca_96628, Pytul_G005233. Genomes sequenced and annotated at Broad (P. infestans, P. parasitica), already have gene ids prefixed by a 4 letter code such as PITG and we have left those as is e.g; PITG_05520 .

Scaffolds:

Since most of these genomes are released as draft assemblies, their genome fasta files are named as follows: a 5 letter organism prefix + strain name _ ‘SC’XXXX, where SC stands for super Contig and XXXX are the 4 digit representation of the scaffold number. Example: PytulBR144_SC1841, PhysoP6497_SC0001, PhyraPr102_SC0008, PhycaLT1534_SC0024, PhyinT30-4_SC0007, HyaarEmoy2_SC0165 etc.

Search Options:

Figure - 1

In the main FungiDB page, there are three search columns available(Figure 1). The search results of these columns return different result types e.g; the first column almost always returns gene records page, the second column returns other data types and the third column has links to tools such as blast, genome browser etc. On the top right hand corner of the page, there are quick search options available for search by gene_id and by gene product name.

Few of the Search strategies are listed below:

1. Curate all infection related genes in Phytophthora sojae as described in the science paper (Tyler et al, 2006) with key word searches such as serine proteases, Metalloproteases, Cysteine proteases, glycosyl hydrolases, pectinesterases, pectate lyases, cutinaes, chitinases, lipases, phospholipases, protease inhibitors, NPP family, PcF family, Six Cys Family, Eight Cys family, Crn Family, nonribosomal peptide synthetases, Polyketide synthases, Cytochrome P450s, CYP51 clan, ABC transporters, PDR, ABCG-half, MDR, MRP, elicitins, Avhs, Crinklers.

Hint: Start your search by Text (product name, Gene Id) search for proteases. After the results are displayed, run another Text(Product name, Gene Id) search and use the next item in your gene list e.g; glycosyl hydrolases. Choose union of both searches and keep repeating this until you are done with your list. Following is the strategy for curating all the infection related genes in Phytophthora sojae.

http://fungidb.org/fungidb/im.do?s=5e276dc38a69ece3

2. Find all the Oomycetes proteins that have an RXLR motif within first 20-50 residues , followed by a dEER motif within 20,60 residues of the first motif. Also see if they have predicted signal peptide leaders and transmembrane domains. (start with similarity-> protein motif search and use the following pattern ^.{20,50}R.LR.{20,60}DEER. Then add a step to it on whether the protein has a signal peptide followed. Then add another step on if the secrted proteins have a transmembrane domain. The stratergy can be found here:

http://fungidb.org/fungidb/im.do?s=53a175bca66b25a7

3. Genomic locations: All the Oomycetes genomes are draft sequences, so searching genomic locations can be little tricky. Click on genome location search (as in figure 2 A) and this will open a page where the pull down menu only displays the genomes that have the complete chromosomes available. Instead, choose genomic sequence id (as in figure 2 B) and fill up the start and stop locations you are interested in. Run the query and go to add step. If you are interested in seeing how many of these genes have orthologs in other oomycetes, just click on add step and click on Eolution->orthology, Phylogenetic profile and choose Oomycetes genomes (Figure 2 C). This results in 32 genes belonging to 25 ortholog groups. Clicking on the ortholog groups and it will display a cascade of glycosyl Hydrolases (Figure 2 D) belonging to different Oomycetes. From the gene_ids, it appears that the genes occur in a syntenic block.

Strategy can be found here: http://fungidb.org/fungidb/im.do?s=dce0a3dacc71bb06

HMMs what they are good at....

2012-06-25T12:47:00.000-07:00

As biologists most of us feel little bit awkward with mathematics and statistics, and more so with the application of mathematics and statistics in interpreting biological data. Recently we see a meteoric rise in availability of biological sequences. However, interpreting them properly still is a bottleneck.
HMMs are there for a long time now, very valuable in speech recognition, weather prediction etc. But recently, we find a surge in application of HMMs in interpreting biological sequences such as gene prediction, profile searches, domain prediction etc. A beginner however gets awed by the complexity and finds it difficult to understand completely how this amazing method/tool can help interpret sequences. Here is my attempt to simplify things a bit.

What is HMM?

Profile HMMs are based on position specific scoring matrices. Lets take the example of the following diagram where we have multiple sequences aligned.
The following are the characteristics:
1. Columns in a multiple sequence alignment is called as a state.
2. Number of alphabets taken together is called as order of the markov chain e.g; zero order, first order and so on.
3. PSSM(Position specific scoring matrix) is calculated for calculating the match states probabilities.

In a HMM, there is a transition probability and an emission probability and there is a state. Lets look at the simplest HMM heuristic:

Suppose say we have a multiple sequence alignment, as represented in Figure - 1.

Figure 1

[ state -> is the column in a multiple sequence alignment.
Transition probabilities are determined using deletions and insertions. Emission probabilities at each state are determined by counting the occurrence of each character in each column. ]

So, lets try to generate the Emission and Transition probabilities from this alignment from the figure.

Emission Probability:
                  A                 T                   G                  C
Pos1          1.0               0                    0                   0

Pos2          0                  0                    0.2                0.8
Pos3          0.2               0                    0.4                0
Pos4          0.6               0                    0.2                0

Transition Probability:( There are 3 transition states, match, deletion, insertion. So, there can be 9 + 2 transition states:
                     m->m     m->i     m->d     i->m    i->i    i->d    d->m      d->i    d->d      s->m       m->e
Trans 1         0             0           0            0         0        0         0             0         0           1             0
Trans 2         1             0           0            0         0        0         0             0         0           0             0
Trans 3         0.6    0 0.4 0 0    0 0 0 0 0 0
Trans 4         0.8          0           0.2         0         0        0         0             0         0           0             0
Trans 5         0             0           0            0         0        0         0             0         0           0             1
[In this table m stands for match, i stands for insertion and d for deletion, s for start and e for end]

Constructing the HMM:

There are multiple algorithms that can be used to build and train HMMs.

i) EM, or Baum-Welsch, algorithm.

This algorithm uses a random MSA, such as the one constructed by ClustalW. The number of match states is defined as the number of conserved columns. This algorithm uses a two stage, iterative process; first, it uses the aligned residues to build a model, and then it realigns the sequences to the model. These two steps are repeated until the model no longer changes.

ii)Viterbi algorithm:

Uses dynamic programming to calculate the best estimate of the emission probabilities.

iii)“Surgery” algorithm:

Uses dynamic adjustments of the HMM length to calculate the best models.

HMMER software comes with a nice user guide. I recommend you read this guide from cover to cover and if you still did not get it, then continue reading my post.

Use of HMM in aligning against another sequence or profile:
After we create the HMM, we can use it to search for sequences that are similar to the ones used to create the HMM. So, we use the model to estimate the probability that a sequence belongs to the model. This probability is often expressed as negative log likelihood score and is dependent on the length of the sequence and the length of the model. Once we have a score, we need a way to assess its significance, just as we did for pairwise sequence alignments.

HMMs have some distinct advantages. They are built on a formal probabilistic basis. They also allow for more sensitive searching than pairwise alignments do. The probabilistic theory can be used to guide the scoring parameters; the theory also allows us to train an HMM using unaligned sequences if we don’t know or trust the alignment. There is a consistent theory behind gap and insertion penalties. Finally, less skill and human intervention is needed to train a good HMM than to train a hand-constructed profile. Because of this, we can create libraries of hundreds of profile HMMs and apply them on a large scale.

Its one major drawback is that its probabilistic basis does not capture higher-order correlations between sequence positions. It assumes that the identity of a particular position is independent of the identity of all of the other positions in the sequence, and we know that not to be true.

Publicly available program HMMER3 just does this.

Components of Hmmer3 Program:
phmmer: Makes a blastp like search.
jackhmmer: Iteratively searches a sequence against a sequence database (something similar to Psi-blast).
hmmbuild: Builds a profile HMM from multiple sequence alignment.
hmmsearch: Searches a profile hmm against a sequence database(Query: HMM, target: AA).
hmmscan: Searches a sequence against a profile hmm database(Query: AA; target: HMM e.g; pfam).
hmmalign: Make a multiple sequence alignment of many sequences to a existing profile.

hmmconvert: Converts profiles into hmm format and vice versa.
hmmemit: Generate sample sequences from profiles.
hmmfetch: Get a profile HMM by name or accession number from the HMM database
hmmpress: Format an hmm database into a binary format to be used by hmmscan.
hmmstat: Show summary statistics for each profile in an hmm database.

In the next post, we will discuss about Applications of Profile search and command line options of Hmmer3 with examples.

Job opportunity as an environment intern!!!

2012-06-07T10:18:00.001-07:00

I got this email to my official email account. Thought will share it here..

www.TerraScapes.org

Environmental Organizer Intern

Job Description:

Do you want to incorporate the environment into your everyday life through environmental activism, advocacy, social justice and event planning? Interning for Terra-Scapes Environmental Consulting is a unique opportunity for a student with a passion for the environment. Engaging the community through educational initiatives about current environmental concerns will be your goal. The Environmental Organizer will work on several pertinent environmental campaigns ranging from protecting our oceans and fisheries to the prohibition of additional chicken farms on the Eastern Shore.

The Environmental Organizer will have the opportunity to gain experience and forge connections in the environmental and green business communities. The position entails working with the following prominent environmental organizations: Pew Environment Group, Sierra Club, Chesapeake Bay Foundation, Lynnhaven River Now, Nature Conservancy, and more.

Responsibilities:

· Research, organize, and execute logistics for environmental workshops and events.

· Recruit and communicate with new and existing volunteers.

· Strengthen and build coalitions with new and existing groups, scientists, students, elected officials, and businesses.

· Work with the media on press releases, advisories, letters to the editor, and opinion editorials.

· Schedule and participate in meetings with local, state, and federal legislators regarding environmental policy initiatives.

· Database organization and entry using Google docs and Excel spreadsheets.

· Research and preparation of reports, documents, and memorandums using Microsoft Office.

· Assist in creation of materials for Terra-Scapes and its coalition partners.

· Develop plan to reach goals with each individual campaign.

Qualifications:

· Grasp office-related computer (Mac and PC) and internet applications.

· Have a passion for the environment and the growth of the green economy.

· Ability to prioritize, work independently, and follow through on concurrent projects and goals.

· Superior organizational skills and ability to meet deadlines required.

· Ability to communicate easily with a diverse range of people.

· Demonstrated experience working independently.

· Owns a vehicle and willing to travel throughout VA for specific workshops or conferences.

· Willingness to meet with Congressional members in Washington, D.C.

· Ability to work from your home/dorm.

Location: Primarily in Hampton Roads, but requires some travel throughout Virginia, North Carolina, Maryland, and Washington, D.C.

To apply, please email resume and cover letter to terra@terrascapes.org

Thanks again!

Terra Pascarosa Duff

TerraScapes Environmental Director

VA Pew Environment Group Representative

Chesapeake Bay Group Sierra Club Chair

www.TerraScapes.org

http://www.pewenvironment.org/

p: 757-291-4736 | f: 866-543-6696 | e: terra@terrascapes.org

"You must be the change you wish to see in the world."
— Mahatma Gandhi

Extracting Fastq files from bam/sam files

2012-05-17T09:11:00.000-07:00

Earlier I have blogged about converting colorspace files into fastq format, which is pretty easy [here]. However, this method does not go without the risk of getting completely wrong reads if there is an error in the read number in the middle.

There is another approach to it which is reliable - that is to extract the fastq files from the sam/bam alignment files. There are quite a few different ways you can achieve this. One of the easiest is to use samtools to extract matching and unmatching reads separately and merging them together. Something like this:

# Extracting mapping and non-mapping reads from bam files
samtools view -f4 input.bam > /data/bowtieOutput/unmapped.sam
# This is for extracting mapping reads from bam files
samtools view -F4 input.bam > /data/bowtieOutput/mapped.sam

Remember, this method extracts files in sam format, so you may have to get into another step to extract fastq files from sam files

But the only problem with this approach is, you end up getting reads that are mapped to the reference sequence multiple times in multiple numbers. For unmapped reads this problem does not exist.

There is another easier approach is also to use bamtofastq . This is a cool software available here.

Simply use command:

# Extracts aligned reads
./bam2fastq --no-unaligned -o out.fastq input.bam
#Extracts un-aligned reads
./bam2fastq --no-aligned -o out.fastq1 input.bam

This is one of the fastest and easiest methods to extract mapping and unmapping reads from bam files into fastq format. However, this also has similar problem like the sam extractor - that is, it extracts reads mapped to multiple locations multiple times. In order to get around this the fastest way would be to make a uniq set from the aligned.fastq file. Here are the steps that needs to be followed to get it done faster.

1. Sort your bam file before running bam2fastq on them. This way the output fastq file will have all the redundant reads in a sorted order and hence easier to remove redundant.
2. Run bam2fastq on them as described above.
3. Write a small perl script to fish out the uniq reads from the fastq output file as below:

open FH, $ARGV[0] or die "Can't open file for reading $!\n";

my $prev="";

while(<FH>){

        if(/\@(\d+_\d+_\d+_)/){

                if($prev eq $1){
                        # Skips reading 3 lines
                        for(0..2){
                         <FH>
                        }
                next;
                }
        $prev = $1;
        }

        print $_;
}

My Results:

                         Reads recovered from bam2fastq               uniq reads
                         ----------------------------------------               -------------
expt 1                  23921306                                          8909517
expt 2                  22685196                                        8924548
expt 3                  40484226                                          15248005
expt 4                  13659968                                          5141253
expt 5                  17185935                                               6862160
expt 6              19855362    8978295
expt 7                  13189432                                                     7067910
expt 8                  14423653 6713690

This is quite consistent with what I wanted to achieve!!

How to code position fixed content on a browser.

2012-04-30T08:54:00.000-07:00

While representing sequence text alignment on a browser, it is very useful to have a static genome reference on the top of the page, that descends as you scroll down. So, you are actually not loosing track of the reference but at the same time you see the read alignments and their position with respect to the genome. The genome references can be long and the browser width is not enough to contain it in a single view. So, there will be a lot of room for horizontal scrolling. In this case you want the reference on the top behave normally as with any horizontal scroll, but vertically it should be position fixed. Here is an example of what I am talking about:

I have created this application after realizing its importance in dealing with nextgen sequence alignment viewer. Here are the few fixed libraries, you may find incredibly useful, so better install/download these in your local server.

A large part of the application have been downloaded from mootools site. The following are the files I have downloaded from mootools site:

mootools-yui-compressed.js

From http://code.google.com/p/joomla-scripts-optimized/downloads/detail?name=mootools-yui-compressed.js&can=2&q=

DemosMain.php from mootools web site: (If you don't find it in the web site, the source is copy pasted below)

Style.php from mootools site (Source code is copy pasted below)

Then bundle up these 3 scripts into a directory say textStyle.

Now create a Php file and write functions that will print your object of interest Some like this:

 open FH, $fileName or die "Can't open file for reading $! \n";  
 my $count=0;  
 my $start;  
 &printHTMLcode;  
 print "<body>";  
 print "<pre>";  
 while(<FH>){  
     if(/(\d+)\t(\d+)\t(\S+)/){  
         if($count == 0){  
             if($genomeStart >= $2){  
             print "<div id=\"DemoItem2\">";  
                 for(my $i=0; $i<=$genomeStart - $2; $i++){  
                     print "&nbsp;";  
                 }  
             &colorCode($genomeSeq);  
             print "</div><br><br>";  
             $start = $2;  
             }  
             else{  
             print "<div id=\"DemoItem2\">";  
             &colorCode($genomeSeq);  
             $start = $genomeStart+1;  
             print "</div><br><br>";  
             }  
         }  
         for(my $i=0; $i<=($2-$start);$i++){  
             print "&nbsp;";  
         }  
     &colorCode($3);  
     }  
 print "\n";  
 $count++;  
 }  
 print "<script src=\"http://www.eumicrobedb.org/develop/textStyle/mootools-yui-c  
 ompressed.js\" type=\"text/javascript\">   \  
     </script>                  \  
     <script src=\"http://www.eumicrobedb.org/develop/textStyle/demosMain.php  
 \" type=\"text/javascript\">                     \  
     </script>                  \  
 ";  
 print "</pre>";  
 print "</body>";  
 close(FH);  
 }  
 sub printHTMLcode  
 {  
 print "<html>";  
 print "<head>";  
 print "<link rel=\"stylesheet\" href=\"http://www.eumicrobedb.org/develop/textSt  
 yle/style.php\" type=\"text/css\" media=\"screen\">";  
 print "<style type=\"text/css\">  
     #DemoItem2{  
         position: fixed;  
         top   : 10px;  
         left  : 6px;  
         width  : 200px;  
         height : 50px;  
     }  
     #DemoItem2{  
         top: 10px;  
     }  
     </style>";  
 print " <script text=\"text/javascript\">  
         function loadDemo(){  
          var demoTwo = new ScrollSpy({  
                  min: 0,  
                  mode: 'horizontal',  
                  onEnter: function(position,enters) {  
                  },  
                  onLeave: function(position,leaves) {  
                  },  
                  onTick: function(position,state,enters,leaves)  
 {  
                 \$(\"DemoItem2\").style.left = -position.x+ 6+\"  
 px\";  
                  },  
                  container: window  
                 });  
         }  
     </script>  
     </head>  
     <body>  
 \n";  
 }

I will explain now what all these means.

The file it opens in the while loop is a data file that has sequence and their starting point info. So ignore that part.

There is a function call to printHTML code (marked in red). This function is a style sheet text wrapped. You don't have to wrap it in a function, you could just write it on the top if you like. But it looks messy. Please note that inside printHTML code, I have defined demoItem2 div class that actually does the text/content to scroll down. The main text that scrolls down is within &colorcode function, that actually colors and prints the genome reference that is suppose to be static on the screen.

In the end I am calling the moo tools js files(marked in yellow) that actually lets the content scroll horizontally.

here are the source code for style.php and demomain.php, if you did not find it from the mootools:

Source for style.php

 <?php  
 print "<br>  
 *{<br>  
     margin: 0px;<br>  
     padding: 0px;<br>  
 }<br>  
 <br>  
     body, html{<br>  
         width: 100%;<br>  
         height: 100%;<br>  
 <br>  
     }<br>  
 <br>  
 body{<br>  
     font-family:Gill Sans,sans-serif;<br>  
     background: #2a2a2a;  <br>  
 }<br>  
 <br>  
 #wrap{<br>  
     width: 100%;<br>  
     height: 100%;<br>  
     min-width: 1000px;<br>  
 }<br>  
 <br>  
 #head{<br>  
 height:156px;<br>  
 position:fixed;<br>  
 width:100%;<br>  
 top: 0px;<br>  
 }<br>  
 <br>  
 #main{<br>  
     padding-top:148px;<br>  
 }<br>  
 <br>  
 #logo{<br>  
     background-image:url(images/logo.png);<br>  
     background-position:0 0;<br>  
     height:108px;<br>  
     left:50px;<br>  
     position:absolute;<br>  
     top:26px;<br>  
     width:90px;<br>  
 }<br>  
 <br>  
 #return{<br>  
     color:#FAC500;<br>  
     font-size:13px;<br>  
     position:absolute;<br>  
     right:50px;<br>  
     text-decoration:none;<br>  
     top:122px;<br>  
 } <br>  
 #return:hover{<br>  
     color: #fff;<br>  
 }<br>  
 <br>  
 img, a{<br>  
     outline: none;<br>  
 }<br>  
 <br>  
 h1{<br>  
     color:#FFFFFF;<br>  
     font-family:Helvetica,\"Helvetica Neue\",Arial,sans-serif;<br>  
     font-size:18px;<br>  
     font-weight:normal;<br>  
     letter-spacing:-1px;<br>  
     position:absolute;<br>  
     right:50px;<br>  
     text-shadow:2px 2px 0 #111111;<br>  
     top:95px;<br>  
 }<br>  
 <br>  
 .hr{<br>  
     border-top:1px dashed #555555;<br>  
     height:1px;<br>  
     position:absolute;<br>  
     top:153px;<br>  
     width:100%;<br>  
 }<br>  
 <br>  
 <br>  
 .TabHolder{<br>  
     height: 26px;<br>  
     width: auto;<br>  
     position: fixed;<br>  
     z-index: 100000;<br>  
 }<br>  
 <br>  
 .leftTab, .rightTab{<br>  
     display: inline-block;<br>  
     width: 17px;<br>  
     background-image: url(\"images/leftTabDarker.png\");<br>  
     height: 26px;<br>  
     margin: 0px;<br>  
     float: left;<br>  
 }<br>  
 <br>  
 .rightTab{<br>  
         background-image: url(\"images/rightTabDarker.png\");<br>  
 }<br>  
 <br>  
 .midTab{<br>  
     display: inline-block;<br>  
     width: 17px;<br>  
     background-image: url(\"images/midTabDarker.png\");<br>  
     height: 26px;<br>  
     margin: 0px;<br>  
     min-width: 200px;<br>  
     width: auto;<br>  
     float: left;<br>  
     font-size: 11px;<br>  
     font-family: 'Lucida Sans Unicode', 'Lucida Grande', sans-serif;<br>  
     line-height: 23px;<br>  
     text-shadow: 1px 1px 0px #fff;<br>  
     padding: 0px 2px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol{<br>  
     text-decoration: none;<br>  
     color: #444;<br>  
     border-right: 1px solid #999;<br>  
     border-left: 1px solid #fff;<br>  
     padding: 0px 10px;<br>  
     display: inline-block;<br>  
     line-height: 12px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol:first-child{<br>  
     border-left: 0px;<br>  
 }<br>  
 <br>  
 .midTab .midTabcol:last-child, .midTab .lastmidTabcol{<br>  
     border-right: 0px;<br>  
 }<br>  
 <br>  
 .midTab a:hover{<br>  
      color: #111;<br>  
 }<br>  
 <br>  
 .midTab a:focus{<br>  
      position: relative;<br>  
      top: 1px;<br>  
 }<br>  
 <br>  
 .goodBrowsers{<br>  
         display:inline-block;<br>  
         position: relative;<br>  
         top: 3px;<br>  
 }<br>  
 <br>  
 #footerTab{<br>  
     bottom: 0px;<br>  
     right: 100px;<br>  
     bottom: -27px;<br>  
 }<br>  
 <br>  
 <br>  
 ";  
 ?>

Source code for DemosMain.php

 var Page = {  
     initialize: function() {  
          // nice one David Walsh  
          var ss = new ScrollSpy({  
                  min: 0,  
                  mode: 'vertical',  
                  onEnter: function(position,enters) {  
                  },  
                  onLeave: function(position,leaves) {  
                  },  
                  onTick: function(position,state,enters,leaves)  
 {  
                     $("head").style.top = -position.y+"px";  
                  },  
                  container: window  
                 });  
                 if($("footerTab")){  
                     $("footerTab").set('morph', {  
                     duration: 500,  
                     transition: 'Sine:out'  
                 });  
                     setTimeout(function(){  
                             $("footerTab").morph({  
                     'bottom': ['-27px','0px']  
                     });  
                     },750);  
                 }  
                 loadDemo();  
     }  
  }  
  window.addEvent("domready", Page.initialize);  
 var ScrollSpy = new Class({  
     /* implements */  
     Implements: [Options,Events],  
     /* options */  
     options: {  
         min: 0,  
         mode: 'vertical',  
         max: 0,  
         container: window,  
         onEnter: $empty,  
         onLeave: $empty,  
         onTick: $empty  
     },  
     /* initialization */  
     initialize: function(options) {  
         /* set options */  
         this.setOptions(options);  
         this.container = $(this.options.container);  
         this.enters = this.leaves = 0;  
         this.max = this.options.max;  
         /* fix max */  
         if(this.max == 0)  
         {  
             var ss = this.container.getScrollSize();  
             this.max = this.options.mode == 'vertical' ? ss.y : ss.x  
 ;  
         }  
         /* make it happen */  
         this.addListener();  
     },  
     /* a method that does whatever you want */  
     addListener: function() {  
         /* state trackers */  
         this.inside = false;  
         this.container.addEvent('scroll',function() {  
             /* if it has reached the level */  
             var position = this.container.getScroll();  
             var xy = this.options.mode == 'vertical' ? position.y :  
 position.x;  
             /* if we reach the minimum and are still below the max..  
 . */  
             if(xy >= this.options.min && xy <= this.max) {  
                     /* trigger Enter event if necessary */  
                     if(!this.inside) {  
                         /* record as inside */  
                         this.inside = true;  
                         this.enters++;  
                         /* fire enter event */  
                         this.fireEvent('enter',[position  
 ,this.enters]);  
                     }  
                     /* trigger the "tick", always */  
                     this.fireEvent('tick',[position,this.ins  
 ide,this.enters,this.leaves]);  
             }  
             else {  
                 /* trigger leave */  
                 if(this.inside)  
                 {  
                     this.inside = false;  
                     this.leaves++;  
                     this.fireEvent('leave',[position,this.le  
 aves]);  
                 }  
             }  
         }.bind(this));  
     }  
 });

Making 300 dpi images from powerpoint slides

2012-04-17T08:58:00.001-07:00

Publication quality images are usually high resolution images with 300 dpi density. For those using screen shots, copy pasting images on powerpoint slides and later saving them into JPEG or any other file format, usually get the image saved as a 96 dpi. I use gimp(GNU Image manipulating program) as an image editing software for my papers. Gimp has options to increase the resolution. You can go to file -> open -> <your image>

Once the image file opened, go to Image and select scaleImage. Change resolution with more pixel density and then say 'scale'. Finally save the image. Now you will notice, the image file is bigger with more pixel density, but it has not made the image any clearer. When you open the file, you may see unwanted black dots here and there, making the image looks more messy. There is a very cool way to this by editing the windows registry. NOTE: before you start doing this, make sure, you have windows service pack-2 installed.

For that open a command line option in windows something like this:

Then from the left menu, click on HKEY_CURRENT_USER -> Software -> microsoft -> office -> 12.0 -> powerpoint -> options . At this stage there will be some binary values already on the right pane. Then go to Edit -> new -> DWORD (32 bit value) and type in any name such as ExportBitMapResolution. Once created, this will appear on your right hand side binary lists. Then click on the ExportBitMatResolution and then set the value as 300 and select decimal. Then come out of it.

Now whenever you try to save an image from a powerpoint slide, it will by default save it as 300 dpi
as against 96 dpi.