Pseudogenes

Human gene annotations were acquired from GENCODE v17 (31 (link)). Protein-coding transcripts were defined as those with ‘protein_coding’ gene biotype and ‘protein_coding’ transcript biotype. The lncRNAs transcripts were defined as those with ‘processed_transcript’, ‘lincRNA’, ‘3prime_overlapping_ncrna’, ‘antisense’, ‘non_coding’, ‘sense_intronic’ or ‘sense_overlapping’ gene biotype. Small non-coding RNA (sncRNA) transcripts were defined as those with ‘snRNA’, ‘snoRNA’, ‘rRNA’, ‘Mt_tRNA’, ‘Mt_rRNA’, ‘misc_RNA’ or ‘miRNA’ gene biotype. Pseudogene transcripts were defined as those with ‘polymorphic_pseudogene’, ‘pseudogene’, ‘IG_C_pseudogene’, ‘IG_J_pseudogene’, ‘IG_V_pseudogene’, ‘TR_V_pseudogene’ or ‘TR_J_pseudogene’ gene biotype.
Mouse and Caenorhabditis elegans gene annotations were extracted from Ensembl Gene Release 72 and LiftOver to mm9/mm10 and ce6/ce10, respectively. Protein-coding, lncRNAs, sncRNAs and pseudogenes were classified using a similar method. Human, mouse and C. elegans circRNA annotations were downloaded from circBase v0.1 (6 (link)).
These transcripts were scanned to find conserved miRNAs target sites using miRanda v3.3a with the ‘-strict’ parameter. The target sites that overlap with any entry of the aforementioned AGO CLIP clusters were considered as the CLIP-supported target sites.

The GENCODE gene set is created by merging the results of manual and computational gene annotation methods. Manual gene annotation has two major modes of operation: clone-by-clone and targeted annotation. ‘Clone-by-clone’ annotation involves ‘walking’ across a genomic region, investigating the sequence, aligned expression data and computational predictions for each BAC clone. In doing so, an expert annotator investigates all possible genic features and considers all possible annotations and biotypes simultaneously. We believe this approach carries substantial advantages. For example, the decision to annotate a locus as protein-coding or pseudogenic benefits from being able to weigh both possibilities in light of all available evidence. This process helps prevent false positive and false negative misclassifications. Targeted annotation is designed to answer specific questions such as ‘is there an unannotated protein-coding gene in this position?’ Ranked target lists are generated by computational analysis based, for example, on transcriptomic data, shotgun proteomic data or conservation measures. Over the last two years mouse annotation has been dominated by the clone-by-clone approach while the human genome has been refined entirely via targeted reannotation except for the annotation of human assembly patches and haplotypes released by the Genome Reference Consortium (15 (link)), which take a clone-by-clone approach.
Over the last two years, we have focused on two broad areas: completing the first pass manual annotation across the entire mouse reference genome and a dedicated effort to improve the annotation of protein-coding genes in human and mouse.
We have completed the annotation of novel protein-coding genes, lncRNAs and pseudogenes, plus QC and updating previous annotation where necessary for mouse chromosomes 9, 10, 11, 12, 13, 14, 15, 16 and 17. These updates bring the fraction of the mouse genome with completed first pass manual annotation to approximately 97%. In addition, we have continued to work with the NCBI and Mouse Genome Informatics project at the Jackson Laboratory to resolve annotation differences for protein-coding, pseudogene and lncRNA loci. For protein-coding genes this is under the umbrella of the Consensus Coding Sequence (CCDS) project (16 (link)).
We have also manually investigated unannotated regions of high protein-coding potential identified by whole genome analysis using PhyloCSF (17 (link)) (a tool described in more detail below). In human, this led to the addition of 144 novel protein-coding genes and 271 pseudogenes (of which 42 were unitary pseudogenes). In mouse, we annotated orthologous loci for all but 11 of the 144 human protein-coding genes. We have also revisited the annotation of all olfactory receptor loci in both human and mouse, using RNAseq data to define 5′ and 3′ UTR sequences for ∼1400 loci. In human we have also targeted a ‘deep dive’ manual reannotation of genes on clinical panels for paediatric neurological disorders to identify missing functional alternative splicing. Incorporating second and third generation transcriptomic data, we reannotated ∼190 genes and added more than 3600 alternatively spliced transcripts, including ∼1400 entirely novel exons and an additional ∼30kb of CDS. We have also completed an effort to capture all recently described unannotated microexons (18 (link)) into GENCODE, and further added an additional 146 novel microexons mined from public SLRseq data (19 (link)).
As part of the CCDS collaboration with RefSeq, we have checked a large subset of human loci where there was disagreement over gene biotype. Similarly, we have checked all UniProt manually annotated and reviewed (i.e. Swiss-Prot) accessions that lack an equivalent in GENCODE. As a result, we added 32 novel protein-coding loci to GENCODE and rejected more than 200 putative coding loci. Finally, we are manually reviewing genes previously annotated as protein-coding, but with weak or no support based on a method incorporating UniProt, APPRIS, PhyloCSF, Ensembl comparative genomics, RNA-seq, mass spectrometry and variation data (20 (link),21 (link)). Of the 821 loci investigated to date, 54 have had their coding status removed while a further 110 potentially dubious cases remain under review.
The approach taken reflects in the kinds of updates captured in the annotation. For example, the targeted reannotation in human leads to the annotation of few novel protein-coding loci but many novel transcripts at updated protein-coding and lncRNA loci. Conversely, in mouse the emphasis on clone-by-clone annotation identifies many more novel loci and transcripts across a broader range of biotypes (Figure 1).

Genomic DNA was extracted from bacteria using a Nucleospin Microbial DNA kit (Thermo Fischer Scientific) as per manufacturer instructions. Further, the QIASeq FX DNA Library Kit (Qiagen) prepared genomic libraries for sequencing. L. amnigena PTJIIT1005 genome was sequenced on NGS (Next Generation Sequencing) Illumina NovaSeq6000 Platform by Redcliffe Lifetech, Noida. A total of 9,365,132 raw reads were obtained; 8,596,940 Illumina reads were de novo assembled using Unicycler (version 0.4.4). The assembled genome sequence was annotated by the tool Prokka 1.12. The complete genome sequence was submitted to NCBI.
Average Nucleotide Identity (ANI) [17 ] measures nucleotide-level genome similarity between the coding regions of two genomes. The complete genome sequence was submitted in FASTA format as an input file. This tool gives the similarity index percentage [18 (link)]. ANI is computed using the formula [19 ]: $$\underset{\mathbf{G}\mathbf{1}\to \mathbf{G}\mathbf{2}}{\mathbf{gANI}}=\frac{\underset{\text{bbh}}{\mathrm{\Sigma }}\left(\text{Percent},,\text{Identity},{}_{.},{}^{*},\text{A},\text{lignment},,\text{length}\right)}{\text{lengths}\text{of}\text{BBH}\text{genes}}$$
Genome annotation of Lelliottia amnigena was done by RAST (Rapid Annotation using Subsystems Technology), PATRIC (The Pathosystems Resource Integration Center), and PGAP (Prokaryotic Genome Annotation Pipeline). Assembled genome sequence was submitted in RAST in FASTA format as input files, assigned functions to the genes. It also predicted the subsystems which were represented in the genome. By using this information, it reconstructs the metabolic network and makes the output file easily downloadable. Similarly, contigs were submitted in PATRIC as input files which provided annotation, subsystem summary, phylogenetic tree, and pathways. NCBI PGAP was used to annotate the bacterial genome where the complete genomic sequence was submitted in FASTA format as an input file, and it predicted the protein-coding regions and functional genome units like tRNAs, rRNA, pseudogenes, transposons, and mobile elements.

Insertion events are known to be caused by various mechanisms and have various consequences [26 (link)]. To characterize and investigate the origins of the detected insertions, we decomposed them into TRs, TEs, tandem duplications (TDs), satellite sequences, dispersed duplications, processed pseudogenes, alternative sequences, “deletions” in GRCh38, and nuclear mitochondrial DNA sequences (NUMTs).
We first applied Tandem Repeats Finder (TRF) [27 ] to all inserted sequences and defined TRs as having (1) element lengths < 50 bp and (2) covering more than 50% of an inserted sequence. After filtering TRs, we identified TEs using RepeatMasker [28 ] if (1) an inserted sequence covered a TE > 50%, (2) the inserted sequence was covered by the TE > 50% (reciprocal overlap), and (3) the total substitutions and indels were < 50% (matching condition).
Previous studies have reported that TDs are understudied but widespread [26 (link), 29 (link)]. After detecting TRs and TEs, we manually reviewed the remaining insertions and found that they contained TDs derived from non-repetitive regions in the reference. We considered these insertions as TDs. To identify this class of insertions, we aligned all insertions except TRs to GRCh38 using BLAT [30 (link)]. We then collected insertions mapped to original breakpoints within 5 bp with > 90% in BLAT identity and defined them as TDs. In this process, missing TRs with long repeat elements were found. Therefore, they were added to the TR callset if (1) an inserted sequence aligned within 500 bp from the insertion breakpoint and (2) the ratio of the total number of matching bases to the insertion length was > 0.5.
To understand the remaining insertions, we manually checked their features by aligning them to the reference using BLAT [30 (link)]. We identified insertions that were aligned from end to end to different chromosomal regions with high identity (> 90%). We defined these insertions as dispersed duplications. Next, we detected insertions aligned to a series of exons and untranslated regions (UTRs) of coding genes with high identity (> 90%) and classified them as processed pseudogenes. We also found other insertions aligned to the alternative sequences (e.g., “alt” or “fix” sequences) on BLAT with high identity (> 90%). We classified them as alternative sequences. Some of the insertions left at this point were thought to have arisen by deletion events in GRCh38 because they were securely aligned to the chimpanzee reference genome (panTro6), although they were classified as insertions when compared with GRCh38 [3 (link)]. We aligned the remaining insertions to the panTro6 assembly and categorized the insertions that lifted over panTro6 with high accuracy (> 90%) within 100 bp of the inserted position on GRCh38 as "deletions” in GRCh38. After this, the remaining insertions were manually reviewed, and features of the genomic regions (segmental duplications or self-chain) were examined.

Sequencing reads were adaptor-trimmed using TrimGalore (v0.0.6, Babraham Bioinformatics) and aligned to the mm10 mouse reference genome, the hg38 human reference genome, or the dm6 drosophila genome (6.27.1) using bowtie2 (v2.4.2⁶⁹ (link)). Only uniquely mapping reads were kept using samtools (v1.12⁷⁰ (link)) and duplicates removed using picard MarkDuplicates (v2.26.10, Broad Institute). Read numbers of all datasets are listed in Table S2. The generated data sets were used to define peak regions using MACS2 (v2.2.6⁷¹ (link)) with parameters --nomodel, read extension by model option (extension by 250 bp), using the pooled Inputs as control file and an FDR value of q = 0.05 with the broad option -b and a cutoff of q=0.1. External datasets were processed identically to own generated datasets. In case external datasets were sequenced paired-end, the paired-end function of MACS2 (v2.1.1) was used (FDR value of q= 0.05) using broad option -b and a cutoff of q=0.1.
External datasets:-

RING1B, SUZ12: Peak set (BED-file) taken from GSE132752⁶⁸ (link)

-

PCGF1 (GSE132752: GSM3891362-GSM3891364, Inputs: GSM3891380-GSM3891382)

-

PCGF2 (GSE132752: GSM3891368-GSM3891370, Inputs: GSM3891380-GSM3891382)

-

PCGF6 (GSE132752: GSM3891374-GSM3891376, Inputs: GSM3891380-GSM3891382)

-

H3K4me3 (ENCODE project [https://doi.org/10.17989/ENCSR000CGO], datasets: ENCFF997CAQ, ENCFF425ZMW)

-

RNAP-Ser5p (GSE128643: GSM3681611-GSM3681612²⁷ (link))

-

CpG Islands: Peak set (BED-file) taken from GSE21442⁶⁷ (link) and converted to mm10 using UCSC LiftOver

-

Enhancers: Enhancer set (BED-file) taken from the FANTOM consortium et al.⁶⁶ (link)

To obtain the final peak dataset, all called peaks of the individual replicates (n= 2-3) were concatenated, sorted using the BEDtools (v2.30.0,⁷³ (link)) SortBed function and overlapping peaks within 500 bp distance merged by the BEDtools Merge function. From this dataset only peaks that overlap the peak set of all individual replicates were kept using the BEDtools Intersect Interval function with -a specified and overlapping blacklisted regions were excluded. The resulting peak sets were used for subsequent analysis. H2BK120ub1-decorated genes were defined by first filtering the Ensembl gene list of mm10 (n=32,025) for not overlapping blacklisted regions and a gene length over 3000 bp to exclude short pseudogenes (n=20,380). 8,103 genes had at least one H2BK120ub1 peak in their transcribed region and 2,332 genes thereof were longer than 50kB. For HCT116 data, the Ensembl list of hg38 was used with the same filtering parameters resulting in a list of n=11,120 (>3kB, H2BK120ub1-positive), n=4,337 (>50kB, H2BK120ub1-positive) genes.