Commit f6f7b098 authored by Andreas Tille's avatar Andreas Tille
Browse files

New upstream version 1.0.3

parent 64b022e1

MANUAL

+77 −1
Original line number Diff line number Diff line
@@ -198,7 +198,83 @@ can be generated from a GI taxid dump: 


### Custom database



TODO: Add toy example for nodes.dmp, names.dmp and seqid2taxid.map

To build a custom database, you need the provide the follwing four files to `centrifuge-build`:



  - `--conversion-table`: tab-separated file mapping sequence IDs to taxonomy IDs. Sequence IDs are the header up to the first space or second pipe (`|`).  

  - `--taxonomy-tree`: `\t|\t`-separated file mapping taxonomy IDs to their parents and rank, up to the root of the tree. When using NCBI taxonomy IDs, this will be the `nodes.dmp` from `ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz`.

  - `--name-table`: '\t|\t'-separated file mapping taxonomy IDs to a name. A further column (typically column 4) must specify `scientific name`. When using NCBI taxonomy IDs, `names.dmp` is the appropriate file.

  - reference sequences: The ID of the sequences are the header up to the first space or second pipe (`|`)



When using custom taxonomy IDs, use only positive integers greater-equal to `1` and use `1` for the root of the tree.



#### More info on `--taxonomy-tree` and `--name-table`



The format of these files are based on `nodes.dmp` and `names.dmp` from the NCBI taxonomy database dump. 



- Field terminator is `\t|\t`

- Row terminator is `\t|\n`



The `taxonomy-tree` / nodes.dmp file consists of taxonomy nodes. The description for each node includes the following

fields:



    tax_id                  -- node id in GenBank taxonomy database

    parent tax_id           -- parent node id in GenBank taxonomy database

    rank                    -- rank of this node (superkingdom, kingdom, ..., no rank)



Further fields are ignored.



The `name-table` / names.dmp is the taxonomy names file:



    tax_id                  -- the id of node associated with this name

    name_txt                -- name itself

    unique name             -- the unique variant of this name if name not unique

    name class              -- (scientific name, synonym, common name, ...)



`name class` **has** to be `scientific name` to be included in the build. All other lines are ignored



#### Example



*Conversion table `ex.conv`*: 

    

    Seq1	11

    Seq2	12

    Seq3	13

    Seq4	11



*Taxonomy tree `ex.tree`*: 



    1	|	1	|	root

    10	|	1	|	kingdom

    11	|	10	|	species

    12	|	10	|	species

    13	|	1	|	species



*Name table `ex.name`*:



    1	|	root	|		|	scientific name	|

    10	|	Bacteria	|		|	scientific name	|

    11	|	Bacterium A	|		|	scientific name	|

    12	|	Bacterium B	|		|	scientific name	|

    12	|	Some other species	|		|	scientific name	|



*Reference sequences `ex.fa`*:



    >Seq1

    AAAACGTACGA.....

    >Seq2

    AAAACGTACGA.....

    >Seq3

    AAAACGTACGA.....

    >Seq4

    AAAACGTACGA.....



To build the database, call



    centrifuge-build --conversion-table ex.conv \

                     --taxonomy-tree ex.tree --name-table ex.name \ 

                     ex.fa ex



which results in three index files named `ex.1.cf`, `ex.2.cf` and `ex.3.cf`.



### Centrifuge classification output

MANUAL.markdown

+78 −1
Original line number Diff line number Diff line
@@ -211,7 +211,84 @@ can be generated from a GI taxid dump: 


### Custom database



TODO: Add toy example for nodes.dmp, names.dmp and seqid2taxid.map

To build a custom database, you need the provide the follwing four files to `centrifuge-build`:



  - `--conversion-table`: tab-separated file mapping sequence IDs to taxonomy IDs. Sequence IDs are the header up to the first space or second pipe (`|`).  

  - `--taxonomy-tree`: `\t|\t`-separated file mapping taxonomy IDs to their parents and rank, up to the root of the tree. When using NCBI taxonomy IDs, this will be the `nodes.dmp` from `ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz`.

  - `--name-table`: '\t|\t'-separated file mapping taxonomy IDs to a name. A further column (typically column 4) must specify `scientific name`. When using NCBI taxonomy IDs, `names.dmp` is the appropriate file.

  - reference sequences: The ID of the sequences are the header up to the first space or second pipe (`|`)



When using custom taxonomy IDs, use only positive integers greater-equal to `1` and use `1` for the root of the tree.



#### More info on `--taxonomy-tree` and `--name-table`



The format of these files are based on `nodes.dmp` and `names.dmp` from the NCBI taxonomy database dump. 



- Field terminator is `\t|\t`

- Row terminator is `\t|\n`



The `taxonomy-tree` / nodes.dmp file consists of taxonomy nodes. The description for each node includes the following

fields:



    tax_id                  -- node id in GenBank taxonomy database

    parent tax_id           -- parent node id in GenBank taxonomy database

    rank                    -- rank of this node (superkingdom, kingdom, ..., no rank)



Further fields are ignored.



The `name-table` / names.dmp is the taxonomy names file:



    tax_id                  -- the id of node associated with this name

    name_txt                -- name itself

    unique name             -- the unique variant of this name if name not unique

    name class              -- (scientific name, synonym, common name, ...)



`name class` **has** to be `scientific name` to be included in the build. All other lines are ignored



#### Example



*Conversion table `ex.conv`*: 

    

    Seq1	11

    Seq2	12

    Seq3	13

    Seq4	11





*Taxonomy tree `ex.tree`*: 



    1	|	1	|	root

    10	|	1	|	kingdom

    11	|	10	|	species

    12	|	10	|	species

    13	|	1	|	species



*Name table `ex.name`*:



    1	|	root	|		|	scientific name	|

    10	|	Bacteria	|		|	scientific name	|

    11	|	Bacterium A	|		|	scientific name	|

    12	|	Bacterium B	|		|	scientific name	|

    12	|	Some other species	|		|	scientific name	|



*Reference sequences `ex.fa`*:



    >Seq1

    AAAACGTACGA.....

    >Seq2

    AAAACGTACGA.....

    >Seq3

    AAAACGTACGA.....

    >Seq4

    AAAACGTACGA.....



To build the database, call



    centrifuge-build --conversion-table ex.conv \

                     --taxonomy-tree ex.tree --name-table ex.name \ 

                     ex.fa ex



which results in three index files named `ex.1.cf`, `ex.2.cf` and `ex.3.cf`.





### Centrifuge classification output

Makefile

+1 −0
Original line number Diff line number Diff line
@@ -191,6 +191,7 @@ CENTRIFUGE_SCRIPT_LIST = centrifuge \ 
	centrifuge-build \

	centrifuge-inspect \

	centrifuge-download \

	centrifuge-kreport \

	$(wildcard centrifuge-*.pl)

aln_sink.h

+2 −2
Original line number Diff line number Diff line
@@ -2305,9 +2305,9 @@ void AlnSinkSam<index_t>::appendMate( 
			case READ_ID:      appendReadID(o, rd.name); break;

			case SEQ_ID:       appendSeqID(o, rs, ebwt.tree()); break;

			case SEQ:          o.append((string(rd.patFw.toZBuf()) + 

										(rdo == NULL? "" : "N" + string(rdo->patFw.toZBuf()))).c_str()); break;

										(rdo == NULL? "" : "_" + string(rdo->patFw.toZBuf()))).c_str()); break;

			case QUAL:         o.append((string(rd.qual.toZBuf()) + 

										(rdo == NULL? "" : "I" + string(rdo->qual.toZBuf()))).c_str()); break;

										(rdo == NULL? "" : "_" + string(rdo->qual.toZBuf()))).c_str()); break;



			case SEQ1:         o.append(rd.patFw.toZBuf()); break;

			case QUAL1:        o.append(rd.qual.toZBuf()); break;

centrifuge

+169 −22
Original line number Diff line number Diff line
@@ -116,6 +116,13 @@ my %read_compress = (); 
my $cap_out = undef;       # Filename for passthrough

my $no_unal = 0;

my $large_idx = 0;



# Variables handling the output format

my $outputFmtSam = 0 ;

my $tabFmtOptIdx = 0 ;

my $needReadSeq = 0 ;

my $removeSeqCols = 0 ;



# Remove whitespace

for my $i (0..$#bt2_args) {

	$bt2_args[$i]=~ s/^\s+//; $bt2_args[$i] =~ s/\s+$//;
@@ -179,6 +186,7 @@ for(my $i = 0; $i < scalar(@bt2_args); $i++) { 
	}

	for my $rarg ("un-conc", "al-conc", "un", "al") {

		if($arg =~ /^--${rarg}$/ || $arg =~ /^--${rarg}-gz$/ || $arg =~ /^--${rarg}-bz2$/) {

			$needReadSeq = 1 ;

			$bt2_args[$i] = undef;

			if(scalar(@args) > 1 && $args[1] ne "") {

				$read_fns{$rarg} = $args[1];
@@ -193,7 +201,57 @@ for(my $i = 0; $i < scalar(@bt2_args); $i++) { 
			last;

		}

	}

	if ($arg eq "--out-fmt" )

	{

		$i < scalar(@bt2_args)-1 || Fail("Argument expected in next token!\n");

		$i++;

		if ( $bt2_args[$i] eq "sam" ) 

		{

			$outputFmtSam = 1 ;

		}

		#$bt2_args[$i] = undef;



	}



	if ( $arg eq "--tab-fmt-cols" )

	{

		$i < scalar(@bt2_args)-1 || Fail("Argument expected in next token!\n");

		$tabFmtOptIdx = $i + 1 ;

	}

}



# Determine whether we need to add two extra columns for seq and qual to out-fmt

if ( $needReadSeq == 1 && ( $tabFmtOptIdx == 0 || $outputFmtSam == 1 ) )

{

	my $i ;

	my $needAdd = 1 ;

	if ( $tabFmtOptIdx != 0 )

	{

		my @cols = split /,/, $bt2_args[ $tabFmtOptIdx ] ;

		foreach my $f (@cols)

		{

			if ( $f eq "readSeq" )

			{

				$needAdd = 0 ;		

				last ;

			}

		}



	}

	else

	{

		push @bt2_args, "--tab-fmt-cols" ;

		push @bt2_args, "readID,seqID,taxID,score,2ndBestScore,hitLength,queryLength,numMatches" ;

		$tabFmtOptIdx = scalar( @bt2_args ) - 1 ;

	}

	

	if ( $needAdd )

	{

		$removeSeqCols = 1 ;

		$bt2_args[ $tabFmtOptIdx ] .= ",readSeq,readQual" ; 

	}

}



# If the user asked us to redirect some reads to files, or to suppress

# unaligned reads, then we need to capture the output from Centrifuge and pass it

# through this wrapper.
@@ -423,6 +481,17 @@ my $cmd = "$align_prog$debug_str --wrapper basic-0 ".join(" ", @bt2_args); 
# Possibly add read input on an anonymous pipe

$cmd = "$readpipe $cmd" if defined($readpipe);



# The function removes the two extra columns that we added to get the read seq and qual

sub RemoveSeqCols

{

	my $line = $_[0] ;

	my @cols = split /\t/, $line ;

	pop @cols ;

	pop @cols ;

	my $tab = "\t" ;

	return join( $tab, @cols ) ;

}



Info("$cmd\n");

my $ret;

if(defined($cap_out)) {
@@ -485,30 +554,87 @@ if(defined($cap_out)) { 
			}

		}

	}



	my $seqIndex = -1 ;

	my $qualIndex = -1 ;

	my $readIdIndex = -1 ;

	if ( $outputFmtSam == 0 )

	{

		my $outputHeader = <BT> ;

		my @cols = split /\t/, $outputHeader ;

		for ( my $i = 0 ; $i < scalar( @cols ) ; ++$i )

		{

			if ( $cols[$i] =~ /readSeq/ )

			{

				$seqIndex = $i ;

			}

			elsif ( $cols[$i] =~ /readQual/ )

			{

				$qualIndex = $i ;

			}

			elsif ( $cols[$i] =~ /readID/ )

			{

				$readIdIndex = $i ;

			}

		}

		if ( $seqIndex == -1 && scalar( keys %read_fhs) == 0 )

		{

			Error( "Must use readSeq in --tabFmtCols in order to output unaligned reads." )  ;

		}

		

		$outputHeader = RemoveSeqCols( $outputHeader )."\n" if ( $removeSeqCols == 1 ) ;

		print {$ofh} $outputHeader ;

	}

	else

	{

		$seqIndex = 9 ;

		$qualIndex = 10 ;

		$readIdIndex = 0 ;

	}



	while(<BT>) {

		chomp;

		my $filt = 0;

		unless(substr($_, 0, 1) eq "@") {

			# If we are supposed to output certain reads to files...

			my $tab1_i = index($_, "\t") + 1;

			my $tab2_i = index($_, "\t", $tab1_i);

			my $fl = substr($_, $tab1_i, $tab2_i - $tab1_i);

			my $unal = ($fl & 4) != 0;

			#my $tab1_i = index($_, "\t") + 1;

			#my $tab2_i = index($_, "\t", $tab1_i);

			#my $fl = substr($_, $tab1_i, $tab2_i - $tab1_i);

			my $unal = 0 ;

			if ( /unclassified/ )

			{

				$unal = 1 ;

			}

			$filt = 1 if $no_unal && $unal;

			if($passthru) {

				if(scalar(keys %read_fhs) == 0) {

				if(scalar(keys %read_fhs) == 0 || $seqIndex == -1 ) {

					# Next line is read with some whitespace escaped

					my $l = <BT>;

					# my $l = <BT>;

				} else {

					my $mate1 = (($fl &  64) != 0);

					my $mate2 = (($fl & 128) != 0);

					my $unp = !$mate1 && !$mate2;

					my $pair = !$unp;

					my @cols = split /\t/ ;

					my $isPaired = 0 ;

					my $pair = 0 ;

					if (  $cols[$seqIndex] =~ /_/ )

					{

						$pair = 1 ;

					}

					my $unp = !$pair ;



					# Next line is read with some whitespace escaped

					my $l = <BT>;

					chomp($l);

					$l =~ s/%(..)/chr(hex($1))/eg;

					#my $l = <BT>;

					#chomp($l);

					#$l =~ s/%(..)/chr(hex($1))/eg;



					if((defined($read_fhs{un}) || defined($read_fhs{al})) && $unp) {

						my $l ;

						if ( $qualIndex != -1 )

						{

							$l = "@".$cols[ $readIdIndex ]."\n".$cols[$seqIndex]."\n+\n".$cols[$qualIndex]."\n" ;

						}

						else 

						{

							$l = ">".$cols[ $readIdIndex ]."\n".$cols[$seqIndex]."\n" ;

						}

						if($unal) {

							# Failed to align

							print {$read_fhs{un}} $l if defined($read_fhs{un});
@@ -517,21 +643,42 @@ if(defined($cap_out)) { 
							print {$read_fhs{al}} $l if defined($read_fhs{al});

						}

					}

					my $warnedAboutLength = 0 ;

					if((defined($read_fhs{"un-conc"}) || defined($read_fhs{"al-conc"})) && $pair) {

						my $conc  = (($fl &   2) != 0);

						if     ($conc && $mate1) {

							print {$read_fhs{"al-conc"}{1}} $l if defined($read_fhs{"al-conc"});

						} elsif($conc && $mate2) {

							print {$read_fhs{"al-conc"}{2}} $l if defined($read_fhs{"al-conc"});

						} elsif(!$conc && $mate1) {

							print {$read_fhs{"un-conc"}{1}} $l if defined($read_fhs{"un-conc"});

						} elsif(!$conc && $mate2) {

							print {$read_fhs{"un-conc"}{2}} $l if defined($read_fhs{"un-conc"});

						my @seq = split /_/, $cols[$seqIndex] ;

						my @qual = ( substr( $cols[$qualIndex], 0, length( $seq[0] ) ), substr( $cols[$qualIndex], length( $seq[0] ) + 1 ) ) ;

						

						my $l1 ;

						my $l2 ;

						if ( $qualIndex != -1 )

						{

							$l1 = "@".$cols[ $readIdIndex ]."\n".$seq[0]."\n+\n".$qual[0]."\n" ;

						}

						else 

						{

							$l1 = ">".$cols[ $readIdIndex ]."\n".$seq[0]."\n" ;

						}

						if ( $qualIndex != -1 )

						{

							$l2 = "@".$cols[ $readIdIndex ]."\n".$seq[1]."\n+\n".$qual[1]."\n" ;

						}

						else 

						{

							$l2 = ">".$cols[ $readIdIndex ]."\n".$seq[1]."\n" ;

						}



						if     ( !$unal) {

							print {$read_fhs{"al-conc"}{1}} $l1 if defined($read_fhs{"al-conc"});

							print {$read_fhs{"al-conc"}{2}} $l2 if defined($read_fhs{"al-conc"});

						} else {

							print {$read_fhs{"un-conc"}{1}} $l1 if defined($read_fhs{"un-conc"});

							print {$read_fhs{"un-conc"}{2}} $l2 if defined($read_fhs{"un-conc"});

						}

					}

				}

			}

		}

		$_ = RemoveSeqCols( $_ ) if ( $removeSeqCols == 1 ) ;

		print {$ofh} "$_\n" if !$filt;

	}

	for my $k (@fhs_to_close) { close($k); }