Long non-coding RNA detection workflow

Given a raw RNA-seq data, a typical workflow for the identification and annotation of lncRNAs is composed of two major parts – a part which is common to processing all RNA-seq data and a part dedicated to lncRNA.

In the common part, the first step (1a) involves generation of a control report for sequence quality, which is typically done using the FastQC software (Andrews et al., Reference Andrews2010). In the second step (2a), raw sequence reads are pre-processed by filtering out contaminations from sequencing adapters and by removing or trimming low-quality reads. The minimum threshold for read quality score is usually set to 20 (Wu et al., Reference Wu, Liu, Guo, Cheng, Ren, Chen, Li, Duan, Sun and Yang2018). The Trimmomatic software (Bolger et al., Reference Bolger, Lohse and Usadel2014) is a popular tool for pre-processing of raw RNA-seq data. In the next step (3a), cleaned sequences are mapped to the reference genome with the most commonly used software tools being: Tophat (Trapnell et al., Reference Trapnell, Pachter and Salzberg2009), Tophat2 (Kim et al., Reference Kim, Pertea, Trapnell, Pimentel, Kelley and Salzberg2013), Bowtie2 (Langmead and Salzberg Reference Langmead and Salzberg2012) or HISAT2 (Kim et al., Reference Kim, Langmead and Salzberg2015). Finally in the last step (4a), assembling of the uniquely mapped sequence reads into transcripts is most often done by Cufflinks (Trapnell et al., Reference Trapnell, Williams, Pertea, Mortazavi, Kwan, van Baren and Pachter2010) or StringTie (Pertea et al., Reference Pertea, Pertea, Antonescu, Chang, Mendell and Salzberg2015).

In order to proceed with the analysis specific to lncRNA, several transcript filtration steps are to be pre-imposed on the annotated RNA-seq data. The first step (1b) is to identify only novel transcripts, which do not correspond to the protein-coding part of the genome and do not represent previously annotated lncRNAs. It is most often performed using the Cuffcompare function of Cufflinks, which classifies all the available assembled transcripts based on their annotation to the pre-specified reference genome in the GFF format. Transcripts classified as ‘A transfrag falling entirely within a reference intron’, ‘Unknown, intergenic transcript’ or ‘Exonic overlap with reference on the opposite strand’ are then selected for downstream analysis. The genomic annotation of lncRNA in livestock is still very scarce as compared to human or mouse genomes; therefore, for livestock species, a much larger number of novel lncRNA transcripts is expected than for humans. Such potential candidate sequences for lncRNA are then subjected to several filtering steps (2b), which can differ depending on the analytic approach. The most important filtering steps comprise (i) filtering by length – removing transcripts shorter than 200 bp and transcripts longer than 10 000 bp containing a single exon, (ii) filtering by sequence content – removing transcripts overlapping with repeat or low complexity regions defined in the reference genome assembly (note, that for livestock genomes, this information is still limited), (iii) filtering by expression level – removing extremely high and extremely low expressed transcripts using the FPKM measure (fragments per kilo base of transcript per million mapped reads) to quantify the expression level; removal thresholds can either be arbitrarily chosen or estimated dynamically from the available data and (iv) filtering by protein-coding potential – various approaches involving removing transcripts containing known protein-coding domains using, for example, Transeq (El-Gebali et al., Reference El-Gebali, Mistry, Bateman, Eddy, Luciani, Potter, Qureshi, Richardson, Salazar, Smart, Sonnhammer, Hirsh, Paladin, Piovesan, Tosatto and Rober Finn2019) or HMMER (Eddy et al., Reference Eddy2011; Finn et al., Reference Finn, Clements and Eddy2011) software; removing transcripts with a significant hit in the Pfam database using, for example, PfamScan software (Bateman et al., Reference Bateman, Birney, Cerruti, Durbin, Etwiller, Eddy, Griffin-Jones, Howe, Marshall and Sonnhammer2002; Finn et al., Reference Finn, Bateman, Clements, Coggill, Eberhardt, Eddy, Heger, Hetherington, Holm, Mistry and Sonnhammer2014); removing transcripts, which products show similarity to known proteins from the RefSeq non-redundant protein database or the UniRef90 database (Suzek et al., Reference Suzek, Wang, Huang, McGarvey and Wu2015) using BLASTX (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990); removing transcripts based on their protein-coding potential level estimated, for example, by Coding Potential Calculator (CPC) (Kong et al., Reference Kong, Zhang, Ye, Liu, Zhao, Wei and Gao2007), Coding Potential Assessment Tool (CPAT) (Wang et al., Reference Wang, Park, Dasari, Wang and Kocher2013), Coding Non-Coding Index (CNCI) software (Sun et al., Reference Sun, Luo, Bu, Zhao, Yu, Zhang, Liu, Chen and Zhao2013) or predictor of Long non-coding RNAs and mEssenger RNAs based on an improved K-mer scheme (PLEK) (Li et al., Reference Li, Zhang and Zhou2014). The next step (3b) comprises merging of the RNA-seq data corresponding to novel lncRNAs, defined by the above workflow, and the known lncRNAs with positions defined in databases such as the ALDB (Li et al., Reference Li, Ma, Zeng, Wang, Geng, Yang and Cui2015) or the NONCODE. Note that various authors apply different thresholds for the percent of sequence identity and the percent of length of the aligned sequence to call lncRNAs. The final downstream analysis (4b) of the combined data sets depends on the experimental hypothesis and on the underlying experimental design. Most typically, it involves (i) the comparison of lncRNA expression levels between experimental conditions using, for example, DESeq2 software (Love et al., Reference Love, Huber and Anders2014), (ii) identification of target genes of differentially expressed lncRNAs, which can either be done by a dedicated software, for example, LncTar (Li et al., Reference Li, Ma, Zeng, Wang, Geng, Yang and Cui2015), by considering the physical proximity between a lncRNA and a protein-coding gene, or by considering high correlations between the expression level of a lncRNA and an messenger RNA (mRNA), and (iii) the functional annotation of target genes to metabolic pathways and/or gene ontologies using, for example, DAVID (Huang et al., Reference Huang, Sherman, Zheng, Yang, Imamichi, Stephens and Lempicki2009), KOBAS (Xie et al., Reference Xie, Mao, Huang, Ding, Wu, Dong, Kong, Gao, Li and Wei2011) programs, or the in-house tools provided by annotation databases, for example, the Gene Ontology database enrichment analysis tool (geneontology.org; Eilbeck et al., Reference Eilbeck, Lewis, Mungall, Yandell, Stein, Durbin and Ashburner2005) or the Reactome database analysis tool (reactome.org).

Function of long non-coding RNA

As it can be seen from the number of lncRNA genes or transcripts listed in Tables 1 to 3, there are more lncRNAs than mRNA molecules transcribed from a DNA template. Moreover, lncRNAs can be identified within various cell compartments, such as nucleus, nucleolus, cytoplasm and mitochondria. This reflects the variety of functions which they exhibit on cell metabolism. Technically, lncRNAs functions are related to (i) transcription through either chromatin modifications by interaction with enzymes or through interactions with transcriptional machinery proteins and miRNAs; (ii) post-transcriptional regulations of mRNA molecules, such as capping, alternative splicing, editing, transport, translation, degradation and stability; and (iii) epigenetic modifications manifested by the regulation of imprinting (Bhat and Jones Reference Bhat and Jones2016; Fernandes et al., Reference Fernandes, Guaman, Vasconcellos, Heinemann, Picardeau and Nascimento2019). On the organism level, lncRNAs are known to be abnormally expressed in many diseases with the most predominant influence on cancer and viral infections thereby playing a role of biomarkers. A practical example of lncRNA being a biomarker in human cancers is HOTAIR, whose overexpression results in development and metastases of several cancer types in humans (Lorenzi et al., Reference Lorenzi, Avila Cobos, Decock, Everaert, Helsmoortel, Lefever, Verboom, Volders, Speleman, Vandesompele and Mestdagh2019). In the healthy physiological state, lncRNAs play a role in organ differentiation during embryogenesis (Grote and Herrmann, Reference Grote and Herrmann2015) as well as in the process of aging (Xing et al., Reference Xing, Bai, Guo, Chen, Hua, Zhang, Ma, Ren, Wang and Wang2017). Most of the applications related to livestock investigate the functional annotation of lncRNAs, manifested by Gene Ontologies (GO) and KEGG pathways assigned to their target genes. These functions however are strongly related to the experimental design applied in each particular study and thus not of a universal nature. For example, in a recent study on pigs, Chen et al. (Reference Chen, Shi, Chen, Li, Li, Zou, Ch and Ch2019) applied the DAVID software to constructed clusters composed of GO terms and KEGG pathways characteristic for genes, which were targets of lncRNAs differentially expressed in relation to growth performance. In addition, You et al. (Reference You, Zhang, Liu, Song, Yang and Lian2019) applied the DAVID software for functional clustering of GO term and KEGG pathways related to target genes of lncRNAs differentially expressed in white leghorn chicken infected by the Marek’s disease virus and in a control, healthy group. In cattle, an example of functional annotation is the study of Gao et al. (Reference Gao, Li, Lai, Zhou, Wu, Huang, Lan, Lei, Chen and Dang2019). Using KOBAS, the authors tested a functional enrichment in GO terms and KEGG pathways of genes targeted by lncRNAs differentially expressed in two developmental stages of testis.

Table 1 Recent (since 2016) studies on long non-coding RNA (lncRNA), including long intergenic non-coding RNA (lincRNA), in Sus scrofa on a genome-wide scale

Table 2 Recent (since 2016) studies on long non-coding RNA (lncRNA), including long intergenic non-coding RNA (lincRNA), in Bos taurus on a genome-wide scale

Table 3 Recent (since 2016) studies on long non-coding RNA (lncRNA), including long intergenic non-coding RNA (lincRNA), in Gallus gallus on a genome-wide scale

Long non-coding RNA in livestock

Genomic annotation of long non-coding RNAs

One of the most complete and therefore most widely used databases that store lncRNAs is the NONCODE. It provides the collection and annotation of ncRNAs, especially lncRNAs, in 17 species, including livestock. The current version of the NONCODE (5.0) contains 548 640 transcripts, identified either by the RNA-seq, expression microarrays, or based on the literature. Yet, the database is far from complete while most of the studies listed in Tables 1 to 3 report a very large number of novel lncRNAs. For instance, Kern et al. (Reference Kern, Wang, Chitwood, Korf, Delany, Cheng, Medrano, Van Eenennaam, Ernst, Ross and Zhou2018) reported that only 18.3% of pig, 1.7% of cattle and 5.7% of poultry lncRNA transcripts from the NONCODE database overlapped with transcripts detected in their analysis.

We annotated lncRNAs of five species from the NONCODE database (accessed on 20 November 2019) using the Variant Effect Predictor software (McLaren et al., Reference McLaren, Pritchard, Rios, Chen, Flicek and Cunningham2010). 51 453 lncRNA genes were annotated to the Sscrofa11.1 (Warr et al., Reference Warr, Affara, Aken, Bickhart, Beiki, Billis, Chow, Eory, Heather, Flicek, Girón, Griffin, Hall, Hannum, Hourlier, Howe, Hume, Izuogu, Kim, Koren, Liu, Manchanda, Martin, Nonneman, O’Connor, Phillippy, Gary, Rosen, Rund, Sargent, Schook, Schroeder, Schwartz, Skinner, Talbot, Tseng, Tuggle, Watson, Smith and Archibald2019; accessed from www.ensembl.org/Sus_scrofa on 20 November 2019 with the corresponding Ensembl ID GCA_000003025.6) reference genome, 25 683 lncRNA genes were annotated to the ARS-UCD1.2 reference genome of Bos taurus (Shamimuzzaman et al., Reference Shamimuzzaman, Le, Unni, Diesh, Triant, Walsh, Tayal, Connat, Hagen and Elsik2020; accessed from www.ensembl.org/Bos_taurus on 20 November 2019 with the corresponding Ensembl ID GCA_002263795.2) and 22 843 lncRNA genes were annotated to the GRCg6a reference genome of Gallus gallus (accessed from www.ensembl.org/Gallus_gallus on 20 November 2019 with the corresponding Ensembl ID GCA_000002315.5). These livestock annotations were compared to the well-annotated mouse GRCm38 (318 287 lncRNA genes) and human GRCh38 (616 532 lncRNA genes) genomes. Since genomic coordinates of lncRNAs were defined in relation to different versions of reference genomes, we converted them to the most current assembly of each species using the liftOver software (Kent et al., Reference Kent, Sugnet, Furey, Roskin, Pringle, Zahler and Haussler2002). Species-specific annotations are visualized in Figure 1. The distribution pattern of lncRNA for Homo sapiens and Mus musculus is very similar. The observed difference between those two species and the livestock species demonstrates the incompleteness of the annotation of lncRNA in livestock. Precisely, a more complete annotation of the lncRNA in humans and mice is manifested by the fact that the categories representing the actual proper annotation for lncRNA (non-coding transcript exon variant, non-coding transcript variant and regulatory region variant) make up 39.35% of all annotations in Homo sapiens and 43.84% in Mus musculus, but only 13.90% in Sus scrofa, 12.54% in Gallus gallus and as little as 3.14% in Bos taurus. In livestock, the missing ncRNA annotation seems to be predominantly assigned to intergenic and intron sequences.

Figure 1 Genomic annotation of long non-coding RNA genes from the NONCODE database for (a) human (Homo sapiens), (b) mouse (Mus musculus), (c) cattle (Bos taurus), (d) poultry (Gallus gallus) and (e) pig (Sus scrofa). UTR = Untranslated Transcribed Region.