Document how to get hg38 12-column BED file from UCSC

[sreichl]

• consider generating it from GFF/GTF file or similar that was already used in an upstream module e.g. rnaseq-pipeline or atacseq_pipline ie something simpler and potentially more correct as it comes from the same source/version as the read mapping ( or that can be easily retrieved.). (Generate 12-column BED file based on GFF from RNA pipeline. )

Instructions to get hg38 12-column BED file from UCSC

Okay, here are the complete step-by-step instructions to download the NCBI RefSeq annotations for the hg38 genome assembly from UCSC and generate a 12-column, gzipped BED file (.bed.gz) where the 4th column contains the gene symbol.

This process involves downloading the data in genePred format, rearranging columns using awk so the gene symbol is first, and then converting to BED format using the genePredToBed utility.

• Prerequisites:*
• Web browser

• Command-line terminal (like Terminal on macOS/Linux, or WSL/Git Bash/Cygwin on Windows)
• Command-line tools: gunzip, gzip, awk (usually pre-installed on Linux/macOS)
• The genePredToBed utility from UCSC (you'll download this if you don't have it).

• Instructions:*

1. Go to the UCSC Table Browser:
   Open your web browser and navigate to: https://genome.ucsc.edu/cgi-bin/hgTables

2. Set Genome and Assembly:

   • clade: Mammal
   • genome: Human
   • assembly: Dec. 2013 (GRCh38/hg38)

3. Select Track and Table:

   • group: Genes and Gene Predictions
   • track: NCBI RefSeq
   • table: ncbiRefSeq

4. Choose Output Format:

   • output format: Select all fields from selected table (This will output data in a format compatible with genePred).

5. Specify Output File:

   • In the output file field, enter a filename. Let's use: hg38_refseq_raw.gp.gz
   • Output field separator: tsv (tab-separated)
   • Make sure the gzip compressed checkbox is checked.

6. Download the Data:

   • Click the get output button. The file hg38_refseq_raw.gp.gz will be downloaded.

7. Open Your Command-Line Terminal:

   • Navigate to the directory where you downloaded the file.

8. Unzip the Downloaded File:

   gunzip hg38_refseq_raw.gp.gz

   • You should now have a file named hg38_refseq_raw.gp.

9. Verify Column Order (Important!):

   • The ncbiRefSeq table usually has the transcript ID (e.g., NM_...) in column 1 and the gene symbol (e.g., TP53) in column 12. Let's check the first line to confirm:

     head -n 1 hg38_refseq_raw.gp | awk '{print "Col2:", $2, " Col13:", $13}'

   • Confirm that Column 2 looks like a transcript ID and Column 13 looks like a gene symbol. If they are in different columns, adjust the column numbers ($2, $13) in the next step accordingly.

10. Swap Columns to Put Gene Symbol First and remove header:

    • Use awk to create a new file where the original 13th column (gene symbol) is moved to the 2nd column, and the original 2nd column (transcript ID) is moved to the 13th column.

awk 'BEGIN{FS=OFS="\t"} /^#/ {next} {print $13, $3, $4, $5, $6, $7, $8, $9, $10, $11}' hg38_refseq_raw.gp > hg38_refseq_genePred_for_conversion.gp
```
* This reads hg38_refseq_raw.gp, swaps columns 2 and 13 for each line using a tab (`\t`) as the separator, and writes the result to `hg38_refseq_swapped_noheader.gp`.

11. Get the genePredToBed Utility:
    • If you don't already have the UCSC Kent command-line utilities installed:
      • Go to: http://hgdownload.soe.ucsc.edu/admin/exe/
      • Navigate into the directory corresponding to your operating system (e.g., linux.x86_64/, macOSX.x86_64/).
      • Download the genePredToBed executable file.
      • Make the file executable: chmod +x genePredToBed
      • You can either place this file in your current directory (and run it using ./genePredToBed) or move it to a directory in your system's PATH.

or using conda

• conda create -n genepredtobed
• conda install bioconda::ucsc-genepredtobed

2. Create the Correctly Formatted Input for genePredToBed:
   • Run this awk command. It selects the necessary columns from hg38_refseq_raw.gp and outputs them in the order genePredToBed expects, with the gene symbol ($2 from the original file) as the first column.

   awk 'BEGIN{FS=OFS="\t"} /^#/ {next} {print $13, $3, $4, $5, $6, $7, $8, $9, $10, $11}' hg38_refseq_raw.gp > hg38_refseq_genePred_for_conversion.gp

   • Explanation:
     • BEGIN{FS=OFS="\t"}: Use tabs as separators.
     • /^#/ {next}: Skip the header line.
     • {print $13, $3, $4, $5, $6, $7, $8, $9, $10, $11}: Print the original columns:
       • $13 (name2 -> gene symbol) -> becomes input col 1 (for BED name)
       • $3 (chrom) -> becomes input col 2 (for BED chrom)
       • $4 (strand) -> becomes input col 3 (for BED strand)
       • $5 (txStart) -> becomes input col 4 (for BED chromStart)
       • $6 (txEnd) -> becomes input col 5 (for BED chromEnd)
       • $7 (cdsStart) -> becomes input col 6 (for BED thickStart)
       • $8 (cdsEnd) -> becomes input col 7 (for BED thickEnd)
       • $9 (exonCount) -> becomes input col 8 (for BED blockCount)
       • $10 (exonStarts) -> becomes input col 9 (for BED blockStarts calculation)
       • $11 (exonEnds) -> becomes input col 10 (for BED blockSizes calculation)
   • This creates a new file hg38_refseq_genePred_for_conversion.gp specifically formatted for genePredToBed.

Context: genePredToBed expects its input file to have a specific 10-column structure derived from the standard genePred format:

1. name (This will become the 4th column / name field in the output BED file)
2. chrom
3. strand
4. txStart
5. txEnd
6. cdsStart
7. cdsEnd
8. exonCount
9. exonStarts (comma-separated list)
10. exonEnds (comma-separated list)

Your current hg38_refseq_swapped_noheader.gp file starts with 0 XM_011541469.2 chr1 - ..., which doesn't match this required input structure. The 0 is the bin field, and the XM_... is the transcript ID. Crucially, the first column is not the gene symbol, which is why the previous attempt failed and why genePredToBed is likely encountering unexpected data (like the - for strand in a position where it might expect a coordinate if the columns were misinterpreted).

3. Run genePredToBed:

   • Now, run the conversion tool on this new, correctly formatted file:

   ./genePredToBed hg38_refseq_genePred_for_conversion.gp hg38_refseq_geneSymbol.bed

   • This should now correctly interpret the columns and generate the BED12 file. The first column of the input ($13 from the original file, the gene symbol) will be placed into the 4th column (name field) of hg38_refseq_geneSymbol.bed.

4. Compress the Final BED File:

   gzip hg38_refseq_geneSymbol.bed

Result:*

You will now have a file named hg38_refseq_geneSymbol.bed.gz. This is the 12-column BED annotation file for hg38 NCBI RefSeq, compressed with gzip, and with the gene symbol correctly placed in the 4th column.

[sreichl]

from gtracks manual

--genes <{path/to/genes.bed.gz,GRCh37,GRCh38,hg19,hg38,Sp9512}>
compressed 6-column BED file or 12-column BED12 file containing gene annotations. Alternatively,
providing a genome identifier will use one of the included gene tracks. (default: GRCh37)

[sreichl]

Gemini Code

#!/bin/bash

# --- Configuration ---
# Input GTF file from Ensembl
GTF_FILE="Homo_sapiens.GRCh38.112.gtf"
# Prefix for output files
OUTPUT_PREFIX="Homo_sapiens.GRCh38.112.genes"

# --- Workflow ---
echo "Step 1: Converting GTF to genePred..."
gtfToGenePred \
  -geneNameAsName \
  -ignoreGroupsWithoutExons \
  "${GTF_FILE}" \
  "${OUTPUT_PREFIX}.genePred"

echo "Step 2: Converting genePred to BED12..."
genePredToBed \
  "${OUTPUT_PREFIX}.genePred" \
  "${OUTPUT_PREFIX}.tmp.bed12"

echo "Step 3: Sorting, compressing, and indexing..."
# Sort the BED file
sort -k1,1 -k2,2n "${OUTPUT_PREFIX}.tmp.bed12" > "${OUTPUT_PREFIX}.bed12"

# Compress with bgzip
bgzip "${OUTPUT_PREFIX}.bed12"

# Index with tabix
tabix -p bed "${OUTPUT_PREFIX}.bed12.gz"

# --- Cleanup ---
echo "Cleaning up intermediate files..."
rm "${OUTPUT_PREFIX}.genePred" "${OUTPUT_PREFIX}.tmp.bed12"

echo "Done! Your final file is ${OUTPUT_PREFIX}.bed12.gz"

[sreichl]

ChatGPT Code

#!/usr/bin/env bash
set -euo pipefail

# Usage check
if [[ $# -ne 2 ]]; then
  echo "Usage: $0 <input.gtf[.gz]> <output.bed12>"
  exit 1
fi

INPUT_GTF=$1
OUTPUT_BED=$2

# temp files
TMP_GTF=$(mktemp --suffix=.gtf)
TMP_GP=$(mktemp)
TMP_SWAP=$(mktemp)
trap 'rm -f "$TMP_GTF" "$TMP_GP" "$TMP_SWAP"' EXIT

# 1. Decompress GTF (if needed)
if [[ $INPUT_GTF == *.gz ]]; then
  zcat "$INPUT_GTF" > "$TMP_GTF"
else
  cp "$INPUT_GTF" "$TMP_GTF"
fi

# 2. Convert GTF → genePredExt, use gene_name as name2
gtfToGenePred \
  -genePredExt \
  -geneNameAsName2 \
  "$TMP_GTF" \
  "$TMP_GP"

# 3. Swap col1 (transcript ID) ↔ col12 (gene symbol)
awk 'BEGIN{FS=OFS="\t"} {t=$1; $1=$12; $12=t; print}' \
  "$TMP_GP" > "$TMP_SWAP"

# 4. Convert genePredExt → BED12 (col1 now gene symbol)
genePredToBed \
  "$TMP_SWAP" \
  "$OUTPUT_BED"

echo "→ BED12 file written to $OUTPUT_BED"

[sreichl]

conda install -c bioconda ucsc-bed
or
conda install -c bioconda ucsc-gtftogenepred ucsc-genepredtobed

better check docs

[sreichl]

# Download the GTF file (check Ensembl for the latest release number)
wget https://ftp.ensembl.org/pub/release-112/gtf/homo_sapiens/Homo_sapiens.GRCh38.112.gtf.gz
wget ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz 

# Unzip it
gunzip Homo_sapiens.GRCh38.112.gtf.gz
gunzip Homo_sapiens.GRCh38.110.gtf.gz

[sreichl]

Think about downloading the GTF file and creating the needed BED file within the pipeline. This could be an extra feature, but it's not required. If it's not provided, we need a genome and version to download from Ensembl. The code is part of RNA sequencing analysis, with one step converting GTF to BED, which I already have a draft for. This would be a major update, making it more reproducible and needing fewer resources and preparations.

Also, when creating the BED file, consider making it only for the required genes. This makes it smaller and more efficient but less versatile. So, think carefully about this option. It might not be a good idea.

[bednarsky]

I think this issue is important - a simple way of getting this kept me from using this module, since it's "only a sanity check anyways"

[sreichl]

Yes, agreed. Hence, its priority is set to high.