Maybe I should simplify my question. I need to tweak the first file (#), make it look like #2, so that it is compatible with the tool that I am using to assign taxonomic ranks to my blast outputs.

Both of those formats look OK enough for BLAST to work, it seems like there might be a more subtle problem.  

Is it possible that the identifier in #1 is not unique, and there are repeats in the file?  This can happen if there is a space within the identifier, or paired-end sequences where the end information is lost, or only part of the identifier is being carried forward - does the pipeline place a limit on the length of identifiers?

One idea is to replace all of the header information with a simple incrementing identifier like:

>A1

...

>A2

...

You could try this for a few sequences manually first, to see whether this is OK for the pipeline.

@Jonathan David Moore, thanks a lot for your detailed answer! I have run megablast for 9 of my datasets, and they all seem to work fine. This is the megablast output I am getting for one of them:

M01409:7:000000000-A45GT:1:1101:18413:1651_1:N:0:2 gi|557880854:10-609 98.84 172 2 0 34 205 600 429 2e-82 307

Your suggestion is to tweak the subset of the data with this simple identifier (A1), and run the megablast, if I am not wrong? Or are you suggesting me to to run my taxonomical classification tool on the modified megablast outputs? I just wanted to make sure. I will surely try this.

The pipeline has been created by my colleague and I do not think it does place a limit on the length of identifiers. Space within the identifier could potentially cause a problem. However, there is no space in #1, that is the part I am getting confused..

You could try to adjust the megablast output, that might be quicker just to test whether uniqueness is the problem.

I did adjust it for a subset initially, only keeping A1 and A2 did not work, apparently taxonomic classification tool is using part of the information in the identifier.

Two separate single end files differ from each other with the last section of the identifier (end1 has _1:N:0:2, end2 has _2:N:0:2 ), and I figured out that it was causing the problem in the pipeline. Once I removed these parts from the two single end fasta files, it worked!!

The only issue right now is that I should remove this from more than 10 million reads (using replace option in gedit). I wish I would have known some scripts that could help me remove those through command line. I am hoping that it does not crash somewhere in the middle. Thanks a lot for the very useful suggestion.

It is strange that the pipeline you use is sensitive to the sequence identifier, but hey... if it works :) Here is a short PERL script that should remove the _1:N:0:2 and the _2:N:0:2 endings within the file "input.txt" (put the correct file name here).

To install perl and run this script you simply need to download and install active perl http://www.activestate.com/activeperl , put the script below into a file with a .pl extension, put the file in the the same folder as the input.txt and run the pl file (command line: perl perlscript.pl).

Please note that this scipt works as intended only if you have only one "_" symbol in your sequence identifier (which seems to be the case, if not: pls message me and I can rework this script for you). Also don't forget to backup your data just in case. 

#strart perlscript; 

use strict;

use warnings;

open INPUT, "input.txt";

open OUTPUT, ">output.txt";

while() {

chomp;

my $string = $_;

if($string =~ /\_/) {

my @s = split(/\_/, $string);

$string = $s[0];

}

print OUTPUT "$string\n";

}

close INPUT;

close OUTPUT;

#end perlscript; 

The script reads the input file line by line, finds lines with "_" symbols and returns them without the "_" and anything to the right from it. Lines without "_" symbols are returned unchanged. 

If you apply this script your input: 

>M01403:7:000000000-A45GT:1:1102:16645:1483_1:N:0:2 

TAATTGATCCGTTAA

Becomes output: 

>M01403:7:000000000-A45GT:1:1102:16645:1483

TAATTGATCCGTTAA

There are also many other ways to do this, but this one is the most simple to explain given the fact that I don't know what operating system you are used to and what tools you are familiar with :)

Hope this helps with your issue.

P.S. Also we recently published a tool (called TUIT) for automatic taxonomic identification of sequences that should not depend on the sequence identifiers or NGS platform http://www.ncbi.nlm.nih.gov/pubmed/24502797

If you allready have the megablast results and if they are xml formatted ( -5 option in Blast) you can run TUIT with "-b output" and it will use your megablast outputs.

Article TUIT, a BLAST-Based tool for taxonomic classification of nuc...

Thanks a lot @Alexander Panchin! It seems to be a very useful script. I did try it with my file, however what this script does is to remove everything after the underscore, hence this is the output I am getting (which you already mentioned):

HWI-ST513:311:C1JADACXX:7:1101:1162:2187

HWI-ST513:311:C1JADACXX:7:1101:1162:2187

....

Means that I am losing all the gi, percentage and length information (an example to megablast output: M01409:7:000000000-A45GT:1:1101:18413:1651_1:N:0:2 gi|557880854:10-609 98.84 172 2 0 34 205 600 429 2e-82 307) without which I cannot process the downstream analysis!

Any suggestions for this issue? I feel that it can be worked out and still function:)

Since I only needed to remove the very last part of the first identifier, i.e _1:N:0:2 and _2:N:0:2, I found a similar script, and just adapted for my file. It works!!!! Script goes like this:

my $replacestring = "_1:N:0:2";

open(my $infile,"

Oh yes. That solution is better. I assumed that your identifier looks like this:  >M01403:7:000000000-A45GT:1:1102:16645:1483_1:N:0:2, didn't notice that more detailed description in the comments. I'm glad that perl helped. It's actually very useful and simple. 

Here is a database/server that you may find useful for future use....

http://www.ezbiocloud.net/eztaxon

Bilgenur, the problematic part of your identifiers are usually indicating paired-end reads. So for MiSeq the naming scheme is, that the forward reads have an additional 1:N:0:2 (for example, the last number depends on the index sequence used) and the reverse reads an additional 2:N:0:2 appended.

For a more detailed description see here:

http://support.illumina.com/help/SequencingAnalysisWorkflow/Content/Vault/Informatics/Sequencing_Analysis/CASAVA/swSEQ_mCA_FASTQFiles.htm

Because both, the forward and reverse reads, are present in one single file (as described above), each identifiers would be doubled. This could be problematic for tools which cannot handle (interleaved) MiSeq paired-end data. Therefore, I assume, the identifiers were modified, as in your case, to the following (adding a "_" between the identifier and the read number and remove the gap/white space):

from

>M01403:7:000000000-A45GT:1:1102:16645:1483 1:N:0:2

to

>M01403:7:000000000-A45GT:1:1102:16645:1483_1:N:0:2

If you now remove the "_1:N:0:2" ore the like, then you will remove the uniqueness of the paired-end read identifiers, which (again) might be problematic for some tools.

If the identifiers are causing the problems in your pipeline (maybe the "_" is the problem?), i would modify the perl script and instead of removing this part change it, from e.g. "_1:N:0:2" to ":1:N:0:2" or to "#1:N:0:2". This way, the IDs remain unique and also the pairing information (is it a forward or reverse read?) stays intact.

- Sebastian

I am now trying to resolve another issue on extracting only the gi numbers. I have an excel file with >20,000 lines containing gi numbers and length information, such as the following:

gi|85726344:1-644

gi|106896103|gb|DQ525233.1

What I would like to do is to only get the gi numbers, for above examples:

gi|85726344

gi|106896103

However, when I use the script I mentioned earlier (see attached), I cannot get rid of the numbers, so for instance when I want to remove (:) I am still ending up with -644, for the first case. Does this script only work with underscores or colon as well?