De novo transcriptome profiling of highly purified human lymphocytes primary cells

To help better understand the role of long noncoding RNAs in the human immune system, we recently generated a comprehensive RNA-seq data set using 63 RNA samples from 13 subsets of T (CD4+ naive, CD4+ TH1, CD4+ TH2, CD4+ TH17, CD4+ Treg, CD4+ TCM, CD4+ TEM, CD8+ TCM, CD8+ TEM, CD8+ naive) and B (B naive, B memory, B CD5+) lymphocytes. There were five biological replicates for each subset except for CD8+ TCM and B CD5+ populations that included 4 replicates. RNA-Seq data were generated by an Illumina HiScanSQ sequencer using the TruSeq v3 Cluster kit. 2.192 billion of paired-ends reads, 2×100 bp, were sequenced and after filtering a total of about 1.7 billion reads were mapped. Using different de novo transcriptome reconstruction techniques over 500 previously unknown lincRNAs were identified. The current data set could be exploited to drive the functional characterization of lincRNAs, identify novel genes and regulatory networks associated with specific cells subsets of the human immune system.

To help better understand the role of long noncoding RNAs in the human immune system, we recently generated a comprehensive RNA-seq data set using 63 RNA samples from 13 subsets of T (CD4 + naive, CD4 + T H 1, CD4 + T H 2, CD4 + T H 17, 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. There were five biological replicates for each subset except for CD8 + T CM and B CD5 + populations that included 4 replicates. RNA-Seq data were generated by an Illumina HiScanSQ sequencer using the TruSeq v3 Cluster kit. 2.192 billion of paired-ends reads, 2 × 100 bp, were sequenced and after filtering a total of about 1.7 billion reads were mapped. Using different de novo transcriptome reconstruction techniques over 500 previously unknown lincRNAs were identified. The current data set could be exploited to drive the functional characterization of lincRNAs, identify novel genes and regulatory networks associated with specific cells subsets of the human immune system.

Background & Summary
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 1 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 2 . 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 3 , in the establishment of the cardiovascular lineage and in the control of somatic tissue differentiation 4 . 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 5 and during naive CD4 + T cells differentiation into distinct helper T cell lineages [6][7][8] . Our results on human CD4 + T lymphocytes specific long intergenic non-coding RNAs (lincRNAs) 7 are in agreement with the findings in mice. In this work 63 RNA samples from 13 subsets of T (CD4 + naive, CD4 + T H 1, CD4 + T H 2, CD4 + T H 17, 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 in Fig. 1a as well as the number of biological replicates for each cell population. After RNA-seq sequencing we exploited different de novo transcriptome reconstruction approaches that led to the identification of over 500 previously unknown lincRNAs 9 . The general experimental design is shown in Fig. 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.

Methods
Purification of primary immunological cell subsets These methods are expanded from our previous article 9 . 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 (see Table 1). 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 10 (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 11 (version 1.4.1) and STAR 12 (version 2.2.0). The reference genome was indexed using Bowtie 13 (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 in Table 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.

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 a de 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 15 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 4 to 5 × 10 4 previously unknown transcripts for each lymphocyte population. The third approach was based on the Genome-guided Trinity 16 pipeline (see Supplementary File 1: example of command lines and Code Availability 9) (release 2012-10-05, http://trinityrnaseq.github.io/#genome_guided) that generates de novo transcripts by local assembly on previously mapped reads from specific locations. We used STAR instead of the Trinity's default aligner GSNAP 17 , 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 18 (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 three de 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 (see Data Citation 3 and Data 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 19 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 20 database (see Data 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 22 (cloned from https://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) (see http://genome.ucsc.edu/FAQ/ FAQformat.html#format5 and Data 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 Cabili et al. 1 , 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 each de 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 (see Supplementary 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'. The de novo annotation has been integrated concatenating it to public reference annotation and the resulting one was used for downstream analyses 7 .
The de novo annotation comprises 563 novel lincRNAs genes and 1,797 novel transcripts, published in our previous work 7 is available in Data Citation 1.

Data Records
In this study we deposited 1 dataset, which contains the RNA-Seq raw reads in FASTQ format (see Data  Citation 1 and Supplementary 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 (see Data Citation 1: the new 563 annotated lincRNAs 7 ) 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. RNAseq data generated were initially analyzed with FastQC and a summary plot with the data from all samples is depicted in Fig. 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 in Fig. 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). 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.

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 3 de 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 24 (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 (see Data 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 25 , so definition of a FPKM threshold would contribute to discriminate low abundant functional transcripts from technical or biological noise. In a recent study 26 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. In Fig. 3c  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.

Usage Notes
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 (see Data 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 from https://github.com/trinityrnaseq/trinityrnaseq.  For the evaluation of the coding potential of de novo transcripts we suggest to use other recently developed software that perform the classification more efficiently than PhyloCSF, such as iSeeRNA, CNCI 27 and CPAT 28 . 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 or in 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.