微科盟全转录组数据分析(两文库版)结题报告

1. 背景介绍

非编码RNA（ non-coding RNA，ncRNA ）是细胞中不进行蛋白翻译的 RNA ，ncRNA 包括 rRNA 、 tRNA 、 snRNA 、 snoRNA 、 lncRNA等,这些 RNA 可以与蛋白配合发挥功能也可以独立发挥催化功能。从发现 tRNA 、 rRNA 等开始，具有生物学作用的 ncRNA 已有多年历史。

这些 ncRNA 包含丰富的调节和功能单元，其中lncRNA占比较大，即大多数转录但不编码蛋白的序列均与lncRNA相关，长链非编码 RNA（ long non-coding RNA，lncRNA ）是指长度超过 200nt 的 RNA ，其本身不编码蛋白，而是以 RNA 的形式形成多层面调控基因的表达。

目前，已从单细胞真核生物以及人类中鉴定出数以万计的 lncRNA 基因座，并依据其相对于蛋白质编码基因的相对位置，将其分为基因间 lncRNA ，内含子 lncRNA ，正义 lncRNA ，反义 lncRNA ，双向 lncRNA 五类。

lncRNA 在表观遗传、顺式或反式转录及转录后水平上参与多层次基因表达调控，参与 X 染色体沉默、基因组印记以及染色质修饰、转录激活、转录干扰、核内运输等生物学进程，但经实验鉴定的与疾病相关的 lncRNA 数量仅为已鉴定基因座的1％不到，其生物学功能有待进一步挖掘。

miRNA 是一类小的调节 RNA，在细胞增殖，细胞死亡，细胞发育和分化，病毒感染，造血，肿瘤发生等生物过程中发挥重要作用。对于 miRNA 的研究还有助于了解淋巴瘤，白血病，癌症等多种疾病的发生发展过程。近年来越来越多的研究证实了 miRNA 参与了包括生长发育、细胞分化、细胞凋亡、脂类代谢、激素分泌、信号转导及逆境胁迫反应等在内的多种生物过程。miRNA 在物种进化中相当保守,其高度的保守性与其功能有着密切的关系。随着miRNA 在不同生物中的大量发现，对 miRNA 功能的研究也越来越引起人们的重视。例如lin-4，let-7 等miRNA其功能及靶基因已经明确。

然而对很多的 miRNA，其功能研究还属于初级阶段。早年研究 miRNA 的功能是通过实验的手段对 miRNA 或对应靶点进行突变，过表达miRNA，来观察个体表型的变化。但其具有没有办法广谱的鉴定其潜在的功能，或者表型没有发生显著的变化。因此利用表达谱及生物信息学的方法研究 miRNA 其及其靶基因成为研究 miRNA 功能的重要方法。

全转录组关联分析是基于非编码RNA分子参与调控mRNA的表达以及内源RNA竞争假说。该假说的内容是可以与miRNA结合的各种类型RNA可以与miRNA构成RNA诱导沉默复合体来mRNA的表达水平进行调节。因此，全转录关联分析的目的是构建各类非编码RNA与mRNA间的靶向调控关系以及非编码RNA之间的互作网络。

2. 项目流程

从 RNA 样品到最终数据获得，样品检测、建库、测序每一个环节都会对数据质量和数量产生影响，而数据质量又会直接影响后续信息分析的结果。为了从源头上保证测序数据的准确性、可靠性，我们对样品检测、建库、测序每一个生产步骤都严格把控，从根本上确保了高质量数据的产出：

图2.1 项目流程图

2.1 总 RNA 样本提取与检测

对 RNA 样品的检测主要包括 4 种方法：

  (1)琼脂糖凝胶电泳分析 RNA 完整性。

  (2)Nanodrop 检测 RNA 的纯度（OD260/280 比值）。

  (3)Qubit 对 RNA 浓度进行精确定量。

  (4)Agilent 2100 精确检测 RNA 的完整。

2.2 文库构建

2.2.1 lncRNA文库构建

RNA 样本检测合格后，使用 Ribo-Zero kit 从总 RNA 样品中去除 rRNA（lncRNA 有些具有与 mRNA 相同的 polyA 尾结构，有些不含有，用去除 rRNA 的方法能够最大程度地保留含有 polyA 尾的 lncRNA ）。向富集得到的 RNA 中加入 fragmentation buffer 将 RNA 打断成小片段。然后以片段化的 RNA 为模板，加入 6bp 随机引物（ random hexamers ）反转录合成 cDNA 第一条链，再加入缓冲液、 dNTPs（ dNTP 中的 dTTP 用 dUTP 取代）、 DNA polymerase I 和 RNase H 合成 cDNA 第二条链。合成的双链 cRNA 经纯化、末端修复、加 A 、连接测序接头处理后，用 USER 酶降解含有 U 的 cDNA 第二链并进行 PCR 富集，最后用 AMPure XP beads 纯化 PCR 产物，得到最终的链特异性文库。文库构建完成后，先使用 Qubit 2.0 进行初步定量，稀释文库，随后用 Agilent 2100 对文库的插入片段大小进行检测，插入片段符合预期后，使用 qPCR 方法对文库的有效浓度进行准确定量，以保证文库质量。文库构建原理图如下：

图2.2.1 lncRNA 文库构建原理图

2.2.2 miRNA文库构建

样品检测合格后，使用 Small RNA Sample Pre Kit 构建文库，利用 Small RNA 的 3’及 5’端特殊结构（5’端有完整的磷酸基团，3’端有羟基），以 total RNA 为起始样品，直接将 Small RNA 两端加上接头，然后反转录合成 cDNA。随后经过 PCR 扩增,PAGE胶电泳分离目标 DNA 片段，切胶回收得到的即为 cDNA 文库。 构建原理图如下：

图2.2.2 miRNA文库构建原理图

2.3 文库构建与质检

文库构建完成后，先使用 Qubit2.0 进行初步定量，稀释文库至 1ng/ul ，随后使用 Agilent 2100 对文库的插入片段大小进行检测，并使用 qPCR 对文库的有效浓度进行准确定量（文库有效浓度 > 2nM ），以保证文库质量。

2.4 上机测序

库检合格后，按照有效浓度及目标下机数据量的需求将不同文库混合在一起， cBOT 成簇后使用 Illumina 高通量测序平台 NovaSeq 6000 进行测序。

3. 生物信息分析流程

lncRNA生物信息分析流程

分析流程示意图如下:

流程图1 lncRNA生物信息分析流程图

miRNA生物信息分析流程

分析流程示意图如下:

流程图2 miRNA生物信息分析流程图

全转录组关联分析流程

全转录组关联分析分为WTS和ceRNA两种分析流程。两种流程都基于全转录组测序和分析数据。WTS分析是基于差异分析比较组合之间差异表达的RNA及预测的靶分子，结果中只对差异表达RNA分子进行分析；ceRNA分析直接依照ceRNA假说，以miRNA为中心，分析样本间分子表达相关性，不依赖于差异分析。最终构建一个关联ncRNA和mRNA的细胞内调控网络。

分析流程示意图如下:

流程图3 关联分析生物信息分析流程图

3.1 原始与过滤数据质量控制

数据质控步骤结果文件夹结构：包含原始数据与过滤数据质量控制结果，sRNA(small RNA)和lncRNA的文件夹结构相同：

--fastqc

----[样本编号]_fastqc.html【fastqc单样本质控结果网页版】

----[样本编号]_fastqc.zip【fastqc单样本质控结果图片

--multiqc

----multiqc_data【multiqc多样本质控结果汇总】

------multiqc各项文件

----multiqc_report.html【multiqc多样本质控结果汇总网页版】

--fastq_stat

----fastq_stat

------[样本编号]_[1/2].fq_stat.txt【各样本质控统计结果表格】

------bases_content_lineplot.[pdf/png/svg]【测序碱基比例图】

------error_barplot.[pdf/png/svg]【测序碱基错误率分布图】

------filter_pieplot.[pdf/png/svg]【各样本过滤结果饼图】

--fastq_stat.csv【各样本质控统计结果汇总表】

--fastq_stat.html【各样本质控统计结果汇总表网页版】

3.1.1 测序数据说明

全转录组测序包含Paired-End测序结果与Single-End测序结果，两者分别是对lncRNA(三文库包含ceRNA)的双端测序结果和sRNA(small RNA)单端测序结果。

测序得到的原始数据转化为序列数据，以 fastq 格式保存， FASTQ 文件是用于存储测序的核苷酸序列和其质量信息的测序生成的标准格式文件，它是序列格式中常见的一种。其核苷酸序列及其测序质量信息都是以 ASCII 码编码的，最初是由一代测序 sanger 开发的，目的是为了将 FASTA 序列与其质量信息合并为一个文件方便研究者查阅，现如今已成为高通量测序数据结果的标准。 FASTQ 文件是以四行为一个基本单元，并对应待测样本的其中一条测序信息， 如下所示：

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT

NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT

+

#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

图3.1.1 测序 fastq 文件格式

上述文件中第一行以“ @ ”开头，随后为 Illumina 测序标识符( Squence Identifiers )和描述文字；第二行是测序片段的碱基序列；第三行以“ + ”开头，随后为 Illumina 测序标识符(也可为空)；第四行是测序片段每个碱基相对应的测序质量值，该行每个字符对应的 ASCII 值减去 33 ，即为该碱基的测序质量值。

3.1.2 测序质量评估

测序过程本身存在机器错误的可能性，测序错误率分布检查可以反映测序数据的质量，序列信息中每个碱基的测序质量值保存在 fastq 文件中。如果测序错误率用 e 表示， Illumina 的碱基质量值用 Qphred 表示，则有： Qphred = -10log10(e) 。 Illumina Casava 1.8 版本碱基识别与 Phred 分值之间的简明对应关系见下表。

表3.1.1 Illumina Casava 1.8 版本碱基识别与Phred分值之间的简明对应关系
Phred分值 碱基错误识别概率 碱基正确识别概率 Q-Score
10 1/10 90% Q10
20 1/100 99% Q20
30 1/1000 99.9% Q30
40 1/10000 99.99% Q40

各样本碱基错误率如下图。其中横轴为碱基位置，纵轴为对应碱基错误率。

图3.1.2 测序碱基错误率分布图

各样本对应位置碱基比例如下图。其中横轴为碱基位置，纵轴为对应位置碱基比例。

图3.1.3 测序碱基比例图

测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于 RNA-seq 技术，测序错误率分布具有两个特点：

(1) 测序错误率会随着测序序列（ Sequenced Reads ）长度的增加而升高，这是由于测序过程中化学试剂的消耗而导致的，并且为 Illumina 高通量测序平台都具有的特征。

(2) 前几个碱基的位置也会发生较高的测序错误率，是因为在 lncRNA-seq 建库过程中反转录需要随机引物，推测前几个碱基测序错误率较高的原因是随机引物碱基和 RNA 模版的不完全结合。

测序错误率分布检查用于检测在测序长度范围内，有无测序错误率异常的碱基位置。一般情况下，每个碱基位置的测序错误率都应该低于 0.5% 。

3.1.3 各样本原始测序数据质量评估

我们对各样本原始数据的序列数目、碱基数目、 Q20 、 Q30 、 GC 占比等指标进行统计，从而直观的对原始数据大小与质量进行描述并汇总在表中。

  1. Sample: 样本编号
  2. File: 数据文件名称
  3. Reads: 序列数目
  4. Bases: 碱基总数
  5. Q20: 符合 Q20 标准碱基数目占总碱基数的比例
  6. Q30: 符合 Q30 标准碱基数目占总碱基数的比例
  7. gc%: GC 百分比

3.1.4 各样本质控后测序数据质量评估

lncRNA数据的质量控制条件如下：

用fastp软件对每一个样本的测序数据做质控处理(参数：--detect_adapter_for_pe)，得到clean data用于下游分析，质控规则如下：

对整条 read 做保留或弃去处理:

1.当一条 read 中 N 碱基的个数达到一定的数量（默认为 5 个）时，弃去整一条 read 。

2.当一条 read 中低质量碱基（默认质量值低于 15 的碱基为低质量碱基）达到一定的比例（默认这个比例值为 40% ）时，弃去整一条 read 。

3.当一条 read 的碱基个数少于一定的数量（默认为 15 个碱基）时，弃去整一条 read 。

4. adapter 序列裁剪:默认情况下，对 pair-end 的 read1 和 read2 做两序列比对，查看比对结果中， read1 的右端是否会超出 read2 的右端， read2 的左端是否会超出 read1 的左端，如果超出的话，表明是测到了 adapter 序列，则会把超出的部分给剪掉。

sRNA(small RNA)数据质量控制条件如下：

用cutadapt软件对每一个样本的测序数据做质控处理(参数：-a "AAAAAAA" -a "GGGGGGG" -a "TTTTTTT" -a "CCCCCCC" -a AGATCGGAAGAGCACACGTCT -m 18 --quality-base 33 --discard-untrimmed)，得到clean data用于下游分析，质控规则如下：

1.去除含有超过连续7个相同碱基的序列。

2.adapter序列裁剪:去除序列中的测序adapter。

3.当序列的碱基个数少于18时，弃去整一条序列。

4.去除简直质量小于33的序列。

我们对各样本过滤后数据的序列数目、碱基数目、 Q20 、 Q30 、 GC 占比等指标进行统计，从而直观的原始数据大小与质量进行描述并汇总在表1.3中。

  1. Sample: 样本编号
  2. File: 数据文件名称
  3. Reads: 序列数目
  4. Bases: 碱基总数
  5. Q20: 符合 Q20 标准碱基数目占总碱基数的比例
  6. Q30: 符合 Q30 标准碱基数目占总碱基数的比例
  7. gc%: GC 百分比

3.1.5 各样本过滤结果统计

各样本过滤结果如下图。饼图各部分为因各原因过滤的序列数目与比例。

图3.1.4 各样本过滤结果饼图

3.2 参考序列比对分析

lncRNA数据比对与统计文件夹结构：

--[样本编号]

----css

----images_qualimapReport【各样本的qualimap软件比对结果图】

----qualimapReport.html【各样本的qualimap软件比对结果分析报告网页版】

----raw_data_qualimapReport【各样本的qualimap软件转录本覆盖度统计结果】

----rnaseq_qc_results.txt【各样本的qualimap软件统计结果统计表】

--map_pieplot.[pdf/png/svg]【样本比对区域统计饼图】

--mapping_stat.[pdf/png/svg]【样本比对统计柱状图】

--mapping_stat.csv【样本比对统计表】

--[样本编号]_hisat2.txt【Hisat2软件比对原始结果】

--map_stat

----[样本编号]_map_region.txt【样本比对区域统计表】

3.2.1 序列比对统计

测序片段是由mRNA随机打断而来，为了确定这些片段对应基因组上的哪些区域，我们使用HISAT2软件将质控后的序列比对到参考基因组上，生成.sam或.bam文件，里面包含过滤后序列与基因组的比对记录。

随后用Qualimap rnaseq软件流程评估比对情况(参数: --paired)

我们对序列比对到基因组上的数目进行统计，并对比对率进行统计。

  1. Sample: 样本编号
  2. Overall alignment rate: 总比对率
  3. Aligned 0 times: 未比对到基因组序列数目与比例
  4. Aligned exactly 1 time: 单一比对到基因组序列数目与比例
  5. Aligned >1 times: 多次比对到基因组序列数目与比例

各样本比对统计柱状图如下。其中横轴比对数目比例，纵轴为样本名，不同颜色代表比对到基因组的比对次数类型。

图3.2.1 样本比对统计柱状图

1.2 比对基因组区域统计

我们对序列比对到基因组上的外显子，内含子，基因间区的数目和比对率进行统计。

  1. Sample: 样本编号
  2. Exonic: 比对到外显子的数目和比例
  3. Intronic: 比对到内含子的数目和比例
  4. Intergenic: 比对到基因间区的数目和比例
  5. Overlappingexon: 比对到多个外显子区的数目与比例

各样本比对区域统计饼图如下。其中各部分为比对到基因组不同区域序列数目与比例

图3.2.2 样本比对区域统计饼图

3.3 转录本拼接与注释

lncRNA组装与筛选文件夹结构：

--[样本编号].gtf.gz【各样本组装结果的gtf注释】

--all_gene.fa.gz【全部组装基因的fasta文件】

--all_transcripts.gtf.gz【全部组装转录本的gtf注释】

--lncRNAFilter

----lncRNA_stat.[csv/html]【全部组装转录本的gtf注释】

----FeatureSignature

------lncRNA_exon_density.[pdf/png]【转录本外显子数目密度分布图】

------lncRNA_length_density[pdf/png]【转录本长度密度分布图】

------lncRNA_orf_density[pdf/png]【转录本ORF密度分布图】

----FilterStat

------Unclassified.gtf.gz【未分类的转录本gtf注释】

------Novel_mRNA_information.[csv/html]【新鉴定出的mRNA信息表】

------Novel_lncRNA_information.[csv/html]【新鉴定出的lncRNA信息表】

------Novel_mRNA.[gtf/fa].gz【新鉴定出的mRNA的gtf注释与fasta序列】

------Known_mRNA.[gtf/fa].gz【已知的mRNA的gtf注释与fasta序列】

------Novel_lncRNA.[gtf/fa].gz【新鉴定出的lncRNA的gtf注释与fasta序列】

------Known_lncRNA.[gtf/fa].gz【已知的lncRNA的gtf注释与fasta序列】

------Novel_lncRNA_classification.[txt/png/pdf]【新鉴定的lncRNA类型分布图】

------Coding_potential_filter.result.[csv/html/png/pdf]【新鉴定的转录本编码潜能结果】

3.3.1 lncRNA组装

Stringtie可准确较为精准的拼接转录本以及实现转录本定量，我们使用Stringtie软件在基于比对到基因组上的结果，将reads拼接成转录本并进行定量[1]。样本组装结果位于以样本名称命名的gtf文件中。转录本序列与注释结果位于all_gene.fa.gz和all_transcripts.gtf.gz文件中。

图3.3.1 Stringtie拼接示意图

表1.1.1 Stringtie 拼接结果展示（示例）
Seqname Source Feature Start End Score Strand Frame Attributes
1 StringTie transcript 3582 3813 1000 - . gene_id "STRG.1"; transcript_id "STRG.1.1"; cov "9.047414"; FPKM "0.886686"; TPM "3.803272";
1 StringTie exon 3582 3813 1000 - . gene_id "STRG.1"; transcript_id "STRG.1.1"; exon_number "1"; cov "9.047414";
1 StringTie transcript 139827 140124 1000 - . gene_id "STRG.2"; transcript_id "STRG.2.1"; cov "2.919463"; FPKM "0.286120"; TPM "1.227258";
1 StringTie exon 139827 140124 1000 - . gene_id "STRG.2"; transcript_id "STRG.2.1"; exon_number "1"; cov "2.919463";
1 StringTie transcript 179469 179894 1000 - . gene_id "STRG.3"; transcript_id "STRG.3.1"; cov "2.816901"; FPKM "0.276069"; TPM "1.184144";
1 StringTie exon 179469 179894 1000 - . gene_id "STRG.3"; transcript_id "STRG.3.1"; exon_number "1"; cov "2.816901";
... ... ... ... ... ... ... ... ...
  1. Seqname：染色体名称;
  2. Source：来源
  3. Feature：序列类型
  4. Start：起始位点
  5. End：终止位点
  6. Score：minor FPKM/major FPKM
  7. Strand：正负链信息
  8. Frame：编码起始位点
  9. Attributes: 序列的其他描述信息

3.3.2 lncRNA过滤

我们将组装转录本与已知数据库比较，将转录本分成已知转录本和未知转录本两组。

未知转录本中，我们根据lncRNA的结构特点：1.长度>200nt，2.不编码蛋白的特点设置了严格的筛选条件，将未知转录本分成Novel_lncRNA和Novel_mRNA两组。

1：外显子数目：选择外显子数目大于等于2的转录本；

2：转录本长度：选择转录本长度大于200bp的转录本；

3：已知注释：使用gffcompare软件,去除与数据库注释外显子区域有重叠的转录本；

4：编码潜能：我们使用CPC2，CNCI，PFAM三款软件对编码潜能进行分析，将结果中没有编码潜能的转录本的交集作为Novel_lncRNA，有编码潜能的转录本的交集Novel_mRNA；

5：根据HGNC进行最终筛选并命名Novel_lncRNA。

lncRNA过滤结果位于lncRNAFilter文件夹中，我们对已知和未知lncRNA数目进行统计。

  1. Types: lncRNA种类，分成Known和Novel。
  2. Total: lncRNA数目

3.3.3 编码潜能分析

我们使用CPC2，CNCI，PFAM三款软件对编码潜能进行分析，我们对转录本进行统计。绘制维恩图，直观的展示各个方法预测出的lncRNA共有和特有的数目，每款软件所代表的圆圈数字之和代表鉴定出的lncRNA数目，交叠的部分表示软件之间共同鉴定出的lncRNA。我们取三款软件预测结果的交集，作为预测得到的lncRNA数据集。

图3.3.2 软件预测lncRNA数目维恩图

  1. transcript_id:转录本ID
  2. CPC:CPC软件鉴定结果
  3. CNCI:CNCI软件鉴定结果
  4. PFAM:PFAM软件鉴定结果

3.3.4 新鉴定lncRNA类型分析

我们将新鉴定的lncRNA根据HGNC数据库及其与已知的mRNA的位置关系分为4种类型：

1.antisense：转录本与相反链的蛋白质编码基因的外显子重叠，或已知反向调控编码区基因；

2.lincRNA：位于基因间区的转录本；

3.sense overlapping：转录本与相同链的蛋白质编码基因的外显子重叠；

4.sense intronic：位于编码基因的内含子区且不与同链上的外显子重叠。

  1. 第1列:lncRNA类型
  2. 第2列:lncRNA数目

饼图展示了四张类型lncRNA的占比情况。其中lincRNA为位于基因间区的lncRNA，antisense为反义链的lncRNA，sense_overlapping为与外显子重叠的lncRNA，sense_intronic为与内含子重叠的lncRNA。

图3.3.3 lncRNA 类型分布图

3.3.5 lncRNA在转录本长度，外显子数目，orf长度特征

我们将组装结果与物种已知的基因注释结果进行比较，筛选出Annotated mRNA和Annotated lncRNA。我们再将新鉴定出的转录本根据编码潜能分成Novel mRNA和Novel lncRNA。我们对lncRNA与mRNA在转录本长度、外显子数目以及ORF长度的对比，从lncRNA与mRNA的差别可以验证我们预测得到的lncRNA是否符合lncRNA的特征。ORF是基因的正常核苷酸序列，从起始密码子到终止密码子的阅读框不允许存在使翻译中断的终止密码子。我们对lncRNA与mRNA的orf进行比较分析,将得到的orf序列转化为蛋白序列，得到长度分布。

我们对已知与未知lncRNA和mRNA的长度、外显子数目、orf长度进行统计，绘制密度分布图，横轴为转录本长度，外显子和orf长度，纵坐标为密度，颜色代表不同类型的转录本。

图3.3.4 已知与未知lncRNA和mRNA长度、外显子数目、orf长度分布图

3.3.6 新鉴定出的lncRNA和mRNA信息汇总

我们对新鉴定出来的lncRNA与mRNA的注释信息进行汇总。

1.新鉴定lncRNA信息汇总

  1. Novel_lncRNA_ID:新鉴定lncRNA ID
  2. Novel_lncRNA_Gene_ID:新鉴定lncRNA所属基因ID
  3. Novel_lncRNA_Gene_Name:新鉴定lncRNA所属基因名称
  4. Gene_Type:基因类型
  5. Status:转录本类型
  6. Chromosome:染色体名称
  7. Start:起始位点
  8. End:终止位点
  9. Strand:链特异性
  10. Exon_number:外显子数目
  11. Length:转录本长度
  12. ORF:ORF长度

1.新鉴定mRNA信息汇总

  1. Novel_mRNA_ID:新鉴定RNA ID
  2. Novel_mRNA_Gene_ID:新鉴定mRNA所属基因ID
  3. Novel_mRNA_Gene_Name:新鉴定mRNA所属基因名称
  4. Gene_Type:基因类型
  5. Status:转录本类型
  6. Chromosome:染色体名称
  7. Start:起始位点
  8. End:终止位点
  9. Strand:链特异性
  10. Exon_number:外显子数目
  11. Length:转录本长度
  12. ORF:ORF长度

3.4 lncRNA定量分析

lncRNA定量分析文件夹结构：

--1.Count

----anno_genes.FPKM.csv【基因表达FPKM表】

----mRNA_transcripts.readcount.csv【mRNA原始序列数目表】

----genes.readcount.csv【基因原始序列数目表】

----lnc_transcripts.readcount.csv【lncRNA原始序列数目表】

----mRNA_genes.readcount.csv【mRNA所在的基因原始序列数目表】

----anno_transcripts.FPKM.csv【转录本表达FPKM表】

----lnc_transcripts.FPKM.csv【lncRNA表达FPKM表】

----mRNA_gene.FPKM.csv【mRNA所在的基因FPKM表】

----mRNA_transcripts.FPKM.csv【mRNA转录本FPKM表】

----transcripts.readcount.csv【转录本原始序列数目表】

--2.Correlation

----lncRNA

------correlation.[csv/pdf/png/svg]【根据lncRNA表达样本相关性热图】

------PCA2D.[pdf/png]【根据lncRNA表达样本二维降维结果图】

------PCA3D.[pdf/png]【根据lncRNA表达样本三维降维结果图】

----mRNA

------correlation.[csv/pdf/png/svg]【根据mRNA表达样本相关性热图】

------PCA2D.[pdf/png]【根据mRNA表达样本二维降维结果图】

------PCA3D.[pdf/png]【根据mRNA表达样本三维降维结果图】

--3.Distribution

----Sample_compare.boxplot.[pdf/png]【FPKM分布箱线图】

----Sample_compare.density.[pdf/png]【FPKM分布密度图】

----Sample_compare.violin.[pdf/png]【FPKM分布提琴图】

3.4.1 lncRNA定量分析

我们使用StringTie软件对已知转录本和预测得到的转录本进行定量，StringTie应用网络流算法以及可选的从头组装来拼接转录本。

为了去除测序深度与转录本长度导致的表达偏差，使基因表达水平具有可比性，我们计算FPKM值，即每百万片段中来自某转录本每千碱基的片段数目。FPKM 归一化了测序深度和转录本长度对片段数目的影响，是目前较为常用表达水平估算方法。

图3.4.1 Stringtie拼接与定量步骤示意图

3.4.2 转录本与基因定量结果

1.基因水平Raw counts原始序列数目统计，genes.readcount.csv文件包括lncRNA和mRNA的原始序列数目的记录。

  1. geneID:基因ID
  2. [样本名称]:对应样本的原始序列数目

2.转录本水平Raw counts原始序列数目统计，transcripts.readcount.csv文件包括lncRNA和mRNA的原始序列数目的记录。

  1. transcriptID:转录本ID
  2. [样本名称]:对应样本的原始序列数目

3.基因水平FPKM

  1. gene_id:基因ID
  2. gene_name:基因名称
  3. [样本名称]:对应样本的FPKM表达量

4.转录本水平FPKM

  1. transcript_id:转录本ID
  2. transcript_name:转录本名称
  3. gene_id:基因ID
  4. gene_name:基因名称
  5. [样本名称]:对应样本的FPKM表达量

转录本可以分成lncRNA和mRNA两种，其对应的raw counts和FPKM值矩阵是上述矩阵的子集。

3.4.2 基因表达水平分布

我们统计了各样本中基因的FPKM分布，通过箱线图、提琴图和密度图展示表达水平的分布情况。箱线图、提琴图的横轴为样本，纵轴为对数化的FPKM值，不同的颜色代表不同的样本，盒形图对应最大值、上四分位数、中值、下四分位数和最小值。密度图横轴为对数化的FPKM值，纵坐标为密度，颜色代表不同的样本。

图3.4.2 样本表达分布箱线图、提琴图和密度图

3.4.3 相关性分析

3.4.3.1 相关性分析

生物学重复通常是所有生物学实验所必须的，目前主流期刊基本都会要求有生物学重复。生物学重复主要具有证明实验操作是可重复的与保证下游分析结果是可靠的两种作用。样本间的基因表达水平相关性是检验实验或样本选择可靠性的重要指标，同一分组的重复之间相似性应大于不同分组间的样本相似性。

我们通过基因表达量计算了样本之间的R2值，并通过R2值代表样本之间的相似性。我们对lncRNA和mRNA分别进行相关性分析。两种RNA的结果格式是相同的。

  1. 横纵轴代表样本名称，数值代表相关性R2值

我们通过热图展示样本间的相关性，横纵轴代表样本名称，颜色与数值代表样本相关性R2值。

图3.4.3 样本相关性热图

3.4.3.2 主成分分析

主成分分析（Principal Components Analysis ,PCA）也常用于评估组间差异及组内样本重复情况，PCA采用线性降维的方法，对数以万计的基因变量进行降维并提取主成分，从而可以很好的反映样本之间的关系。我们对所有样本的基因表达值进行PCA分析，并使用第1，2，3主成分绘制PCA图。理想条件下，组内重复应更为相似，在图中会聚在一起，而组间的样本的样本相似度不如组内那么高，倾向于不聚在一起。

我们分别在二维和三维层面上绘制PCA结果。二维PCA图展示三个主成分维度两两组合的点图，三维PCA图在三维展示样本聚类情况，不同形状和颜色代表不同的分组与样本。

图3.4.4 PCA图

3.5 lncRNA文库可变剪接分析

lncRNA文库可变剪接分析文件夹结构：

--1.daslist

----比较条件

------[比较条件]_[SE/RI/A3SS/A5SS/MXE].txt【各比较条件下五种可变剪接类型统计结果表】

------[比较条件]_[SE/RI/A3SS/A5SS/MXE]_significant.txt【显著变化的各比较条件下五种可变剪接类型统计结果表】

--2.dasplot

----比较条件

------SE

--------[SE可变剪接事件].[pdf/png]【显著变化的SE类型各样本sashimi图】

------RI

--------[RI可变剪接事件].[pdf/png]【显著变化的RI类型各样本sashimi图】

------MXE

--------[MXE可变剪接事件].[pdf/png]【显著变化的MXE类型各样本sashimi图】

------A3SS

--------[A3SS可变剪接事件].[pdf/png]【显著变化的A3SS类型各样本sashimi图】

------A5SS

--------[A5SS可变剪接事件].[pdf/png]【显著变化的A5SS类型各样本sashimi图】

3.5.1 可变剪接事件鉴定

可变剪切是pre-RNA经过剪接后形成不同的mRNA的过程，是基因多样性的重要机制。对生命活动过程有重要影响，因此近些年来可变剪切分析逐渐成为热点。

可变剪切差异分析包括可变剪切事件定量及表达差异显著性分析。我们使用rMATS进行可变剪切分析。该软件可以鉴定五种可变剪切事件，分别为SE（外显子跳跃）、RI（内含子滞留）、MXE（互斥外显子）、A5SS（差异5’剪切位点）、A3SS（差异3’剪切位点）。

图3.5.1 rMATS分析可变剪切的五种类型

SNP和INDEL分析是结构分析的重要内容，对肿瘤等研究具有重要意义。我们通过picard-tools等工具对比对结果进行染色体坐标排序、去重等处理，最后通过GATK分别进行SNP Calling和Indel Calling，并对原始结果进行过滤，得到分析结果。

变异位点主要分为SNP与INDEL，下表展示的是INDEL结果。SNP结果格式类似。

  1. GeneID:可变剪接所在基因编号
  2. geneSymbol:可变剪接所在基因名称
  3. chr:染色体名称
  4. strand:链特异性
  5. exonStart_0base:外显子起始位置
  6. exonEnd:外显子终止位置
  7. upstreamES:上游exon起始位置
  8. upstreamEE:上游exon终止位置
  9. downstreamES:下游exon起始位置
  10. downstreamEE:下游exon终止位置
  11. IJC_SAMPLE_1:Inclusion Isoform表达量
  12. SJC_SAMPLE_1:Skipping Isoform表达量
  13. IJC_SAMPLE_2:nclusion Isoform表达量
  14. SJC_SAMPLE_2:Skipping Isoform表达量
  15. IncFormLen:Inclusion Isoform有效长度
  16. SkipFormLen:Skipping Isoform有效长度
  17. PValue:显著性p值
  18. FDR:显著性FDR值
  19. IncLevel1:比较组1 Inclusion Isoform的占比
  20. IncLevel2:比较组2 Inclusion Isoform的占比
  21. IncLevelDifference:差异

3.5.2 差异剪接可视化

我们对显著差异的可变剪接进行可视化展示，SE类型如图所示：标题为可变剪接三个外显子所在的染色体坐标及正负链信息，跨外显子比对测序数使用连接外显子junction边界的弧线表示。弧线的粗细和比对到junction上的reads数成正比，同时弧线上的数字指出了junction reads的数目，右上方标注了各个样本可变剪接事件所在的基因及IncLevel值。红色和橘黄色分别表示不同的样本，inclusion level值代表exon Inclusion Isoform在两个Isoform总表达量的比值，可以直接看出不同样本中可变剪接事件的表达差异

图3.5.3 剪接可视化

3.6 lncRNA文库SNP和INDEL分析

lncRNA文库SNP和INDEL分析文件夹结构：

--1.snpsite

----all_sample_[INDEL/SNP].[txt/vcf]【全部样本INDEL和SNP突变位点统计表】

----[样本编号]_[INDEL/SNP].[txt/vcf]【各样本INDEL和SNP突变位点统计表】

--2.snpeff

----[样本编号]_[INDEL/SNP]_snpeff.txt【各样本INDEL和SNP突变位点SNPEFF注释表】

--3.stat

----INDEL_[impact/region].[pdf/png/svg/txt]【INDEL突变影响程度与区域柱状图】

----SNP_[function/impact/region].[pdf/png/svg/txt]【SNP突变功能，影响程度与区域柱状图】

3.6.1 SNP和INDEL鉴定与过滤

SNP是基因组上由单核苷酸多态性，其数量极多，多态性丰富。每一个SNP位点发生转换和颠换，二者之比为1:2。SNP在CG序列上出现最为频繁，而且多是C转换为T，CG中的C常为甲基化的，自发地脱氨后即成为胸腺嘧啶。一般而言，SNP是指变异频率大于1%的单核苷酸变异。

Indel是指样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。

SNP和INDEL分析是结构分析的重要内容，对肿瘤等研究具有重要意义。我们通过picard-tools等工具对比对结果进行染色体坐标排序、去重等处理，最后通过GATK分别进行SNP Calling和Indel Calling，并对原始结果进行过滤，得到分析结果。

变异位点主要分为SNP与INDEL，下表展示的是INDEL结果。SNP结果格式类似。

vcf后缀文件为vcf格式变异位点文件。txt后缀文件为变异位点提取的基因型信息。

  1. CHROM:染色体名称
  2. POS:位点位置
  3. ID:变异位点的ID，默认使用.
  4. REF:参考碱基
  5. ALT:突变碱基
  6. GT(C1,C2,C3,Y1,Y2,Y3):各样本的基因型，以冒号隔开
  7. AD(C1,C2,C3,Y1,Y2,Y3):各样本基因型分别支持REF和ALT的测序深度
  8. DP(C1,C2,C3,Y1,Y2,Y3):各样本在变异位点的测序深度

3.6.2 SNP和INDEL功能注释

我们使用snpEff软件对鉴定出来的突变位点进行注释，并评估其功能及突变的影响程度。

  1. CHROM:染色体名称
  2. POS:位点位置
  3. ID:变异位点的ID，默认使用.
  4. REF:参考碱基
  5. ALT:突变碱基
  6. DP:变异位点的测序深度
  7. AD:基因型分别支持REF和ALT的测序深度
  8. GT:基因型，以冒号隔开
  9. GeneID:变异位点所在的基因编号
  10. GeneName:变异位点所在的基因名称
  11. FeatureID:变异位点所在的转录本编号
  12. Biotype:变异位点所在的转录本类型
  13. HGVS_C:HGVS注释
  14. HGVS_P:在蛋白水平的HGVS注释
  15. EFFECT:变异位点的影响
  16. IMPACT:变异位点的影响程度

3.6.3 SNP和INDEL位点统计分析

我们对变异位点功能，变异位点区域和变异位点影响对位点的重要程度进行统计和绘图，变异位点功能是从同义突变、错义突变、无义突变三个角度分析。变异位点区域是从EXON、INTRON、INTERGENIC等基因结构区域角度进行分析。变异位点影响从HIGH、MODERATE、LOW、MODIFIER（自身无表型效应，和别的突变位点同时存在才会产生影响）四个角度进行分析。柱状图横坐标为样本，纵坐标为突变数目，不同的颜色代表不同的突变的类型，位置，严重程度等。

图3.6.1 变异位点功能，区域和影响柱状图

3.7 lncRNA文库差异分析

lncRNA文库差异分析文件夹结构：

--比较条件

----lncRNA

------lncrna_all_results.csv【lncRNA差异比较结果统计表】

------diff_lncrna_results.csv【lncRNA显著差异比较结果统计表】

------diff_lncrna_count.csv【lncRNA显著差异原始序列数目表】

------diff_lncrna_[up/down].csv【显著上调/下调lncRNA结果统计表】

------diff_lncrna_heatmap.[pdf/png/svg]【lncRNA显著差异比较结果热图】

------lncRNA_MA_plot.[pdf/png/svg]【lncRNA MA图】

------lncRNA_volcano_plot.[pdf/png/svg]【lncRNA火山图】

------diff_lncrna_list.txt【lncRNA显著差异表】

----mRNA

------gene_all_results.csv【mRNA差异比较结果统计表】

------diff_gene_results.csv【mRNA显著差异比较结果统计表】

------diff_gene_count.csv【mRNA显著差异原始序列数目表】

------diff_gene_[up/down].csv【显著上调/下调mRNA结果统计表】

------mRNA_diff_gene_heatmap.[pdf/png/svg]【mRNA显著差异比较结果热图】

------mRNA_MA_plot.[pdf/png/svg]【mRNA MA图】

------mRNA_volcano_plot.[pdf/png/svg]【mRNA火山图】

------diff_gene_list.txt【mRNA显著差异表】

3.7.1 差异分析

在全转录组分析中，确定不同样本组间的基因表达量和lncRNA表达量是否有差异是分析的核心内容之一。差异分析主要分成三个步骤。

1.对原始的read count进行归一化处理。

2.通过统计学模型进行假设检验，获得存在显著差异的基因并给出p-value。

3.进行多重假设检验校正，获得经校正的存在显著差异的基因。

我们一般要求每个组中均具有生物学重复，并使用DESeq2软件进行基因表达差异显著性分析。一般来说，如果一个基因在两组样品中表达差异(FC，Fold Change)达到2倍以上，我们认为这样的基因是存在表达差异的。而为了判断这种差异究竟是由于各种误差导致还是本质差异，我们需要对所有基因在这两个样本中的表达量数据进行假设检验。而转录组分析是针对成千上万个基因的，随着分析的基因数的增多，假阳性率也在累积增多，这时需要对假设检验的p-value进行校正（padj值）控制假阳性的比例。差异倍数和padj值的选择可以根据实际情况进行调整。经验值（差异大于2倍且padj小于0.05）并不能完全适用所有的情况，当按照经验值筛选后得到的基因数目过少很有可能导致下游的功能富集分析没有显著性结果，可以适当放宽筛选标准。如果得到的差异基因数目过多，则需使用更加严格的阈值标准进行筛选。

3.7.2 差异基因/lncRNA列表

我们对每个比较组合进行差异分析并生成差异表达基因列表和差异表达lncRNA列表

  1. lncrnaid:lncrna ID
  2. baseMean:平均表达水平
  3. log2FoldChange:log2差异倍数
  4. lfcSE:差异倍数标准差
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. padj:矫正后显著性p值
  8. lncRNA_Gene_Name:lncRNA所属基因名称
  9. [样本编号]:样本原始序列数目

3.7.3 差异基因/lncRNA统计可视化

我们使用MA图和火山图来直观展示每个组合比较的差异基因分布情况，从而展示了基因丰度、变化幅度与统计显著性之间的关系。

一般来说，丰度较大，变化明显且显著性较高的基因可以作为后续分析与实验研究的靶点。MA图展示的是利用两组归一化后表达量的均值和差异倍数之间的关系，靠近右下和左上的基因即为丰度较高且变化幅度较大的基因。火山图展示的是差异倍数与padj值之间的关系，靠近左上角和右上角的基因即为统计显著性较强且变化幅度较大的基因。

MA图横坐标为归一化后表达量的均值，纵坐标为差异倍数(log2FC)，红色和蓝色分别代表显著上调和显著下调的基因，灰色为不显著变化基因。标签标注的点为显著变化排名前20的基因名称。

图3.7.1 差异基因MA图

火山图横坐标为差异倍数，纵坐标为显著性，红色和蓝色分别代表显著上调和显著下调的基因，灰色为不显著变化基因。标签标注的点为显著变化排名前20的基因名称。

图3.7.2 差异基因火山图

3.7.4 差异基因/lncRNA聚类分析

将所有比较组的差异基因取并集之后作为差异基因集。两组以上的实验可以对差异基因集进行聚类分析，将表达模式相近的基因聚在一起，并用热图的形式进行展示，我们按照显著性排序，并按顺序绘制排名前25、全部显著差异基因表达热图.我们采用层次聚类的方法对表达量进行聚类分析，并取每个样本表达量与样本表达量均值之差显示表达量的变化。

热图横坐标为不同样本，纵坐标为不同的基因，取显著性前25的基因绘制热图，图中值为每个样本表达量与样本表达量均值之差，左侧的聚类信息显示了表达模式相似的差异基因。

图3.7.3 显著差异lncRNA聚类热图

3.8 lncRNA 靶基因预测

lncRNA靶基因预测分析文件夹结构：

--比较条件

----diff_lncrna_co_expression_target.txt【lncRNA共表达调控靶基因】

----diff_lncrna_co_location_target.txt【lncRNA邻近调控靶基因结果】

3.8.1 lncRNA靶基因预测

lncRNA在转录水平和转录后水平上广泛参与生物体各类重要生理过程，对靶基因的表达进行调控。lncRNA调控靶基因的机制多种多样，最常见的对靶基因的方式是邻近调控和共表达调控。邻近调控是指lncRNA对附近基因的调控作用，我们将lncRNA上下游100kb内的基因作为邻近靶基因。共表达调控是指lncRNA转录后作用于较远的靶基因的作用方式，我们对多个样本间的表达相关性进行分析，在样本数大于5时选择相关性大于0.95的基因作为靶基因。

邻近调控靶基因结果如下表

  1. lncrna_id:lncrna ID
  2. baseMean:平均表达水平
  3. log2FoldChange:log2差异倍数
  4. lfcSE:差异倍数标准差
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. padj:矫正后显著性p值
  8. Chromosome:靶基因所属染色体编号
  9. lncRNA_Gene_ID:lncRNA所属基因
  10. lncRNA_Gene_Symbol:lncRNA所属基因Symbol名称
  11. lncRNA_Status:已知/新鉴定的lncRNA
  12. lncRNA_Start:lncRNA起始位点
  13. lncRNA_End:lncRNA终止位点
  14. lncRNA_Strand:lncRNA链特异性
  15. Gene_Id:靶基因ID
  16. Gene_Symbol:靶基因symbol名称
  17. Gene_Start:靶基因起始位点
  18. Gene_End:靶基因终止位点
  19. Gene_Strand:靶基因链特异性
  20. Distance:lncRNA与靶基因距离
  21. Location:靶基因与lncRNA位置关系

共表达调控靶基因结果如下表

  1. lncrna_id:lncrna ID
  2. baseMean:平均表达水平
  3. log2FoldChange:log2差异倍数
  4. lfcSE:差异倍数标准差
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. padj:矫正后显著性p值
  8. Chromosome:靶基因所属染色体编号
  9. lncRNA_Gene_ID:lncRNA所属基因
  10. lncRNA_Gene_Symbol:lncRNA所属基因Symbol名称
  11. Status:已知/新鉴定的lncRNA
  12. mRNA_Gene_Id:靶基因ID
  13. mRNA_Gene_Symbol:靶基因symbol名称
  14. Pearson_correlation:lncRNA与靶基因表达相关性
  15. Pvalue:共表达相关显著性

3.9 差异基因功能富集分析

lncRNA文库显著差异基因富集分析文件夹结构：

--比较条件

----Gene_All

------All_diff_GO_enrich_results.txt【全部显著基因GO富集分析结果表】

------All_diff_KEGG_enrich_results.txt【全部显著基因KEGG富集分析结果表】

------All_diff_GOenrich_barplot.[pdf/png/svg]【全部显著基因GO富集分析柱状图】

------All_diff_GOenrich_dotplot.[pdf/png/svg]【全部显著基因GO富集分析点图】

------All_diff_KEGGenrich_barplot.[pdf/png/svg]【全部显著基因KEGG富集分析柱状图】

------All_diff_KEGGenrich_dotplot.[pdf/png/svg]【全部显著基因KEGG富集分析点图】

------All_diff_GOenrich_mix_barplot.[pdf/png/svg]【全部显著基因GO富集分析分组柱状图】

------All_diff_GOenrich_mix_dotplot.[pdf/png/svg]【全部显著基因GO富集分析分组点图】

------All_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目结果表】

------All_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目结果表】

------All_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目柱状图】

------All_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目点图】

------All_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目柱状图】

------All_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目点图】

------All_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组柱状图】

------All_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组点图】

----Gene_Up【上调基因富集分析结果，结构同上】

------Up_diff_GO_enrich_results.txt

------Up_diff_KEGG_enrich_results.txt

------Up_diff_GOenrich_barplot.[pdf/png/svg]

------Up_diff_GOenrich_dotplot.[pdf/png/svg]

------Up_diff_KEGGenrich_barplot.[pdf/png/svg]

------Up_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Up_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Gene_Down【下调基因富集分析结果，结构同上】

------Down_diff_GO_enrich_results.txt

------Down_diff_KEGG_enrich_results.txt

------Down_diff_GOenrich_barplot.[pdf/png/svg]

------Down_diff_GOenrich_dotplot.[pdf/png/svg]

------Down_diff_KEGGenrich_barplot.[pdf/png/svg]

------Down_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Down_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Diff_Gene_Term.txt【差异基因表】

3.9.1 富集分析介绍

我们对差异靶基因进行富集分析。通过富集分析，可以找到不同条件下的差异靶基因与哪些生物学功能或通路显著性相关，从而揭示和理解生物学过程的基本分子机制。

我们采用clusterProfiler软件对差异靶基因集进行GO功能富集分析，KEGG通路富集分析等。富集分析基于超几何分布原理，其中差异靶基因集为差异分析所得基因并注释到GO或KEGG数据库的基因集，背景基因集为所有参与差异分析并注释到GO或KEGG数据库的基因集。富集分析结果是对每个差异比较组合的所有差异基因集进行富集。

图3.9.1 富集分析原理图

3.9.2 GO富集分析

GO(Gene Ontology)是描述基因功能的综合性数据库（http://www.geneontology.org/)，可分为生物过程（Biological Process）和细胞组成（Cellular Component）分子功能（Molecular Function）三个部分。GO 功能显著性富集分析给出与基因集背景相比，在差异表达基因中显著富集的 GO 功能条目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库的各个 term 映射，计算每个 term 的基因数目，然后找出与整个基因集背景相比，在差异表达基因中显著富集。GO富集分析方法为clusterProfiler的enricher函数，该方法使用超几何分布检验的方法获得显著富集的GO Term。统计被显著富集的各个GO term中的基因数，以柱状图的形式展示。GO富集以padj小于0.05为显著富集。

  1. ID:GO条目编号
  2. Description:GO条目名称
  3. GeneRatio:对应条目含有的显著差异基因数/显著差异总基因数
  4. BgRatio:对应条目基因数/总基因数
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. p.adjust:矫正后显著性p值
  8. qvalue:显著性q值
  9. geneID:条目含有的显著差异基因ID
  10. Count:条目含有的显著差异基因数目
  11. geneID_Symbol:条目含有的显著差异基因Symbol名称
  12. Term:GO条目编号
  13. Type:GO条目所属类型

根据GO富集分析结果，在GO的三个分类中（cellular_component:细胞组分；biological_process：生物学过程；molecular_function：分子功能）分别统计每个分类上调基因、下调基因及所有基因的数目，以柱状图和气泡图展示。

柱形图和气泡图是富集分析结果的图形化展示方式。在此图中，富集程度通过GeneRatio、p.adjust和富集到此通路上的基因个数等指标来衡量并可视化。其中GeneRatio指差异表达的基因中位于该pathway条目的基因数目与所有有注释基因中位于该pathway条目的基因总数的比值。p.adjust是做过多重假设检验校正之后的Pvalue，p.adjust的取值范围为[0，1]，越接近于零，表示富集越显著。我们挑选了给类型富集最显著的10条目在该图中进行展示，若富集的条目不足20条，则全部展示。

图3.9.2 GO富集柱状图

图3.9.3 GO富集气泡图

3.9.3 KEGG富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库，其中整合了基因组化学和系统功能信息等众多内容。信息分析时，我们将KEGG数据库按照动物，植物，真菌等进行分类，依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的pathway，KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。

KEGG的富集结果表格形式与GO富集结果相似，此处不再赘述。

图3.9.4 KEGG富集柱状图

图3.9.5 KEGG富集气泡图

3.10 差异 lncRNA 靶基因功能富集分析

lncRNA文库显著差异lncRNA靶基因富集分析文件夹结构：

--比较条件

----Gene_All

------All_diff_GO_enrich_results.txt【全部显著基因GO富集分析结果表】

------All_diff_KEGG_enrich_results.txt【全部显著基因KEGG富集分析结果表】

------All_diff_GOenrich_barplot.[pdf/png/svg]【全部显著基因GO富集分析柱状图】

------All_diff_GOenrich_dotplot.[pdf/png/svg]【全部显著基因GO富集分析点图】

------All_diff_KEGGenrich_barplot.[pdf/png/svg]【全部显著基因KEGG富集分析柱状图】

------All_diff_KEGGenrich_dotplot.[pdf/png/svg]【全部显著基因KEGG富集分析点图】

------All_diff_GOenrich_mix_barplot.[pdf/png/svg]【全部显著基因GO富集分析分组柱状图】

------All_diff_GOenrich_mix_dotplot.[pdf/png/svg]【全部显著基因GO富集分析分组点图】

------All_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目结果表】

------All_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目结果表】

------All_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目柱状图】

------All_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目点图】

------All_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目柱状图】

------All_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目点图】

------All_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组柱状图】

------All_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组点图】

----Gene_Up【上调基因富集分析结果，结构同上】

------Up_diff_GO_enrich_results.txt

------Up_diff_KEGG_enrich_results.txt

------Up_diff_GOenrich_barplot.[pdf/png/svg]

------Up_diff_GOenrich_dotplot.[pdf/png/svg]

------Up_diff_KEGGenrich_barplot.[pdf/png/svg]

------Up_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Up_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Gene_Down【下调基因富集分析结果，结构同上】

------Down_diff_GO_enrich_results.txt

------Down_diff_KEGG_enrich_results.txt

------Down_diff_GOenrich_barplot.[pdf/png/svg]

------Down_diff_GOenrich_dotplot.[pdf/png/svg]

------Down_diff_KEGGenrich_barplot.[pdf/png/svg]

------Down_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Down_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Diff_Gene_Term.txt【差异基因表】

3.10.1 富集分析介绍

介绍信息与3.9 差异基因功能富集分析相同

3.10.2 GO富集分析

GO(Gene Ontology)是描述基因功能的综合性数据库（http://www.geneontology.org/)，可分为生物过程（Biological Process）和细胞组成（Cellular Component）分子功能（Molecular Function）三个部分。GO 功能显著性富集分析给出与基因集背景相比，在差异表达基因中显著富集的 GO 功能条目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库的各个 term 映射，计算每个 term 的基因数目，然后找出与整个基因集背景相比，在差异表达基因中显著富集。GO富集分析方法为clusterProfiler的enricher函数，该方法使用超几何分布检验的方法获得显著富集的GO Term。统计被显著富集的各个GO term中的基因数，以柱状图的形式展示。GO富集以padj小于0.05为显著富集。

  1. ID:GO条目编号
  2. Description:GO条目名称
  3. GeneRatio:对应条目含有的显著差异基因数/显著差异总基因数
  4. BgRatio:对应条目基因数/总基因数
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. p.adjust:矫正后显著性p值
  8. qvalue:显著性q值
  9. geneID:条目含有的显著差异基因ID
  10. Count:条目含有的显著差异基因数目
  11. geneID_Symbol:条目含有的显著差异基因Symbol名称
  12. Term:GO条目编号
  13. Type:GO条目所属类型

根据GO富集分析结果，在GO的三个分类中（cellular_component:细胞组分；biological_process：生物学过程；molecular_function：分子功能）分别统计每个分类上调基因、下调基因及所有基因的数目，以柱状图和气泡图展示。

柱形图和气泡图是富集分析结果的图形化展示方式。在此图中，富集程度通过GeneRatio、p.adjust和富集到此通路上的基因个数等指标来衡量并可视化。其中GeneRatio指差异表达的基因中位于该pathway条目的基因数目与所有有注释基因中位于该pathway条目的基因总数的比值。p.adjust是做过多重假设检验校正之后的Pvalue，p.adjust的取值范围为[0，1]，越接近于零，表示富集越显著。我们挑选了给类型富集最显著的10条目在该图中进行展示，若富集的条目不足20条，则全部展示。

图3.10.2 GO富集柱状图

图3.10.3 GO富集气泡图

3.10.3 KEGG富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库，其中整合了基因组化学和系统功能信息等众多内容。信息分析时，我们将KEGG数据库按照动物，植物，真菌等进行分类，依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的pathway，KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。

KEGG的富集结果表格形式与GO富集结果相似，此处不再赘述。

图3.10.4 KEGG富集柱状图

图3.10.5 KEGG富集气泡图

3.11 sRNA(small RNA)文库序列长度分析

sRNA(small RNA)文库序列长度分析文件夹结构：

--[样本编号]Read_Length_distribution.[pdf/png/svg]【sRNA(small RNA)片段的长度分布统计图】

--[样本编号]_readlength.txt.gz【sRNA(small RNA)文库序列片段的长度分布结果】

--[样本编号]_readlength_dis.txt【sRNA(small RNA)片段的长度分布统计表】

3.11.1 sRNA(small RNA)长度筛选

一般来说，动物sRNA(small RNA)的长度区间为 18~35nt，植物sRNA(small RNA)的长度区间为 18~30nt，长度分布的峰能帮助我们判断sRNA(small RNA) 的种类，如miRNA集中在22nt，siRNA集中在24nt等。长度分布结果如显示多数样本的序列在21~22nt出现了峰值，表明数据较为可靠。对各样品的过滤后序列的长度分布统计图如下。其中横坐标为序列最大碱基位置， 纵坐标为对应序列个数。

图3.11.1 sRNA(small RNA)片段的长度分布统计图

3.12 sRNA(small RNA)文库序列来源分析

sRNA(small RNA)文库序列来源分析文件夹结构：

--RNA_type_number.[pdf/png/svg]【sRNA(small RNA)来源统计图】

--sRNA(small RNA)_source_summary.txt【sRNA(small RNA)来源统计表】

--[样本编号]_filtered_collapsed.fasta【各样本去重后fasta序列文件】

--[样本编号]_rfam_filtered.fasta【各样本去重去除已知RNA后fasta序列文件】

3.12.1 sRNA(small RNA)来源分类注释

因为测序文库中的 sRNA(small RNA) 来源于不同类型的分子，因此需要弄清楚从细胞中不同类型的分子获得的读数的分布，例如mRNA，miRNA，tRNA，rRNA，snoRNA，repeats，siRNA和piRNA。在将读数与不同类型的分子比对后，通常重要的是计算读数和独特序列的长度分布。sRNA(small RNA)读数和独特序列的正确长度分布是sRNA(small RNA)测序文库的良好反映。如果在测序文库中改变sRNA(small RNA)的长度分布，则可能是由于sRNA(small RNA)生成途径中的关键基因突变或测序文库的不良品质引起的。下表是sRNA(small RNA)测序文库中的读段和独特序列的长度分布，以及映射到不同类型分子的读段和独特序列的分布。

  1. rna_type: RNA种类，包括lncrna，ribozyme，rrna，snorna，snrna，srna，trna，unannotated。
  2. [样本编号]: 样本编号

下图是sRNA(small RNA)测序文库中的sRNA来源统计图。横轴为样本，纵轴为不同类型RNA的序列数目。

图3.12.1 sRNA(small RNA)来源统计图

3.13 miRNA鉴定，分布，二级结构预测与碱基偏好性分析

miRNA鉴定，分布，二级结构预测与碱基偏好性分析文件夹结构：

--Known_position

----Known_miRNA-fst-bias.[pdf/png/svg]【已知miRNA首位碱基偏好柱状图】

----Known_miRNA-NT-bias.[pdf/png/svg]【已知miRNA碱基分布统计柱状图】

----[样本编号].position【已知miRNA碱基分布统计】

----[样本编号].firstbase【已知miRNA首位碱基偏好统计】

--Novel_position

----Novel_miRNA-fst-bias.[pdf/png/svg]【新鉴定出miRNA首位碱基偏好柱状图】

----Novel_miRNA-NT-bias.[pdf/png/svg]【新鉴定出miRNA碱基分布统计柱状图】

----[样本编号].position【新鉴定出miRNA碱基分布统计】

----[样本编号].firstbase【新鉴定出miRNA首位碱基偏好统计】

--Structure_plot

----[miRNA名称].[jpg/pdf]【miRNA结构图】

--hairpin_Known_Novel.fa【已知与新鉴定出的miRNA发夹结构fasta序列】

--hairpin_mature_Known_Novel.fa【已知与新鉴定出的成熟miRNA fasta序列】

--hairpin_mature_Known_Novel.pairs【已知与新鉴定出的成熟miRNA碱基配对结果】

--known_miRNA.map_stat.[txt/html]【已知miRNA比对统计结果】

--novel_miRNA.map_stat.[txt/html]【新鉴定出的miRNA比对统计结果】

--Known_Novel_star_readcount.txt【已知与新鉴定出的miRNA定量结果】

--mature_Known_Novel.fa【已知与新鉴定出的miRNA成熟fasta序列】

--miRBase_Known_Novel.mrd【miRBase已知与新鉴定的miRNA信息】

--miRNA.map_Known_Novel.fas【已知与新鉴定出的miRNA fasta序列】

--Novel_hairpin.pos【新鉴定出的miRNA发夹结构位置】

--Novel_hairpin.str【新鉴定出的miRNA发夹结构】

--Novel_star.fa【新鉴定出的miRNA star fasta序列】

3.13.1 miRNA样本比对统计

在鉴定 miRNA 时，我们将序列与miRBase中的注释进行比对，得到各样品匹配上的sRNA的详细情况。miRNA具有保守性，有些保守miRNA可能在一个物种中具有多个旁系同源物，利用miRNA的保守性和数据库注释，来分析miRNA。

  1. Types: 比对上的miRNA类型。Mapped mature: 比对上的miRNA成熟体，Mapped hairpin: 比对上的miRNA前体，Mapped uniq sRNA: 比对上的sRNA的种类。Mapped total sRNA: 比对上的sRNA的数目。
  2. Total: 总共匹配上的各类型数目。
  3. [样本编号]: 各样本匹配上的各类型数目。

3.13.2 miRNA基因表达水平分布

我们统计了各样本中miRNA的TPM分布，通过箱线图和密度图展示表达水平的分布情况。箱线图横轴为样本，纵轴为对数化的TPM值，不同的颜色代表不同的样本，盒形图对应最大值、上四分位数、中值、下四分位数和最小值。密度图横轴为对数化的TPM值，纵坐标为密度，颜色代表不同的样本。

下图为miRNA茎环结构图，红色字母代表成熟miRNA。短线是碱基的匹配，突起的环则代表碱基不匹配。报告中给出的为结构示意图，具体结果位于Structure_plot文件夹中。

图3.13.1 miRNA茎环结构图（示意图）

我们对miRNA的细节进行描述，如下表：

miRNA annotation
>mmu-let-7a-1
total read count 430285
mmu-let-7a-5p read count 429223
mmu-let-7a-1-3p read count 1026
remaining read count 36
exp ffffffffffff5555555555555555555555fffffffffffffffffffffffffffff3333333333333333333333fffffffff
pri_seq uucacugugggaugagguaguagguuguauaguuuuagggucacacccaccacugggagauaacuauacaaucuacugucuuuccuaaggugau
pri_struct .(((((.(((((.(((((((((((((((((((((.....(((...((((....)))).))))))))))))))))))))))))))))).))))). #MM +/-
C2_479540_x2 ...........augagguaguagguuguauag.............................................................. 0 +
C3_292339_x2 ...........augagguaguagguuguauagu............................................................. 0 +
  1. miRNA名称
  2. total readcount: 比对上的miRNA数目
  3. miRNA名称-5p: 比对上5p的miRNA数目
  4. miRNA名称-3p: 比对上3p的miRNA数目
  5. remaining read count: 比对上前体但是未比对上成熟体的数目
  6. exp:miRNA前体注释，5代表5'区域，3代表3'区域，f代表前体上的其他区域
  7. pri_seq: 前体序列
  8. pri_struct: miRNA前体的二级结构
  9. 样本名字_序列编号_xn，n是该序列出现的次数

3.13.3 miRNA样本表达统计

我们对已知，新鉴定出来的和新鉴定出来的前体miRNA在各样本中的表达量进行统计。

  1. miRNA: miRNA编号
  2. [样本编号]: 各样本miRNA的表达序列数目。

3.13.4 已知miRNA和新鉴定出的miRNA首位碱基偏好及前22位碱基偏好分析

miRNA在由前体发育为成熟体时其过程是由Dicer酶切完成的，酶切位点的特异性使得miRNA成熟体序列首位碱基具有很强的偏向性，因此我们对不同长度miRNA的首位点碱基分布和各位点碱基分布统计。miRNA的5’端第一个碱基对U有很强的偏好性，对G则有抗性;第二到第四个碱基缺乏比其他位置碱基通常缺乏C(第四个碱基有时例外)。统计miRNA首位碱基偏好是检验数据质量的依据。

已知miRNA和新鉴定出的miRNA序列首位碱基偏好性图横轴为序列长度，纵轴为各碱基比例。

图3.13.2 已知miRNA首位碱基偏好柱状图

图3.13.3 新鉴定出miRNA首位碱基偏好柱状图

已知miRNA和新鉴定出的miRNA各位点碱基分布统计柱状图横轴为位置，纵轴为各碱基比例。

图3.13.4 已知miRNA碱基分布统计柱状图

图3.13.5 新鉴定出miRNA碱基分布统计柱状图

3.14 miRNA定量与降维分析

miRNA定量与降维分析文件夹结构：

--Quant

----Raw_counts.txt【miRNA原始序列数定量结果表】

----TPM.txt【miRNA TPM定量结果表】

----TPM_interval.txt【miRNA TPM定量结果表】

--cor_pearson.[txt/png/pdf/svg]【样本皮尔森相关系数表与图】

--miRNA_PCA.[pdf/png/svg]【样本PCA图】

--TPM_boxplot.[pdf/png/svg]【miRNA TPM箱线图】

--TPM_density_distribution.[pdf/png/svg]【miRNA TPM密度图】

3.14.1 miRNA定量分析

我们对已知，新鉴定出来的和新鉴定出来的前体miRNA在各样本中的表达量进行统计。

  1. miRNA: miRNA编号
  2. [样本编号]: 各样本miRNA的表达序列数目。

在原始数据的基础上，我们计算了TPM表达量，并对各样本不同表达水平的miRNA数目和比例进行统计。

  1. miRNA: miRNA编号
  2. [样本编号]: 各样本miRNA的TPM表达量。
  1. Category: miRNA表达量分组
  2. [样本编号]: 各组miRNA在各样本的数目和比例。

3.14.2 miRNA表达水平分布

我们统计了各样本中miRNA的TPM分布，通过箱线图和密度图展示表达水平的分布情况。箱线图的横轴为样本，纵轴为对数化的TPM值，不同的颜色代表不同的样本，盒形图对应最大值、上四分位数、中值、下四分位数和最小值。密度图横轴为对数化的TPM值，纵坐标为密度，颜色代表不同的样本。

图3.14.2 样本表达分布箱线图和密度图

3.14.3 相关性分析

3.14.3.1 相关性分析

生物学重复通常是所有生物学实验所必须的，目前主流期刊基本都会要求有生物学重复。生物学重复主要具有证明实验操作是可重复的与保证下游分析结果是可靠的两种作用。样本间的基因表达水平相关性是检验实验或样本选择可靠性的重要指标，同一分组的重复之间相似性应大于不同分组间的样本相似性。

我们通过基因表达量计算了样本之间的R2值，并通过R2值代表样本之间的相似性。我们对miRNA进行相关性分析。

  1. 横纵轴代表样本名称，数值代表相关性R2值

我们通过热图展示样本间的相关性，横纵轴代表样本名称，颜色与数值代表样本相关性R2值。

图3.14.3 样本相关性热图

3.14.3.2 主成分分析

主成分分析（Principal Components Analysis ,PCA）也常用于评估组间差异及组内样本重复情况，PCA采用线性降维的方法，对数以万计的基因变量进行降维并提取主成分，从而可以很好的反映样本之间的关系。我们对所有样本的基因表达值进行PCA分析，并使用第1，2，3主成分绘制PCA图。理想条件下，组内重复应更为相似，在图中会聚在一起，而组间的样本的样本相似度不如组内那么高，倾向于不聚在一起。

我们分别在二维PCA图展示三个主成分维度两两组合的点图，不同颜色代表不同的分组与样本。

图3.14.4 PCA图

3.15 miRNA差异分析

miRNA差异分析文件夹结构：

--比较条件

----miRNA_all_results.csv【miRNA差异分析结果表】

----diff_miRNA_results.csv【显著差异miRNA结果表】

----diff_miRNA_count.csv【显著差异miRNA原始序列数目表】

----diff_miRNA_[up/down].csv【显著上调/下调miRNA结果表】

----diff_miRNA_heatmap.[pdf/png/svg]【显著差异miRNA热图】

----miRNA_MA_plot.[pdf/png/svg]【miRNA差异分析MA图】

----miRNA_volcano_plot.[pdf/png/svg]【miRNA差异分析火山图】

3.15.1 差异分析

在全转录组分析中，确定不同样本组间的miRNA表达量是否有差异是分析的核心内容之一。差异分析主要分成三个步骤。

1.对原始的read count进行归一化处理。

2.通过统计学模型进行假设检验，获得存在显著差异的基因并给出p-value。

3.进行多重假设检验校正，获得经校正的存在显著差异的基因。

我们一般要求每个组中均具有生物学重复，并使用DESeq2软件进行基因表达差异显著性分析。一般来说，如果一个基因在两组样品中表达差异(FC，Fold Change)达到2倍以上，我们认为这样的基因是存在表达差异的。而为了判断这种差异究竟是由于各种误差导致还是本质差异，我们需要对所有基因在这两个样本中的表达量数据进行假设检验。而转录组分析是针对成千上万个基因的，随着分析的基因数的增多，假阳性率也在累积增多，这时需要对假设检验的p-value进行校正（padj值）控制假阳性的比例。差异倍数和padj值的选择可以根据实际情况进行调整。经验值（差异大于2倍且padj小于0.05）并不能完全适用所有的情况，当按照经验值筛选后得到的基因数目过少很有可能导致下游的功能富集分析没有显著性结果，可以适当放宽筛选标准。如果得到的差异基因数目过多，则需使用更加严格的阈值标准进行筛选。

3.15.2 差异基miRNA列表

我们对每个比较组合进行差异分析并生成差异miRNA列表

  1. miRNA:miRNA ID
  2. baseMean:平均表达水平
  3. log2FoldChange:log2差异倍数
  4. lfcSE:差异倍数标准差
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. padj:矫正后显著性p值
  8. symbol:miRNA ID
  9. [样本编号]:样本原始序列数目

3.15.3 差异miRNA统计可视化

我们使用MA图和火山图来直观展示每个组合比较的差异基因分布情况，从而展示了基因丰度、变化幅度与统计显著性之间的关系。

一般来说，丰度较大，变化明显且显著性较高的基因可以作为后续分析与实验研究的靶点。MA图展示的是利用两组归一化后表达量的均值和差异倍数之间的关系，靠近右下和左上的基因即为丰度较高且变化幅度较大的基因。火山图展示的是差异倍数与padj值之间的关系，靠近左上角和右上角的基因即为统计显著性较强且变化幅度较大的基因。

MA图横坐标为归一化后表达量的均值，纵坐标为差异倍数(log2FC)，红色和蓝色分别代表显著上调和显著下调的基因，灰色为不显著变化基因。标签标注的点为显著变化排名前20的基因名称。

图3.15.1 差异miRNA MA图

火山图横坐标为差异倍数，纵坐标为Deseq2软件给出的差异显著性值，红色和蓝色分别代表显著上调和显著下调的基因，灰色为不显著变化基因。标签标注的点为显著变化排名前20的基因名称。

图3.15.2 差异miRNA火山图

3.15.4 差异miRNA聚类分析

将所有比较组的差异基因取并集之后作为差异基因集。两组以上的实验可以对差异基因集进行聚类分析，将表达模式相近的基因聚在一起，并用热图的形式进行展示，我们按照显著性排序，并按顺序绘制排名前25、全部显著差异基因表达热图.我们采用层次聚类的方法对表达量进行聚类分析，并取每个样本表达量与样本表达量均值之差显示表达量的变化。

热图横坐标为不同样本，纵坐标为不同的基因，取显著性前25的基因绘制热图，图中值为每个样本表达量与样本表达量均值之差，左侧的聚类信息显示了表达模式相似的差异基因。

图3.15.3 显著差异miRNA聚类热图

3.16 miRNA靶基因预测分析

miRNA靶基因预测分析文件夹结构：

--比较条件

----diff_miRNA_target.txt【差异miRNA靶基因表】

--target_gene_description.txt【差异miRNA靶基因注释】

3.16.1 miRNA靶基因预测

动物中大部分miRNA与靶基因间完全或几乎完全互补配对，并在互补位点的第10或11位核苷酸上发生切割作用以调控靶基因的表达。切割后会产生2个片段，5'剪切片段和3'剪切片段。其中，5'剪切片段包含5'帽子结构和3'羟基;3'剪切片段则包含自由的5'单磷酸和3'polyA尾巴。有参动物miRNA靶基因预测为miRanda、RNAhybrid两个软件的交集，得到miRNA和靶基因间的对应关系。

  1. miRNA_id:miRNA ID
  2. baseMean:平均表达水平
  3. log2FoldChange:log2差异倍数
  4. lfcSE:差异倍数标准差
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. padj:矫正后显著性p值
  8. target_gene:miRNA靶基因

3.17 显著差异miRNA靶基因富集分析

显著差异miRNA靶基因富集分析文件夹结构：

--比较条件

----Gene_All

------All_diff_GO_enrich_results.txt【全部显著基因GO富集分析结果表】

------All_diff_KEGG_enrich_results.txt【全部显著基因KEGG富集分析结果表】

------All_diff_GOenrich_barplot.[pdf/png/svg]【全部显著基因GO富集分析柱状图】

------All_diff_GOenrich_dotplot.[pdf/png/svg]【全部显著基因GO富集分析点图】

------All_diff_KEGGenrich_barplot.[pdf/png/svg]【全部显著基因KEGG富集分析柱状图】

------All_diff_KEGGenrich_dotplot.[pdf/png/svg]【全部显著基因KEGG富集分析点图】

------All_diff_GOenrich_mix_barplot.[pdf/png/svg]【全部显著基因GO富集分析分组柱状图】

------All_diff_GOenrich_mix_dotplot.[pdf/png/svg]【全部显著基因GO富集分析分组点图】

------All_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目结果表】

------All_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目结果表】

------All_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目柱状图】

------All_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目点图】

------All_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目柱状图】

------All_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因KEGG富集分析显著条目点图】

------All_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组柱状图】

------All_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)【全部显著基因GO富集分析显著条目分组点图】

----Gene_Up【上调基因富集分析结果，结构同上】

------Up_diff_GO_enrich_results.txt

------Up_diff_KEGG_enrich_results.txt

------Up_diff_GOenrich_barplot.[pdf/png/svg]

------Up_diff_GOenrich_dotplot.[pdf/png/svg]

------Up_diff_KEGGenrich_barplot.[pdf/png/svg]

------Up_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Up_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Up_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Up_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Gene_Down【下调基因富集分析结果，结构同上】

------Down_diff_GO_enrich_results.txt

------Down_diff_KEGG_enrich_results.txt

------Down_diff_GOenrich_barplot.[pdf/png/svg]

------Down_diff_GOenrich_dotplot.[pdf/png/svg]

------Down_diff_KEGGenrich_barplot.[pdf/png/svg]

------Down_diff_KEGGenrich_dotplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_barplot.[pdf/png/svg]

------Down_diff_GOenrich_mix_dotplot.[pdf/png/svg]

------Down_diff_signGO_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGG_enrich_results.txt(如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signKEGGenrich_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_barplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

------Down_diff_signGOenrich_mix_dotplot.[pdf/png/svg](如果没有显著富集的通路则结果中没有该文件)

----Diff_Gene_Term.txt【差异基因表】

3.17.1 富集分析介绍

介绍信息与3.9 差异基因功能富集分析相同

3.17.2 GO富集分析

GO(Gene Ontology)是描述基因功能的综合性数据库（http://www.geneontology.org/)，可分为生物过程（Biological Process）和细胞组成（Cellular Component）分子功能（Molecular Function）三个部分。GO 功能显著性富集分析给出与基因集背景相比，在差异表达基因中显著富集的 GO 功能条目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库的各个 term 映射，计算每个 term 的基因数目，然后找出与整个基因集背景相比，在差异表达基因中显著富集。GO富集分析方法为clusterProfiler的enricher函数，该方法使用超几何分布检验的方法获得显著富集的GO Term。统计被显著富集的各个GO term中的基因数，以柱状图的形式展示。GO富集以padj小于0.05为显著富集。

  1. ID:GO条目编号
  2. Description:GO条目名称
  3. GeneRatio:对应条目含有的显著差异基因数/显著差异总基因数
  4. BgRatio:对应条目基因数/总基因数
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. p.adjust:矫正后显著性p值
  8. qvalue:显著性q值
  9. geneID:条目含有的显著差异基因ID
  10. Count:条目含有的显著差异基因数目
  11. geneID_Symbol:条目含有的显著差异基因Symbol名称
  12. Term:GO条目编号
  13. Type:GO条目所属类型

根据GO富集分析结果，在GO的三个分类中（cellular_component:细胞组分；biological_process：生物学过程；molecular_function：分子功能）分别统计每个分类上调基因、下调基因及所有基因的数目，以柱状图和气泡图展示。

柱形图和气泡图是富集分析结果的图形化展示方式。在此图中，富集程度通过GeneRatio、p.adjust和富集到此通路上的基因个数等指标来衡量并可视化。其中GeneRatio指差异表达的基因中位于该pathway条目的基因数目与所有有注释基因中位于该pathway条目的基因总数的比值。p.adjust是做过多重假设检验校正之后的Pvalue，p.adjust的取值范围为[0，1]，越接近于零，表示富集越显著。我们挑选了给类型富集最显著的10条目在该图中进行展示，若富集的条目不足20条，则全部展示。

图3.17.2 GO富集柱状图

图3.17.3 GO富集气泡图

3.17.3 KEGG富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库，其中整合了基因组化学和系统功能信息等众多内容。信息分析时，我们将KEGG数据库按照动物，植物，真菌等进行分类，依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的pathway，KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。

KEGG的富集结果表格形式与GO富集结果相似，此处不再赘述。

图3.17.4 KEGG富集柱状图

图3.17.5 KEGG富集气泡图

3.18 miRNA,lncRNA,mRNA关联分析

非编码RNA分子可以参与调控蛋白编码基因的表达。根据ceRNA (competing endogenous RNA) 假说，细胞内的各类含有miRNA结合位点的RNA通过竞争性地与miRNA结合构成RISC达到调控基因表达的目的。因此，全转录关联分析是对各类ncRNA与mRNA间的调控关系以及ncRNA之间的互作网络进行构建与分析。

图3.18.1 miRNA调控基因表达原理示意图

3.18.1 LncRNA与mRNA关联分析

LncRNA通过共定位或共表达来调控靶基因表达。我们对差异表达 的lncRNA靶基因与差异表达的基因进行上下调分组分析，来定义其调控关系。当显著差异表达的lncRNA的靶基因同时也是显著差异表达的mRNA时，则该mRNA可能受到lncRNA直接或间接调控。

图3.18.2 lncRNA靶基因与差异mRNA共表达venn图

3.18.2 miRNA与lncRNA关联分析

部分miRNA来源于lncRNA，在分析miRNA的靶向 lncRNA时，我们需要将这些潜在的前体lncRNA过滤掉。使用Blast软件通过同源分析来寻找这一类lncRNA。

  1. miRNA:miRNA 前体编号ID
  2. hairpin_start:hairpin起始位置
  3. hairpin_end:hairpin终止位置
  4. lnc_id:lncRNA编号
  5. lnc_start:lncRNA起始位置
  6. lnc_end:lncRNA终止位置

lncRNA上带有miRNA结合位点时，miRNA可以与lncRNA结合导致lncRNA的降解。我们通过miRanda软件预测差异表达的miRNA的靶向lncRNA。

  1. miRNA:miRNA 编号ID
  2. lncRNA:miRNA靶基因

3.18.3 miRNA与mRNA关联分析

mRNA上带有miRNA结合位点时，miRNA可以与mRNA结合导致mRNA的降解。我们通过miRanda软件预测差异表达的miRNA的靶向mRNA。

  1. miRNA:miRNA 编号ID
  2. mRNA:mRNA靶基因

3.18.4 GO富集分析

  1. ID:GO条目编号
  2. Description:GO条目名称
  3. GeneRatio:对应条目含有的显著差异基因数/显著差异总基因数
  4. BgRatio:对应条目基因数/总基因数
  5. stat:显著性统计量
  6. pvalue:显著性p值
  7. p.adjust:矫正后显著性p值
  8. qvalue:显著性q值
  9. geneID:条目含有的显著差异基因ID
  10. Count:条目含有的显著差异基因数目
  11. geneID_Symbol:条目含有的显著差异基因Symbol名称
  12. Term:GO条目编号
  13. Type:GO条目所属类型

根据GO富集分析结果，在GO的三个分类中（cellular_component:细胞组分；biological_process：生物学过程；molecular_function：分子功能）分别统计每个分类上调基因、下调基因及所有基因的数目，以柱状图和气泡图展示。

柱形图和气泡图是富集分析结果的图形化展示方式。在此图中，富集程度通过GeneRatio、p.adjust和富集到此通路上的基因个数等指标来衡量并可视化。其中GeneRatio指差异表达的基因中位于该pathway条目的基因数目与所有有注释基因中位于该pathway条目的基因总数的比值。p.adjust是做过多重假设检验校正之后的Pvalue，p.adjust的取值范围为[0，1]，越接近于零，表示富集越显著。我们挑选了给类型富集最显著的10条目在该图中进行展示，若富集的条目不足20条，则全部展示。

图3.18.3 GO富集柱状图

图3.18.4 GO富集气泡图

3.17.3 KEGG富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库，其中整合了基因组化学和系统功能信息等众多内容。信息分析时，我们将KEGG数据库按照动物，植物，真菌等进行分类，依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出差异基因相对于所有有注释的基因显著富集的pathway，KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。

KEGG的富集结果表格形式与GO富集结果相似，此处不再赘述。

图3.18.5 KEGG富集柱状图

图3.18.6 KEGG富集气泡图

4. 参考文献

[1]Nath， J. and A. Nath， A Modified Algorithm for Quantifying of Pre-mature MiRNAs using Some Fractal Parameters[J]. International Journal of Advanced Computer Research 2012. 2:202-207.

[2]. Bartel DP.MicroRNAs: genomics， biogenesis， mechanism， and function[J].Cell.2004，116(2):281-97.

[3] LONG D， LEE R， WILLIAMS P， et al. Potent effect of target structure on microRNA function [J]. Nature Structural & Molecular Biology， 2007， 14(4): 287-94.

[4].郭韬，李广林，魏强，梁永宏.植物MicroRNA功能的研究进展[J].西北植物学报.2011，31(11):2347-54.

[5]Cock， P.J.A.， Fields， C.J.， Goto， N.， Heuer， M.L.， and Rice， P.M. (2010). The Sanger FASTQ file format for sequences with quality scores， and the Solexa/Illumina FASTQ variants. Nucleic acids research 38， 1767-1771.

[6]Erlich， Y.， and Mitra， P.P. (2008). Alta-Cyclic: a self-optimizing base caller for next-generation sequencing. Nature methods 5， 679-682

[7]Jiang， L.， Schlesinger， F.， Davis， C.A.， Zhang， Y.， Li， R.， Salit， M.， Gingeras， T.R.， and Oliver， B.(2011). Synthetic spike-in standards for RNA-seq experiments. Genome research 21， 1543-1551.

[8]Meyers， B.C.， et al.， Criteria for annotation of plant MicroRNAs[J]. Plant Cell， 2008. 20(12): 3186-90.

[9]Axtell， M.J. and B.C. Meyers， Revisiting Criteria for Plant MicroRNA Annotation in the Eraof Big Data[J]. Plant Cell， 2018. 30(2): 272-284.

[10]Michael I Love，Wolfgang Huber，Simon Anders.(2014).Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.Genome Biology，DOI:10.1186/s13059-014-0550-8.

[11] German， M.A.， et al.， Construction of Parallel Analysis of RNA Ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome[J]. Nat Protoc， 2009. 4(3):356-62.

[12]Young， M.D.， Wakefield， M.J.， Smyth， G.K.， and Oshlack， A. goseq: Gene Ontology testing for RNA-seq datasets.

[13]Kanehisa M， Araki M， Goto S， Hattori M， Hirakawa M， et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids research36:D480–484.


5. 分析软件列表

分析步骤 软件名称 版本
测序数据质控 fastqc v0.12.1
测序数据质控 multiqc 1.17
测序数据质控 fastp 0.23.4
测序数据质控 cutadapt 4.3
数据比对 hisat2 2.2.1
比对数据统计 qualimap 2.2.1
序列组装 stringtie 2.2.1
转录本鉴定 gffcompare v0.12.6
lncRNA预测 CPC2.py 0.1
lncRNA预测 CNCI CNCI version 2 Feb 28, 2014
lncRNA预测 PFAM 1.0
基因表达定量 stringtie 2.2.1
统计与绘图 R 4.3.2
可变剪接分析 rmats v4.2.0
可变剪接分析结果绘图 rmats2sashimiplot v1.0
突变分析 GATK v4.1.1.0
突变注释 snpEff 4.3t
差异分析 DESeq2 1.42.0
富集分析 clusterProfiler 4.6.2
RNA类型注释 blast 2.14.0+
miRNA数据比对，定量与结构分析 mirdeep2 0.1.3
miRNA靶基因预测 miRNada v3.3a
miRNA靶基因预测 RNAhybrid -

6. 联系我们

地址：广东省深圳市南山区留仙大道1201号大学城创客小镇

邮编：518055