<p>Several years ago I created rsidx to address this particular use case. <a href="https://github.com/bioforensics/rsidx" rel="noreferrer">https://github.com/bioforensics/rsidx</a></p>
<pre><code># VCF should be sorted by genomic coordinates and indexed by tabix
rsidx index dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx
rsidx search dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx rs3114908 rs10756819
</code></pre>

CC BY-SA 4.0
2023-07-10T13:14:09.930Z
2023-07-10T13:14:09.930Z
Daniel Standage
<p>rsids_chr1.txt contains a list of rsids, one per line.</p>
<pre><code>bcftools view --include ID==@rsids_chr1.txt dbSNP_v156_GRCh38.p14.vcf.gz -Ou | bcftools query -f '%CHROM  %ID %POS  %REF  %ALT{0}\n' &gt; dbSNP_v156_GRCh38.p14_AF.vcf
</code></pre>

CC BY-SA 4.0
2023-07-10T11:02:00.933Z
2023-07-10T11:02:00.933Z
Cristian Riccio
<p>This was part of a custom liftover function which handled multiple alleles and reverse compliments. I do not know if it will be useful. I believe there was a separate thing somewhere with lookup that I may add. It handled apparent non-1:1 mappings in addition to platform-specific &quot;rsIDs.&quot;</p>
<p>converter = get_lifter(current, target)</p>
<pre><code>def convert_coordinates(input_file, output_file, converter):
        line_count = 0
        skipped = 0
        strand_flips = 0
        with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
            for line in infile:
                # Skip headers
                if line.startswith(&quot;#&quot;):
                    outfile.write(line)
                    continue

                # VCF files are tab-separated
                parts = line.strip().split('\t')

                # Get the chromosome and position
                chrom = parts[0]
                pos = int(parts[1])

                # Convert the coordinates
                try:
                    new_coords = converter.query(chrom, pos)
                except KeyError:
                    #print(f&quot;Chromosome {chrom} not found in the liftover chain file.&quot;)
                    skipped += 1
                    continue

                # If there's no new coordinates, we skip this variant
                if not new_coords:
                    skipped += 1
                    continue

                # Get the new chromosome, position and strand
                new_chrom, new_pos, strand = new_coords[0]

                # If the new strand is negative, reverse complement the REF and ALT alleles
                if strand == '-':
                    parts[3] = reverse_complement(parts[3])
                    parts[4] = reverse_complement(parts[4])
                    strand_flips += 1

                # Replace the old chromosome and position with the new ones
                parts[0] = new_chrom
                parts[1] = str(new_pos)

                # Write the line back into the file
                outfile.write('\t'.join(parts))
                outfile.write('\n')

                # Print a progress statement every 1000 lines
                line_count += 1
                if line_count % 100000 == 0:
                    print(f&quot;Processed {line_count} lines&quot;)

        # Calculate and print the percentage of variants skipped
        percentage_skipped = (skipped / line_count) * 100
        percentage_flips = (strand_flips / line_count) * 100
        print(f&quot;\nFinished processing. Out of {line_count} total lines:&quot;)
        print(f&quot;{percentage_skipped}% of variants were skipped due to missing chromosomal data.&quot;)
        print(f&quot;{percentage_flips}% of variants required a strand flip.&quot;)
</code></pre>

CC BY-SA 4.0
2023-10-07T05:17:27.253Z
2023-10-07T05:39:56.463Z
2023-10-07T05:39:56.463Z
BigMistake
<p>I have a list of about 100,000 rsids and I want to get their genomic coordinates on the GRCh38 genome build. Is there a command line tool that allows me to do this? If yes, which one?</p>
<p>I have tried bcftools but, given the error message I got, I believe that it does not support search by rsid:</p>
<pre><code>bcftools view -R rsids_chr1_subset.txt -f PASS -v snps dbSNP_v156_GRCh38.p14.vcf.gz &gt; rsids_chr1.vcf
[E::bcf_sr_regions_init] Could not parse 1-th line of file rsids_chr1_subset.txt, using the columns 1,2[,-1]
Failed to read the regions: rsids_chr1_subset.txt
</code></pre>

CC BY-SA 4.0
2023-10-07T05:13:04.103Z
The tricky thing is doing this fast. If I remember correctly, I have successfully used an implementation written in C a long time ago that did just this. It also had a Python version. The name escapes me at the moment. However, this just reminded me that I ended up writing my own implementation because there was an issue. I will make an answer.
BigMistake
<bcftools><rsid>
How can I retrieve the GRCh38 coordinates of a list of rsids?
# How can I retrieve the GRCh38 coordinates of a list of rsids?

(Asked by: username_0 on 2023-07-10)

I have a list of about 100,000 rsids and I want to get their genomic coordinates on the GRCh38 genome build. Is there a command line tool that allows me to do this? If yes, which one?
I have tried bcftools but, given the error message I got, I believe that it does not support search by rsid:
bcftools view -R rsids_chr1_subset.txt -f PASS -v snps dbSNP_v156_GRCh38.p14.vcf.gz > rsids_chr1.vcf
[E::bcf_sr_regions_init] Could not parse 1-th line of file rsids_chr1_subset.txt, using the columns 1,2[,-1]
Failed to read the regions: rsids_chr1_subset.txt

**Tags:** bcftools, rsid
> **Comments:**
> **username_1**: The tricky thing is doing this fast. If I remember correctly, I have successfully used an implementation written in C a long time ago that did just this. It also had a Python version. The name escapes me at the moment. However, this just reminded me that I ended up writing my own implementation because there was an issue. I will make an answer.  

---

## 3 answers

### Answer by username_2 on 2023-07-10. Score: 5

Several years ago I created rsidx to address this particular use case. https://github.com/bioforensics/rsidx
# VCF should be sorted by genomic coordinates and indexed by tabix
rsidx index dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx
rsidx search dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx rs3114908 rs10756819


---

### Answer by username_0 on 2023-07-10. Score: 2

rsids_chr1.txt contains a list of rsids, one per line.
bcftools view --include ID==@rsids_chr1.txt dbSNP_v156_GRCh38.p14.vcf.gz -Ou | bcftools query -f '%CHROM  %ID %POS  %REF  %ALT{0}\n' > dbSNP_v156_GRCh38.p14_AF.vcf


---

### Answer by username_1 on 2023-10-07. Score: 1

This was part of a custom liftover function which handled multiple alleles and reverse compliments. I do not know if it will be useful. I believe there was a separate thing somewhere with lookup that I may add. It handled apparent non-1:1 mappings in addition to platform-specific "rsIDs."
converter = get_lifter(current, target)
def convert_coordinates(input_file, output_file, converter):
        line_count = 0
        skipped = 0
        strand_flips = 0
        with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
            for line in infile:
                # Skip headers
                if line.startswith("#"):
                    outfile.write(line)
                    continue

                # VCF files are tab-separated
                parts = line.strip().split('\t')

                # Get the chromosome and position
                chrom = parts[0]
                pos = int(parts[1])

                # Convert the coordinates
                try:
                    new_coords = converter.query(chrom, pos)
                except KeyError:
                    #print(f"Chromosome {chrom} not found in the liftover chain file.")
                    skipped += 1
                    continue

                # If there's no new coordinates, we skip this variant
                if not new_coords:
                    skipped += 1
                    continue

                # Get the new chromosome, position and strand
                new_chrom, new_pos, strand = new_coords[0]

                # If the new strand is negative, reverse complement the REF and ALT alleles
                if strand == '-':
                    parts[3] = reverse_complement(parts[3])
                    parts[4] = reverse_complement(parts[4])
                    strand_flips += 1

                # Replace the old chromosome and position with the new ones
                parts[0] = new_chrom
                parts[1] = str(new_pos)

                # Write the line back into the file
                outfile.write('\t'.join(parts))
                outfile.write('\n')

                # Print a progress statement every 1000 lines
                line_count += 1
                if line_count % 100000 == 0:
                    print(f"Processed {line_count} lines")

        # Calculate and print the percentage of variants skipped
        percentage_skipped = (skipped / line_count) * 100
        percentage_flips = (strand_flips / line_count) * 100
        print(f"\nFinished processing. Out of {line_count} total lines:")
        print(f"{percentage_skipped}% of variants were skipped due to missing chromosomal data.")
        print(f"{percentage_flips}% of variants required a strand flip.")


---

