Output

resquiggle

f5c resquiggle aligns raw signals to basecalled reads. f5c resquiggle output is explained below.

resquiggle TSV output

The default output is an intuitive TSV format with the following columns.

If a corresponding base has no corresponding signal samples, a . will be printed.

resquiggle PAF output

The above TSV output is bulky. Specifying -c will produce a condensed output in PAF-like format (inspired by UNCALLED) with the following columns: The query is the raw-signal and the target is the basecalled-read.

For DNA reads, column 8 is a closed coordinate (inclusive) and column 9 is an open coordinate (not-inclusive); and, column 8 coordinate is always smaller than column 9. However, this is different in direct-RNA reads because sequencing of direct-RNA happens in the reverse direction (3’->5’) and therefore the raw signal is also in the reverse direction. The basecalled read output by basecallers are however in the correct direction (5’->3’). Thus, For RNA reads, column 8 coordinate will be larger than that on column 9. column 9 is a closed coordinate (inclusive) while column 8 is an open coordinate (not-inclusive) in this case, in contrary to DNA.

10,11 and 12 are yet to be decided. column 12 is always 255 for now.

Following optional tags are present, which are subsequently described in detail after the table:

The sh and sc tag values can be used to scale the raw signal to the pore model. This can be done as: scaled pA current values = (pA - sh) / sc where pA = (raw_signal + offset) * range / digitisation)).

ss tag

ss string is a custom encoding that compacts the signal-base alignment. It can be thought of as an extended version CIGAR string that accommodates for signal alignment needs. This ss string was inspired by the –sam option in Nanopolish for eventalign.

Consider the example 8,5,4,8I4,3D4,5, for DNA. This means 8 signal samples map to the starting base of the sequence; the next 5 samples to the next base, the next 4 samples to the next base; 8 next samples are missing a mapping in the basecalled read (insertion to reference); 4 samples map to the next base; 3 bases in the basecalled read have no corresponding signal samples (deletion); 4 samples map to the next base; and 5 samples map to the next base. Note that the start indexes of read and signal are the absolute values in columns 8 and 3 above. the ss string is relative to this.

The ,, D and I can be though of as three different operations. The , after a number means step one base in the basecall reference, while step number of samples preceding the , in the raw signal. The D after a number means step number of bases preceding the D in the basecalled read and no stepping in the raw signal. The I after a number means step number of samples preceding the I in the raw signal and no stepping in the basecalled read.

DNA Example

To make things further clear, given below is an illustration for an alignment of DNA raw-signal to a basecalled read.

Each dot in the top sequence in the illustration represent a k-mer and the number above each dot is the corresponding k-mer index. If the basecalled read is ACGGTAACTATAC and assuming the k-mer size in the k-mer model is 6, the 0th k-mer is ACGGTA, 1st k-mer is CGGTAA … and the 7th k-mer is CTATAC. Each dot in the bottom sequence in the illustration represent a raw-signal sample and the number below each dot is the corresponding signal index.

The TSV output from resquiggle will look like below:

The PAF output from resquiggle will look like below (the header is not present in the actual output):

direct-RNA Example

Now see the illustration below for direct-RNA.

Note that the RNA is sequenced 3’->5’ end, so the raw signal is 3’->5’ direction. As the basecaller outputs the basecalled read in 5’->3’ direction, the basecalled read is reversed to be 3’->5’ in the illustration (note: indices in illustration denote the actual index in the basecalled read in 5’->3’ direction). If the basecalled read in 5’->3’ direction is ACGGUAACUAUA and assuming the k-mer size in the k-mer model is 5, the 0th k-mer is ACGGU, 1st k-mer is CGGUA … and the 7th k-mer is CUAUA.

The TSV output will look like below:

The PAF output will look like below:

C code snippet to parse ss tag

A C code snippet that converts the value in the ss tag (a readable code which is not optimised) is given below:

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>

#define MAX_LEN_KMER 20000

int main(){

    //paf fields
    char *ss="2,3,2,2I1,2D2,3,";
    int start_raw=2; int end_raw=17; int len_raw_signal=20;
    int start_kmer=0; int end_kmer=8; int len_kmer=8;

    // Raw signal start index for the corresponding k-mer and Raw signal end index for the corresponding k-mer
    int st_raw_idx[MAX_LEN_KMER]; int end_raw_idx[MAX_LEN_KMER];

    //intialise to -1
    for(int i=0; i<MAX_LEN_KMER; i++){ st_raw_idx[i]=end_raw_idx[i]=-1; }

    int st_k = start_kmer; int end_k = end_kmer; //if DNA, start k-kmer index is start_kmer column in paf and end k-kmer index is end_kmer column in paf
    int8_t rna = start_kmer > end_kmer ? 1 : 0; //if RNA start_kmer>end_kmer in paf
    if(rna){ st_k = end_kmer; end_k = start_kmer; } //if RNA, start k-kmer index is end_kmer column in paf and end k-kmer index is start_kmer column in paf

    int i_k = st_k; int i_raw = start_raw; //current k-mer index and current raw signal index

    //buffer for storing digits preceding each operation and its index
    char buff[11]; int i_buff=0;

    while(*ss){

        if(*ss==',' || *ss=='I' || *ss=='D'){
            if(i_buff <= 0){ fprintf(stderr,"Bad ss: Preceding digit missing\n"); exit(1); }//if nothing in buff

            buff[i_buff]=0; //null terminate buff
            int num = atoi(buff);
            if(num < 0){ fprintf(stderr,"Bad ss: Cannot have negative numbers\n"); exit(1); }
            i_buff=0; buff[0]=0; //reset buff

            if(*ss=='I'){ //if an insertion, current raw signal index is incremented by num
                i_raw += num;
            } else if(*ss=='D'){ //if an deletion, current k-mer index is incremented by num
                i_k += num;
            } else if (*ss==','){ //if a mapping, increment accordingly and set raw signal indices for the current k-mer
                end_raw_idx[i_k] = i_raw; i_raw += num;
                st_raw_idx[i_k] = i_raw; i_k++;
            }
        } else {
            if(!isdigit(*ss)){ fprintf(stderr,"Bad ss: A non-digit found when expected a digit\n"); exit(1); }
            buff[i_buff++]=*ss;
        }
        ss++;
    }

    if(i_raw!=end_raw){ fprintf(stderr,"Bad ss: Signal end mismatch\n"); exit(1); } //current raw signal index should be equal to end_raw
    if(i_k!=end_k){ fprintf(stderr,"Bad ss: Kmer end mismatch\n"); exit(1); } //current k-mer index should be equal to end_k

    for(int i=st_k; i<end_k; i++){
        if(end_raw_idx[i]==-1){
            if(st_raw_idx[i] != -1) { fprintf(stderr,"Bad ss: This shoud not have happened\n"); exit(1); }//if st_raw_idx[i] is -1, then end_raw_idx[i] should also be -1
            printf("%d\t.\t.\n", rna ? len_kmer-i-1 : i);
        }else {
            printf("%d\t%d\t%d\n", rna ? len_kmer-i-1 : i, end_raw_idx[i], st_raw_idx[i]);
        }
    }
}

Note that the beginning of the alignment could be crude due to the adaptor and barcode in the raw signal. In future versions, an option for this trimming may be introduced. Until then, if there is a huge step in the alignment in the beginning, that part could be trimmed by the user.

eventalign

eventalign TSV output

The default TSV output is as below:

Command line options can be used to modify default columns or print additional columns. --samples will print two additional columns, namely start_idx and end_idx.

Following optional columns are controlled via command line options:

eventalign PAF output

Specifying -c will generate the output in a PAF-like format (much smaller than the TSV output, f5c v1.3 or higher). This is similar to the PAF output format in f5c resquiggle with major difference being that the “basecalled read” in resquiggle is now the “refererence sequence” in eventalign. Assumming that the reader is well familiarised with the PAF output explained in resquiggle, that information is not repeated here, instead only a summary is given. Unlike in resquiggle, the strand column (column 5) can be now both ‘+’ and ‘-‘. The query is the raw-signal and the target is the reference.

For DNA reads, column 8 is a closed coordinate (inclusive) and column 9 is an open coordinate (not-inclusive); and, column 8 coordinate is always smaller than column 9. However, this is different in direct-RNA reads because sequencing of direct-RNA happens in the reverse direction (3’->5’) and therefore the raw signal is also in the reverse direction. The sequences in the reference transcriptome are however in the correct direction (5’->3’). Thus, For RNA reads, column 8 coordinate will be larger than that on column 9. column 9 is a closed coordinate (inclusive) while column 8 is an open coordinate (not-inclusive) in this case, in contrary to DNA.

10,11 and 12 are yet to be decided. column 12 is always 255 for now.

Following optional tags are present, which are subsequently described in detail after the table:

The sh and sc tag values can be used to scale the raw signal to the pore model. This can be done as: scaled pA current values = (pA - sh) / sc where pA = (raw_signal + offset) * range / digitisation)). Please refer to the explanations in f5c resquiggle about the ss tag. Some examples are given below.

DNA examples

Positive strand

Assume we have a read signal named rid0 of 1000 signal samples, mapped to a reference contig named ctg0 of 35 bases Assume a k-mer size of 6. We have a total of 30 k-mers in the reference. Assume the signal-reference alignment looks like in teh figure below. Assume that the 12-24th bases (0-index; bed-like) inthis contig are TTGATGGTGGAA. Thus, 12th kmer is TTGATG, 13th k-mer is TGATGG, 14th k-mer is GATGGT, .. and the 18th k-mer is GTGGAA.

The TSV output from resquiggle will look like below (assume --print-read-name and --signal-index are provided):

The PAF output from eventalign will look like below (the header is not present in the actual output):

Negative strand

Assume we have a read signal named rid1 of 1000 signal samples, mapped to a reference contig named ctg0 of 35 bases Assume a k-mer size of 6. We have a total of 30 k-mers in the reference. Assume the signal-reference alignment looks like in the figure below (note: indices in illustration denote the actual index in the + strand of the reference genome in 5’->3’ direction).

Assume that the 11-24th bases (0-index; bed-like) in this contig are ATTGATGGTGGAA. Thus, 11th kmer is ATTGAT, 12th k-mer is TTGATG, 13th k-mer is TGATGG, .. 17th k-mer is GGTGGA and the 18th k-mer is GTGGAA.

Th negative strand is like:

5' ATTGATGGTGGAA 3' + strand
   |||||||||||||
3' TAACTACCACCTT 5' - strand

The Reverse complement is thus TTCCACCATCAAT. The 11th k-mer ATTGAT in the + strand relates to ATCAAT in the - strand, 12th k-mer TTGATG relates to CATCAA, 13th k-mer TGATGG relates to CCATC ,… , 17th k-mer GGTGGA relates to TCCACC and 18th k-mer GTGGAA relates to TTCCAC.

The TSV output from resquiggle will look like below (assume --print-read-name and --signal-index are provided):

The PAF output from eventalign will look like below (the header is not present in the actual output):

RNA examples

Positive strand

Assume we have a read signal named rid0 of 3000 signal samples, mapped to a reference transcript (or can be a ctg in the reference genome) named trn0 of 65 bases. Assume a k-mer size of 5. We have a total of 61 k-mers in the reference.

Assume the signal-reference alignment looks like in the figure below. Note that the RNA is sequenced 3’->5’ end, so the raw signal is 3’->5’ direction. However, as transcipts in the reference are in 5’->3’ direction, the transcript is reversed to be 3’->5’ in the illustration (note: indices in illustration denote the actual index in the transcriptd in 5’->3’ direction).

Assume that the 45-56th bases (0-index; bed-like) in this transcript in 5’->3’ direction is GAGAGCCCTGA. Then, 45th kmer is GAGAG, 46th k-mer is AGAGC, 47th k-mer is GAGCC, .. and the 51st k-mer is CCTGA.

The TSV output from resquiggle will look like below (assume --print-read-name and --signal-index are provided):

The PAF output from eventalign will look like below (the header is not present in the actual output):

Note that start_kmer and end_kmer are otherway round compared to DNA.

Negative strand

Assume we have a read signal named rid1 of 500 signal samples, mapped to a reference contig named ctg1 of 20 bases. Assume a k-mer size of 5. We have a total of 24 k-mers in the reference. Assume the signal-reference alignment looks like in the figure below (note: indices in illustration denote the actual index in the + strand of the reference genome in 5’->3’ direction).

Assume that the 4-13th bases (0-index; bed-like) in this contig are AAATGGCTGA. Thus, 4th kmer is AAATG, 5th k-mer is AATGG, .. 8th k-mer is GGCTG and the 9th k-mer is GCTGA.

Th negative strand is like:

5' AAATGGCTGA 3' + strand
   ||||||||||
3' TTTACCGACT 5' - strand

The Reverse complement is thus TCAGCCATTT. The 4th k-mer AAATG in the + strand relates to CATTT in the - strand, 5th k-mer AATGG relates to CCATT ,… , 8th k-mer GGCTG relates to CAGCC and 9th k-mer GCTGA relates to TCAGC.

The TSV output from resquiggle will look like below (assume --print-read-name and --signal-index are provided):

The PAF output from eventalign will look like below (the header is not present in the actual output):

eventalign SAM output

Specifying -a will generate the output in SAM format (much smaller than the TSV output, f5c v1.3 or higher). This output SAM file will contain the following two additional tags added to the S/BAM file given as the input to f5c:

si tag contains four comma separated values start_raw, end_raw, start_kmer and end_kmer, respectively. Those values have the same as the columns 3,4,8 and 9 in the PAF format explained above when specified along with -c.

call-methylation

Same as Nanopolish output.