【原】中国荷斯坦奶牛新的lncRNA全基因组鉴定

咱们《生信技能树》的B站有一个lncRNA数据分析实战，缺乏配套笔记，所以我们安排了100个lncRNA组装案例文献分享，以及这个流程会用到的100个软件的实战笔记教程！

下面是100个lncRNA组装案例文献分享

标题：Genome Wide Identification of Novel Long Non-coding RNAs and Their Potential Associations With Milk Proteins in Chinese Holstein Cows

标题：中国荷斯坦奶牛新的lncRNA全基因组鉴定及其与乳蛋白的关系

杂志：Frontiers in Genetics（2018）

通讯作者：Shengli Zhang

机构：中国农业大学动物科学技术学院，农业部动物遗传繁育重点实验室，动物繁育国家工程实验室

文章链接:https:///10.3389/fgene.2018.00281

摘要：

本研究采用全转录组RNA测序技术，对3头乳蛋白率极高和3头乳蛋白含量低的中国荷斯坦奶牛的乳腺组织样本进行了LncRNA转录组分析。

在这项研究中，通过5个严格的步骤和编码潜力的筛选，共鉴定出6450个lncRNA转录本。总共鉴定出31个lncRNAs和18个新基因在高乳蛋白样品(HP)和低乳蛋白样品(LP)中存在差异表达。**通过生物信息学分析选择差异表达的LncRNA预测目标基因，然后整合差异表达的mRNA数据、基因功能、基因本体(GO)和途径、全基因组关联研究(GWAS)和数量性状位点(QTL)信息，以及网络分析以进一步描述潜在的相互作用。**有几个LncRNA(如XLOC_059976)可作为预测乳蛋白含量的候选标记。

这是第一个对与奶牛乳蛋白特性相关的lncRNAs和mRNAs进行全球表达谱分析的研究。这些结果为牛奶蛋白质的合成提供了重要的信息和见解，也为未来牛奶品质的改善提供了潜在的目标。

关键词：长非编码RNA，乳腺，转录组，牛奶蛋白，综合研究

背景知识：乳蛋白是人体最重要的营养物质之一。牛奶中蛋白质的数量和组成在很大程度上由奶牛的遗传因素决定，虽然已通过QTL定位、候选基因分析、GWAS或NGS技术在奶牛中鉴定出一些影响产奶量和成分的致病基因和突变，但乳蛋白的合成和分泌涉及复杂的过程，需要进行彻底的检测。已证明产奶性状受到强大的表观遗传调控。表观遗传基因调控机制是通过对染色质结构的调节来实现的，既可以抑制基因表达，也可以增强基因表达。lncRNAs与家畜的发育、代谢和免疫调节以及复杂性状的适应和表型变异有关。

数据和方法

选择产后60天(泌乳高峰期) 一胎 6只的牛奶蛋白含量极其高和低的各三头(高≥3.5%和低≤3.0%)健康，无乳腺炎的荷斯坦奶。测序数据公开可以获取，编号是：PRJNA416150

1、候选转录本组装流程

• 测序数据：Illumina Hiseq 2500；125 bp paired-end reads.；去除rRNA建库，链特异性；PRJNA416150

• 质控：Fastqc

• 过滤：Fastx_toolkit (0.0.13)

• 比对：TopHat2   bovine genome(UMD 3.1),

• 组装：

• Cufflinks（ “min-frags-per-transfrag = 0”, “–library-type = fr-firststrand” and “–mask-file = ncRNA.gtf”）
• Scripture,（pairedEnd）
• Stringtie （ (version 1.0.1)   (-f 0.01 -c 0.01).）
• Transcomb   (V1.0)  “-s first”, “-l 200” and “-e 50”

2、lncRNA过滤筛选和鉴定

• a.转录本类别   i  u  x
• b.长度≥200nt，外显子大于2
• c.FPKM≥0.3。过低表达量被认为是转录噪音。
• d.开放阅读框小于120aa  ;预测软件：TransDecoder (3.0.1)
• e.转录本编码能力评估CPC <0 ,PLEK <0 ;CNCI 大于<0
• f.转录本转换成 氨基酸序列评估编码能力 ：HMMER  Pfam30.0

这个流程在文章有一个示意图，见文末！

3.保守分析

使用PhyloFIT计算物种间保守区和非保守区的系统发育模型。模型和HMM转换参数被输入到phastCons中，以计算一组lncRNAs和编码基因的保守分数

4.差异分析

P value <0.05

5.富集分析

DAVID

GO 显著富集：p-value <0.05

IPA检验lncRNA相关基因富集  。-log2(p-value)>1.3为显著富集

6.LncRNAs和miRNAs之间的相互作用

利用BLASTN找到已知miRNA前体，使用RNAfold （R包）对lncRNA转录本的二级结构进行预测。预测miRNA与lncRNA相关性使用MiRanda  score>160,energy <-15。

7.靶基因预测和功能分析

使用perl脚本检查lncRNA上下游10kb和100kb的编码基因。用Pearson法计算lncRNA与编码基因的表达相关性，p值<0.05。从AnimalQTLdb中提取了乳蛋白性状的QTL信息。3从12个GWAS研究中收集到972个与乳蛋白性状相关的显著SNPs

结果

1.高通量测序

• 比对率：

• 91.55% (range:91.10–92.02%) of the reads were mapped to the bovine genome (Ensembl UMD3.1)

• 83.22% (range: 82.01–84.41%) were uniquely mapped reads

• 8.33%(range: 7.17–9.44%) were multi-mapped reads

2.全基因组lncRNA鉴定及特征

总共6,450 lncRNA 转录本  5,256 lncRNA 位点，可以分成：

• 4,292 lincRNAs  （基因间区lncRNA）
• 810  ilncRNAs（内含子lncRNA）
• 1,348  lncNATs（反义lncRNA）

然后就是标准流程，看组装到的新lncRNA的一些特性，外显子数量和长度，表达量情况等等：

A. lncRNA表达水平低。

B. LncRNAs的大小明显小于蛋白质编码转录本、新基因和假基因。

C. lncRNA外显子主要集中在两个

D. lncRNA明显短于其他

F. lncRNAs的保守性比蛋白质编码区低

3.差异分析

在高蛋白和低蛋白之间，q-value<0.05，共有31个差异表达的lncRNA，其中15个上调，16个下调。

其中18个是 新lncRNA，6个上调，12个下调。

这个差异分析比较容易复现，基本上看我六年前的表达芯片的公共数据库挖掘系列推文即可；

• 解读GEO数据存放规律及下载，一文就够
• 解读SRA数据库规律一文就够
• 从GEO数据库下载得到表达矩阵 一文就够
• GSEA分析一文就够（单机版+R语言版）
• 根据分组信息做差异分析- 这个一文不够的
• 差异分析得到的结果注释一文就够

4.IncRNAs和miRNAs之间关系的预测

• 为了确定lncRNA是否真的是miRNA前体。作者比较了从miRBase 获得的lncRNA序列和miRNA序列，发现13个lncRNA含有8个完整的miRNA前体。对lncRNA转录本二级结构的预测表明，一些lncRNA含有稳定的miRNA前体发夹结构。

为了研究鉴定的lncRNA是否被miRNAs靶向，作者用Miranda软件分析了6450个lncRNA转录本。共有4972个LncRNA转录本被预测为788个牛miRNAs的靶标。

• 其中，206个lncRNA被miR-15a、miR-486、miR-135、miR-101a、miR-152和miR-139作为靶标，据报道它们与乳蛋白合成有关。一个差异表达的LncRNA(XLOC_059976)被预测为miR-139和miR-152的靶标，这意味着XLOC_059976可能作为乳蛋白合成过程的调节因子。

5.LncRNA与蛋白编码基因表达水平的相关性

5,251个lncRNA与18,227个mRNAs显著相关(p值<0.05)，其中前面的差异分析得到的统计学显著的31个差异表达的lncRNA与11,161个mRNAs显著相关(p值<0.05）。详见：TABLE S11 | Pearson correlations between protein-coding genes and differentially expressed lncRNAs.

候选LncRNAs、mRNAs和通路的网络图。黄色三角形、绿色圆圈和粉红色菱形分别代表lncRNA、mRNA和通路。

6.LncRNAs功能的预测

大量lncRNA的功能是未知的，但是它们主要是cis-regulators，所以可以根据它们临近的蛋白编码基因功能来近似推断，然后表达量的相关性也可以类推到。

• 根据位置关系推断 使用bedtools等工具！

• 表达量的相关性， 比如杂志Cancer Medicine, 2020的文章《 Genome-wide DNA methylation analysis by MethylRad and the transcriptome profiles reveal the potential cancer-related lncRNAs in colon cancer》，在进行结直肠癌相关lncRNA的功能富集分析，就是采用LncRN2Target v2.0和StarBase分析与15个lncRNA共表达的蛋白编码基因，其中lncRNA HULC和ZNF667-AS1分别鉴定到28个、9个共表达的蛋白编码基因！

预测了30个LncRNAs，调控34个影响乳蛋白合成的基因。

三个DEG(FOS、IRF2和SOCS2)被发现参与PRL信号通路(Li C.et al.，2016)。据报道，PRL对蛋白激酶C的刺激可能与PRL刺激FOS mRNA的积累有关，而蛋白激酶C的激活对PRL在乳制品合成和有丝分裂中的所有作用都是必不可少的。

XLOC_2427074、XLOC_2500996、XLOC_2938516和XLOC_593004的表达与fos基因显著相关，表明它们在乳蛋白合成中起重要作用。

IGFBP2与XLOC_1186672、XLOC_1243232、XLOC_1284424和XLOC_2273208相关，IgFBP5与XLOC_051681、XLOC_059976和XLOC_2280878相关。据此推测，lncRNA-igfbp2/igfbp5对可能参与胰岛素样生长因子途径，影响乳蛋白合成。

XLOC_051681与16个基因的表达相关。

总结：

本研究使用四个软件进行组装转录本。选取至少在两个样本，两个软件中都鉴定到的转录本，进行转录本的筛选。

对差异lncRNA和mRNA表达的综合解释表明，30个lncRNA可能调控34个影响乳蛋白合成的基因。

推测的LncRNA XLOC_059976可能是预测乳蛋白组成表型的关键候选生物标志物。

文末友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

• 数据挖掘（GEO,TCGA,单细胞）2021第4期
• 生信爆款入门-2021第4期