I read this article liked it a lot. So, copy pasting it for future reference:
http://www.nature.com/scitable/topicpage/Small-Non-coding-RNA-and-Gene-Expression-1078
Friday, April 03, 2009
Monday, March 30, 2009
Fungus farmers show way to new drugs
[I came across this interesting story in nature news, thought I will put it in my blog:]
Ant colonies could be key to advances in biofuels and antibiotics.
-Erika Check Hayden
Studies of bacteria on leaf-cutting ants could yield new antibiotics.Studies of bacteria on leaf-cutting ants could yield new antibiotics.M. MOFFETT/FLPA
In a mutually beneficial symbiosis, leaf-cutting ants cultivate fungus gardens, providing both a safe home for the fungi and a food source for the ants. But this 50-million-year-old relationship also includes microbes that new research shows could help speed the quest to develop better antibiotics and biofuels.
Ten years ago, Cameron Currie, a microbial ecologist then at the University of Toronto in Ontario, Canada, discovered that leaf-cutting ants carry colonies of actinomycete bacteria on their bodies (C. R. Currie et al. Nature 398, 701–704; 1999). The bacteria churn out an antibiotic that protects the ants' fungal crops from associated parasitic fungi (such as Escovopsis). On 29 March, Currie, Jon Clardy at the Harvard Medical School in Boston and their colleagues reported that they had isolated and purified one of these antifungals, which they named dentigerumycin, and that it is a chemical that has never been previously reported (D.-C. Oh et al. Nature Chem. Bio. doi: 10.1038/nchembio.159; 2009). The antifungal slowed the growth of a drug-resistant strain of the fungus Candida albicans, which causes yeast infections in people.
“These ants are walking pharmaceutical factories.”
Because distinct ant species cultivate different fungal crops, which in turn fall prey to specialized parasites, researchers hope that they will learn how to make better antibiotics by studying how the bacteria have adapted to fight the parasite in an ancient evolutionary arms race. "These ants are walking pharmaceutical factories," says Currie, now at the University of Wisconsin, Madison.
That's not the end to the possible applications. The ant colonies are also miniature biofuel reactors, Currie reported on 25 March at the Genomics of Energy & Environment meeting at the Joint Genome Institute in Walnut Creek, California. Each year, ants from a single colony harvest up to 400 kilograms of leaves to feed their fungal partners. But no one has worked out how the fungi digest the leaves, because samples of fungus grown in petri dishes can't break down cellulose, a tough molecule found in plant cells. Researchers are keenly interested in better ways to break down cellulose, because it might allow them to make more efficient biofuels than those made from sugary foods, such as maize (corn).
So Currie and his colleagues sequenced small segments of DNA from bacteria and other organisms living in fungus gardens in three Panamanian leaf-cutting ant colonies. They then compared the DNA against databases to identify what species were living in the fungus gardens, and what genes they contained.
This 'metagenomics' approach found that there are many species of bacteria in the fungus gardens that are capable of breaking down cellulose. The team also detected the genetic signatures of fungal enzymes that can break down cellulose, which raises the question of why the fungi can't break down cellulose in the laboratory.
Currie suggests that the newfound bacterial and fungal enzymes might be efficient at digesting cellulose because they have evolved for centuries along with the ant-fungal symbiosis. This could mean that the fungus can only break down cellulose in its natural context, or that the enzymes Currie detected are brought into the colony from outside. "The idea is that the ants' long evolutionary history may help us in our own attempts to break down plant biomass," he says.
Other researchers call Currie's findings interesting, but say they wanted to see a more thorough analysis of the data. "It's interesting that he found these fungal enzymes in the gardens that he didn't expect [based on] what the fungus was capable of doing by itself," says John Taylor, a mycologist at the University of California, Berkeley.
Taylor says that Currie's continued scrutiny of the lives of ants provides insights into the web of interactions necessary for the survival of any single species. "I think the coolest thing about this is that you start with one organism, and then you find more and more organisms involved in the relationship," he says. It may take a village to raise a child; it seems it also takes a village to break down cellulose.
Ant colonies could be key to advances in biofuels and antibiotics.
-Erika Check Hayden
Studies of bacteria on leaf-cutting ants could yield new antibiotics.Studies of bacteria on leaf-cutting ants could yield new antibiotics.M. MOFFETT/FLPA
In a mutually beneficial symbiosis, leaf-cutting ants cultivate fungus gardens, providing both a safe home for the fungi and a food source for the ants. But this 50-million-year-old relationship also includes microbes that new research shows could help speed the quest to develop better antibiotics and biofuels.
Ten years ago, Cameron Currie, a microbial ecologist then at the University of Toronto in Ontario, Canada, discovered that leaf-cutting ants carry colonies of actinomycete bacteria on their bodies (C. R. Currie et al. Nature 398, 701–704; 1999). The bacteria churn out an antibiotic that protects the ants' fungal crops from associated parasitic fungi (such as Escovopsis). On 29 March, Currie, Jon Clardy at the Harvard Medical School in Boston and their colleagues reported that they had isolated and purified one of these antifungals, which they named dentigerumycin, and that it is a chemical that has never been previously reported (D.-C. Oh et al. Nature Chem. Bio. doi: 10.1038/nchembio.159; 2009). The antifungal slowed the growth of a drug-resistant strain of the fungus Candida albicans, which causes yeast infections in people.
“These ants are walking pharmaceutical factories.”
Because distinct ant species cultivate different fungal crops, which in turn fall prey to specialized parasites, researchers hope that they will learn how to make better antibiotics by studying how the bacteria have adapted to fight the parasite in an ancient evolutionary arms race. "These ants are walking pharmaceutical factories," says Currie, now at the University of Wisconsin, Madison.
That's not the end to the possible applications. The ant colonies are also miniature biofuel reactors, Currie reported on 25 March at the Genomics of Energy & Environment meeting at the Joint Genome Institute in Walnut Creek, California. Each year, ants from a single colony harvest up to 400 kilograms of leaves to feed their fungal partners. But no one has worked out how the fungi digest the leaves, because samples of fungus grown in petri dishes can't break down cellulose, a tough molecule found in plant cells. Researchers are keenly interested in better ways to break down cellulose, because it might allow them to make more efficient biofuels than those made from sugary foods, such as maize (corn).
So Currie and his colleagues sequenced small segments of DNA from bacteria and other organisms living in fungus gardens in three Panamanian leaf-cutting ant colonies. They then compared the DNA against databases to identify what species were living in the fungus gardens, and what genes they contained.
This 'metagenomics' approach found that there are many species of bacteria in the fungus gardens that are capable of breaking down cellulose. The team also detected the genetic signatures of fungal enzymes that can break down cellulose, which raises the question of why the fungi can't break down cellulose in the laboratory.
Currie suggests that the newfound bacterial and fungal enzymes might be efficient at digesting cellulose because they have evolved for centuries along with the ant-fungal symbiosis. This could mean that the fungus can only break down cellulose in its natural context, or that the enzymes Currie detected are brought into the colony from outside. "The idea is that the ants' long evolutionary history may help us in our own attempts to break down plant biomass," he says.
Other researchers call Currie's findings interesting, but say they wanted to see a more thorough analysis of the data. "It's interesting that he found these fungal enzymes in the gardens that he didn't expect [based on] what the fungus was capable of doing by itself," says John Taylor, a mycologist at the University of California, Berkeley.
Taylor says that Currie's continued scrutiny of the lives of ants provides insights into the web of interactions necessary for the survival of any single species. "I think the coolest thing about this is that you start with one organism, and then you find more and more organisms involved in the relationship," he says. It may take a village to raise a child; it seems it also takes a village to break down cellulose.
Thursday, March 26, 2009
Using galaxy server
I recently came across galaxy server for genome sequence manipulation. While I do large scale computing, I myself have thought several times to create a web based tool for large scale sequence manipulation, that will be able to handle huge data files. My requirements would often be to subtract sequences from one file to another, fetch sequences from huge fasta file containing names from a list file, find duplicate triplicate sequences, find sequences of a certain length etc. I think the galaxy server provides many of that through a web server. The link is here
creating a wiki page for community portal
I have been toying with the idea of creating a wiki for our VMD web portal. The rationale is to create a user document in a wiki, which will also be used for bug tracking versioning etc. I am now looking at freely available tools/softwares for wiki creation( a compiled version is available here ).
For me documentation in wiki could be a great idea, since users can edit this page at ease and complement each others observation, instead of one person writing a whole document. Docuwiki seem to be the simplest one without the intervention of database. It just handles files and so I would assume is much safer option for our servers.
For me documentation in wiki could be a great idea, since users can edit this page at ease and complement each others observation, instead of one person writing a whole document. Docuwiki seem to be the simplest one without the intervention of database. It just handles files and so I would assume is much safer option for our servers.
Monday, March 09, 2009
Find pattern using grep (egrep)
One of the quickest and easiest ways to find patterns in a huge file without writing any scripts is through unix grep. I find it extremely useful to fetch strings having a certain pattern.
Looking for Tabs:
Many manuals don't say a lot about finding a tab character in a pattern. If you have a pattern "gene TAB exon", then the way to look for the file would be:
grep 'gene[[:blank:]]exon'
Looking for Tabs:
Many manuals don't say a lot about finding a tab character in a pattern. If you have a pattern "gene TAB exon", then the way to look for the file would be:
grep 'gene[[:blank:]]exon'
Monday, March 02, 2009
Implementing binary search in searching a list of sequences from a complex file
In sequence analysis work, it becomes imperative to search a list of sequences from a sequence file or a file having complex sets of data. One such complex dataset I can think of is an ACE file. An ACE file has multiple components:
Suppose you have a list of sequences in a file and you got to search which sequence belongs to which Contig from the ACE file. Then probably, you have to read the ACE file first as an hash of arrays, something like this:
After reading ACE into hash of arrays, need to iterate over the hash and get the sequence of arrays searched against each entry in the file. If the sequence file size is hugh, its good to implement a binary search(O(log2N)) instead of a serial search(O(N)). But knowing the fact that binary search can be implemented on a sorted array or list, it is more probable to sort the list file instead of sorting the array of hash. Because once sorted the list file stays sorted and does not need any sorting.
Implementation:
my @list = sort(@list); # sorting the sequence names in list file
AS contigs reads
CO contig_name bases reads segments compl (CAP3: segments=0)
sequence
BQ base_qualities
AF read1 compl padded_start_consensus (negatives meaning?)
AF read2 ..
BS segments
RD read1 bases info_items info_tags (latter two set to 0 by CAP3)
sequence
QA read1 qual_start qual_end align_start align_end
DS (phred header? left empty by CAP3)
RD read2 ...
Suppose you have a list of sequences in a file and you got to search which sequence belongs to which Contig from the ACE file. Then probably, you have to read the ACE file first as an hash of arrays, something like this:
while(){
if(/^CO Contig/){
$key = $_;
$flag = 1;
}
if($flag && /^[ATGC][ATGC][ATGC]/){
chomp;
$seq.=$_;
}
if(/^BQ/){
push(@arr,$seq);
$seq = '';
}
if(/^AF\s+(\S+)/){
push(@arr,$1);
}
if(/^BS/ && $flag){
$HOA{$key} = [ @arr ];
$flag = 0;
@arr = ();
}
}
After reading ACE into hash of arrays, need to iterate over the hash and get the sequence of arrays searched against each entry in the file. If the sequence file size is hugh, its good to implement a binary search(O(log2N)) instead of a serial search(O(N)). But knowing the fact that binary search can be implemented on a sorted array or list, it is more probable to sort the list file instead of sorting the array of hash. Because once sorted the list file stays sorted and does not need any sorting.
Implementation:
my @list = sort(@list); # sorting the sequence names in list file
# The following implementation is a linear search and takes
# much longer time
foreach my $key(keys %HOA){
my @arr = @{ $HOA{$key}} ;
my @found;
for(my $i=1;$i<=$#arr;$i++){
for(my $j=0;$j<=$#list;$j++){
chomp($list[$j]);
if($arr[$i] =~ $list[$j]){
push(@found, $arr[$i]);
splice(@list,$j,1); # removing the
# element from the list
last;
}
}
}
Binary search:
for(my $i=1;$i<=$#arr;$i++){
push(@found, bsearch($arr[$i],\@list);
}
sub bsearch {
my ($x, $a) = @_; # search for x in array a
my ($l, $u) = (0, @$a - 1); # lower, upper end of search interval
my $i; # index of probe
while ($l <= $u) {
$i = int(($l + $u)/2);
print($i, "\n");
if ($a->[$i] < $x) {
$l = $i+1;
}
elsif ($a->[$i] > $x) {
$u = $i-1;
}
else {
return $i; # found
}
}
return -1; # not found
}
Thursday, February 19, 2009
Cool perl tips
While dealing with very huge files, it is often necessary not to read the entire file into an array. This will be a huge RAM cost to the machine. Iterating through a file handle is often a easy solution. But often we need to go back and forth the file as we find a match. The easiest way to achieve this will be to use perl 'seek' function.
Similarly, Redo operator can also be explored to move back in loop.
Searching patterns in 2 data files:
In bioinformatics work, it is often necessary to search a list of names from a file from another file, that either has a list of sequences or blast output or anything else. So, the trick to reduce the search time will be keep the larger file as the outer loop(mostly as a while loop), the inside loop could be a while loop or a for loop, depending on the file size. If there is a for loop for the inner loop, then it is possible to get rid of the element that has just had a match to the outer element. This way with every iteration, the inner loop gets shorter and the time spent in searching gets shorter. In perl it is quite easy to remove an array element by simply using a delete function.
delete($array[index]) will do the job. But just delete removes the element while leaving $array[index] space empty. In order to remove the element along with index use splice function instead.
splice(@array,index,offset) --> will remove offset number of elements and indices starting from index. i.e; splice(@array,2,7) will remove 7 indices starting from 2.
The whole construct may look something like this:
If you have a subroutine and you need to pass more than one array, then don't pass by value. If passed by value then the list is flattened and the first and second element of the array is passed as 2 arrays. Example: @arr1 = [1,2,3];
@arr2 = [6,7,8];
func(@arr1,@arr2);
will print: Array1 is 1 and array2 is 2
sub func
{
my @a = shift;
my @b = shift;
print "Array1 is @a and array2 is @b\n";
}
In order to get the values intact pass by reference:
func(\@arr1,\@arr2);
in
sub func{
$a = shift;
$b = shift;
@arr1 = @$a; # de-referencing array
@arr2 = @$b;
Here is a time and space saver if you want to sort a hash by value and iterate over the sorted list:
# sorting tmp hash by value
foreach my $sorted(sort {$tmp{$b} <=> $tmp{$a}} keys %tmp){
my %geneHash;
for my $gene(keys %{$HOH1{$sorted}}){
my @coord = sort{$a <=> $b} @{ $HOH1{$sorted}{$gene}{'pos'}};
my $strand1 = $HOH1{$sorted}{$gene}{'strand'};
$geneHash{$gene} = $coord[0];
}
foreach my $sortedgene(sort {$geneHash{$a} <=> $geneHash{$b}} keys %geneHash ){
print FH "$sortedgene\t$serielNumber\n";
print FH1 ">
$serielNumber\n";
print FH1 "$transcript{$sortedgene}\n";
`sed 's/\t$sortedgene/$serielNumber/' output.gff`;
$serielNumber++;
}
}
# One variation if you want to sort the sequence file in descending order of their length. step -1:
Read fasta file into a hash called as %hash
# sort hash by the length of the value
foreach my $scaffold(sort {length($hash{$b}) <=> length($hash{$a}) } keys %hash) {
print ">scaffold_$num\n$hash{$scaffold}\n";
print FH "$scaffold\tscaffold_$num\t".length($hash{$scaffold})."\n"; $num++;
}
To make it to the ascending order simply change this to:
foreach my $scaffold(sort {length($hash{$a}) <=> length($hash{$b}) } keys %hash)
while(){ some pattern.. seek(FH, - length($_), 1); # will point the File handle # one line back }
NOTE: Just reading one line of a file at a time can be done by simply;
$_= ; # where FH stands for file handle.
Next call to $_= ; will read the next line. This is a great way to read
2 files together especially if they both are representing sequence and qual files
and the number of lines and values are in same order in both the files. Simply iterate
one file with a while loop and call $_= inside the while loop for the second
file.
Similarly, Redo operator can also be explored to move back in loop.
Searching patterns in 2 data files:
In bioinformatics work, it is often necessary to search a list of names from a file from another file, that either has a list of sequences or blast output or anything else. So, the trick to reduce the search time will be keep the larger file as the outer loop(mostly as a while loop), the inside loop could be a while loop or a for loop, depending on the file size. If there is a for loop for the inner loop, then it is possible to get rid of the element that has just had a match to the outer element. This way with every iteration, the inner loop gets shorter and the time spent in searching gets shorter. In perl it is quite easy to remove an array element by simply using a delete function.
delete($array[index]) will do the job. But just delete removes the element while leaving $array[index] space empty. In order to remove the element along with index use splice function instead.
splice(@array,index,offset) --> will remove offset number of elements and indices starting from index. i.e; splice(@array,2,7) will remove 7 indices starting from 2.
The whole construct may look something like this:
while(Passing Array by Reference){ if(/Query= (\S+)/){ $flag = 0; for(my $i; $i<=$#arr;$i++){ my @tmp = split(/\t/,$arr[$i]); if($1 eq $tmp[0]){ print "$_"; $flag = 1; delete($arr[$i]); # or splice(@arr,$i,1); last; } } next; } if($flag){ print "$_"; } }
If you have a subroutine and you need to pass more than one array, then don't pass by value. If passed by value then the list is flattened and the first and second element of the array is passed as 2 arrays. Example: @arr1 = [1,2,3];
@arr2 = [6,7,8];
func(@arr1,@arr2);
will print: Array1 is 1 and array2 is 2
sub func
{
my @a = shift;
my @b = shift;
print "Array1 is @a and array2 is @b\n";
}
In order to get the values intact pass by reference:
func(\@arr1,\@arr2);
in
sub func{
$a = shift;
$b = shift;
@arr1 = @$a; # de-referencing array
@arr2 = @$b;
Here is a time and space saver if you want to sort a hash by value and iterate over the sorted list:
# sorting tmp hash by value
foreach my $sorted(sort {$tmp{$b} <=> $tmp{$a}} keys %tmp){
my %geneHash;
for my $gene(keys %{$HOH1{$sorted}}){
my @coord = sort{$a <=> $b} @{ $HOH1{$sorted}{$gene}{'pos'}};
my $strand1 = $HOH1{$sorted}{$gene}{'strand'};
$geneHash{$gene} = $coord[0];
}
foreach my $sortedgene(sort {$geneHash{$a} <=> $geneHash{$b}} keys %geneHash ){
print FH "$sortedgene\t$serielNumber\n";
print FH1 ">
$serielNumber\n";
print FH1 "$transcript{$sortedgene}\n";
`sed 's/\t$sortedgene/$serielNumber/' output.gff`;
$serielNumber++;
}
}
# One variation if you want to sort the sequence file in descending order of their length. step -1:
Read fasta file into a hash called as %hash
# sort hash by the length of the value
foreach my $scaffold(sort {length($hash{$b}) <=> length($hash{$a}) } keys %hash) {
print ">scaffold_$num\n$hash{$scaffold}\n";
print FH "$scaffold\tscaffold_$num\t".length($hash{$scaffold})."\n"; $num++;
}
To make it to the ascending order simply change this to:
foreach my $scaffold(sort {length($hash{$a}) <=> length($hash{$b}) } keys %hash)
Thursday, February 12, 2009
Finding overlap
In sequence analysis, it is a common practice to find overlaps between 2 sets of features. This may be a coding sequence with repeats or coding sequences from 2 different prediction programs or 2 different coding sequences and so on. If data is present in gff files and are read into hashes with gene being the key holding the positions in an array, then the comparison can be made the following ways:
Lets find overlap between a gene and a repeat element.
geneLeft, repeatLeft be the first element of gene and repeat array respectively and geneRight, repeatRight be the last element of the gene and repeat array respectively.
#if repeat elements are occurring towards the left side of the gene element:
# and overlaps at the left hand side
repeatLeft <= geneLeft and repeatRight >= geneLeft # Note here Not repeatRight <= geneRight -> This condition will be true even if the repeat is occuring way upstream of the gene element and even if they don't overlap
# If repeat occurs at the right han side of the gene and overlaps at the right side
repeatLeft <= geneRight and repeatRight >= geneRight
# If repeat element lies within the gene element
geneLeft <= repeatLeft and repeatRight <= geneRight
#If gene element occurs inside the repeat element
repeatLeft <= geneLeft and geneRight <= repeatRight
If any of these conditions become true then it will find overlap(Join them with 'or')
if( ($coord1[$i] <= $coord[0] && $coord1[$i+1] >= $coord[-1]) ||
($coord1[$i] >= $coord[0] && $coord1[$i+1] <= $coord[-1]) ||
($coord1[$i] <= $coord[0] && $coord1[$i+1] >= $coord[0]) ||
($coord1[$i] >= $coord[0] && $coord1[$i] <= $coord[-1] )){
Tuesday, February 10, 2009
Using multidimensional hashes in storing sequence objects
One of the most important functions of multidimensional perl hashes has been their usage in reading sequence related data. The most common data type one reads from a genome sequencing prospective is a gff file. A gff file has the following attributes that are most likely captured in a data structure:
1. Name of the scaffold or chromosome on which the gene resides.
2. Name of the gene
3. Gene positions
4. Strand
Other attributes are less important, so lets ignore them at the moment.
So, now a scaffold will have multiple genes and a gene will have multiple locations. The ideal data structure that can store all these information will be a has of hashes and one of the hash holds arrays. In other words:
scaffold_1 => {
gene1 => { pos => [positions], strand => + or minus},
gene2 => { pos => [positions], strand => + or minus},
.......
},
scaffold_2 => {
gene1 => {pos => [positions], strand => + or minus},
gene2 => {pos => [positions], strand => + or minus},
.......
},
.....
Populating this hash from gff file:
sub test_gff{
my $fileName = shift;
my %HOH;
open PRED, $fileName or die "Can't open file $!\n";
while(){
chomp;
my @line = split(/\t/,$_);
push(@{$HOH{$line[0]}->{$line[8]}->{'pos'}}, $line[3], $line[4]);
$HOH{$line[0]}->{line[8]}->{'strand'}, $line[5];
}
close PRED;
return %HOH;
}
Now Printing the values:
my %HOH1 = &utility::test_gff($ARGV[0]); # Read the file
for my $scaffold(keys %HOH1){
for my $gene(keys %{$HOH1{$scaffold}}){
my @coord = @{ $HOH1{$scaffold}{$gene}{'pos'}};
my $strand = $HOH1{$scaffold}{$gene}{'strand'};
print "$scaffold -> $gene -> @coord -> $strand\n";
}
}
# If you want to sort the co-ordinates,Then that can be done as below:
my @coord = sort{$a <=> $b} @{ $HOH1{$scaffold}{$gene}{'pos'}};
# Values will be printed as :
Contig632 -> 713969 -> 1687 3391 -> +
Contig1122 -> 707832 -> 3109 3183 3194 3200 3200 3232 3232 3251 -> +
Contig1294 -> 714344 -> 3744 4473 -> +
Contig1294 -> 714343 -> 24 1652 -> -
Contig772 -> 707541 -> 14 64 126 147 147 182 -> +
Contig772 -> 713832 -> 4 81 126 147 147 194 -> +
Contig772 -> 711369 -> 1 64 130 147 147 204 -> +
Cool...
Adding one more dimension:
Mostly in draft genome sequences the genome assembly is not in great shape, so most often same genes can be found in multiple locations(Not necessarily through gene duplication). Sometimes the same gene could be present in the same scaffold in more than one location. In that case to stop over writing or appending the positions for the same gene, another check can be implemented, that is the gene serial number. So, you can just add position as another attribute with seriel number being the key.
push(@{$HOH{$line[0]}->{$line[8]}->{$serialNum}->{'pos'}}, $line[3], $line[4]);
$HOH{$line[0]}->{line[8]}->{$serialNum}->{'strand'}, $line[5];
And while iterating add one more for loop:
for my $scaffold(keys %HOH1){
for my $gene(keys %{$HOH1{$scaffold}}){
for my $serialNum(keys %{%HOH1{$scaffold}{$gene}}){
my @coord = @{ $HOH1{$scaffold}{$gene}{$serialNum}{'pos'}};
my $strand = $HOH1{$scaffold}{$gene}{$serialNum}{'strand'};
print "$scaffold -> $gene -> @coord -> $strand\n";
}
}
The ease at which we can use perl data structures make it so desirable. The coder does not need to worry about memory allocation stuff. Earlier I used to use linked list in C to store sequence data. Then I will be stuck in debugging stuff a lot. Now after I discovered this, I have almost stopped coding in C. Happy perl coding..
Thursday, January 29, 2009
Editing publication quality images with GIMP
GIMP stand for GNU Image Manipulating Program. I have this software in my computer on and off for the last few years. I was adamant to not budge for paid softwares like Adobe photo shop and like. Gimp used to sit there in my computer and annoy me with its small icon, so I used to uninstall it often. When I needed to manipulate some figures for publication, I would install it back on. Not until very recently I had no clue how powerful and useful gimp can be. It can even beat the best image editing softwares in market..
Challenges:
Recently, I faced a challenge merging 3 pathway maps. They can be as complicated as a circuit diagram!! I had to basically merge Image 1 and Image 2 and make a composite image12. Then I had to merge this with image 3 in such a way that, imae3 - image12 will be in a different color and image12 - image3 will be in a different color.
Open image 1 as layer(file --> open as layer), then open image 2 as layer as well. Then go to layer menu --> Merge down --> save image as image12. In image12, you will see the extra portions of image 1 as little lighter objects, but image 2 appears as it is. Now open image12 as layer then open image 3 as layer. Then go to layer menu --> merge down --> save as image123.
Now in image123 you see the extra portions in 12 that is not present in 3 as lighter objects. Well enough..
Now do the reverse. Open image 3 as layer first followed by image 12 as layer. Then merge as before and name it as image 312.
What you will see is in image 312 --> All that is present extra in 3 are not present in image 12 will appear as solid bordered transparent object. Distinct from the upper layer. In Image 123 you will see just the reverse. The ones that are present extra in image 12 and are not present in image 3. Now after all these, you can happily open image 312 and start playing around with colors. Open image 312 and as shown in figure select the shape from the toolbox menu. Select the area neatly that need to be colored differently. Then choose color. Then select color menu and just click on the selected area in the image. Viola.. that changes the color so well....
Friday, January 16, 2009
The power of vim
vim is a very powerful text editor having many great features that can prove to be a time saver if you are editing a large file.
I use the substitution operation in vim more often than not. One can a do a number of things with just very little practice.
1. Make a newline:
%s/$/\r/g # will make the end of the line into a new line character.
2. Replace annoying ^M characters that often happens from a DOS to unix transition
%s/ctrlVM/BY anything you like/ # Remember here ctrl-V makes a ^ and M makes an M
3. If you want to substitute something from line 20,30 here is how you can do
:20,30s/old/new/ -> remember here no '%' before 's'
4. If you want to change the order of two words from Hypothetical protein -> protein Hypothetical
:%s/\(Hypothetical\) \(protein\)/\2 \1/ -> This will change all occurrences of Hypothetical protein to protein Hypothetical.
Some More Tricks in vim
Delete lines that contain pattern:
:g/pattern/d
Delete all empty lines:
:g/^$/d
Delete lines in range that contain pattern:
:20,30/pattern/d
or with marks a and b:
:'a,'b/pattern/d
Substitute all lines for first occurance of pattern:
:%s/pattern/new/
:1,$s/pattern/new/
Substitute all lines for pattern globally (more than once on the line):
:%s/pattern/new/g
:1,$s/pattern/new/g
Find all lines containing pattern and then append -new to the end of each line:
:%s/\(.*pattern.*\)/\1-new/g
Substitute range:
:20,30s/pattern/new/g
with marks a and b:
:'a,'bs/pattern/new/g
Swap two patterns on a line:
:s/\(pattern1\)\(pattern2\)/\2\1/
Capitalize the first lowercase character on a line:
:%s/\([a-z]\)/\u\1/
more concisely:
:%s/[a-z]/\u&/
Capitalize all lowercase characters on a line:
:%s/\([a-z]\)/\u\1/g
more concisely:
:s/[a-z]/\u&/g
Capitalize all characters on a line:
:s/\(.*\)/\U\1\E/
Example: In a sequence file, if you have multiple headers and you are particularly
interested in substituting line Hp_ENSC_10a23 -> Hp_ENSC_10A23 without touching any
other line this is for you..
:%s/ENSC_\(.*\)\([a-z]\)/ENSC_\1\u\2/g
Regular Expressions:
Regular expressions in unix is more or less similar like perl - ya perl borrowed regexp from Unix..
But using substitution, one need to take care of adding '\' (esc) character in order to match a value.
For example if a line with super_2[TAB]5997239[TAB]5997668 needs to change to super_2:5997239-5997668, then instead of using '\S+' as the first non space patter use \(\S\+\). So, in other words the syntax will be
:%s/^\(\S\+\)\t\(\S\+\)\t/\1:\2-/g
Capitalize the first character of all words on a line:
:s/\<[a-z]/\u&/g
Uncapitalize the first character of all words on a line:
:s/\<[A-Z]/\l&/g
Change case of character under cursor:
~
Change case of all characters on line:
g~~
Change case of remaining word from cursor:
g~w
Increment the number under the cursor:
Decrement the number under the cursor:
redraw:
Turn on line numbering:
:set nu
Turn it off:
:set nonu
Number lines (filter the file through a unix command and replace with output):
:%!cat -n
Sort lines:
:%!sort
Sort and uniq:
:%!sort -u
Read output of command into buffer:
:r !ls -l
Refresh file from version on disk:
:e!
Open a new window:
n
Open a new window with the same file (split):
s
Split window vertically:
v
Close current window:
c
:q
Make current window the only window:
o
Cycle to next window:
w
Move to window below current window:
j
Move to window above current window:
k
Move to window left of current window:
h
Move to window right of current window:
l
Set textwidth for automatic line-wrapping as you type:
:set textwidth=80
Turn on syntax highlighting
:syn on
Turn it off:
:syn off
Force the filetype for syntax highlighting:
:set filetype=python
:set filetype=c
:set filetype=php
Use lighter coloring scheme for a dark background:
:set background=dark
Htmlize a file using the current syntax highlighting:
:so $VIMRUNTIME/syntax/2html.vim
Or, htmlize from a command prompt:
in 2html.sh put:
#!/bin/sh
vim -n -c ':so $VIMRUNTIME/syntax/2html.vim' -c ':wqa' $1 > /dev/null 2> /dev/null
Now just run: shell> 2html.sh foo.py
Getting rid of tab while copy pasting from a vim
edited file:
By default when you try to copy paste a block from a file edited with vim, it automatically appends character in the beginning of each line. In order to get rid of that in your .vimrc file put a line saying "set paste". This will do the job.
File Recovery:
Many times while editing a file the connection to the server may be lost. This is what I face all the time. In that case you can recover all the changes you made earlier by saying vim -r .file_name.swp
Executing Shell Commands
You can execute a command with the shell, for example, you can compile you C code when editing it, :!gcc kernel.c, or :!gcc % to compile the current editing file. Executing other commands is the same.
Keyword Matching
If you are lazy as me, you will find this function wonderful :) In insert mode, you can type a few characters of a word, e.g., if there is a variable called phosphatidylinositolphosphate, you may just type pho then C-p. Vim will find out the last word in the file started with characters pho, if it is not you want, you can re-type C-p to match the further previous ones. Similarly, C-n can do that for finding the next matches. Therefore you do not worry about the mistyping of the long variables, or uncommon-used words.
Insert/Delete an indent
C-T/C-D insert/delete an indent in the current line, no matter what column the cursor is in. There are handy hot keys when you are editing language with ambiguous block boundary such as Python, which you have to delete the indent yourself when necessary.
Settings for .vimrc
Some useful setting and mapping for .vimrc are listed as follows:
set nocp " nocompatible
set sw=4 " shiftwidth
set et " expandtab
set wm=8 " wrapmargin
set bs=2 " backspace
set ru " ruler
set ic " ignorecase
set is " incsearch
set scs " smartcase: override the 'ic' when searching
" if search pattern contains uppercase char
set vb t_vb= " set visual bell and disable screen flash
set backup
NOTE
Sometimes find replace and delete lines with a particular pattern can very quickly done using sed. Yes, vim takes longer time to complete a request if the file is large. So, instead try sed.
sed '/^--/d' lane7.1 > tmp # Will delete all lines with pattern --
# from file lane7.1 and will print in tmp
# -- comes as a pattern delimiter from 'grep'
sed 's/^@/>/g' lane7.1 > tmp # Will replace all lines beginning with @ to >
# This is useful in converting fastq files to fasta
I use the substitution operation in vim more often than not. One can a do a number of things with just very little practice.
1. Make a newline:
%s/$/\r/g # will make the end of the line into a new line character.
2. Replace annoying ^M characters that often happens from a DOS to unix transition
%s/ctrlVM/BY anything you like/ # Remember here ctrl-V makes a ^ and M makes an M
3. If you want to substitute something from line 20,30 here is how you can do
:20,30s/old/new/ -> remember here no '%' before 's'
4. If you want to change the order of two words from Hypothetical protein -> protein Hypothetical
:%s/\(Hypothetical\) \(protein\)/\2 \1/ -> This will change all occurrences of Hypothetical protein to protein Hypothetical.
Some More Tricks in vim
Delete lines that contain pattern:
:g/pattern/d
Delete all empty lines:
:g/^$/d
Delete lines in range that contain pattern:
:20,30/pattern/d
or with marks a and b:
:'a,'b/pattern/d
Substitute all lines for first occurance of pattern:
:%s/pattern/new/
:1,$s/pattern/new/
Substitute all lines for pattern globally (more than once on the line):
:%s/pattern/new/g
:1,$s/pattern/new/g
Find all lines containing pattern and then append -new to the end of each line:
:%s/\(.*pattern.*\)/\1-new/g
Substitute range:
:20,30s/pattern/new/g
with marks a and b:
:'a,'bs/pattern/new/g
Swap two patterns on a line:
:s/\(pattern1\)\(pattern2\)/\2\1/
Capitalize the first lowercase character on a line:
:%s/\([a-z]\)/\u\1/
more concisely:
:%s/[a-z]/\u&/
Capitalize all lowercase characters on a line:
:%s/\([a-z]\)/\u\1/g
more concisely:
:s/[a-z]/\u&/g
Capitalize all characters on a line:
:s/\(.*\)/\U\1\E/
Example: In a sequence file, if you have multiple headers and you are particularly
interested in substituting line Hp_ENSC_10a23 -> Hp_ENSC_10A23 without touching any
other line this is for you..
:%s/ENSC_\(.*\)\([a-z]\)/ENSC_\1\u\2/g
Regular Expressions:
Regular expressions in unix is more or less similar like perl - ya perl borrowed regexp from Unix..
But using substitution, one need to take care of adding '\' (esc) character in order to match a value.
For example if a line with super_2[TAB]5997239[TAB]5997668 needs to change to super_2:5997239-5997668, then instead of using '\S+' as the first non space patter use \(\S\+\). So, in other words the syntax will be
:%s/^\(\S\+\)\t\(\S\+\)\t/\1:\2-/g
Capitalize the first character of all words on a line:
:s/\<[a-z]/\u&/g
Uncapitalize the first character of all words on a line:
:s/\<[A-Z]/\l&/g
Change case of character under cursor:
~
Change case of all characters on line:
g~~
Change case of remaining word from cursor:
g~w
Increment the number under the cursor:
Decrement the number under the cursor:
redraw:
Turn on line numbering:
:set nu
Turn it off:
:set nonu
Number lines (filter the file through a unix command and replace with output):
:%!cat -n
Sort lines:
:%!sort
Sort and uniq:
:%!sort -u
Read output of command into buffer:
:r !ls -l
Refresh file from version on disk:
:e!
Open a new window:
Open a new window with the same file (split):
Split window vertically:
Close current window:
:q
Make current window the only window:
Cycle to next window:
Move to window below current window:
Move to window above current window:
Move to window left of current window:
Move to window right of current window:
Set textwidth for automatic line-wrapping as you type:
:set textwidth=80
Turn on syntax highlighting
:syn on
Turn it off:
:syn off
Force the filetype for syntax highlighting:
:set filetype=python
:set filetype=c
:set filetype=php
Use lighter coloring scheme for a dark background:
:set background=dark
Htmlize a file using the current syntax highlighting:
:so $VIMRUNTIME/syntax/2html.vim
Or, htmlize from a command prompt:
in 2html.sh put:
#!/bin/sh
vim -n -c ':so $VIMRUNTIME/syntax/2html.vim' -c ':wqa' $1 > /dev/null 2> /dev/null
Now just run: shell> 2html.sh foo.py
Getting rid of tab while copy pasting from a vim
edited file:
By default when you try to copy paste a block from a file edited with vim, it automatically appends
File Recovery:
Many times while editing a file the connection to the server may be lost. This is what I face all the time. In that case you can recover all the changes you made earlier by saying vim -r .file_name.swp
Executing Shell Commands
You can execute a command with the shell, for example, you can compile you C code when editing it, :!gcc kernel.c, or :!gcc % to compile the current editing file. Executing other commands is the same.
Keyword Matching
If you are lazy as me, you will find this function wonderful :) In insert mode, you can type a few characters of a word, e.g., if there is a variable called phosphatidylinositolphosphate, you may just type pho then C-p. Vim will find out the last word in the file started with characters pho, if it is not you want, you can re-type C-p to match the further previous ones. Similarly, C-n can do that for finding the next matches. Therefore you do not worry about the mistyping of the long variables, or uncommon-used words.
Insert/Delete an indent
C-T/C-D insert/delete an indent in the current line, no matter what column the cursor is in. There are handy hot keys when you are editing language with ambiguous block boundary such as Python, which you have to delete the indent yourself when necessary.
Settings for .vimrc
Some useful setting and mapping for .vimrc are listed as follows:
set nocp " nocompatible
set sw=4 " shiftwidth
set et " expandtab
set wm=8 " wrapmargin
set bs=2 " backspace
set ru " ruler
set ic " ignorecase
set is " incsearch
set scs " smartcase: override the 'ic' when searching
" if search pattern contains uppercase char
set vb t_vb= " set visual bell and disable screen flash
set backup
NOTE
Sometimes find replace and delete lines with a particular pattern can very quickly done using sed. Yes, vim takes longer time to complete a request if the file is large. So, instead try sed.
sed '/^--/d' lane7.1 > tmp # Will delete all lines with pattern --
# from file lane7.1 and will print in tmp
# -- comes as a pattern delimiter from 'grep'
sed 's/^@/>/g' lane7.1 > tmp # Will replace all lines beginning with @ to >
# This is useful in converting fastq files to fasta
Wednesday, January 14, 2009
Plotting Histograms and boxplot with R
Plotting huge files in R can be a cake walk. Here is how you can proceed:
Suppose say you have a huge file with 2 million entries in[TAB] format.
Then go ahead and read the file in the following format:
>data <- read.table(file="fileName",head=TRUE,sep=",")
# read.table is a file reading function
# Note here, if you have no header then say head=FALSE and if you data is separated by
# \t then say sep="\t"
#You may like to see summary of your data by saying
> summary(data)
V1 V2
HPAA-aaa57a02.b1: 1 Min. :0.000000
HPAA-aaa57a02.g1: 1 1st Qu.:0.009985
HPAA-aaa57a03.b1: 1 Median :0.020565
HPAA-aaa57a03.g1: 1 Mean :0.021674
HPAA-aaa57a04.b1: 1 3rd Qu.:0.033898
HPAA-aaa57a04.g1: 1 Max. :0.065789
(Other) :1140845
# Pretty Cool..
If you have not written down the column names in the data file, don't worry R provides one for you.
>names(data)
>[1] "V1" "V2"
# So, V1 and V2 are your column names.
# To Plot a Histogram simply say
>hist(data$V2)
# If you did not select any output device by default it will plot as a Rplots.ps file
# on your working directory
# You can however copy this to a file of your interest
>dev.copy(postscript,"myfile.ps",height=7,width=7)
# Here Height and width are in inches
# Sometimes it does not allow you to do a dev.copy when your dev.list() is NULL.
# You can populate it by first plotting something or by doing a hist(data$V2)
# or just a plot(data$V2)
# Then do
>hist(data$V2)
# You can check how many devices are attached to your screen by doing a
> dev.list()
postscript postscript postscript
2 3 4
# You can switch off your current device by doing a
>dev.off()
postscript
2
> dev.list()
postscript postscript
2 3
# I Plotted my histogram using the following command:
5
> hist(data$V2,breaks=100,main="Trace File Blast hit Frequency plot, cutoff=30"
,xlab="Y=M/(R+L - 2O)",col="red",las=1,freq=TRUE)
> lines(density(data$V2))
This will save the file in the current file "myfile.ps". Every subsequent call to a plot function will keep plotting in different pages of the same file.
Alternatively a truehist() can be called to plot histograms. But before that do a
>require(MASS)
> truehist(data$V2,nbins=100,main="Trace File Blast hit Frequency plot, cutoff=
30",xlab="Y=M/(R+L - 2O)",col="blue",las=1,prob=TRUE)
Making a boxplot:
Plotting a boxplot is equally easier using R. Follow the following commands:
data <- read.table(file="bgcorrect.out",head=FALSE,sep="\t")
summary(data)
#boxplot(data$vals,main='normalized data after synthetic normalization', ylab='G
ene Expression values');
boxplot(log2(data))
help(boxplot)
SAVING OUTPUT AS TEXT FILE
#In case you are just interested to save the histogram output file into a text file
> h <- hist(data$V2,breaks=100,plot=FALSE) # Plot= FALSE will not plot, default true
> save(h,ascii=TRUE,file="data.txt") # Will save data in ascii format, default=FALSE
# If you want to read the content of h
>h
# You will get the following
# $breaks -> break points
# $count -> frequency
# $density -> plot density
# $intensities -> plot intensity
# $mids -> mid points of data
So, just printing a single variable will be
>h$mids # Will print the mid point..
Suppose say you have a huge file with 2 million entries in
Then go ahead and read the file in the following format:
>data <- read.table(file="fileName",head=TRUE,sep=",")
# read.table is a file reading function
# Note here, if you have no header then say head=FALSE and if you data is separated by
# \t then say sep="\t"
#You may like to see summary of your data by saying
> summary(data)
V1 V2
HPAA-aaa57a02.b1: 1 Min. :0.000000
HPAA-aaa57a02.g1: 1 1st Qu.:0.009985
HPAA-aaa57a03.b1: 1 Median :0.020565
HPAA-aaa57a03.g1: 1 Mean :0.021674
HPAA-aaa57a04.b1: 1 3rd Qu.:0.033898
HPAA-aaa57a04.g1: 1 Max. :0.065789
(Other) :1140845
# Pretty Cool..
If you have not written down the column names in the data file, don't worry R provides one for you.
>names(data)
>[1] "V1" "V2"
# So, V1 and V2 are your column names.
# To Plot a Histogram simply say
>hist(data$V2)
# If you did not select any output device by default it will plot as a Rplots.ps file
# on your working directory
# You can however copy this to a file of your interest
>dev.copy(postscript,"myfile.ps",height=7,width=7)
# Here Height and width are in inches
# Sometimes it does not allow you to do a dev.copy when your dev.list() is NULL.
# You can populate it by first plotting something or by doing a hist(data$V2)
# or just a plot(data$V2)
# Then do
>hist(data$V2)
# You can check how many devices are attached to your screen by doing a
> dev.list()
postscript postscript postscript
2 3 4
# You can switch off your current device by doing a
>dev.off()
postscript
2
> dev.list()
postscript postscript
2 3
# I Plotted my histogram using the following command:
5
> hist(data$V2,breaks=100,main="Trace File Blast hit Frequency plot, cutoff=30"
,xlab="Y=M/(R+L - 2O)",col="red",las=1,freq=TRUE)
> lines(density(data$V2))
This will save the file in the current file "myfile.ps". Every subsequent call to a plot function will keep plotting in different pages of the same file.
Alternatively a truehist() can be called to plot histograms. But before that do a
>require(MASS)
> truehist(data$V2,nbins=100,main="Trace File Blast hit Frequency plot, cutoff=
30",xlab="Y=M/(R+L - 2O)",col="blue",las=1,prob=TRUE)
summary(data)
#boxplot(data$vals,main='normalized data after synthetic normalization', ylab='G
ene Expression values');
boxplot(log2(data))
help(boxplot)
SAVING OUTPUT AS TEXT FILE
#In case you are just interested to save the histogram output file into a text file
> h <- hist(data$V2,breaks=100,plot=FALSE) # Plot= FALSE will not plot, default true
> save(h,ascii=TRUE,file="data.txt") # Will save data in ascii format, default=FALSE
# If you want to read the content of h
>h
# You will get the following
# $breaks -> break points
# $count -> frequency
# $density -> plot density
# $intensities -> plot intensity
# $mids -> mid points of data
So, just printing a single variable will be
>h$mids # Will print the mid point..
Friday, January 09, 2009
RepeatModeler Woes...
Continued from my last post..
Folks, if you are tempted to run a de novo repeat finder such as RepeatModeler, then this is for you.
I happened to get copies of all the dependencies along with RepeatModeler and got them installed in my linux machine. Here is what is not documented in any installation manual.
Installation Caveats:
1. If you installed RepeatMasker without getting a copy of RepBase, it will not complain and installation will be fine. But when you complete installation of RepeatModeler(If you could), then it will complain big time - "you did not install RepBase on repeatMasker", so go back and re-configure RepeatMasker all over again.
SO, LESSON NO 1 ALWAYS GET REPBASE BEFORE INSTALLING REPEATMASKER
2. TRFinder is named as trfLinux.. when you obtain it from the download site. The installation program looks for just "trf" program in the folder. So, RENAME the executable to trf before installation.
3. The biggest trouble comes with installation of RepeatModeler. It will run smooth till it asks you the path for wublast installation. WUBLAST is currently named as ABBLast and is no longer available. You only get the 1998 version from the web site. If you are lucky and you have a copy of wublast, then the installation complains about xdformat.
ERROR:
/home/sutripa/wublast/xdformat is too old. Must have one dated 3/16/2005 or new
er. Install a newer version of wublast and re-run configure.
In order to fix this, you can do little bit of hacking:
Go to the "configure" file in repeatmodeler package and comment the block
#unless ( $result =~ /option requires an argument -- m/ ) {
# die "$wuLocation/xdformat is too old. Must have one dated 3/16
/2005 "
# . "or newer. Install a newer version of wublast and re-run
"
# . "configure.\n";
#}
(This is located between line number 247 - 255)
4. While giving the path for RECON, say $RECON_PATH/bin. Otherwise configuration script will not find the executables.
Once all these are taken care of installation of repeatmodeler can be done.
RUNNING REPEATMODELER:
Again this is not any lesser challenge than installation.
Goto RepeatModeler file and on line 234 do the following:
Replace line
#my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $genomeDB`;
with
my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -n -r $genomeDB`;
AND ALSO Line number 500
Replace line
#my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $sampleFa
staFile`;
with
my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -n -r $sampleFastaF
ile`;
Use the README file for creating a database using the xdformat program from wublast.
/BuildXDFDatabase -name elephant elephant.fa
IMPORTANT:
[Before Running this, check your fasta file. Make the header of fasta as >string[SPACE]number format, where the number could ideally be the length of the sequence. Leave nothing else trailing on the header.]
If you fail to do this, while running wublast the program is going to complain
either of these Two. Either sequence not found OR Database is empty.
After doing all this, you may like to run the program as
RepeatModeler -database Ha.test >& Ha.out
Alas!!You thought the problem is finally solved and go check the run file. It still has error like:
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-3_1-40000) and possibly more! Ignoring all entries on this line.
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-1_1280001-1320000) and possibly more! Ignoring all entries on this line.
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-1_520001-560000) and possibly more! Ignoring all entries on this line.
This is beacuse there is a silly bug in line 1590.
Go check one of the directories RM_12849.TueJan131331092009/round-1
you will see the header is printed like this: >>seq-3_1-40000
Now go to line 1590 in RepeatModeler and Replace
#print FASTA ">$seqID" . "_$start-$end\n";
With
print FASTA "$seqID" . "_$start-$end\n";
Now you can run the program.
Even after doing all these, I still have no luck in running this script. The run file ended up giving me a "unix broken pipe" error which is way beyond my control to fix.
Argh...
Folks, if you are tempted to run a de novo repeat finder such as RepeatModeler, then this is for you.
I happened to get copies of all the dependencies along with RepeatModeler and got them installed in my linux machine. Here is what is not documented in any installation manual.
Installation Caveats:
1. If you installed RepeatMasker without getting a copy of RepBase, it will not complain and installation will be fine. But when you complete installation of RepeatModeler(If you could), then it will complain big time - "you did not install RepBase on repeatMasker", so go back and re-configure RepeatMasker all over again.
SO, LESSON NO 1 ALWAYS GET REPBASE BEFORE INSTALLING REPEATMASKER
2. TRFinder is named as trfLinux.. when you obtain it from the download site. The installation program looks for just "trf" program in the folder. So, RENAME the executable to trf before installation.
3. The biggest trouble comes with installation of RepeatModeler. It will run smooth till it asks you the path for wublast installation. WUBLAST is currently named as ABBLast and is no longer available. You only get the 1998 version from the web site. If you are lucky and you have a copy of wublast, then the installation complains about xdformat.
ERROR:
/home/sutripa/wublast/xdformat is too old. Must have one dated 3/16/2005 or new
er. Install a newer version of wublast and re-run configure.
In order to fix this, you can do little bit of hacking:
Go to the "configure" file in repeatmodeler package and comment the block
#unless ( $result =~ /option requires an argument -- m/ ) {
# die "$wuLocation/xdformat is too old. Must have one dated 3/16
/2005 "
# . "or newer. Install a newer version of wublast and re-run
"
# . "configure.\n";
#}
(This is located between line number 247 - 255)
4. While giving the path for RECON, say $RECON_PATH/bin. Otherwise configuration script will not find the executables.
Once all these are taken care of installation of repeatmodeler can be done.
RUNNING REPEATMODELER:
Again this is not any lesser challenge than installation.
Goto RepeatModeler file and on line 234 do the following:
Replace line
#my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $genomeDB`;
with
my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -n -r $genomeDB`;
AND ALSO Line number 500
Replace line
#my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $sampleFa
staFile`;
with
my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -n -r $sampleFastaF
ile`;
Use the README file for creating a database using the xdformat program from wublast.
IMPORTANT:
[Before Running this, check your fasta file. Make the header of fasta as >string[SPACE]number format, where the number could ideally be the length of the sequence. Leave nothing else trailing on the header.]
If you fail to do this, while running wublast the program is going to complain
either of these Two. Either sequence not found OR Database is empty.
After doing all this, you may like to run the program as
RepeatModeler -database Ha.test >& Ha.out
Alas!!You thought the problem is finally solved and go check the run file. It still has error like:
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-3_1-40000) and possibly more! Ignoring all entries on this line.
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-1_1280001-1320000) and possibly more! Ignoring all entries on this line.
FastaDB:_cleanIndexAndCompact - ERROR: Multiple fasta seqs appear on one line (>
seq-1_520001-560000) and possibly more! Ignoring all entries on this line.
This is beacuse there is a silly bug in line 1590.
Go check one of the directories RM_12849.TueJan131331092009/round-1
you will see the header is printed like this: >>seq-3_1-40000
Now go to line 1590 in RepeatModeler and Replace
#print FASTA ">$seqID" . "_$start-$end\n";
With
print FASTA "$seqID" . "_$start-$end\n";
Now you can run the program.
Even after doing all these, I still have no luck in running this script. The run file ended up giving me a "unix broken pipe" error which is way beyond my control to fix.
Argh...
RepeatMasker blues..
I read rave reviews on RepeatMasker for repeat finding, so decided to get hold of a copy and get it used. I used this software before and did not like it much, but after getting some leading genome research institutes using it for their newly sequenced genomes, I got tempted. I used to use TRFinder from Boston university for finding direct repeats earlier.
I went to the web page of RepeatMasker that directs me to an array of other dependencies to be installed first. RepeatMasker depends on a repeat library RepBase which can be obtained from http://www.girinst.org. You go there and have to get yourself registered. Your registration need to be accepted by the group and then they send you an email carrying your login id and password. This may take upto 2 working days. Only after this you can download the RepBase database.
However, the trouble with RepBase is, it is a predefined library and will not detect any new repeat that may be there in your genome. So, the best option is to run a "De Novo" repeat finder. From the same web site you can get a link to RepeatModeler that acts as a wrapper around some de novo repeat finders such as RepeatScout and RECON.
Fair enough, then you go to install RepeatModeler...
This web page indicates you have to have a number of dependencies:
-- Of course you need to have perl 5.8 or above installed first.
-- you will need the Trfinder that I was mentioning earlier.
-- then you will need to download and install repeatmasker. This will not install unless you have RepBase in place.
--For RepeatMasker, you will also need phred/phrap/cross_match package. This can be done by writing to the authors and getting a copy by email. This may take a few business days.
--After you are done with repeatmasker you will need RECON to be installed
--Then RepeatScout need to be downloaded and installed.
--You will need WUBLAST also to be installed. The recommended site from the repeatmodeler page points to WUBLAST site at http://blast.wustl.edu/. When you go there you get redirected to http://www.advbiocomp.com . In advbiocomp site, they say that wublast is no longer supported and you can get an older version of WUblast for free. Then you get the binaries, inflate and keep them in place.
Now the final part: Installation of RepeatModeler. You do a ./configure. It asks you perl path, repeatMasker path etc. When it comes to WuBlast path, it exits saying that you xdformat program is older than 2005. What the hell, you back and check in wublast site, that says the program is as old as 1998. There are no newer version. There you get stuck.
Out of frustartion, I went and looked at the configure file in RepeatModeler. In line 237 there is a statement:
else {
my $result = `$wuLocation/xdformat -m 2>&1`;
# Two responses are expected:
# Good Version:
# /blah/wublast/xdformat: option requires an argument -- m
# Old Version:
# /blah/wublast/xdformat: invalid option -- m
unless ( $result =~ /option requires an argument -- m/ ) {
die "$wuLocation/xdformat is too old. Must have one dated 3/16/
5 "
. "or newer. Install a newer version of wublast and re-run
. "configure.\n";
}
I went back and asked the authors of xdformat, what this -m option means, they said to ignore it. Now I can change this, if only I find what is the substitute of this. Because my search in repeatModeler directory says me -m2 option is used:
gus@tyler-lab:~/RepeatModeler$ grep --regex=-m -r .
./configure: my $result = `$wuLocation/xdformat -m 2>&1`;
./Refiner: $malign->serializeOUT( "$wrkDir/$consName-malign-$round.ser"
)
./RepeatClassifier:# wublastx the simple-masked consensi vs the transposable ele
ment
./RepeatModeler:my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $genomeD
B`;
./RepeatModeler: my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $s
ampleFastaFile`;
./RepeatModeler: # . " -mincount 3 -tandemdist 500 -output $fastaFile.lfreq
"
./RepeatModeler:"$workDir/sampleDB-$round.fa $workDir/consensi.fa -warnings -kap
-wordmask=dust wordmask=seg maskextra=10 -hspmax=20 V=0 B=250 Q=25 R=5 W=7 S=25
0 gapS2=250 S2=125 X=250 gapX=500 -matrix=comparison.matrix"
I am still waiting for a response...
I went to the web page of RepeatMasker that directs me to an array of other dependencies to be installed first. RepeatMasker depends on a repeat library RepBase which can be obtained from http://www.girinst.org. You go there and have to get yourself registered. Your registration need to be accepted by the group and then they send you an email carrying your login id and password. This may take upto 2 working days. Only after this you can download the RepBase database.
However, the trouble with RepBase is, it is a predefined library and will not detect any new repeat that may be there in your genome. So, the best option is to run a "De Novo" repeat finder. From the same web site you can get a link to RepeatModeler that acts as a wrapper around some de novo repeat finders such as RepeatScout and RECON.
Fair enough, then you go to install RepeatModeler...
This web page indicates you have to have a number of dependencies:
-- Of course you need to have perl 5.8 or above installed first.
-- you will need the Trfinder that I was mentioning earlier.
-- then you will need to download and install repeatmasker. This will not install unless you have RepBase in place.
--For RepeatMasker, you will also need phred/phrap/cross_match package. This can be done by writing to the authors and getting a copy by email. This may take a few business days.
--After you are done with repeatmasker you will need RECON to be installed
--Then RepeatScout need to be downloaded and installed.
--You will need WUBLAST also to be installed. The recommended site from the repeatmodeler page points to WUBLAST site at http://blast.wustl.edu/. When you go there you get redirected to http://www.advbiocomp.com . In advbiocomp site, they say that wublast is no longer supported and you can get an older version of WUblast for free. Then you get the binaries, inflate and keep them in place.
Now the final part: Installation of RepeatModeler. You do a ./configure. It asks you perl path, repeatMasker path etc. When it comes to WuBlast path, it exits saying that you xdformat program is older than 2005. What the hell, you back and check in wublast site, that says the program is as old as 1998. There are no newer version. There you get stuck.
Out of frustartion, I went and looked at the configure file in RepeatModeler. In line 237 there is a statement:
else {
my $result = `$wuLocation/xdformat -m 2>&1`;
# Two responses are expected:
# Good Version:
# /blah/wublast/xdformat: option requires an argument -- m
# Old Version:
# /blah/wublast/xdformat: invalid option -- m
unless ( $result =~ /option requires an argument -- m/ ) {
die "$wuLocation/xdformat is too old. Must have one dated 3/16/
5 "
. "or newer. Install a newer version of wublast and re-run
. "configure.\n";
}
I went back and asked the authors of xdformat, what this -m option means, they said to ignore it. Now I can change this, if only I find what is the substitute of this. Because my search in repeatModeler directory says me -m2 option is used:
gus@tyler-lab:~/RepeatModeler$ grep --regex=-m -r .
./configure: my $result = `$wuLocation/xdformat -m 2>&1`;
./Refiner: $malign->serializeOUT( "$wrkDir/$consName-malign-$round.ser"
)
./RepeatClassifier:# wublastx the simple-masked consensi vs the transposable ele
ment
./RepeatModeler:my $dbIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $genomeD
B`;
./RepeatModeler: my $sampleIndex = `$RepModelConfig::XDFORMAT_PRGM -m2 -n -r $s
ampleFastaFile`;
./RepeatModeler: # . " -mincount 3 -tandemdist 500 -output $fastaFile.lfreq
"
./RepeatModeler:"$workDir/sampleDB-$round.fa $workDir/consensi.fa -warnings -kap
-wordmask=dust wordmask=seg maskextra=10 -hspmax=20 V=0 B=250 Q=25 R=5 W=7 S=25
0 gapS2=250 S2=125 X=250 gapX=500 -matrix=comparison.matrix"
I am still waiting for a response...
Sunday, January 04, 2009
Author List
Recently there was a post by physioprof on the co-first author that appears second in the list. While I don't agree to with views that second listed first author to be completely ignored, but it does raise some question about authorship in papers.
While the younger lot generally go for quantity over quality and try to get themselves in any paper they find getting out of the group, many older un-productive ones that have been lying dormant for a long time get into papers through their experience. In many instances there are at least 30% of the people in a multi author paper get a free ride. There are also few others who generally get a ride through their influence or by having a good relationship with the PI. I know atleast a bunch of people who see their names in the paper by periodically checking their previous PI's website. They have no idea how and when the paper get published. Now this raises a serious question on genuine authorships in papers. In many journals there is a section to mention which author did what to the paper and I think that should be the norm
While the younger lot generally go for quantity over quality and try to get themselves in any paper they find getting out of the group, many older un-productive ones that have been lying dormant for a long time get into papers through their experience. In many instances there are at least 30% of the people in a multi author paper get a free ride. There are also few others who generally get a ride through their influence or by having a good relationship with the PI. I know atleast a bunch of people who see their names in the paper by periodically checking their previous PI's website. They have no idea how and when the paper get published. Now this raises a serious question on genuine authorships in papers. In many journals there is a section to mention which author did what to the paper and I think that should be the norm
Friday, January 02, 2009
Author List and Order
Recently there was a post by physioprof on the co-first author that appears second in the list. While I don't agree to his views that second listed first author to be completely ignored, but it does raise some question about authorship in papers.
While the younger lot generally go for quantity over quality and try to get themselves in any paper they find getting out of the group, many older un-productive ones that have been lying dormant for a long time
While the younger lot generally go for quantity over quality and try to get themselves in any paper they find getting out of the group, many older un-productive ones that have been lying dormant for a long time
Thursday, October 30, 2008
Protein of the month
Every month protein spotlight from Expasy highlights a protein.The archive can be found here.
This month is the month for 2 proteins that has ability to make blood group A, B, AB to universal donor blood groups. There are in fact 2 molecules one for blood group A and one for blood group B. These are 2 bacterial glycosidases – α-N-acetylgalactosaminidase and α-1,3-galactosidase A.
As the blood groups are identified by the presence of an antigen on its surface that makes it either A, B or AB. The O group does not have any known antigen. This particular antigen is determined by a carbohydrate structure on the end of oligosaccharide chains on glycoproteins and glycolipids lodged in the cell’s membrane. Blood group A is defined by a terminal α-1,3-linked N-acetylgalactosamine (GalNAc). Blood group B is defined by a terminal α-1,3-linked galactose (Gal). And blood group O (which should really be pronounced ‘zero’) lacks both of these monosaccharides but instead is defined by an α-1,2-linked fucose. If there was a way of wiping off the monosaccharide entities (GalNAc and Gal) on A cells and B cells respectively, then blood groups A and B could be converted into the universal blood group O. What scientists then needed to find was something very small which could snap off the monosaccharide tips: an enzyme of some sort.
α-N-acetylgalactosaminidase and α-1,3-galactosidase A – which have the power to cleave and discard the GalNAc and Gal entities thus converting a red blood cell into a cell with neither a recognisable A nor B blood group antigen. In other words, they can convert A and B blood groups into a universal O blood group. Alpha-N-acetylgalactosaminidase is known in greater molecular detail than its homolog α-1,3-galactosidase A. Surprisingly, unlike other glycosidases of the same family, it needs the help of a cofactor (NAD+), which is nestled in the depths of a narrow tunnel, to cleave the GalNAc monosaccharide. Another part of the enzyme forms a crater which is large enough to accommodate an antigen from which the monosaccharide tip can be cleaved.
This month is the month for 2 proteins that has ability to make blood group A, B, AB to universal donor blood groups. There are in fact 2 molecules one for blood group A and one for blood group B. These are 2 bacterial glycosidases – α-N-acetylgalactosaminidase and α-1,3-galactosidase A.
As the blood groups are identified by the presence of an antigen on its surface that makes it either A, B or AB. The O group does not have any known antigen. This particular antigen is determined by a carbohydrate structure on the end of oligosaccharide chains on glycoproteins and glycolipids lodged in the cell’s membrane. Blood group A is defined by a terminal α-1,3-linked N-acetylgalactosamine (GalNAc). Blood group B is defined by a terminal α-1,3-linked galactose (Gal). And blood group O (which should really be pronounced ‘zero’) lacks both of these monosaccharides but instead is defined by an α-1,2-linked fucose. If there was a way of wiping off the monosaccharide entities (GalNAc and Gal) on A cells and B cells respectively, then blood groups A and B could be converted into the universal blood group O. What scientists then needed to find was something very small which could snap off the monosaccharide tips: an enzyme of some sort.
α-N-acetylgalactosaminidase and α-1,3-galactosidase A – which have the power to cleave and discard the GalNAc and Gal entities thus converting a red blood cell into a cell with neither a recognisable A nor B blood group antigen. In other words, they can convert A and B blood groups into a universal O blood group. Alpha-N-acetylgalactosaminidase is known in greater molecular detail than its homolog α-1,3-galactosidase A. Surprisingly, unlike other glycosidases of the same family, it needs the help of a cofactor (NAD+), which is nestled in the depths of a narrow tunnel, to cleave the GalNAc monosaccharide. Another part of the enzyme forms a crater which is large enough to accommodate an antigen from which the monosaccharide tip can be cleaved.
Tuesday, October 28, 2008
When under attack, plants can signal microbial friends for help
A new research finding at University Delaware proves that plants are not just sitting ducks, they react in many different ways to an external attack. The full story can be found here
Wednesday, October 22, 2008
Substantial biases in ultra-short read data sets from high-throughput DNA sequencing
This paper appeared in September issue of Nucleic Acid Research(September; 36(16): e105.)
They generated and analyzed two Illumina 1G ultra-short read data sets, i.e. 2.8 million 27mer reads from a Beta vulgaris genomic clone and 12.3 million 36mers from the Helicobacter acinonychis genome. They found that error rates range from 0.3% at the beginning of reads to 3.8% at the end of reads. Wrong base calls are frequently preceded by base G. Base substitution error frequencies vary by 10- to 11-fold, with A > C transversion being among the most frequent and C > G transversions among the least frequent substitution errors.
Sequencing Stratergies
To achieve high throughput, the new approaches apply different strategies. 454 Life Sciences has adapted pyrosequencing to a microbead format to sequence 400 000 DNA fragments simultaneously, resulting in a per-run dataset of 100 Mbp with reads averaging 250 bp. SOLiD(Applied Biosystems’ Sequencing by Oligonucleotide Ligation and Detection) sequencing also uses templates immobilized onto microbeads. Here, the sequence of the template DNA is decoded by ligation assays involving oligonucleotides labeled with different fluorophores. The SOLiD read length is currently 25–35 bases, and 2–3 Gbp of data can be collected during an 8-day run. Solexa sequencing is based on amplifying single molecules attached to the surface of a flow cell to generate clusters of identical molecules, followed by sequencing using fluorophore-labeled reversible chain terminators. Solexa sequencing proceeds a base at a time and read length depends on the number of sequencing cycles. Current Illumina sequencing instrumentation achieves read lengths of 36 bases. The Solexa flow cell is composed of eight separately loadable lanes. Since each lane has a capacity of about 5 million reads, > 40 million reads can be generated in a run of 3 days, equivalent to > 1.3 Gbp.
This is quite a finding, since many of us are rapidly moving into high throughput sequencing...
They generated and analyzed two Illumina 1G ultra-short read data sets, i.e. 2.8 million 27mer reads from a Beta vulgaris genomic clone and 12.3 million 36mers from the Helicobacter acinonychis genome. They found that error rates range from 0.3% at the beginning of reads to 3.8% at the end of reads. Wrong base calls are frequently preceded by base G. Base substitution error frequencies vary by 10- to 11-fold, with A > C transversion being among the most frequent and C > G transversions among the least frequent substitution errors.
Sequencing Stratergies
To achieve high throughput, the new approaches apply different strategies. 454 Life Sciences has adapted pyrosequencing to a microbead format to sequence 400 000 DNA fragments simultaneously, resulting in a per-run dataset of 100 Mbp with reads averaging 250 bp. SOLiD(Applied Biosystems’ Sequencing by Oligonucleotide Ligation and Detection) sequencing also uses templates immobilized onto microbeads. Here, the sequence of the template DNA is decoded by ligation assays involving oligonucleotides labeled with different fluorophores. The SOLiD read length is currently 25–35 bases, and 2–3 Gbp of data can be collected during an 8-day run. Solexa sequencing is based on amplifying single molecules attached to the surface of a flow cell to generate clusters of identical molecules, followed by sequencing using fluorophore-labeled reversible chain terminators. Solexa sequencing proceeds a base at a time and read length depends on the number of sequencing cycles. Current Illumina sequencing instrumentation achieves read lengths of 36 bases. The Solexa flow cell is composed of eight separately loadable lanes. Since each lane has a capacity of about 5 million reads, > 40 million reads can be generated in a run of 3 days, equivalent to > 1.3 Gbp.
This is quite a finding, since many of us are rapidly moving into high throughput sequencing...
Wednesday, October 15, 2008
Working with a 64 bit machine
I thought I will write this in my note book, but instead I ended up putting it on my blog site! This problem has taken a lot of my working time to solve. Recently I have been moving our entire web portal from a 32 bit machine to a 64 bit machine. The data transfer was quite smooth except for some small troubles here and there - like with PHP versions. In older PHP versions, I did not need to write 'php' near the opening tag(
I have written a C program to calculate the Fickett score as well as the coding potential of a given sequence. Then plot the values on fly using perl GD graphics module. To my surprise, the same GD code plotted like Fig-1 in 32 bit machine where as in 64 bit machine the output was a straight line(Fig-2). I tried to compare the computed text output in 2 machines. Although they were not identical but they had quite good similarity. What the heck. Even I got the numerical output and plotted on 32 bit machine still got figures like the first one. Then on a quick look I found that the newer bioperl module(perl -MGD -e 'print $GD::VERSION'; 2.41) is unable to plot anything when there is a number like '-inf'. This is handled very well with older(2.35) GD packages. Apart from this I was also wondering how 'atof' function works in C. Earlier, the codon usage table was written where there was no '0' preceeding the decimal point (like 0.02344). The array that reads this value is just an array of size 5 and I used fscanf function to read this.My syntax was:
(fscanf(fp,"%s",store[i])!=EOF) . The file was something like this :
GCA .0200 GCC .019 GCG .0238 GCT .0255
TGC .0104 TGT .008 GAC .0356 GAT .022
GAA .0234 GAG .0427 TTC .017 TTT .0213
GGA .014 GGC .022 GGG .009 GGT .0178
No matter whatever I did, fscanf always mis-behaved when it is reading a float as a string (with syntax %s).So, instead of reading .0200 and GCC as 2 different variables it sometimes read .0200GCC as one variable and GCC as another. It practically did not ignore the space in between. So, when I converted the string .0200GCC to a float using atof, it produced a different result than just .0200. As a result of which I would get values like -inf when the log likelihood values are calculated, where it should have been some negative integer. Now having learn t this lesson, I would suggest anyone reading a mixed file of characters and float or any other data types to read them any other way then to read it as string. For this particular case something like
(fscanf(fp,"%s%f",tmp,&val) !=EOF) would be more than appropriate..
Subscribe to:
Posts (Atom)