3

Starting with a chromosome # and position, I am trying to get chromStart and chromEnd values for a .bed file, but I am not sure how to calculate chromEnd when I have a variant that is an insertion or deletion.

My file is not a VCF file but it gives the chromosome number, chromosomal position (pos), reference (ref) allele, and alternate (alt) allele for each genetic variant, like you would find in a VCF file. The tricky part is that not all of my variants are single nucleotide base changes--I also have insertions and deletions that can be multiple bases long.

For SNPs I am setting chromStart = pos - 1, and chromEnd = start + 1 = pos. (Note, .bed format requires 0-based indexing and the end position is exclusive. See https://genome.ucsc.edu/FAQ/FAQformat.html)

I am unsure of how to handle indels though (insertions and deletions). Currently I am adding the length of the reference allele to the starting position to get the end position (chromEnd = start + len(ref)). I am not sure if this is right though.

I will be using this .bed file with a .bed file from Gencode and bedtools. Does anyone know the standard way the research community (or at least Gencode) treats indels when determining chromosomal position for .bed files? Thanks!

Sarah
  • 486
  • 1
  • 4
  • 18
  • 2
    Usually people don't convert such things to BED files. bedtools understands VCF format, so why not convert to that, since it seems most of the information is already there? – Devon Ryan Oct 04 '17 at 06:53
  • 4
    Could you explain why you are trying to make this into a bed file? Is it for displaying it in a genome browser? I don't really see how a bed file could even show an indel, unless it were by using the coordinates of the reference genome. – terdon Oct 04 '17 at 07:24
  • Hi, @terdon I am trying to make this into a BED file so that I can find which variants overlap exonic areas in the genome. I have a BED file from Gencode listing protein coding regions, and I want to know which of my variants overlap with the exome. – Sarah Oct 04 '17 at 17:54
  • Hi @Devon Ryan, if the best thing to do would be to convert it to VCF format, will the bedtools intersect command be able to compare a VCF file with a BED file, or would I need to convert the Genecode BED to VCF? – Sarah Oct 04 '17 at 17:55
  • bedtools is able to compare BED and VCF files, you don't need to convert anything at that point. – Devon Ryan Oct 05 '17 at 06:57
  • The best way to determine which variants affect exomes is to run through an annotation tool like the VEP (http://www.ensembl.org/info/docs/tools/vep/index.html) for which you need a VCF. – Emily_Ensembl Oct 06 '17 at 15:30

1 Answers1

1

You can convert to VCF and use BEDOPS vcf2bed --insertions, vcf2bed --deletions and vcf2bed --snvs to segregate the VCF into three subsets.

For convenience, BEDOPS conversion tools deal with index shifts, including vcf2bed. The resulting BED file is 0-indexed, half-open.

Alex Reynolds
  • 3,135
  • 11
  • 27
  • Thanks @Alex Reynolds. Is there any way to do this but without segregating the VCF into subsets? – Sarah Oct 20 '17 at 01:29
  • Coming back year later I can say genomic coordinate systems can be tricker than I ever imagined years ago when I wrote this question, especially when converting between different formats lol. Here is a good article explaining 1 based vs 0 based positions. https://www.biostars.org/p/84686/ . VCF format is generally 1 based and BED format is generally 0 based. But also some programs use 1-ish or 0-ish indexing so it gets confusing (looking at you VarSeq ;) ) – Sarah Mar 28 '24 at 23:28