Reading Input Files

The content on this page is meant for people interested in developing VNTRseek, either external users or users in our lab. Content here may change as we make changes to the underlying code in VNTRseek. This page should serve as a reference for guiding design decisions and for listing assumptions made by the various parts of the pipeline.

All functions responsible for reading input files are found in lib/ProcInputReads.pm.

VNTRseek 1.09 adds support for BAM files, but it can also read FASTA and FASTQ files. Currently this is done by splitting up the job of reading sequences from the processing of those sequences.

A function called fork_trf is called for up to $max_processes times, a user-defined variable which should be at most the number of available cores on the system. Each call to fork_trf is passed a unique value which is the index of the call (ie, 4 for the the 5th call). fork_trf will fork a process and then call the appropriate file reading function.

The file reader is determined using the $input_format argument to fork_trf, which is looked up in a hash with sequence formats as keys and references to file reading functions as values. fork_trf will then open a pipe to TRF, call the file reader repeatedly in a loop, and pass the output to another function, pipe_to_trf, which will process the read as needed and then pass it to the TRF pipe.

File readers are expected to return a 2-element list on every call: a FASTA format header, and a sequence with no newline characters in it. When there are no more reads in the input file, the file reader should return undef.

Important: FASTA headers are expected to be unique. As such, the comment portion of the FASTA/FASTQ header should not be removed if this means that the headers will no longer be unique. SRA reads, for example, add the pair indicators, /1 and /2, after the comment portion so it should be kept.

File readers take the following as required arguments: the directory containing all the input files, the detected compression format, the number of files processed so far (this is also the index of the file to be processed next), a reference to the total number of files to process, and a reference to a list of file names relative to the directory.

Writing file readers

Additional input formats can be supported by simply adding in another reader function which takes the same arguments described above, as long as it produces the expected output: a two-element list consisting of a FASTA format header and a read sequence with no whitespace.

The file reader function can use the supplied arguments in any way necessary to read the input file. Below is an example of how the FASTA/FASTQ parser from VNTRseek 1.09 works:

sub read_fastaq {
    my ( $input_dir, $compression, $files_processed,
        $files_to_process, $filelist )
        = @_;

    # Since we are using seqtk, use pipe open mode
    my $openmode = "-|";

    warn "Processing file " . $filelist->[$files_processed] . "\n";
    my $filename
        = ( ($compression) ? $decompress_cmds{$compression} : "" ) . '"'
        . "$input_dir/"
        . $filelist->[$files_processed] . '"'
        . "| seqtk seq -a -S";

    # $files_processed contains how many files processed so far.
    # Use to index into filelist 
    # warn $filename;
    local $/ = ">";

    # warn "Filename/command = '$filename'\n";
    open my $fasta_fh, $openmode, $filename
        or die "Error opening file " . $filename;

   # Consume first empty record because of the way $/ splits the FASTA format.
    <$fasta_fh>;
    return sub {
        local $/ = ">";
        my $fasta_rec = <$fasta_fh>;
        return () unless ($fasta_rec);
        my ( $header, $seq ) = split( /\n+/, $fasta_rec );
        chomp $header;
        chomp $seq;
        # warn "header: '$header'";
        # my $seq = join( "", @seqlines );

        # warn "seq: '$seq'";
        return ( ">". $header, $seq );
    };
}

This function ignores the compression format when deciding how to open the file. Since it uses seqtk for reading in the files, the file open mode is always a pipe open, -|.

The $files_processed variable here is used as a simple index into the list of files in the directory ($filelist). Multiple instances of the file reading function may be run (as described earlier), each one getting a unique value for $files_processed as they are spawned, guaranteeing that two instances may never process the same input file.

We don't make use of seqtk's ability to read gz files here, and instead use the proper decompression command in order to support other formats. The string for the command is stored in a hash and is looked up using the $compression variable. seqtk is given the -a option to ensure the output is FASTA, and the -S option to remove whitespace (needed after seqtk 1.2). We then read in the FASTA output from seqtk, expecting the header to be on one line and the entire sequence to be on the next (rather than split over multiple lines).

We use ">" as a record separator to simplify reading the output from seqtk, but this means we need to skip the first empty record due to how the record separator splits records.

Finally, the header (with the ">" character added back in) and the sequence are returned to the caller, fork_trf. fork_trf will pass on the header and sequence, along with other arguments, to pipe_to_trf.

For our current BAM file implementation, things are a bit different. Instead of using the $filelist array reference that was passed in, we use samtools index to retrieve a list of the chromosomes in the file, and build a list out of that. This is done in advance of the actual forking, but after the directory has been scanned so that we know the format to expect is BAM, and to replace the $filelist array with a list of regions in the BAM file rather than actual files. The $files_processed variable is still used as an index, only it indexes the list from samtools.

Passing to TRF

The pipe_to_trf function takes in arguments specifying how to process the sequence, if necessary; a file handle to the TRF process pipe; and a FASTA header and sequence as strings. If the argument $reverse_read is 1, this function will also produce the reverse complement.

Important: the reverse complement of the sequence must always follow the forward sequence. This order is assumed by the redund binary called in the next step. Any changes to this function must preserve this order.

The TRF pipe must be a file handle and the caller must ensure that the pipe is open. The file handle does not have to be a pipe to a TRF process; for example, during testing, this can be a pipe to STDOUT or to a file.