摘要
长非编码RNA（Long noncoding RNA, lncRNA）是一类转录本长度大于200 nt、缺乏蛋白质编码能力的RNA分子，可通过与mRNA、蛋白质和miRNA等分子相互作用从而以多种形式在不同水平调控基因表达。近年来，lncRNA在人类疾病发生发展、动物生长发育和植物应对环境胁迫等方面的研究一直是国内外研究热点和持续关注的焦点。目前，关于人类、动植物lncRNA的研究主要采用lncRNA测序（lncRNA-sequencing, lncRNA-seq）的方法。lncRNA-seq数据的分析依赖众多软件，且大部分是在Linux系统通过命令行的方式运行，对于大多数研究者开展数据分析可能存在困难。因此，文中在调研现有的lncRNA-seq数据分析流程的基础上，基于Linux系统构建lncRNA-seq本地分析平台，平台主要包括以下几个步骤：首先采用trimmomatic对原始测序数据进行质量控制，使用软件Hisat2和StringTie进行序列的比对、组装，通过Salmon软件进行转录本的定量，结合CPC2、CNCI、Pfam 和PLEK对新的转录本进行编码能力预测。另外，以2组马铃薯测序数据为例详细介绍数据分析过程，以便同领域研究者在服务器上构建本地化分析平台、灵活地开展数据分析。
Abstract
Long noncoding RNA (lncRNA) is a type of transcript with length longer than 200 nt and lacks protein coding ability. It can regulate the gene expression in various forms at different levels by interacting with mRNA, protein, and miRNA. In recent years, the studies of lncRNA in the occurrence and development of human diseases, animal growth and development, and plant response to environmental stress have attracted considerable attention. At present, the researches on the human, animal and plants lncRNAs are mainly performed based on the lncRNA-sequencing (lncRNA-seq) method. However, the analysis of lncRNA-seq data depends on many software and most of them are run in the Linux system, which might be difficult for most of researchers. Therefore, the paper investigates the existing lncRNA-seq data analysis process and constructs the localized lncRNA-seq analysis platform based on the Linux system. The platform mainly includes the following steps: first, the quality of the original sequencing data is controlled by the trimmomatic software and obtain the clean data; then, the clean data are mapped to the reference genome by hisat2 and the transcripts are assembled by stringtie; further, the quantification of the transcript is performed by salmon; the coding ability of the new transcript is predicted by four software, including CPC2, CNCI, Pfam and PLEK. In addition, two groups of potato sequencing data are taken as an example to introduce the data analysis process in detail, which can help the related researchers build a localized analysis platform on the server and carry out data analysis flexibly.
Keywords
长非编码RNA（Long noncoding RNA，lncRNA）是一类长度为200~10 000 nt 之间的RNA，在调控基因表达和多种生物学过程中发挥重要作用。鉴于lncRNA的功能重要性，lncRNA相关研究为国内外持续研究热点领域[1]。普通转录组测序（RNA-sequencing，RNA-seq）方法基于oligo（dT）捕获技术，主要利用mRNA通常都带有poly（A）尾巴的特性，将mRNA和部分lncRNA与总RNA分离。因此基于公共数据库中海量的普通转录组测序数据，越来越多的动植物lncRNA被鉴定出来。近年来，关于人类、动植物lncRNA的研究主要采用lncRNA测序（lncRNA-sequencing，lncRNA-seq）的方法，该测序方法采用去核糖体链特异性建库技术。与普通RNA-seq技术相比，lncRNA-seq方法检获的lncRNA更为全面，可同时揭示大量无poly（A）结构的lncRNA。但从lncRNA/RNA-seq数据中挖掘lncRNA需依赖众多软件、转件安装复杂，且大部分是在Linux系统通过命令行的方式运行，而大多数实验生物学家可能对Linux系统并不熟悉。
目前，国内外存在部分针对RNA-seq数据分析的本地和可视化分析平台，如Galaxy和BioJupies[2-3]。Galaxy作为生物医学领域最受欢迎的在线生物信息分析工具之一包含多种分析软件，如在tools下选择RNA-seq中的Hisat2和Stringtie可分别进行序列比对和组装。但是该网站仅含单个lncRNA编码能力预测软件，而为了保证数据的可靠性，对lncRNA编码能力的预测往往需要结合多个软件同时分析。BioJupies分析平台可使研究者通过直观的界面自动化挖掘GEO测序数据，也可提交原始测序数据进行分析。但是，BioJupies分析平台要求原始数据文件的大小在5 GB以下，同时不涉及lncRNA相关分析。
目前也存在部分流程可从RNA-seq数据中鉴定lncRNA，如UClncR采用Tophat/Hisat2比对序列、Stringtie进行转录本组装，进一步预测候选lncRNA，但该流程仅针对人类数据[4]。LncPipe流程可完成原始转录组数据到lncRNA鉴定的分析，该流程仅适用于人类、小鼠、果蝇和斑马鱼[5]。可见，目前已存在的从lncRNA/RNA-seq数据中鉴定lncRNA的流程均存在一定的局限性，如仅可应用于特定物种。因此，本论文拟结合应用案例，基于Linux系统构建从lncRNA/RNA-seq数据中挖掘lncRNA的本地化分析平台，能够从根本上解决实验生物学家生物信息学数据分析的难题，帮助同领域研究者在服务器构建本地化分析平台，灵活地开展数据分析。
1 系统设计
1.1 数据质量控制
测序原始数据可能包含低质量序列，影响后续数据分析，因此需采用质量控制软件删除低质量序列，同时删除测序过程中人为添加的接头和引物序列。目前lncRNA/RNA-seq原始数据的质量控制广泛采用Trimmomatic和fastp软件处理[6-7]。其中Trimmomatic软件的引用次数高达25 000余次，本论文采用该软件进行原始数据质控，如图1。
图1lncRNA数据分析流程图
Fig.1Workflow chart of lncRNA analysis
1.2 序列比对和组装
筛选快速准确的比对、组装软件。lncRNA-seq针对有参考基因组的样本进行，因此从lncRNA/RNA-seq数据中挖掘lncRNA需将序列比对到参考基因组，从而确定序列起始位点。目前序列比对主流软件包括TopHat2、STAR和Hisat2。TopHat2软件比较耗时，可能需要花费几天的时间来处理单个测序样本[8]。STAR软件使用后缀数组，处理速度比TopHat2快。然而，与Burrows-Wheeler 算法相比，后缀数组方法对内存要求较高[9]。较之以上两款软件，Hisat2运行速度最快、需内存低，精确度也较好[10]。就组装软件而言，目前StringTie和Cufflinks使用比较多，与Cufflinks相比，StringTie运行速度快、具有很高的准确性和新型异构体的发现能力，对于有多重异构体的物种更擅长[11]。鉴于以上原因，本文选择Hisat2进行序列比对，StringTie进行转录本组装。
1.3 转录本定量和编码能力预测
转录本定量软件包括kallisto、 eXpress、Sailfish和Salmon等[12]。鉴于Salmon软件具有准确度高、运行速度快、所需内存低等优点，该软件广泛应用于lncRNA/RNA-seq数据转录本的定量。
随后需要筛选lncRNA鉴定/编码能力预测软件。近年来研究者们开发出多种lncRNA鉴定软件，但是不同软件所获得的结果不一致。因此，需结合多种软件提高lncRNA鉴定的准确性。目前，广泛使用的转录本编码能力预测软件包括CPC2[13]、CNCI[14]、 Pfam[15]和PLEK[16]等。本文结合该四种软件来预测转录本的编码能力。
2 应用
2.1 分析系统及案例数据来源
分析平台基于Linux系统CentOS 7（Version: 3.10.0-1160.el7.x86_64）操作系统构建。
应用案例采用NCBI数据库中野生型马铃薯和转基因抗病马铃薯分别感染致病疫霉（Phytophthora infestans）24 h后的转录组测序数据（PRJNA318049），数据采用的测序方法为双末端测序，测序平台为Illumina HiSeq 2000。其中，数据SRR3371859、 SRR3371860和SRR3371861为野生型马铃薯叶片感染致病疫霉24 h后的转录组测序数据，数据SRR3371864、SRR3371865和SRR3371866为转基因抗病马铃薯叶片进行相同处理后的数据[17]。马铃薯基因组文件下载链接：http://spuddb.uga.edu/data/DM_1-3_516_R44_potato_genome_assembly.v6.1.fa.gz；马铃薯基因注释文件下载链接：http://spuddb.uga.edu/data/DM_1-3_516_R44_potato.v6.1.hc_gene_models.gff3.gz。
2.2 软件安装
分析平台主要涉及的软件如表1所示，软件的安装主要采用Anaconda，基于Anaconda的bioconda通道可方便快捷地一键安装绝大多数lncRNA-seq生物信息学分析软件，包括hisat2、stringtie、gffcompare、samtools、salmon、gffread、faSomeRecords、 PLEK和trimmomatic等，可极大程度地简化安装过程中软件的下载和环境配置。Anaconda安装的命令在CSDN博客上面有很多详细教程可供学习。Anaconda安装完成后，可先创建RNAseq环境，将上述软件安装在相同环境下，命令如下（括号内为注释内容，$符号代表一条命令的起始）：
表1所用分析软件及其功能和来源
Table1Analysis software used and its functions and sources
$ conda create-n RNAseq（创建名为RNAseq的环境）
$ conda activate RNAseq（激活RNAseq环境）
$ conda install-c bioconda hisat2（安装hisat2软件）
在使用软件之前，需要先激活相应的环境，如使用hisat2软件之前，需要先激活RNAseq环境（命令为：conda activate RNAseq），针对某种软件也可单独创建环境，如：
$ conda create-n stringtie
$ conda activate stringtie
$ conda install-c bioconda stringtie
少部分软件无法通过conda安装，可通过表1中的链接进行下载，然后在服务器进行安装。例如sratoolkit软件的下载与安装如下：
$ mkdir sratoolk（创建名为sratoolk的文件夹）
$ cd sratoolk（进入sratoolk目录）
$ curl-O https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz（下载当前最新Linux版本的sratoolkit软件）
$ tar xvf sratoolkit.current-centos_linux64.tar.gz（文件解压缩）
$ cd sratoolkit.3.0.1-centos_linux64/bin（进入bin目录）
$ echo export PATH=$ PATH:/../sratoolk/sratoolkit.3.0.1-centos_linux64/bin ＞＞~/.bashrc（配置环境，“..”为sratoolk上级详细路径）
$ cd ..（返回上一目录）
$ source~/.bashrc（配置环境）
2.3 数据分析流程
为方便数据分析，可建立一个专用文件夹（如localplatform）用于整个分析流程所产生的数据的存储。
第一步，原始测序数据的准备，基因组和注释文件的下载：
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-14/SRR003/371/SRR3371859.sralite.1（下载原始测序数据）
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-14/SRR003/371/SRR3371860.sralite.1
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-11/SRR003/371/SRR3371861/SRR3371861.sralite.1
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-11/SRR003/371/SRR3371864/SRR3371864.sralite.1
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-11/SRR003/371/SRR3371865/SRR3371865.sralite.1
$ wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos5/sra-pub-zq-14/SRR003/371/SRR3371866.sralite.1
$ wget http://spuddb.uga.edu/data/DM_1-3_516_R44_potato_genome_assembly.v6.1.fa.gz（下载基因组文件）
$ wget http://spuddb.uga.edu/data/DM_1-3_516_R44_potato.v6.1.hc_gene_models.gff3.gz（下载注释文件）
$ gzip-d DM_1-3_516_R44_potato_genome_assembly.v6.1.fa.gz（基因组文件解压缩）
$ mv DM_1-3_516_R44_potato_genome_assembly.v6.1.fa PotatoGenome.fa（将基因组文件重命名为PotatoGenome.fa）
$ gzip-d DM_1-3_516_R44_potato.v6.1.hc_gene_models.gff3.gz（注释文件解压缩）
$ mv DM_1-3_516_R44_potato.v6.1.hc_gene_models.gff3 Potatogenome.gff3（注释文件重命名）
$ gffread Potatogenome.gff3-T-o Potatogenome.gtf（注释文件格式转换）
第二步，分别采用sratoolkit软件将下载的原始测序数据转化为fastq格式、trimmomatic软件对转化后的双端测序数据进行质量控制，六个数据分别运行以下两条命令：
$ /../sratoolk/sratoolkit.3.0.1-centos_linux64/bin/fastq-dump /localplatform/ SRR3371859.sralite.1--split-3
$ trimmomatic PE-phred33-threads 10 /localplatform/SRR3371859.sralite.1_1.fastq /localplatform/SRR3371859.sralite.1_2.fastq /localplatform/SRR3371859.sralite.1_1.clean.fastq /localplatform/SRR3371859.sralite.1_1.unpaired.fastq /localplatform/SRR3371859.sralite.1_2.clean.fastq /localplatform/SRR3371859.sralite.1_2.unpaired.fastq AVGQUAL:30 1＞sample1.QC.log2＞&1
第三步，构建hisat2索引，采用hisat2软件分别将6个质控后的数据比对到参考基因组；通过samtools软件将文件格式由sam转化为bam：
$ hisat2_extract_splice_sites.py Potatogenome.gtf ＞ Potatogenome.ss（构建索引）
$ hisat2_extract_exons.py Potatogenome.gtf ＞ Potatogenome.exon（构建索引）
$ hisat2-build-p 16--exon Potatogenome.exon--ss Potatogenome.ss PotatoGenome.fa genome_tran（构建索引）
$ hisat2--dta-x /localplatform/genome_tran-1 /localplatform/SRR3371859.sralite.1_1.clean.fastq-2 /localplatform/SRR3371859.sralite.1_2.clean.fastq-S /localplatform/SRR3371859.sam（序列比对）
$ samtools sort-@ 8-o /localplatform/SRR3371859.bam /localplatform/SRR3371859.sam（格式转化）
第四步，采用stringtie软件基于以下命令分别对6个样本进行转录本的组装；然后对所有样本的转录本进行合并，其中mergelist.txt文件为六个样本组装后的gtf文件名称：
$ stringtie-p 8-G /localplatform/Potatogenome.gtf-o /localplatform/SRR3371859.gtf-l SRR3371859 /localplatform/SRR3371859.bam（转录本组装）
$ stringtie--merge-p 16-G /localplatform/Potatogenome.gtf-o /localplatform/localplatformstringtie_merged.gtf /localplatform/mergelist.txt（转录本合并）
第五步，采用gffcompare软件评估组装后的转录本与注释文件的匹配情况：
$ gffcompare-r/localplatform/Potatogenome.gtf-G-o localplatformmerged/localplatform/localplatformstringtie_merged.gtf
第六步，采用salmon软件对六个样本分别进行转录本的定量，通过该步获得的每个转录本在六个样本中的表达量文件（localplatform_Transcript.count）可用于转录本的差异表达分析：
$ gffread-w /localplatform/localplatform_total_ense.fa-g /localplatform/PotatoGenome.fa/localplatform/localplatformmerged.annotated.gtf（获取所有转录本的序列）
$ salmon index-t /localplatform/localplatform_total_ense.fa-i /localplatform/localplatform_transcripts_index（构建索引）
$ salmon quant-i /localplatform/localplatform_transcripts_index--validateMappings-l A-1 /localplatform/SRR3371859.sralite.1_1.clean.fastq-2 /localplatform/SRR3371859.sralite.1_2.clean.fastq-o /localplatform/SRR3371859.quant（定量）
$ salmon quantmerge--quants /localplatform/*.quant--column NumReads-o /localplatform/localplatform_Transcript.count（将6个样本的定量结果合并在一个文件）
第七步，筛选新的转录本，获取其外显子和长度信息：
$ awk′$ 2 == 0 &&$ 3 == 0 &&$ 4 == 0 &&$ 5 == 0 &&$ 6 == 0&&$ 7 == 0 {next} 1′ localplatform_Transcript.count ＞ localplatform_Transcript_count2（删除在六个样本中表达量均为0的转录本）
$ awk′$ 1!~/Soltu/′ localplatform_Transcript_count2 ＞ temp.txt &&mv temp.txt localplatform_Transcript_count3（筛选新的转录本）
$ awk′NR＞1 {print $ 1}′ localplatform_Transcript_count3 ＞ first_column_localplatform_Transcript_count3（获取新的转录本的ID）
$ awk′NR＞1 {print$ 5，$ 6，$ 10}′ localplatformmerged.localplatformstringtie_merged.gtf.tmap＞localplatformstringtie_merged.gtf.tmap2（获取所有转录本的ID、外显子和长度信息）
结合文件first_column_localplatform_Transcript_count3和localplatformstringtie_merged.gtf.tmap2，获取新转录本的外显子、长度信息，保存在文件Novel_transcript_exon_length中，脚本如下：
$ awk′{print $ 1", "}' first_column_localplatform_Transcript_count3 ＞ first_column_localplatform_Transcript_count3_2
$ awk′{print $ 1", ", $ 2, $ 3}′ localplatformstringtie_merged.gtf.tmap2 ＞ localplatformstringtie_merged.gtf.tmap2_2
$ grep-Ff first_column_localplatform_Transcript_count3_2 localplatformstringtie_merged.gtf.tmap2_2 | awk ′{print $ 0}′ ＞ test.txt
$ awk′{$ 1=$ 1; gsub (/, $ /, "", $ 1) ; print}′ test.txt ＞ Novel_transcript_exon_length
第八步，筛选外显子数大于1、长度大于等于200的新转录本（Novel_transcript），并且获取其序列信息（Novel_transcript.fa），进行后续蛋白质编码能力预测：
$ awk′ ($ 2＞1 &&$ 3＞=200) {print$ 0}′ Novel_transcript_exon_length ＞ Novel_transcript_exon_length_2
$ awk′{print $ 1}′ Novel_transcript_exon_length_2 ＞ Novel_transcript
$ faSomeRecords localplatform_total_ense.fa Novel_transcript Novel_transcript.fa
第九步，基于不同软件对外显子数和长度满足要求的新转录本进行蛋白质编码能力预测：
（1）采用CPC2软件进行蛋白质编码能力预测：
$ mkdir cpc2result
$ CPC2.py-i /localplatform/Novel_transcript.fa-o /localplatform/cpc2result/Novel_transcript_cpc2
根据结果文件Novel_transcript_cpc2.txt中的“coding_probability”列可查看每个转录本的编码能力。
（2）采用CNCI软件进行蛋白质编码能力预测：
$ mkdir CNCIresult
$ python CNCI.py-f /localplatform/Novel_transcript.fa-o /localplatform/CNCIresult/CNCI_novel_transcript_result-m pl-p 50
根据结果文件CNCI.index中的“index”列可查看每个转录本的编码能力。注：若分析物种为动物，则指定模型为“ve”，即“-m ve”
（3）采用Pfam软件进行蛋白质编码能力预测：
$ mkdir Pfamresult
$ transeq/localplatform/Novel_transcript.fa /localplatform/Pfamresult/Novel_transcript_AA.fa-frame=1
$ pfam_scan.pl-fasta /localplatform/Pfamresult/Novel_transcript_AA.fa-dir /PfamData/-out /localplatform/Pfamresult/Pfam_novel_transcript_result
结果文件Pfam_novel_transcript_result“seq id”列中存在的转录本具备蛋白质编码能力，因此筛选“seq id”列中不存在的转录本为潜在lncRNA。
（4）采用PLEK软件进行蛋白质编码能力预测：
$ mkdir PLEKresult
$ python PLEK.py-fasta /localplatform/Novel_transcript.fa-out /localplatform/PLEKresult/PLEK_novel_transcript_result-thread 30
结果文件PLEK_novel_transcript_result第一列显示为“Non-coding”则表明对应的转录本为lncRNA。
综上，通过以上九步分析流程可以获得转录本的表达量信息，以及基于四款软件的lncRNA预测结果。基于以上结果，研究者可采用R语言进行转录本差异表达分析，筛选不同样本间差异表达的lncRNA，以及进行lncRNA共表达靶基因预测等。基于R语言的下游数据分析流程较为简单，不依赖服务器即可完成，因此本文不再进行展示。
3 讨论与展望
lncRNA作为人、动植物体内重要的调控因子，广泛地参与不同的生物学过程。近年来关于lncRNA的研究依然为持续研究热点，通过lncRNA相关关键词在PubMed数据库中检索共检获58 433篇文章（截止到2024.7），近五年来的年文章数量均超过6 000余篇。目前，通过lncRNA/RNA-seq方法鉴定、筛选差异表达的lncRNA为揭示重要科学问题提供研究思路，但是大多数研究者对lncRNA分析流程并不熟悉。
文中整合已有软件，基于Linux系统构建了一套可以从lncRNA/RNA-seq数据中挖掘lncRNA的本地化分析平台。平台采用CPC2[13]、CNCI[14]、Pfam[15]和PLEK[16]四种软件进行转录本编码能力预测，提高lncRNA鉴定的准确性。该四种软件均适用于动植物lncRNA的鉴定，因此本平台可同时应用于动植物lncRNA的研究。另外，以马铃薯转录组数据为应用案例，详细展示数据分析流程，包括质量过滤、序列比对、组装、注释、定量和转录本编码能力预测等每一步所采用软件的使用方法以及详细的运行脚本。同时，数据分析过程中涉及到文本的处理，比如提取转录本ID等，因此本文在数据分析模块也展示了常见的文本处理命令的使用方法（如awk文本处理工具），可帮助相关研究者在本地服务器构建数据分析平台，灵活开展数据分析。
文中所构建的从lncRNA/RNA-seq数据鉴定lncRNA的本地化分析平台也存在一定的局限性。主要体现在需要研究者具备初级的Linux操作基础，能够在服务器进行简单的脚本运行。
