1
$\begingroup$

I am using a pipeline from nanopore to detect structural variants (SVs) in a human sample with long-reads sequencing. The first steps of the pipeline are:

  1. index the reference genome with minimap2
  2. map the long reads to a reference genome using same tool as above
  3. Get BED coordinates from the .bam file.

While the first two steps uses minimap, the third one uses a python script written by the developers of the pipeline, which can be found here.

After the pipeline calls bamref2bed i.e. the python script above, my target.bed looks like the following:

chr1    0   249250621
chr2    0   243199373
chr3    0   198022430
chr4    0   191154276
chr5    0   180915260
chr6    0   171115067
chr7    0   159138663
chr8    0   146364022
chr9    0   141213431
chr10   0   135534747
chr11   0   135006516
chr12   0   133851895
chr13   0   115169878
chr14   0   107349540
chr15   0   102531392
chr16   0   90354753
chr17   0   81195210
chr18   0   78077248
chr19   0   59128983
chr20   0   63025520
chr21   0   48129895
chr22   0   51304566
chrX    0   155270560
chrY    0   59373566

As far as I am concerned the first three columns of a bed file refer to (I) chr name, (II) start position (III) end position. The main thing I am wondering is how come all chromosomes have a start position of 0. Also, reading about converting bam to bed files I guess the most important feature would be the coverage across the whole chromosomes, which in this file is missing.

But my main question would be: is it normal to have a start position of 0 for all chromosomes or is there something wrong here?

$\endgroup$
1
  • 2
    $\begingroup$ BED intervals are usually zero-based, half-open indexed. Issues aside from whether these particular intervals are meaningful results, a starting position of 0 is entirely normal for a BED element. Cite: genome.ucsc.edu/blog/… $\endgroup$ Oct 22 '19 at 9:12
4
$\begingroup$

I cannot access the pipeline from nanopore, but looking at the python script you provided:

with pysam.AlignmentFile(bam_file, "r") as bam:
    for tid in range(0, bam.nreferences):
        ref_name = bam.get_reference_name(tid)
        if len(list(filter(lambda x: x in ref_name, filter_names))) > 0:
            continue

        print(ref_name,
              "0",
              str(bam.get_reference_length(ref_name)),
              sep='\t')

So the script gets the length of the chromosome and by default, outputs the start as 0.Zero based coordinate system is sometimes used in UCSC systems, so I am only guessing the authors of the pipeline use this convention somewhere further down.

As for the coverage you mentioned, that is usually output as a score in the 5th column (bed format). In those situations, usually the coordinates of the alignments in the bam file are converted into bed format. In this pipeline, i think the authors are simply extracting the chromosome information out and putting it into a tabular form.

$\endgroup$
1
  • $\begingroup$ Thanks, that's useful information, I was not aware of the zero based coordinate system. In fact I managed to find the other .bed file with the coverage per chromosome, so I guess the .bed file that I mentioned is only displaying the chromosomes length. $\endgroup$
    – BCArg
    Oct 22 '19 at 9:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.