cannin/biostars_qa huggingface dataset:
Question: Hi everybody,

I have some nucleotide sequence (fasta format) that contains many sequences with invalid characters, like 'e', 'q', 'i' and 'l', there is also a item named "unavailable sequence" for some header. These characters were reported me as a error during working with MEME-chip. Could anybody please let me know how to find these characters and also the item "unavailable sequence" and remove them? Any suggestion and command warmly welcomed.

Answer: cleanup fasta:

sed -e '/^[^>]/s/[^ATGCatgc]/N/g' file.fa
Question: Dear developers

I have issue downloading the kissplice from website http://kissplice.prabi.fr/tools/kiss2rt/
KisSplice2RefTranscriptome version 1.3.2

I am not able to open this page.

could you please resolve it.

Kind regards

Archana

Answer: Dear Archie,

I guess it was a downtime on the FTP server, since now it seems to be working. I was able to download KisSplice2RefTranscriptome version 1.3.2 when writing this message. Could you please retry and check if it is working for you now?

Kind regards.
Question: The term HSP has been used in many research papers while describing BLAST output. What are High Scoring Pairs (HSP)? I'm not very successful in Googling.

Answer: HSP corresponds to the matching region between the query sequence and the database hit sequence.
Question: Hey everybody. I have written a custom python script (below) but it is too slow and I need pointers as to where the code can be improved.

What this script is doing:
I have two files, the first is a large fastq file (single-end reads) and the second is a text file with all of the header lines in the fastq file that I want to keep. The text file and fastq file are in same order with some of the fastq reads not having a corresponding header line in the text file. This script walks down both files and writes the desired reads to a new file. This script works but it is too slow, my large fastq file is ~20GB (~100millionreads) and my current estimate puts me at ~20days for it to finish (I ran it on a small sample). As someone with not much experience in bigdata, I'm sure there is something obvious that I could be doing to speed this up, any help is appreciated.

#!/usr/bin/python3

import sys
from itertools import islice
# first input (sys.argv[1]) is the original fq file name
# second input (sys.argv[2]) is the quality threshold used (e.g. '30')

print("\nentering python script")

header_file_name = sys.argv[2] + "_good_read_headers.txt"
newfile_name = sys.argv[2] + "_good_reads.fastq"

print ("...making new fastq file for: " + goodreads)

file_iterator = 0
head_iterator = 0

while True:
four_line_fq, header = [], []

with open(sys.argv[1], 'r') as name_original: # opening original fq file
lines_gen = islice(name_original, int(file_iterator * 4), int(file_iterator * 4 + 4))
for i in lines_gen:
four_line_fq.append(i)

with open(header_file_name, 'r') as name_header: # opening good header names only file
headers_gen = islice(name_header, head_iterator, head_iterator + 1)
for i in headers_gen:
header.append(i)

if len(four_line_fq) == 0: # if at end of original fq file
print("...end of original fastq file reached, done writing")
break

if four_line_fq[0] == header[0]:
with open(newfile_name, 'a+') as new:
for i in four_line_fq:
new.write(i)
head_iterator += 1
file_iterator += 1

print("...exiting python script\n")

Answer: Using Biopython as advised. Assuming you have a header.txt file with headers in one column (separate by "\n")

from Bio import SeqIO

###Open your result file
filtered_fastq = open("your_filtered_fastq.fastq", "a")
###Create an array to keep headers
header_array=[]
###Stock headers in header_array
with open('header.txt', 'r') as f:
for line in f:
###Remove the \n of your header
header_array.append(line[:-1])
###Read your fastq file
###SeqIO read the file record by record. Record is (Name of read, sequence, "+", quality)
for record in SeqIO.parse('your_fastq.fastq','fastq'):
###For each read if fastq header exist in your header_array, record will be write in your result file
if record.id in header_array:
###Write the record
SeqIO.write(record, filtered_fastq,'fastq')
###Close file
filtered_fastq.close()

I'm interested by the final calculation time, if you could write it in comment, Thanks !
Question: I am trying to edit an assembly and I have a .fa file with the sequence I am looking at.
I am ultimately trying to create a .vfc file to look at errors in my sequence. I am using the script:

bwa index draft1_assembly.fa

bwa mem draft1_assembly.fa reads_1.fastq reads_2.fastq> draft1.sam

samtools view -b -o draft1.bam -S draft1.sam

samtools sort draft1.bam draft1.sorted

samtools index draft1.sorted.bam

samtools faidx draft1_assembly.fa

samtools tview draft1.bam draft1_assembly.fa

samtools mpileup -uf draft1_assembly.fa draft1.bam | bcftools call -vc –O v - >draft1_snps_indels.vcf

I have successfully created a .sam file using:

$bwa index H_vit_Rb.fa

$bwa mem H_vit_Rb.fa reads_1.fastq reads_2.fastq> H_vit_Rb.sam

I then successfully created a .bam file using:

$samtools view -b -o H_vit_Rb.bam -S H_vit_Rb.sam

But when I try to create a .sorted file using:

$samtools sort H_vit_Rb.bam H_vit_Rb.sorted

This message appears, not letting me create a .sorted file

[bam_sort] Use -T PREFIX / -o FILE to specify temporary and final output files

Usage: samtools sort [options...] [in.bam]

Options:

-l INT Set compression level, from 0 (uncompressed) to 9 (best)

-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]

-n Sort by read name

-o FILE Write final output to FILE rather than standard output

-T PREFIX Write temporary files to PREFIX.nnnn.bam

-@, --threads INT

Set number of sorting and compression threads [1]

--input-fmt-option OPT[=VAL]

Specify a single input file format option in the form of OPTION or OPTION=VALUE

-O, --output-fmt FORMAT[,OPT[=VAL]]...

Specify output format (SAM, BAM, CRAM)

--output-fmt-option OPT[=VAL]

Specify a single output file format option in the form of OPTION or OPTION=VALUE

--reference FILE

Reference sequence FASTA FILE [null]

What does this mean and how do i fix it?

Answer: Based on the usage in the error it's quite clear:

You use `$samtools sort H_vit_Rb.bam H_vit_Rb.sorted` and you should use `samtools sort H_vit_Rb.bam > H_vit_Rb.sorted` since the result is by default written to stdout. So add the redirection...
Question: Hi All,

How can I create a lgen file which has only one individual genotyped for 300,000 SNPs. I mean in the context of sample ids, what shall I use? I created the lgen, fam files below. Then I used plink 1.9 to generate ped file (log file shown belos). The problem is that genotypes in ped file mismatch genotype in lgen file.

lgen_file:

200638550003_R01C01 200638550003_R01C01 rs1000002 A G
200638550003_R01C01 200638550003_R01C01 rs1000003 A G
200638550003_R01C01 200638550003_R01C01 rs10000030 A G
200638550003_R01C01 200638550003_R01C01 rs10000037 G G

fam_file:

200638550003_R01C01 200638550003_R01C01 0 0 0 0

**map_file:**

3 rs1000002 0 183635768
3 rs1000003 0 98342907
4 rs10000030 0 103374154
4 rs10000037 0 38924330

ped_file:

200638550003_R01C01 200638550003_R01C01 0 0 0 -9 A A G A G A G A G G G A G A A A A A G G

log_file:

PLINK v1.90b4.6 64-bit (15 Aug 2017) Options in effect: --lfile test_1 --no-fid --no-parents --no-pheno --no-sex --out test_1_1.9 --recode
Hostname: quser11 Working directory: /projects/b1042/BurridgeLab Start time: Wed Sep 13 15:23:58 2017
Random number seed: 1505334238
128715 MB RAM detected; reserving 64357 MB for main workspace.
Processing .lgen file... done.
--lfile: test_1_1.9-temporary.bed + test_1_1.9-temporary.bim +
test_1_1.9-temporary.fam written.
10 variants loaded from .bim file.
1 person (0 males, 0 females, 1 ambiguous) loaded from .fam.
Ambiguous sex ID written to test_1_1.9.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1 founder and 0 nonfounders present.
Calculating allele frequencies... done.
10 variants and 1 person pass filters and QC.
Note: No phenotypes present.
--recode ped to test_1_1.9.ped + test_1_1.9.map ... done.

Thanks Tarek

Answer: Wrote a response earlier, but it seems to have disappeared.

The main issue is that your .map file is unsorted; plink will sort it by position before generating the new .ped + .map fileset, and I'm guessing that's why there appears to be a discrepancy.

Also, it doesn't look like you posted the entire .lgen or .map input files, since they have only 4 variants but the .log file states that 10 were loaded.
Question: I'm running IGV in batch mode following the steps described in the [IGV website](https://software.broadinstitute.org/software/igv/batch).

When I enter a command like the following one in the terminal, IGV opens, runs as expected, and saves a snapshot. Everything OK.

bash igv.sh --batch=script.txt

However, this is done by opening IGV including its graphical user interface (GUI). I think this is unnecessary and probably slows the program, but so far I have not been able to run IGV in batch mode without a GUI being opened automatically, nor when I remove the `new` header command in the script, which is supposed to "open a new session" .

Any ideas on how to run IGV purely in CLI mode?

Answer: Yes. Check out the script here, which can create and run your IGV snapshots on a .bam file with a .bed file as target regions to snapshot:

https://github.com/stevekm/IGV-snapshot-automator

It uses xvfb in an automated fashion as described by Pierre.

If you just want the CLI command to use, it would be something like this:

xvfb-run --auto-servernum --server-num=1 java -Xmx4000m -jar bin/IGV_2.3.81/igv.jar -b igv_script.bat

Note that `xvfb-run` is a wrapper script that usually comes with your system's `xvfb` installation, but may be platform-specific.
Question: Hello,

I have a big newick tree (7000+ nodes) I want to pinpoint a specific list of nodes (change color, highlight or change font size ... ) is there any tool that allows me to do it by command line or python code ? if yes I need the title of the option in the tool that allows to do this

Answer: I think the [ETE toolkit][1] is exactly what you describe

Tree visualization is fully programmatic, check [this part of the tutorial][2].

[1]: http://etetoolkit.org
[2]: http://pythonhosted.org/ete2/tutorial/tutorial_drawing.html
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: I have to give a shout out for the "R Graph Gallery". When learning to make graphs/visualizations, this resource is indispensable. Sometimes when I'm trying to think of a way to visualize some data, I just browse through the gallery for inspiration. Then, if you find something you like, all the code to get you started is right there.

http://addictedtor.free.fr/graphiques/ (now dead)

UPDATED LINKS:

http://gallery.r-enthusiasts.com/allgraph.php (new home of "R Graph Gallery")

http://www.sr.bham.ac.uk/~ajrs/R/r-gallery.html

http://research.stowers-institute.org/efg/R/
Question: Hi,

Anybody has the list of virulent and non-virulent strains name in Staphylococcus aureus?
Pls respond

Answer: I would use [PATRIC database][1] (Bioinformatics Database and Analysis Resource Center).

There is a list of all *Staphylococcus aureus* strains (currently: [8863 genomes][4]) [*last update: 30-Jan-2017*].

You can then **filter the list of organisms according to**:

1. Disease (e.g. *pneumonia*, *gastroenteritis*, *fever* etc.)
2. Host name (e.g. *human*)
3. Antibiotic resistance (i.e. *resistant*, *susceptible*, *intermediate*)

Finally, you can download the table as *Text* or *Excel* files.

**Note:** PATRIC articles: [2014][2], [2017][3].

[1]: https://www.patricbrc.org
[2]: https://doi.org/10.1093/nar/gkt1099
[3]: https://doi.org/10.1093/nar/gkw1017
[4]: https://www.patricbrc.org/view/Taxonomy/1280#view_tab=genomes
Question: Hi,

I am getting error while using StringTie (1.3.3b) on aligned reads from both HISAT2 (2.1.0) and STAR (2.5.3a). Please help me to rectify the problem. The error is related to **missing headers in the alignment files** but when I print header using samtools I can see the headers in those alignment files (hisat2 and star aligned).

**Steps used:**

**HISAT2 alignment:**

>hisat2 --dta --known-splicesite-infile spliceSites.txt --phred33 --qc-filter -x ht2_index/hg38 -U sample1.fastq.gz | samtools view -bS -o sample1_hisat2.bam - ;

>>samtools sort -o sample1_hisat2_sorted.bam -O sample1_hisat2.bam ;

>> samtools view -H sample1_hisat2_sorted.bam | head

@HD	VN:1.0	SO:coordinate
@SQ	SN:1	LN:248956422
@SQ	SN:10	LN:133797422
@SQ	SN:11	LN:135086622
@SQ	SN:12	LN:133275309
@SQ	SN:13	LN:114364328
@SQ	SN:14	LN:107043718
@SQ	SN:15	LN:101991189
@SQ	SN:16	LN:90338345
@SQ	SN:17	LN:83257441

> stringtie sample1_hisat2_sorted.bam -G hg38.GTF -B -e -x -C hg38.GTF -p 8 -A hg38_gene_abund.tab -o hg38_assembled.gtf

**Error:**

[samopen] no @SQ lines in the header.
[sam_read1] missing header? Abort!

**STAR alignment:**

>STAR --readFilesCommand zcat --outSAMmapqUnique 255 --outSAMattributes All --outSAMattrIHstart 0 --outSAMtype BAM SortedByCoordinate --outSAMunmapped None --outFileNamePrefix path/sample1_star --runMode alignReads --genomeDir genome_dir/STAR_index/ --readFilesIn sample1.fastq.gz ;

>> samtools view -H sample1_starAligned.sortedByCoord.out.bam | head

@HD	VN:1.4	SO:coordinate
@SQ	SN:1	LN:248956422
@SQ	SN:10	LN:133797422
@SQ	SN:11	LN:135086622
@SQ	SN:12	LN:133275309
@SQ	SN:13	LN:114364328
@SQ	SN:14	LN:107043718
@SQ	SN:15	LN:101991189
@SQ	SN:16	LN:90338345
@SQ	SN:17	LN:83257441

> stringtie sample1_starAligned.sortedByCoord.out.bam -G hg38.GTF -B -e -x -C hg38.GTF -p 8 -A hg38_gene_abund.tab -o hg38_assembled.gtf

**Error:**

[samopen] no @SQ lines in the header.
[sam_read1] missing header? Abort!

Answer: I think the issue is not in your file but in your `stringtie` command. It looks like it could use a bit of tweaking. The `-x` option requires a parameter and is interpreting the `-C hg38.GTF` as input, resulting in the error message. Try removing the `-x` (or specifying a file) and seeing if that helps.
Question: Hi All,

Sorry for technical question.

I am trying to make a boxplot of gene expression data for some samples, and wanted to highlight specific points (expression value) based on gene name "TP53". I am able to highlight color in scatter plot but couldn't in boxplot, my code for scatter plot is

plot(datamtx, col=ifelse(rownames(datamtx)=="TP53","red","black"))

I tried similar function for Boxplot but doesn't work..

boxplot(datamtx, labels=ifelse(rownames(datamtx)=="TP53","red","black"))

and

boxplot(datamtx, col=ifelse(rownames(datamtx)=="TP53","red","black"))

please help.

Thanks

Answer: # Simulate a two column (tumor and normal) data with 10 genes. Gene names as row names
df=data.frame(normal=rnorm(10,60,5), tumor=rnorm(10,70,4))
row.names(df)=paste0("gene_",seq(1:10))
# Draw box plot
boxplot(df)
# Highlight the gene "gene_4" on box plot. Gene names are represented by row names
stripchart(df["gene_4",], vertical = T, method = "jitter", add=T, pch=20, col=as.factor(colnames(df)))
Question: Hi guys, is there a way to remove from Seurat3 object (all slots) a set of cells listed in a vector?

I tried:

> my_object_without_cells = subset(object = my_dataset, cells = cells_to_be_removed)

but the cells are still there. I checked the n of cells and the n of genes and it is unchanged.

Thank you in advance

Answer: ## Example data:
pbmc_raw <- read.table(
file = system.file('extdata', 'pbmc_raw.txt', package = 'Seurat'),
as.is = TRUE
)
pbmc_small <- CreateSeuratObject(counts = pbmc_raw)

length(colnames(pbmc_small)) ## 80

## list the cells (=colnames) that you want to remove:
toRemove <- c("ATGCCAGAACGACT", "CATGGCCTGTGCAT")

## filter them out:
pbmc_filtered <- pbmc_small[,!colnames(pbmc_small) %in% toRemove]

length(colnames(pbmc_filtered)) ## 78
Question: Hi,
This might be an obvious question...
I have an SAM file looking something like this, my question is: How to you find the contig name (if it is possible). I have looked up the definitions of the header but do not quite get if there is any contig names, is it the same as the QNAME and if so do you need to decode it like the Flags?

QNAME	FLAG	RNAME	POS	MAPQ	CIGAR	MRNM	MPOS	ISIZE	SEQ	QUAL	OPT
@HD	VN:1.0	SO:unsorted
@SQ	SN:1	LN:249250621
@SQ	SN:2	LN:243199373
@SQ	SN:3	LN:198022430
@SQ	SN:4	LN:191154276
..........................................................
@SQ	SN:GL000245.1	LN:36651
@SQ	SN:GL000246.1	LN:38154
@SQ	SN:GL000247.1	LN:36422
@SQ	SN:GL000248.1	LN:39786
@SQ	SN:GL000249.1	LN:38502
@RG	ID:S0143_RC	SM:S0143	PU:no_index	LB:exome	CN:Lab	PL:ILLUMINA
@PG	ID:bowtie2	PN:bowtie2	VN:2.0.2
HISEQ:179:H9A7VADXX:1:1101:2961:2239	83	11	58711034	42	101M	=	58710857	-278	TCTGAAGTGGAGCTTCTAGTATCCCCAGGAGCGCGAAGTGAACACGGAAGGTACCTGCAGGATCCAATTGTGTCCATTGATCTCTCAGAGTGGCTGAGGNT	BFFFFFFFFFFFFFFFFFFBFFFFFFFFFFFFIFIIIIFIIFFFIIIIIIFFFIIIIIIIIIIIIIIIIIIIIIIIIIIIFIIIIIIIFFFFFFFFFB0#B	AS:i:-1	XS:i:-36	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:99A1	YS:i:0	YT:Z:CP	RG:Z:S0143_RC
HISEQ:179:H9A7VADXX:1:1101:2961:2239	163	11	58710857	42	99M	=	58711034	278	CACTTTACCTTTTTTGTCTATAAATTCATTTTGACCACGAGGCACCCCCGGAGCCTCGGTGAATCTGCTGTGATTTTGTAGGCTGCCGGATTCACAAAT	BBBFFFFFFFFFFIIIFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIFFFFFFFFBBBFBFFFFFFFFFBFFFFBBBBBFFFFFFFBFFFFFFFFB	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:99	YS:i:-1	YT:Z:CP	RG:Z:S0143_RC
HISEQ:179:H9A7VADXX:1:1101:2867:2247	83	17	67173782	42	101M	=	67173769	-114	AGTGAGAACTAAATTAAATAATGTATAGTGAGGCCAGATGTGGTGGCTCACACCTGTAATCCCAGCACTTTGAGAGGCTGAGGCGGGTGGATCACCTGAGG	BBFBFFBBBBFFFFFFFFBFBFBBBBBBFFFFFFFFBFIFFFFFFF

Answer: I assume you want to know to which chromosome/contig your reads are aligned. That is the RNAME. QNAME is the name of your read.
Your header contains lines with SQ, indicating the chromosome names are `1`, `2`, `3`, `4`, `5`, `6`,... Other annotations may use `chr1`, `chr2`, `chr3`, `chr4`, `chr5`, `chr6`,...

From your example: the first two reads (a pair) are aligned to `11`, the next pair to chromosome `17` and the last line to chromosome `9`.
Question: Hello,

I would like to check whether L1 repetitive elements are modulated between my treatment and control via RNA-Seq. I have read several papers that have done so but their methods are not clear enough for me as a biologist to reproduce. I have analyzed my data using the Tuxedo suite and have analyzed the "unique" genes. I am wondering what modifications have to be taken into account to accommodate the repetitive nature of LINEs.
1- I have an understanding that some aligners filter out reads that map to several places in the genome. Are my LINE reads being filtered out by tophat?
2- If so, how do I align them?
3- when using cufflinks, intead of using RefSeq, I am assuming I would have to use a repetitive element model?

Thank you!

Answer: If I understand your question correctly, you want to identify the expression levels of of LINE-1 repeats in your RNA-Seq samples? If that is the case follow these instructions.

1. Make a GTF format file your repeat elements or download them from UCSC/Galaxy and run
2. `tophat -G LINE1-repeats.gtf -o treat-rnaseq yourgenome_ebwt_base treat-rnaseq.fastq`
3. `tophat -G LINE1-repeats.gtf -o control-rnaseq yourgenome_ebwt_base control-rnaseq.fastq`
4. `cuffdiff -G LINE1-repeats.gtf treat-rnaseq.bam control_rnaseq.bam`

Step 2 and 3 do the map the RNA-Seq reads to your repeat elements in the genome.

Step 4 calculates the differential expression of your repeat elements in your treatment and control.
Question: Hello! I have a question.

After use MAFFT to aligne a multifasta file, I generated a consensus sequence using seqinr::consensus in R like this:

">seq1
ACCATTCA-AAT--GAAACAAAAATTT-GACCACA-CGTGTG-TGTTCACACACTGGAC-CCAAACAACCC-GAAC"

I would like to design a pair of primers targeting this sequence. How I am afraid the gaps can interfere in this process.

Can I design primers for sequences containing gaps in the middle of it? If yes, there a maximum percentage of gaps the sequence can have?

Thank you in advance.

Answer: You can do this using degenerate primers. At each gap, you would indicate a mixed base (i.e. more than one nt; see link below). The result will be a mix of primers, each with a different nt at that position. It's best to use as few mixed bases as possible in order to minimize the number of primers in the tube. More primers mean more potential unwanted products, as well as less of the correct primer in your tube. If you know the amino acid sequence, you can possibly deduce what the gap could be using a codon table. This could help you add a mixed base of 2 or 3 (e.g. A/G or A/C/T) instead of all 4. The second and third links below give guidelines for designing degenerate primers.

[Mixed Bases][1]

[Degenerate Primer Design][2]

[Degenerate Primer Guidelines][3]

[1]: https://www.idtdna.com/pages/products/custom-dna-rna/mixed-bases
[2]: https://bitesizebio.com/18992/a-primer-for-designing-degenerate-primers/
[3]: https://www.qiagen.com/us/knowledge-and-support/knowledge-hub/bench-guide/pcr/introduction/guidelines-for-degenerate-primer-design-and-use
Question: library(GenomicRanges)

gr=GRanges(seqnames=c("chr1","chr2","chr2"),
ranges=IRanges(start=c(50,150,200),end=c(100,200,300)),
strand=c("+","-","-")
)

I want to get all intervals between 200-300 on chromosome 2, the minus strand. How do I do that?

gr[seqnames(gr) == "chr2" & start(gr) > 200 & end(gr) < 300 & strand(gr) == "-"]

does not work as it does not find the overlaps with 200-300, but rather the intervals strictly contained in that range.

Answer: Make a query range, then `subsetByOverlaps` (sorry in comments by mistake mentioned `findOverlaps`):

q=GRanges(seqnames="chr2",
ranges=IRanges(start = 200, end = 300),
strand="-")

subsetByOverlaps(gr, q)
Question: I have a very large VCF file where the 'ID' column is a unique ID comprising of 'chr:bp'. I would like to update the 'ID' column to dbSNP IDs.

I have downloaded a bed file [chr, from, to, rsid], which I have sorted and tabix indexed. The bedfile is for hg19, which is correct for my data, chromosomes are formatted with 'chr' and their is no header.

It seems that the BCFtools annotate function does allow 'ID' column to be updated, but I am not clear how. I have tried;

i) `bcftools annotate -a dbsnp.bed.gz -c 'CHROM,POS,-,ID' my.vcf.gz`

ii) `bcftools annotate -a dbsnp.bed.gz -c 'CHROM,FROM,TO,ID' my.vcf.gz`

neither of which updated the ID column. I also tried removing the 'ID' first from the VCF `-R` , and piped the vcf into the two commands above. Perhaps this is not the right tool? Any advice appreciated.

Answer: bcftools annotate -c CHROM,FROM,TO,ID -a my_ids.bed.gz -o output.vcf input.vcf.gz

works fore me.

It did also take me some time to get it to work and it is hard to debug were the mistake is:

Some things to try / check:

1. VCF is 1 based, BED is zero based. POS 10 in VCF is start 9 end 10 in BED. https://genome.ucsc.edu/FAQ/FAQformat#format1

2. Make sure your chromosome names match exactly, Chr1 and chr_1 are not the same for bcftools.

3. Remove the quotes around the column list 'CHROM,FROM,TO,ID' -> CHROM,FROM,TO,ID

4. Test with a very small subset of your VCF and BED file that should produce an annotated VCF file. This makes it faster to debug and test different options / formattings until you get it right.