How to Work With Your tellmeGen Raw Data

Woolf Software
genotyping raw-data vcf bioinformatics imputation
A huge regular grid of glowing pods set in wet tidal sand at night, with gaps being filled by thin creeping filaments.

By the end of this guide you will have taken the .zip that tellmeGen provides and turned it into something you can analyse with confidence. Along the way you will confirm which reference genome build the coordinates refer to. You will convert the file into a valid VCF with a correct reference allele at every site, where VCF is the Variant Call Format used as the standard text format for genomic variants. You will then measure call rate and heterozygosity, check that the genetic sex matches your expectation, and annotate the calls against dbSNP and ClinVar. If you have a reason to, you can also produce imputed dosages and polygenic scores.

You will need the raw data export and roughly 20 GB of free disk. On the software side you need bcftools 1.17 or newer, plink 1.9 and plink2, and samtools for indexing FASTA files. You also need reference genomes for both GRCh37 and GRCh38 with their .fai indexes. Everything below runs on a laptop in under an hour, with the exception of imputation, which runs on a remote server.

1. Extract the export and read the header before anything else

The first step is simply to get the file onto disk and look at it, because a few minutes of reading the header will save you from mistakes that are hard to detect later. tellmeGen delivers raw data from the account area under Settings, then Kits, then Download Raw Data. What arrives is a zip containing a single plain text file. Resist the temptation to open it in a spreadsheet: Excel will silently convert marker identifiers and truncate the file at a million rows.

unzip -l tellmegen_raw.zip
unzip tellmegen_raw.zip -d tmg/
head -30 tmg/*.txt

The header is a block of comment lines beginning with #, and it is worth reading carefully. It usually names the array platform, the genome build the positions are reported against, and the orientation convention for alleles. Below the header you should find four tab-separated columns in the layout used by most consumer exports: marker identifier (mostly rs numbers, with some chip-internal identifiers), chromosome, position, and a two-character genotype such as AG. Verify the column order rather than assuming it, since some vendors emit Allele1 and Allele2 as separate columns.

With the layout confirmed, count what you have.

# number of markers
grep -vc '^#' tmg/raw_data.txt

# markers per chromosome
awk -F'\t' '!/^#/{print $2}' tmg/raw_data.txt | sort | uniq -c | sort -k2,2V

# genotype value distribution, including no-calls
awk -F'\t' '!/^#/{c[$4]++} END{for(g in c) print c[g], g}' tmg/raw_data.txt | sort -rn | head -20

A current Illumina Global Screening Array export lands in the high hundreds of thousands of markers. The genotype distribution tells you two useful things: the no-call encoding, typically --, NN, or 00, and whether the file contains insertion/deletion pseudo-alleles written as I and D. Those indel codes carry no sequence and cannot be represented faithfully in VCF, so we drop them and note how many were lost.

Two framing points are worth holding in mind before you go further. The first concerns scope. A genotyping array interrogates a fixed and pre-selected panel of sites, on the order of a million out of roughly three billion base pairs. The file therefore describes well under a tenth of a percent of your genome and contains nothing at all about the positions it does not probe. Whole-genome approaches were built precisely to remove that selection step 1.

The second point is more encouraging. Reproducibility on common polymorphisms is the strong suit of array genotyping, so the calls in this file are generally consistent between runs and between vendors using the same chip. Where results differ between consumer companies, the cause is usually the interpretation layer applied on top rather than the underlying genotypes.

2. Determine the genome build empirically

Every downstream step depends on knowing which reference build the coordinates refer to, and header text about builds is sometimes wrong or absent. Rather than trust it, test it. The reliable approach is to convert against both references and compare how many sites survive reference checking, because bcftools convert --tsv2vcf looks up each position in the FASTA and discards records where neither observed allele matches the reference base.

# strip indel codes and no-calls, keep the four columns tsv2vcf expects
awk -F'\t' 'BEGIN{OFS="\t"} !/^#/ && $4 !~ /[IDND-]/ && length($4)==2 {print $1,$2,$3,$4}' \
  tmg/raw_data.txt > tmg/clean.tsv

bcftools convert --tsv2vcf tmg/clean.tsv -f GRCh37.fa -s TMG01 -Oz -o tmg/b37.vcf.gz 2> b37.log
bcftools convert --tsv2vcf tmg/clean.tsv -f GRCh38.fa -s TMG01 -Oz -o tmg/b38.vcf.gz 2> b38.log

bcftools index -t tmg/b37.vcf.gz && bcftools index -t tmg/b38.vcf.gz
bcftools index -n tmg/b37.vcf.gz; bcftools index -n tmg/b38.vcf.gz

The result is easy to read. The correct build retains almost every site, with only a small residue of losses from strand-ambiguous or probe-specific markers, while the wrong build retains a fraction consistent with chance matching of a single base. The difference is unmistakable. Consumer arrays have historically reported GRCh37 coordinates, but confirm it yourself and record the answer in a README next to the file.

One related detail is worth checking while you are here: chromosome naming. If your FASTA uses chr1 and the raw file uses 1, --tsv2vcf will fail on every record. Rename with bcftools annotate --rename-chrs, supplying a two-column mapping file, rather than rewriting the FASTA.

3. Normalise, fix ploidy, and sort

Now that the coordinates are right, the alleles need attention. Arrays report alleles relative to the probe design, and a minority of markers will appear on the opposite strand to the reference. bcftools +fixref identifies and repairs these when you supply a reference variant set such as the dbSNP VCF.

bcftools +fixref tmg/b37.vcf.gz -- -f GRCh37.fa -i dbsnp_b151_GRCh37.vcf.gz \
  | bcftools sort -Oz -o tmg/fixed.vcf.gz
bcftools index -t tmg/fixed.vcf.gz
bcftools +fixref tmg/fixed.vcf.gz -- -f GRCh37.fa   # re-run in report mode to confirm

The report mode output should show nearly all sites as ref-match and very few as flipped or unresolved. The remaining task in this step is ploidy. The export writes male X and Y genotypes as single letters or as homozygous diploid calls, and the mitochondrial genome is haploid. Use bcftools +fixploidy with a ploidy definition file appropriate to the sample, and inspect X heterozygosity in the next step rather than trusting the vendor’s sex assignment.

4. Run quality control

Quality control on a single sample comes down to three statistics, and they are quick to produce. Convert to PLINK binary format and compute call rate, autosomal heterozygosity, and X chromosome heterozygosity.

plink --vcf tmg/fixed.vcf.gz --make-bed --out tmg/tmg
plink --bfile tmg/tmg --missing --het --out tmg/qc
plink --bfile tmg/tmg --split-x b37 --make-bed --out tmg/tmg.sx
plink --bfile tmg/tmg.sx --check-sex --out tmg/sexcheck

Each statistic answers a different question. A sample call rate below roughly 98 percent suggests a degraded sample or a failed chip section, and you should look at missingness by chromosome to see whether the loss is uniform or clustered. Autosomal heterozygosity far above the population expectation is the signature of sample contamination or a mixture of two DNA sources. In --check-sex, the inbreeding coefficient F on the non-pseudoautosomal X separates cleanly: values near 1 indicate one X, values near 0 indicate two. An intermediate value is not a finding to interpret on your own, and a result that contradicts your expectation should go to a clinician for confirmatory testing in an accredited laboratory rather than into a blog post.

One further check earns its keep. Merge your genotypes with a public reference panel on the intersecting markers, prune for linkage disequilibrium, and project onto reference principal components with plink2 --pca and --score using the reference loadings. Estimating ancestry this way costs twenty minutes and catches sample swaps that call rate will not.

5. Annotate the calls

Positions alone mean little. The next step is to attach identifiers, allele frequencies, and clinical assertions from curated databases. Keep each source in its own INFO field so that you can always trace where a particular claim came from.

bcftools annotate -a clinvar_GRCh37.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNDN,INFO/CLNREVSTAT \
  -Oz -o tmg/annotated.vcf.gz tmg/fixed.vcf.gz
bcftools index -t tmg/annotated.vcf.gz

bcftools view -i 'INFO/CLNSIG ~ "Pathogenic"' tmg/annotated.vcf.gz \
  | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%INFO/CLNSIG\t%INFO/CLNREVSTAT[\t%GT]\n'

For consequence prediction on transcripts, Ensembl VEP with a local cache is the tool we would use, invoked as vep --cache --offline --assembly GRCh37 --everything --vcf. Bear in mind that transcript models and consequence calls are themselves computational inferences over the genome sequence, an inference problem with a long methodological history and real residual uncertainty, particularly outside well-characterised coding exons 2.

Read the ClinVar hits with discipline. Filter hard on review status, keeping only assertions with multiple submitters and no conflicts, and treat everything else as a lead rather than a result. There is a structural reason for this caution: arrays place a fixed probe at a rare variant site, and rare variants are exactly where probe performance is weakest, because cluster separation is trained on very few carriers. A rare pathogenic call from a consumer array is therefore a hypothesis that requires targeted sequencing in a clinical laboratory, and the interpretation of any confirmed result belongs with a genetic counsellor or physician. What this guide describes is measurement, and diagnosis is a separate undertaking.

6. Impute, if you have a reason to

Imputation uses a reference haplotype panel to infer genotypes at sites the chip did not measure, taking the typed markers to tens of millions of positions with per-variant confidence scores. It is worthwhile if you intend to compute polygenic scores, because score weight files reference many variants absent from any single array. Public servers such as TOPMed or the Michigan Imputation Server expect per-chromosome, sorted, bgzipped VCFs on GRCh38. Lift over first.

picard LiftoverVcf I=tmg/fixed.vcf.gz O=tmg/hg38.vcf.gz \
  CHAIN=hg19ToHg38.over.chain.gz REJECT=tmg/rejected.vcf R=GRCh38.fa

for c in $(seq 1 22); do
  bcftools view -r chr${c} tmg/hg38.vcf.gz -Oz -o tmg/chr${c}.vcf.gz
done

Expect a few thousand markers in rejected.vcf, mostly in segmental duplications and other repetitive sequence where the two assemblies disagree. This is a general and well-characterised difficulty rather than a flaw in your file 3. After imputation returns, filter on the estimated imputation quality and discard variants with R² below 0.8 for any quantitative use. Remember that imputation accuracy degrades for individuals whose ancestry is poorly represented in the reference panel.

Polygenic scores then follow from PGS Catalog harmonised weight files:

plink2 --vcf chr1.dose.vcf.gz dosage=DS \
  --score PGS000123_hmPOS_GRCh38.txt 3 4 6 header cols=+scoresums list-variants \
  --out pgs_chr1

Adjust the three column indices to match the identifier, effect allele, and weight columns of your specific file. Then check list-variants output to see how many weights were matched. That number matters, because a score built from 30 percent of the intended variants is not comparable to the published distribution. A raw score is also meaningless without the ancestry-matched reference distribution it was standardised against.

7. Know what the file cannot contain

The most common analytical error is treating an array export as a genome, so it is worth being explicit about what is missing. Three categories are simply absent, each for a different reason.

  • Rare and private variants are absent because the chip was never designed to probe them.
  • Copy number variants are absent because detecting them requires the raw probe intensity values (log R ratio and B allele frequency), which the text export does not include.
  • Structural rearrangements, large repeat expansions, and phased haplotypes across long distances are absent. Resolving them requires reads long enough to span the region, which is the specific capability third-generation single-molecule platforms were developed to provide 4 5.

Short-read approaches mitigate but do not eliminate the ambiguity in repetitive regions, where mapping becomes underdetermined 3. The practical upshot is that a negative result in this file is not evidence of absence, for anything.

If you want coverage of those categories, the measurement has to change rather than the analysis. The choice among sequencing platforms is a real tradeoff among read length, per-base accuracy, and cost 6 7.

Common problems

A handful of failure modes account for most of the trouble people run into, and each has a specific cause and remedy.

bcftools convert --tsv2vcf drops nearly every record. Either the build is wrong, which step 2 resolves, or the chromosome names do not match the FASTA. In the second case bcftools annotate --rename-chrs resolves it. Check both before suspecting the file itself.

The X chromosome shows heterozygous calls in a male sample. Some of these are genuine pseudoautosomal region markers, which is why --split-x exists. If they persist outside the pseudoautosomal boundaries at appreciable rate, suspect contamination and look at autosomal heterozygosity for corroboration.

Genotypes appear reversed relative to another source, for instance CC where you expected GG. This is a strand orientation difference rather than an error, and bcftools +fixref against dbSNP resolves it. Strand-ambiguous A/T and C/G markers cannot be resolved this way and are safest excluded from any cross-dataset comparison.

Marker counts differ from another vendor’s export of the same person. Different chips probe different panels, so overlap between two consumer arrays is typically partial. Compare genotypes only on the intersecting rsIDs, after fixing strand, and expect concordance well above 99 percent on the shared common markers.

A ClinVar pathogenic hit appears in the output. Do not act on it and do not treat it as a result. Confirm the specific variant by an orthogonal method in an accredited laboratory, and take the confirmed result to a clinician who can interpret it in the context of your family and clinical history.

Finally, a word on the upload services that invite you to hand over this file for free reinterpretation. Your genotype file is a permanent identifier that also describes your relatives, and once it is copied you cannot withdraw it. Everything in this guide runs locally on your own machine, which we think is the correct default.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. K G Skryabin, E B Prokhortchouk, A M Mazur, et al. Combining Two Technologies for Full Genome Sequencing of Human. Acta Naturae, 2009. https://doi.org/10.32607/20758251-2009-1-3-102-107 ↩

  2. Christopher B Burge, Samuel Karlin. Finding the genes in genomic DNA. Current Opinion in Structural Biology, 1998. https://doi.org/10.1016/s0959-440x(98)80069-9 ↩

  3. Todd J. Treangen, Steven L. Salzberg. Repetitive DNA and next-generation sequencing: computational challenges and solutions. Nature Reviews Genetics, 2011. https://doi.org/10.1038/nrg3117 ↩ ↩2

  4. E. E. Schadt, S. Turner, A. Kasarskis. A window into third generation sequencing. Human Molecular Genetics, 2010. https://doi.org/10.1093/hmg/ddq481 ↩

  5. Mohit K. Midha, Mengchu Wu, Kuo-Ping Chiu. Long-read sequencing in deciphering human genetics to a greater depth. Human Genetics, 2019. https://doi.org/10.1007/s00439-019-02064-y ↩

  6. Michael L. Metzker. Emerging technologies in DNA sequencing. Genome Research, 2005. https://doi.org/10.1101/gr.3770505 ↩

  7. Akash Gautam. Applications of DNA Sequencing Technologies for Current Research. Techniques in Life Science and Biomedicine for the Non-Expert, 2022. https://doi.org/10.1007/978-3-030-94230-4_23 ↩