Reference/annotation data

Twitcher Workbench needs the reference genome matching with the one used in the original Twitcher analyses to be able to reconstruct and visualise the TSM events. Gene and repeat annotation files are required for the optional annotation steps.

The large reference files are not stored in Git and are not baked into the Docker image. The repository keeps only the folder structure. Before running the application with bundled resources, place reference genomes, annotation files, and repeat-mask files under the ref/ directory.

The application treats files under ref/ as ready-to-use production resources. The app will not sort, normalize, bgzip, tabix-index, or otherwise rewrite bundled resources at runtime. This keeps the web viewer responsive, makes Docker/HPC usage predictable, and allows ref/ to be mounted read-only.

Important! Chromosome names must match

The reference FASTA, gene annotation GFF3, and repeat-mask GFF3 must use the same chromosome naming convention. Human genome data is provided both in the form 1, 2, … and chr1, chr2, … . Both formats are accepted by Twitcher Workbench, but all the files must confirm to one format.

Expected ref/ structure

The directory names are critical, and the files within them can be named freely as long as their suffixes match (e.g., *.gff3.gz, *.gff3.gz.tbi). Recommended structure for an hg38 setup is:

ref/
├── genomes/
│   └── hg38/
│       ├── hg38.no_alt.fa.gz
│       ├── hg38.no_alt.fa.gz.fai
│       └── hg38.no_alt.fa.gz.gzi
├── annotations/
│   └── hg38/
│       ├── Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz
│       └── Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz.tbi
└── masks/
    └── hg38/
        ├── hg38.no_alt.masking.sorted.gff3.gz
        └── hg38.no_alt.masking.sorted.gff3.gz.tbi

The app scans these folders and shows available bundled resources in dropdown menus. If the folders are empty, the app still works. Users can upload a custom reference genome, gene annotation GFF3, and repeat-mask GFF3 through the UI.


Preparing bundled data

Required tools

Install the following command-line tools before preparing reference resources:

samtools
bgzip
tabix

On Debian/Ubuntu-based systems, samtools, bgzip, and tabix are commonly available from packages such as samtools and tabix.


Reference genome FASTA

Place the FASTA file under:

ref/genomes/<genome_name>/

Example:

ref/genomes/hg38/hg38.no_alt.fa.gz

The FASTA must be indexed:

samtools faidx ref/genomes/hg38/hg38.no_alt.fa.gz

Expected files:

hg38.no_alt.fa.gz
hg38.no_alt.fa.gz.fai
hg38.no_alt.fa.gz.gzi

If .gzi is missing, regenerate the compressed FASTA with bgzip and index again.

Check FASTA chromosome names:

cut -f1 ref/genomes/hg38/hg38.no_alt.fa.gz.fai | head

The output defines the chromosome names that annotation and repeat-mask files must match.

The version of the hg38 genome shown above can be downloaded from the HGSVC2 data release:

wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/HGSVC2/technical/reference/20200513_hg38_NoALT/hg38.no_alt.fa.gz

Gene annotation GFF3

Place gene annotation files under:

ref/annotations/<genome_name>/

Recommended hg38 example:

ref/annotations/hg38/Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz
ref/annotations/hg38/Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz.tbi

The GFF3 file must be:

  • chromosome-name compatible with the FASTA
  • sorted by chromosome and start position
  • compressed with bgzip
  • indexed with tabix -p gff

Ensembl gene annotations

Twitcher Workbench have been tested with Ensembl gene annotations. The latest release of human annotation files are available at the Ensembl FTP site. Annotations for the release 115 can be downloaded with:

wget https://ftp.ensembl.org/pub/release-115/gff3/homo_sapiens/Homo_sapiens.GRCh38.115.chr.gff3.gz

Ensembl uses chromosomes IDs 1, 2, …, while the FASTA file contains chr1, chr2, … . Secondly, the Ensembl gene annotation file is not sorted by position and cannot be indexed.

There are many ways for fixing this, and the commands below provide one way.

zcat Homo_sapiens.GRCh38.115.chr.gff3.gz | \
    awk -F '\t' '$1~/^#/{print $0 > "headers.txt"}$1!~/^#/{print $0}' | \
    sort -k1,1 -k4,4n | awk -F '\t' '{print "chr"$0}' > body.txt
( grep -v "###" headers.txt && cat body.txt ) | \
    bgzip -c > Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz
tabix Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz

After placing the file in the right directory, validate that a known region returns features:

tabix ref/annotations/hg38/Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz \
  chr13:24882135-24885186 | head

If this command returns gene records, JBrowse should also be able to display the gene annotation track in that region.


Repeat-masking GFF3

Place repeat-mask files under:

ref/masks/<genome_name>/

Recommended hg38 example:

ref/masks/hg38/hg38.no_alt.masking.sorted.gff3.gz
ref/masks/hg38/hg38.no_alt.masking.sorted.gff3.gz.tbi

The repeat-mask GFF3 file must also be:

  • chromosome-name compatible with the FASTA
  • sorted by chromosome and start position
  • compressed with bgzip
  • indexed with tabix -p gff

Generation of repeat annotations

There are many methods for repeat annotations. The annotations provided for download were created by combining the outputs of dustmasker, TRF and RepeatMasker. For reproducibility, the commands used to create those files are provided.

Dustmasker was run and annotations were converted to GFF:

dustmasker -in <(zcat hg38.no_alt.fa.gz) -linker 5 -outfmt interval -out hg38.no_alt.dustmasker

cat hg38.no_alt.dustmasker | \
    awk '{OFS="\t";
          if($1~/>chr/){chr=substr($1,2)}
          else{print chr,"dustmasker","low_complexity",$1+1,$3+1,".",".",".","."}
    }' | bgzip -c > hg38.no_alt.dustmasker.gff3.gz

tabix hg38.no_alt.dustmasker.gff3.gz

TRF was run and annotations were converted to GFF:

trf <(zcat hg38.no_alt.fa.gz) 2 7 7 80 10 50 500 -h -d -l 6

cat 63.2.7.7.80.10.50.500.dat | \
    awk 'NR<7{next}
        $1~/Sequence/{seq=$2;next}
        $1~/Parameters/{next}
        length($1)>1{OFS="\t";print seq,"trf","tandem_repeat",$1,$2,$8,".",".",$3";"$4";"$6";"$7";"$13";"$14}' | \
    bedtools sort -i - | bgzip -c > hg38.no_alt.trf.gff3.gz

tabix hg38.no_alt.trf.gff3.gz

RepeatMasker annotations were downloaded and converted to GFF:

wget http://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.fa.out.gz

zcat hg38.fa.out.gz | awk '$5~/chr/{OFS="\t";print $5,"repeatmasker",$11,$6,$7,".",".",".",$10}' | \
    bgzip -c > hg38.no_alt.repeatmasker.gff3.gz

tabix hg38.no_alt.repeatmasker.gff3.gz

The annotations are redundant as each method reports simple repeats. In those cases, TRF was chosen as the primary method, RepeatMasker the secondary method, and dustmasker annotations were only included if neither of the two first methods had reported them.

First, “simple” and “complex” repeats reported by RepeatMasker were separated:

zgrep ")n$" hg38.no_alt.repeatmasker.gff3.gz | bgzip -c > hg38.no_alt.repeatmasker_tr.gff3.gz
zgrep -v ")n$" hg38.no_alt.repeatmasker.gff3.gz | bgzip -c > hg38.no_alt.repeatmasker_other.gff3.gz

tabix hg38.no_alt.repeatmasker_other.gff3.gz

Next, intersecting annotations were discarded:

bedtools intersect -v -f 0.9 -a hg38.no_alt.repeatmasker_tr.gff3.gz -b hg38.no_alt.trf.gff3.gz | \
    bgzip -c > hg38.no_alt.repeatmasker_tr_keep.gff3.gz

tabix hg38.no_alt.repeatmasker_tr_keep.gff3.gz

bedtools intersect -v -f 0.9 -a hg38.no_alt.dustmasker.gff3.gz -b hg38.no_alt.trf.gff3.gz \
    hg38.no_alt.repeatmasker_tr_keep.gff3.gz hg38.no_alt.repeatmasker_other.gff3.gz | \
    bgzip -c > hg38.no_alt.dustmasker_keep.gff3.gz

tabix hg38.no_alt.dustmasker_keep.gff3.gz

bedtools intersect -v -f 0.9 -a hg38.no_alt.trf.gff3.gz -b hg38.no_alt.repeatmasker_other.gff3.gz | \
    bgzip -c > hg38.no_alt.trf_keep.gff3.gz

tabix hg38.no_alt.trf_keep.gff3.gz

Finally, the nonredundant files were combined and sorted:

for i in {1..22} X Y; do echo chr$i; done > chroms.txt

cat chroms.txt | while read c; do
    cat <(tabix hg38.no_alt.repeatmasker_other.gff3.gz $c) <(tabix hg38.no_alt.repeatmasker_tr_keep.gff3.gz $c) \
     <(tabix hg38.no_alt.dustmasker_keep.gff3.gz $c) <(tabix hg38.no_alt.trf_keep.gff3.gz $c) | sort -k1,1 -k4,5n;
done | bgzip -c > hg38.no_alt.masking.sorted.gff3.gz

tabix hg38.no_alt.masking.sorted.gff3.gz

Check that all bundled resources are compatible

Check FASTA contigs:

cut -f1 ref/genomes/hg38/hg38.no_alt.fa.gz.fai | head

Check gene annotation contigs:

tabix -l ref/annotations/hg38/Homo_sapiens.GRCh38.115.ucsc.sorted.gff3.gz | head

Check repeat-mask contigs:

tabix -l ref/masks/hg38/hg38.no_alt.masking.sorted.gff3.gz | head

The names should match. For a chr-style hg38 FASTA, all three should use chr1, chr2, etc.