With emerging technologies it is becoming evident that the vast majority of the genome is transcribed (the so-called ‘dark matter of the genome’) and produces a diverse population of non-protein-coding RNAs (ncRNAs), including long non-coding RNAs (lncRNAs). LncRNAs are transcripts of more than 200 base pairs in length that are often expressed with higher cell-specificity compared to protein-coding genes¹ despite having lower expression levels. LncRNAs fold in functional domains that allow them to interact with other RNA molecules, DNA and proteins exerting a plethora of different functions in the cells, as chromatin remodeling (XIST, HOTAIR), transcriptional activation or repression (DBE-T, NeSR, lincRNA-Cox2), competition with microRNAs (linc-MD1, PTEN ceRNAs), splicing (sno-lncRNAs), RNA trafficking (NRON), mRNA stability (TINCR), imprinting (KCNQ1OT1) and translation (lncRNA-21p), among others². LncRNAs are also frequently expressed only in specific developmental stages, hinting to their involvement in cell fate determination. Moreover, lncRNAs have been implicated in the maintenance of stem cell pluripotency and differentiation³, in the establishment of the cardiovascular lineage and in the control of somatic tissue differentiation⁴. Altogether these findings clearly point out the fundamental role of lncRNAs in the control of cell differentiation and in the maintenance of cell identity. Indeed in the mouse immune system lncRNAs expression changes during naive to memory CD8⁺ T cell differentiation⁵ and during naive CD4⁺ T cells differentiation into distinct helper T cell lineages^6–8. Our results on human CD4⁺ T lymphocytes specific long intergenic non-coding RNAs (lincRNAs)⁷ are in agreement with the findings in mice. In this work 63 RNA samples from 13 subsets of T (CD4⁺ naive, CD4⁺ T_H1, CD4⁺ T_H2, CD4⁺ T_H17, CD4⁺ T_reg, CD4⁺ T_CM, CD4⁺ T_EM, CD8⁺ T_CM, CD8⁺ T_EM, CD8⁺ naive) and B (B naive, B memory, B CD5⁺) lymphocytes were collected. The hierarchy of T and B cells during differentiation of the analyzed subsets is depicted inFig. 1a as well as the number of biological replicates for each cell population. After RNA-seq sequencing we exploited differentde novo transcriptome reconstruction approaches that led to the identification of over 500 previously unknown lincRNAs⁹. The general experimental design is shown inFig. 1b. As recent findings suggest that lncRNAs might contribute to the definition of lymphocytes identity and to the modulation of their functional plasticity, our data set could be used as a resource to guide the validation and functional characterization of lincRNAs and to identify genes and regulatory networks associated with specific cells subsets of the human immune system.

(a) Hierarchical representation of the different cell subset originating from hematopoietic stem cells. In this study 13 human primary lymphocyte subsets were profiled: CD4⁺ naive; CD4⁺ T_h1; CD4⁺ T_h2; CD4⁺ T_h17; CD4⁺ T_reg; CD4⁺ T_EM; CD4⁺ T_CM; CD8⁺ T_CM; CD8⁺ T_EM; CD8⁺ naive; B naive; B memory; B CD5⁺. The number of biological replicates and the expressed genes (FPKM>0.21) for each population is indicated. The total number of samples profiled in this study is 63. (b) General overview of the bioinformatic steps and approaches used for the identification of novel lincRNAs.

Full size image

Purification of primary immunological cell subsets

These methods are expanded from our previous article⁹. Blood buffy coat cells of healthy donors were obtained from Fondazione Istituto di Ricovero e Cura a Carattere Scientifico Ca’Granda Ospedale Maggiore Policlinico in Milan, and peripheral blood mononuclear cells were isolated by ficoll-hypaque density-gradient centrifugation. The ethical committee of Fondazione Istituto di Ricovero e Cura a Carattere Scientifico Ca’Granda Ospedale Maggiore Policlinico approved the use of peripheral blood mononuclear cells from healthy donors for research purposes, and informed consent was obtained from subjects. Human blood primary lymphocyte subsets were purified to a purity of >95% by cell sorting through the use of various combinations of surface markers (seeTable 1).

RNA isolation and RNA sequencing

Total RNA was isolated with a mirVana Isolation Kit (Ambion). Libraries for Illumina sequencing were constructed from 100 ng of total RNA with the Illumina TruSeq RNA Sample Preparation Kit v2 (Set A). The libraries generated were loaded on to the cBot automated clonal amplification system (Illumina) for clustering on a HiSeq Flow Cell v3. The libraries clustered on a HiSeq Flow Cell v3 were then sequenced with a HiScanSQ optical imaging system (Illumina). A paired-end run (with a read length of 100 bases) was performed with an SBS Kit v3 DNA sequencing kit (Illumina). Real-time analysis and base calling was performed with HiSeq Control Software (version 1.5, Illumina). CASAVA (version 1.8.2, Illumina) software was used to demultiplex reads into specific sample and groups, the software was configured to operate with ‘--mismatches='1'’ allowing one mismatch during the identification of the indexes (Data Citation 1).

RNA-seq trimming and mapping

To improve sequence quality, samples data were cleaned by Trimmomatic¹⁰ (version 0.30) using the following parameters (LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36) giving as input the forward and reverse FASTQ sequences for each sample. Only the reads that passed the quality or length threshold on both strands were considered for mapping. The whole data set was aligned to human genome assembly GRCh37 (Genome Reference Consortium Human Build 37) using both TopHat¹¹ (version 1.4.1) and STAR¹² (version 2.2.0). The reference genome was indexed using Bowtie¹³ (version 0.12.9) for TopHat alignment. Both TopHat and STAR were used with default parameters; only for TopHat we specified the mate-inner-dist parameter for each sample of our data set (see the associated Metadata Record). Overall read depth and coverage information of the dataset is reported inTable 2. RNA-seq data from the Illumina Human BodyMap 2.0 project (Data Citation 2) consisting of 16 human tissues were downloaded, processed and mapped using the same criteria.

Public reference annotation

Ensembl database (version 67 from May 2012, seeData Citation 3) annotation was integrated with a previously published catalogue of lincRNAs¹ (seeData Citation 4) using Cuffcompare which is provided by the Cufflinks¹⁴ (version 2.1.1) suite. BioMart was used to categorize Ensembl annotation in different classes by their biotype: ‘lincRNA’ (5,804 genes), protein-coding genes (21,976 genes), receptor-encoding using GO term GO:000487 (2,043 genes encoding molecules with receptor activity function) and the class of genes encoding molecules involved in metabolic processes corresponding to GO term GO:0008152 (7,756 genes). The final public reference annotation consisted of 195,392 transcripts that referred to 62,641 genes, 11,170 of which were non-redundant lincRNA-encoding genes.

De novo genome-based transcripts reconstruction

To identify putative novel genes, not yet annotated and specifically expressed in our datasets, we combined multiple tools and their outputs following ade novo genome-based transcripts reconstruction procedure. Samples were aggregated in meta datasets corresponding to the 13 lymphocyte populations. These meta datasets were aligned to the reference genome using two mappers: TopHat and STAR. The resulting 26 alignments were used as independent inputs for Cufflinks configured to use the RABT¹⁵ assembler for the identification of novel transcripts. The following parameters were used in combination with Cufflinks: ‘-g’ to guide the assembly by the public reference annotation. With these approaches we identified about 3×10⁴ to 5×10⁴ previously unknown transcripts for each lymphocyte population. The third approach was based on the Genome-guided Trinity¹⁶ pipeline (seeSupplementary File 1: example of command lines and Code Availability 9) (release 2012-10-05,http://trinityrnaseq.github.io/#genome_guided) that generatesde novo transcripts by local assembly on previously mapped reads from specific locations. We used STAR instead of the Trinity’s default aligner GSNAP¹⁷, as it performed better in terms of both accuracy and computing time. For the first alignment phase STAR was used with the default parameters. The ‘Genome-guided Trinity’ suite was used with the parameters suggested in the main documentation (default). Each candidate transcript was then processed via the Program to Assemble Spliced Alignments¹⁸ (PASA,http://pasapipeline.github.io/). PASA is a genome annotation tool that reconstructs the complete transcript and gene structures, resolves incongruences derived from transcript misalignments and alternatively splices events, refines the public reference annotation and proposes new transcripts and genes in case no previous annotation can explain the new data. PASA was configured to use STAR as aligner. We recompiled STAR to enable it for handling long reads (putative transcripts); the file ‘IncludeDefine.h’, from the source code, was modified setting the variable ‘MAX_READ_LENGTH’ to a value of ‘100000’. Recompiling the source tree using the GNU ‘make’ utility with the command ‘make STARlong’ generated the desired modified binary version of STAR.

Identification of previously unknown lincRNA-encoding genes

Data generated by the three different approaches, TopHat/Cufflinks; STAR/Cufflinks; STAR/Trinity, were separately processed to identify unknown lincRNA-encoding genes.

The threede novo methods applied to each lymphocyte population, generated transcripts and genes without prior knowledge on their ability to encode for proteins or not. In order to identify only the putative novel lincRNAs, known transcripts and previously unknown isoforms of already annotated genes were filtered out. To perform this filtering we compared the public available reference annotations (seeData Citation 3 andData Citation 4) with the datasets produced by each approach using a custom script (see Code Availability 10). This comparison can be performed using more consolidated tools as the UCSC bedtools¹⁹ or Cuffcompare. Transcriptional noise and low polymerase fidelity can create artifactual transcripts therefore only multi-exonic transcripts longer than 200 bases were retained in our analysis. Protein family domains available from Pfam²⁰ database (seeData Citation 5) were searched in all transcripts using the HMMER3 (ref. 21) algorithm and those transcripts that matched at least one of all six possible frames were discarded. Another criteria commonly accepted to define lincRNA is the evaluation of their coding potential; absence of coding potential is distinctive of putative lincRNA. PhyloCSF²² (cloned fromhttps://github.com/mlin/PhyloCSF on Oct. 2013) is a comparative genomics method (phylogenetic codon substitution frequency) built upon a multiple sequence alignment of 29 mammalian genomes in multi-alignment file format (MAF) (seehttp://genome.ucsc.edu/FAQ/FAQformat.html#format5 andData Citation 6). The entire set of novel transcripts that passed the previous filters was used as input for PhyloCSF. Transcripts scoring more than 100 decibans (PhyloCSF scores were obtained using option--frames=6) were excluded from the final catalog. This threshold was calculated by Cabiliet al.¹, as it corresponds to a false-negative rate of 6% for coding genes (i.e., 6% of coding genes are classified as noncoding) and a false-positive rate of ~10% (i.e., 9.5% of noncoding transcripts are classified as coding). They optimized PhyloCSF specificity and sensitivity threshold for the classification of coding and noncoding transcripts on the RefSeq reference sequence database of the National Center for Biotechnology Information (RefSeq coding and RefSeq lincRNAs).

De novo transcriptome data integration

In order to create a comprehensive and unique annotation of novel lincRNAs identified in lymphocytes, duplicates generated by the three approaches adopted must be resolved. To accomplish this task Cuffcompare was used. For eachde novo reconstruction approach Cuffcompare merged the transcripts generated by all the populations. The result is a set of three distinct annotations corresponding to TopHat/Cufflinks, STAR/Cufflinks, STAR/Trinity/PASA. These three lincRNA sets were further merged to generate a non redundant atlas of lincRNAs in human lymphocytes and only those genes identified by at least two out of the three software programs were considered. After data integration through Cuffcompare, a custom script (seeSupplementary File 2 and Code Availability 11) was used to remove and substitute the internal gene id (XLOCs) and internal transcript id (TCONs) assigned by the software with their original and public names.

New lincRNAs were then uniquely identified with a name that contains the prefix ‘linc-’; the Ensembl gene name of the nearest protein-coding gene (irrespective of the strand); the location of the lincRNA relative to the sense of transcription of the nearest protein-coding gene: ‘up’ or ‘down’; the description of the concordance of the transcription between the lincRNA and its nearest coding gene: ‘sense’ or ‘antisense’; a counter to distinguish between lincRNA that share the same nearest protein-coding gene. An example of template name is ‘linc-geneX-(up|down)-(sense|antisense)_#n’. Thede novo annotation has been integrated concatenating it to public reference annotation and the resulting one was used for downstream analyses⁷.

Thede novo annotation comprises 563 novel lincRNAs genes and 1,797 novel transcripts, published in our previous work⁷ is available inData Citation 1.

Code availability

1. 1

   CASAVA (version 1.8.2, Illumina,https://support.illumina.com/sequencing/sequencing_software/casava.html), using the mismatch option--mismatches='1'

2. 2

   Trimmomatic¹⁰ (version 0.30,http://www.usadellab.org/cms/index.php?page=trimmomatic), PE, -phred33, LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, MINLEN:36

3. 3

   TopHat¹¹ (version 1.4.1,https://ccb.jhu.edu/software/tophat/index.shtml), mate-inner-dist was set for each dataset to the InnerSize filed available fromSupplementary Table 1

4. 4

   STAR¹² (version 2.20,https://github.com/alexdobin/STAR/releases), default

5. 5

   Bowtie¹³ (version 0.12.9,http://bowtie-bio.sourceforge.net/index.shtml), default

6. 6

   Cufflinks¹⁴ (version 2.1.1,http://cole-trapnell-lab.github.io/cufflinks/), option: -g Tells Cufflinks to use the supplied reference annotation a GTF file to guide RABT assembly.

7. 7

   HMMER3 (ref. 21) (version 3.0,http://hmmer.janelia.org/), default

8. 8

   PhyloCSF²² (release Oct. 2013,https://github.com/mlin/PhyloCSF), --frames=6

9. 9

   Custom template for Trinity/Pasa,Supplementary File 1,https://gist.github.com/helios/7cb7a6afbe625749e824 under GPL v3.0

10. 10

    Custom Biogem package, freely available athttps://github.com/ingm-oss/bioruby-genomic-intervals under GPL v3.0

11. 11

    Custom rename_gtf.rb,Supplementary File 2, a ruby script freely available athttps://github.com/ingm-oss/rename_gtf under GPL v3.0

In this study we deposited 1 dataset, which contains the RNA-Seq raw reads in FASTQ format (seeData Citation 1 andSupplementary Table 1),which is a simplified version of the ISA-TAB (see the associated Metadata Record). This dataset contains 63 samples in total, grouped by 13 lymphocyte subsets with 4 or 5 biological replicates each.Supplementary Table 1 is an XLSX with the following header: Source, the original source name used by the lab; Name, assigned by the provider; SubSet, the lymphocyte subset; Antibody, the antibodies used for sorting; InnerSize, the estimated inner size; R1 URI, the forward reads uri for download; R1 MD5SUM, checksum for the forward reads; R2 URI, the reverse reads uri for download; R2 MD5SUM, checksum for the reverse reads. The annotation of the 563 newly described lincRNA (seeData Citation 1: the new 563 annotated lincRNAs⁷) is a General Transfer Format (GTF).

RNA-seq raw data quality

Assessing the quality of the data performing the Quality Control (QC) is crucial to the whole study. RNA-seq data generated were initially analyzed with FastQC and a summary plot with the data from all samples is depicted inFig. 2a. The quality of the reads during the sequencing tends to decrease but it can be further improved using specific software that removes low quality bases reducing the length of the read or directly discard the whole read when its quality is too low. To perform the trimming and filtering Trimmomatic was run on each sample and the data were later on reanalysed with FastQC to confirm the quality improvements. The summary of the resulting data is shown inFig. 2b Another criteria to measure the QC for NGS reads is the % of GC content, which is improved by the filtering (Fig. 2c). Moreover the trimming step did not dramatically decrease the final number of reads (Fig. 2d).

(a) Phred quality score of the average distribution over all reads across all samples in each base before and (b) after trimming. (c) %GC content before and after trimming. (d) Detailed overview of the human lymphocyte subsets profiled: raw reads (black), the reads trimmed and filtered by quality (blue), and (e) the comparison of the mapped reads using TopHat (light green) and STAR (light orange).

Full size image

During the study two mapping software were used, TopHat and STAR. To exclude the possibility of discordance between the two aligners, the mapping results were compared to assess their mapping performance. The alignments with the two software showed a good concordance (96%) with a slight advantage of STAR in terms of mapped reads (Fig. 2e).

RNA-seq biological replicates

Biological replicates are fundamental to guarantee data consistency, in this study the lymphocyte populations profiled have 4 biological replicates for B CD5⁺ and CD8⁺ T_CM and 5 biological replicates for all the other populations. In order to establish the congruency among biological replicates Principal Component Analysis (PCA) (Fig. 3a) and hierarchical clustering (Fig. 3b) were performed. A good separation between B and T cells samples is achieved by PCA on normalized read counts using DESeq2 (ref. 23). Comparable results are obtained using hierarchical clustering on the same data. Moreover, similarity between biological replicates of the same population showed a good consistency and correlation among them.

(a) Principal Component Analysis (PCA) performed using DESeq2 rlog-normalized RNA-seq data. Loadings for principal components 1 (PC1) and PC2 are reported in graph (on x and y-axes). (b) Hierarchical clustering analyses performed using DESeq2 rlog-normalized RNA-seq data. Color code (from white to dark blue) refers to the distance metric used for clustering (dark blue corresponds to the maximum of correlation values). (c) Violin plot of the normalized FPKM values for the newly identified lincRNAs, previously annotated lincRNAs and transcription factors genes. The black line represents the normalized FPKM threshold (0.21 FPKM).

Full size image

De novo transcripts identification

Multiple combinations of software and filters were used for the identification of lincRNAs in the 13 lymphocytes populations. Moreover, we considered only newly described lincRNAs detected in at least 2 out of 3de novo approaches to improve the reliability of the data.

LincRNAs discrimination between coding and non-coding RNA depends on the algorithm used to asses the coding potential, in this study was used PhyloCSF. The final dataset of putative lincRNAs was further processed using iSeeRNA²⁴ (webserver version 1.2.2) in order to verify our results using a different approach based on Support Vector Machines (SVM). The classification we obtained is highly concordant, in fact ~99% of the putative lincRNAs contained in the final catalogue (seeData Citation 1) are classified as 'noncoding' also by iSeeRNA.

Expression threshold definition

As reported in literature, many lincRNAs are expressed at lower levels than protein coding genes²⁵, so definition of a FPKM threshold would contribute to discriminate low abundant functional transcripts from technical or biological noise. In a recent study²⁶ an approach based on the integration of RNA-seq and CHIP-seq data was used for the definition of a sensible FPKM threshold. 17 human cell lines from ENCODE project were analyzed to establish a relationship between gene expression levels and promoters activities. The expression cutoff was set where the fraction of genes associated to active promoters is equal to the fraction of them associated to repressed promoters. We considered a threshold of 0.21 FPKM that is the mean of the data reported for different cell types in the paper. InFig. 3c is shown that newly identified lincRNAs (selected with FPKM expression values >0.21), have higher expression levels compared to previously annotated ones with expression values above the threshold and with no threshold (light green and light blue). We then considered for the downstream analysis only genes whose expression values were at least 0.21 FPKM in one population.

This study was performed on the version 67 from May 2012 of Ensembl GRCh37. In order to access and use the catalogue of newly described lincRNAs generated in this study (seeData Citation 1 ), researchers must update it to the most updated genome version using the liftover software from UCSC (https://genome.ucsc.edu/cgi-bin/hgLiftOver) or the assembly converter from Ensembl (http://www.ensembl.org/Homo_sapiens/Tools/AssemblyConverter).

Software used during this study went through minor and major code base updates. The more notable software suite that has been updated during the time is Trinity and is strongly suggested to use the latest release downloadable fromhttps://github.com/trinityrnaseq/trinityrnaseq.

For the evaluation of the coding potential ofde novo transcripts we suggest to use other recently developed software that perform the classification more efficiently than PhyloCSF, such as iSeeRNA, CNCI²⁷ and CPAT²⁸. It has been demonstrated that these algorithms have a higher level of accuracy, and execution times are considerably faster.

Our transcriptome analysis includes thirteen human primary T and B cells whereas most of the available immune system datasets are limited to mice samples, cell lines orin vitro expanded cells. Therefore this study represents a valuable resource for those researches who need to access and analyze the expression patterns of both coding and non-coding transcripts in human lymphocytes. Moreover the thorough analysis we performed to assess the expression of both novel and previously annotated lincRNAs in these cells sets the grounds for further studies on the still largely uncharacterized function of long non-coding RNA in human lymphocytes subsets.

How to cite this article: Bonnal, R. J. P.et al.De novo transcriptome profiling of highly purified human lymphocytes primary cells.Sci. Data 2:150051 doi: 10.1038/sdata.2015.51 (2015).

M. Moro and M.C. Crosti for technical assistance with cell sorting; B.J. Haas and A. Dobin for help with the integration of genome-guided Trinity with STAR aligner; S. Donizetti for the IT support; the Istituto Nazionale Genetica Molecolare Bioinformatics Facility for support; Supported by Il Consiglio Nazionale delle Ricerche–Il Ministero dell’Istuzione dell’Universita e della Ricerca (Progetto Bandiera EPIGEN), Fondazione Cariplo (grant n. 2013-0955), the Associazione Italiana per la Ricerca sul Cancro (grant n. IG2013-ID14596), the European Research Council (grant n. 269022 to S.A.; grant n. 617978 to M.P.) and Fondazione Romeo ed Enrica Invernizzi.

Authors and Affiliations

1. Istituto Nazionale Genetica Molecolare ‘Romeo ed Enrica Invernizzi’, Via F. Sforza 35, Milan, 20122, Italy

   Raoul J.P. Bonnal, Valeria Ranzani, Alberto Arrigoni, Serena Curti, Ilaria Panzeri, Paola Gruarin, Sergio Abrignani, Grazisa Rossetti & Massimiliano Pagani

2. Department of Clinical Sciences and Community Health, Università degli Studi di Milano, Via Festa del Perdono 7, Milano, 20122, Italy

   Sergio Abrignani

3. Department of Medical Biotechnology and Translational Medicine, Università degli Studi di Milano, Via Festa del Perdono 7, Milano, 20122, Italy

   Massimiliano Pagani

Contributions

R.J.P.B set up the bioinformatics pipelines and wrote the manuscript; V.R. set up the bioinformatics pipelines and contributed to the preparation of the manuscript; A.A. set up the bioinformatics pipelines and contributed to the preparation of the manuscript; G.R. designed and performed the main experiments, analyzed the data and contributed to the preparation of the manuscript; S.C., I.P. and P.G. prepared the library and performed the experiments; S.A. and M.P. designed the study, supervised research and wrote the manuscript; and all authors discussed and interpreted the results.

Corresponding authors

Correspondence toRaoul J.P. Bonnal,Grazisa Rossetti orMassimiliano Pagani.

Competing interests

The authors declare no competing financial interests.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visithttp://creativecommons.org/licenses/by/4.0 Metadata associated with this Data Descriptor is available athttp://www.nature.com/sdata/ and is released under the CC0 waiver to maximize reuse.

Reprints and permissions