\subsection*{Sample selection}
Bulk RNA-seq from a curated set of 1,139 flash frozen colorectal tissue samples was re-processed and analyzed for the presence of unannotated HERV transcripts \citep{Dampier2020}.
For discovery of associations between tissue phenotypes and HERV expression, paired-end sequencing from 833 colon biopsies was used.
Of those biopsies, 462 were from healthy, 62 were from tumor adjacent normal, and 309 were from tumor tissue.
This set of samples is referred to as the discovery cohort.

For validation of discovered associations, two cohorts of samples from independent subjects were used.
The first (independent cohort \#1) included single-end sequencing from 276 colon biopsies of which 12 were healthy, 56 were tumor adjacent, and 208 were tumor.
The second (independent cohort \#2) included paired-end sequencing from 30 pairs of tumor and matched adjacent normal biopsies (i.e. 15 samples of both tissue phenotypes paired by subject).
To permit valid inference from regular linear models, only one tissue phenotype was used per subject in the discovery cohort and independent cohort \#1 (i.e. all samples were from independent subjects).

Samples were originally gathered from publicly available RNA-seq datasets hosted on the Genomic Data Commons and Sequence Read Archive as well as from a recently generated dataset of healthy mucosal biopsies obtained during screening colonoscopies at the Catalan Institute of Oncology in Barcelona, Spain \citep{Dampier2020, DiezObrero2020}.
A complete description of sample selection and cohort construction for the curated dataset is available in Dampier et al. (2020) \citep{Dampier2020}, and preliminary bioinformatic processing is described in the Supplementary Methods.

Although race was not analyzed in this study, we report that the discovery cohort included 37 Asian, 83 Black, and 660 White subjects, while the race of 53 others was not available.
Independent cohort \#1 included 12 Asian, 4 Black, and 17 White subjects, while the race of 243 others was not available.
Independent cohort \#2 included 7 Asian, 1 Black, and 7 White subjects.
Although no human subjects were recruited during this study, all previous studies from which data was used followed experimental protocols approved by institutional review committees and met the guidelines of relevant governmental agencies.

\subsection*{Human endogenous retrovirus annotation}
HERVs were identified using the genomic coordinates (GRCh38) and feature names of spliced transcripts of HERV elements from a recently published human cancer transcriptome annotated with the GENCODE v24 basic gene model \citep{Frankish2018} and a model of genome-wide repeat elements \citep{Attig2019}.
The cancer transcriptome was assembled from RNA-seq reads of 24 subjects from each of 32 cancer types obtained from The Cancer Genome Atlas (TCGA), as previously described \citep{Attig2019}.
Names and coordinates of repeat elements were generated using RepeatMasker \citep{Smit2015} configured with nhmmer \citep{Wheeler2013} in sensitive mode using the Dfam 2.0 library (v150923) \citep{Hubley2015}, as previously described \citep{Attig2017}.
The method relies on profile hidden Markov models representing known human repeat families to classify repeat elements throughout the genome.
Previously unannotated HERV transcripts were selected for the current study in multiple steps described in the Supplementary Methods.

\subsection*{Transcript quantification}
Transcript quantification was performed with Salmon v1.2.1 \citep{Patro2017} in mapping-based mode with the \verb|--gcBias| and \verb|--validateMappings| flags set.
Inputs to the Salmon \verb|quant| command were RNA-seq reads in FASTQ format and a custom Salmon index including unannotated HERV transcripts.
Creation of the Salmon index is described in the Supplementary Methods.

\subsection*{Proviral locus expression analysis}
Associations between HERV locus expression and tissue phenotype were tested using the \emph{Bioconductor} suite of genomics analysis packages \citep{bioc} in the R statistical programming language \citep{R}.
Data processing and analysis steps are described in detail in the Supplementary Methods.
Briefly, transcript expression estimates were converted to gene counts using \emph{tximport} \citep{Soneson2015}, confounding variables were estimated using \emph{sva} \citep{sva}, and expression levels were modeled using \emph{DESeq2} \citep{Love2014}.

Counts of all HERV loci with sufficient expression were modeled with phenotype and latent factors as explanatory variables.
Due to missing information for some samples, effects of demographic and technical factors were not modeled directly but instead captured with latent factors as described in the Supplementary Methods.
Significant associations between locus counts and tissue phenotype were identified by testing whether the absolute value of the effect of phenotype was at least two fold change for a given locus.
Multiple hypothesis testing correction was performed using a false discovery rate of 5\%.
All figures were produced with \emph{ggplot2} \citep{Wickham2016} as described in the Supplementary Methods.

\subsection*{Biomarker performance assessment}
Expression levels of select HERV loci were used as binary classifiers to distinguish tumor from non-tumor samples.
Classification was performed by logistic regression using the \verb|glm| function in R.
Transcript per million abundance estimates from Salmon were used as predictors of tissue phenotype.
The performance of each classifier was assessed by calculating the area under the receiver operating characteristic curve using the \emph{pROC} \citep{pROC} package.

All code is available on GitHub (https://github.com/dampierch/herv).
