Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-10T13:14:09.930Z
Answers.LastActivityDate: 2023-07-10T13:14:09.930Z
Answers.UserDisplayName: Daniel Standage
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-10T11:02:00.933Z
Answers.LastActivityDate: 2023-07-10T11:02:00.933Z
Answers.UserDisplayName: Cristian Riccio
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-07T05:17:27.253Z
Answers.LastActivityDate: 2023-10-07T05:39:56.463Z
Answers.LastEditDate: 2023-10-07T05:39:56.463Z
Answers.UserDisplayName: BigMistake
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-07T05:13:04.103Z
Comments.Text: 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.
Comments.UserDisplayName: BigMistake
Tags: <bcftools><rsid>
Title: How can I retrieve the GRCh38 coordinates of a list of rsids?
ThreadText: # 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.")
  
  ```
  
  
  ---
  
  