测序寻找新的lncRNA并分析，完整的实验就应该这么做！

Posted 2023-04-08

参考技术A Characterization of complete lncRNAs transcriptome reveals the functional and clinical impact of lncRNAs in multiple myeloma

完整lncRNAs转录组的表征揭示了lncRNAs在多发性骨髓瘤中的功能和临床影响

发表期刊：Leukemia

发表日期：2021 Feb 17

影响因子：8.665

DOI:  10.1038/s41375-021-01147-y

一、研究背景

多发性骨髓瘤(MM)是一种以骨髓中浆细胞(PC)不受控制的克隆性增殖为特征的血液学肿瘤。尽管这种疾病的治疗取得了进展，目前中位生存期为7年，但它仍然被认为是一种无法治愈的恶性肿瘤，因为大多数MM患者对治疗产生抗药性导致疾病进展。

虽然约90%的基因组被转录成RNA，但只有1-2%被翻译成蛋白质，这凸显了人类细胞中非编码转录组的规模。lncRNAs表达的失调可以影响不同类型癌症（包括MM）发病和/或进展的相关途径。在MM中，少量lncRNAs的表达改变与患者的进展和生存有关，提示这些元素在疾病的发病机制中起着一定的作用。

二、材料与方法  

  1   数据来源

1）从38名新诊断的未经治疗的MM患者和3名健康供体中获得骨髓抽吸标本（GSE151063）

2）使用IA14发布的多发性骨髓瘤研究基金会（MMRF）CoMMpass研究数据集中的生存数据（n=542）

3）SMILO敲除相关的MARS-seq 数据：GSE134057

  2   分析流程

1）ssRNA-seq（测序）

2）lncRNAs注释：对基因组中的lncRNAs的位置进行注释

3）差异表达分析：limma软件包

4）样本异质性研究和基于FC的基因表达变化选择：采用变异系数(CV)对样本的变异性进行研究；将lncRNAs分类为上调（至少50%的样本中logFC大于1，小于25%的样本中logFC小于-1），下调（至少50%的样本中logFC小于-1，小于25%的样本中logFC大于1）和无变化（其余lncRNAs）；使用R软件包实现的SOM神经网络对lncRNA进行聚簇

5）染色质组蛋白标记分析：定义了89个lncRNAs与MM中染色质标记的de novo增益

6）lncRNA SMILO的研究和表征：跨越SMILO启动子的CpGs的DNA甲基化数据来源于本组之前发表的数据；细胞培养；SMILO敲除实验；RT-qPCR；增殖和凋亡测定；干扰素α治疗实验

7）MARS-seq：使用DEseq2包进行归一化和差异基因表达分析；基因本体（GO）和基因集富集分析（GSEA）

8）生存研究：单cox、多cox

三、结果展示

01 - MM的整个lncRNAs转录组的特征分析

对从38名MM患者骨髓中纯化的PC进行配对端ssRNA-seq。通过长度、低编码潜能和表达水平来筛选这样的转录本，鉴定出40,511个新型lncRNA，它们在38个MM患者样本中至少有3个表达（图1A）。在新的MM患者样本中验证了其中一些新型lncRNA的表达（补充图1A）。

MM中表达的编码基因和lncRNA基因数量的比较，后者包括：（1）以前在Gencode G19中注释的lncRNAs（G19lncRNAs），（2）以前工作中在不同B细胞亚群中鉴定的lncRNAs（BC鉴定的lncRNAs），（3）在MM患者样本中发现的一组新的lncRNAs（MM-identified lncRNAs）。在MM中发现的新型lncRNAs构成了所研究的lncRNAs群体中最大的一组，占MM PC中所有表达基因的56%（图1B）。为了确定MM细胞的特定基因组区域是否与lncRNAs的转录增加有关，分析了这些元件的全基因组分布，观察到编码基因和长非编码基因在染色体中均匀分布（图1C）。

接下来，lncRNAs根据其与编码基因的距离进行分类，显示上游转录物是最常见的类型，其次是下游lncRNAs，以及位于编码基因内部的lncRNAs（图1D）。与之前注释的lncRNAs相比，MM中发现的lncRNAs位于编码基因内部的比例更高（图1D）。此外，编码基因内藏有这种MM识别的lncRNAs的表达明显高于其余没有MM识别的lncRNAs的编码基因(图1E)，这表明特定编码基因的表达增加可能引发MM细胞中lncRNAs子集的调控，或者反之亦然。

这些结果表明，编码基因和lncRNA基因，可能一起并从基因组的相同区域编码，可能是肿瘤发展的关键参与者。

02 - lncRNAs在MM中表达的异质性和特异性

接下来，作者比较了MM和从健康捐献者骨髓（BMPC）中分离的正常PC之间的lncRNAs转录组。尽管MM标本中发现了大量的lncRNAs，但只有571个lncRNAs和78个编码基因有差异。分析了MM PC和BMPCs中lncRNAs和编码基因的CV，检测到MM比BMPCs中所有类型转录物的表达异质性程度更高（补充图1B）。MM样本中lncRNAs的表达异质性明显高于编码基因（图2A；补充图1B），这一发现可能解释了检测到的差异表达lncRNAs数量较少的原因，这表明这些元素可能有助于理解疾病的临床异质性。

为了检测异常表达的lncRNAs以解释这种异质性的方式，单独比较了每个MM患者和BMPCs的表达谱。利用基于FC的基因表达变化标准，在MM患者中发现了10351个过度表达和9535个下调的lncRNAs（图2B）。其中，检测到的lncRNAMALAT1，在以前的MM研究中描述过，还在一系列新的MM患者中验证了一些差异表达的lncRNAs（补充图1C）。

接下来，作者的目的是从先前的分析中确定的B细胞分化背景下在MM-PCs中失调的lncRNAs子集，因为它们可能代表疾病的特定治疗靶点。为此，分析了这19886个lncRNAs在B细胞分化状态的不同正常亚群中的表达，并与MM PC中的表达进行比较。观察到lncRNAs的三种不同的表达模式（图2C）。簇1包含2760个lncRNAs，在B细胞分化过程中有不规则的表达模式，在MM PC中有一致的高表达。簇2包含675个lncRNAs，在整个B细胞分化过程中低表达，在MM PC中略有增加。最后，簇3显示，在整个B细胞分化过程中，989个lncRNAs的表达非常低且均匀，MM样本明显增加。最后一种表达模式表明存在一组几乎完全在MM-PCs中表达的lncRNAs（称为MM特异性lncRNAs）。

03 - MM特异性lncRNAs的调节

为了确定MM中特异性lncRNA的表达是否受表观遗传调控，分析了之前工作中6个定义常见染色质状态的组蛋白标记的ChIP-seq数据（H3K4me3、H3K4me1、H3K27ac、H3K36me3、H3K27me3和H3K9me3）。

与正常B细胞亚群相比，观察到MM中MM特异性lncRNAs位点的活性组蛋白标记在全局范围内增加（图2D；补充图2A），且主要与活性启动子和增强子有关（图2E；补充图2A）。虽然大多数MM特异性lncRNA表现出活性染色质标记的增加，但这些lncRNA中只有一小部分(989个中的89个)呈现出新生的染色质激活，即存在于正常B细胞亚群中的抑制性标记被MM标本中的激活性染色质修饰所取代(图2F；补充图2B)。这89个lncRNAs的表达显示新生的表观遗传激活(图2F)明显高于其他MM特异性lncRNAs(图2G；补充图2B，C)。

04 - MM特异性lncRNA SMILO对MM细胞的生存至关重要

在MM中从新生表观组激活区域表达的89个lncRNA中，作者发现了LINC00582（ENSG00000229228，命名为SMILO）（图3A），以及由两个外显子组成的基因间lncRNA，位于TSNAX和DISC1编码基因之间，转录自染色体带1q42.2的负链，这是MM患者中经常扩增的基因组区域。

SMILO的表达在整个B细胞分化过程中无法检测到，除了在一些BMPCs中的边缘表达水平外（图3B），与BMPCs相比，64%的MM患者的SMILO表达上调。SMILO的表达在1q扩增的患者中显著升高，尽管这种表达的增加并不是这组患者所独有的（补充图3）。与正常PC形成对比，SMILO位点的新的表观基因组激活与MM PC的DNA甲基化丢失有关（图3C；补充图4A）。这些结果表明，除了1q扩增外，表观遗传机制也参与了MM患者中SMILO的激活及其过度表达。

敲除SMILO导致三种MM细胞系增殖率下降，凋亡细胞百分比增加（图3D；补充图4B），表明SMILO过度表达对MM细胞的生存至关重要。SMILO敲除后，KMS-11细胞中的RNA-seq分析显示，分别有84个和110个基因的下调和上调（图3E）。SMILO敲除后下调的编码基因富集在几个调节基因表达的过程中以及MM细胞的相关已知功能和途径（图3F）。抑制SMILO表达后，上调基因富集的首要途径之一是I型干扰素(IFN)信号通路(图3G;补充图4C)，其失调已被证明是MM细胞稳态的关键。此外，敲除SMILO导致几个干扰素刺激基因上调，说明MM中SMILO上调维持了这些编码基因的抑制，从而对MM细胞产生抗凋亡和增殖作用。这些结果在另外两个骨髓瘤细胞系中通过qPCR进一步验证（图3H；补充图4D）。

通过在MM.1S、MM.1R和KMS-11mm细胞系中添加不同浓度的IFNα，证明了IFN途径参与MM细胞的死亡。IFNα的使用引发了细胞凋亡的增加、细胞增殖的减少和不同ISG的上调（补充图4E-G）。此外，内源性逆转录病毒（ERVs）的表达，在抑制SMILO后上调（图3I；补充图4H），表明这些元件可能负责IFN途径的激活。

总之，数据表明，SMILO过度表达是MM细胞存活所必需的，其抑制可能触发ERVs的过度表达和IFN途径的激活，最终导致诱导细胞自主死亡，可能通过免疫原性细胞死亡（图3J）。

05 - MM特异性lncRNAs的预后价值

本研究最终旨在确定MM特异性lncRNAs的表达是否对MM患者具有预后价值，为此使用了来自IA14 CoMMpass研究中542名患者的RNA序列数据。由于CoMMpass中包含的RNA-seq数据只能提供关于先前注释的lncRNAs，将分析限制在Gencode中注释的89个MM特异性lncRNAs中的7个。在CoMMpass研究中包括的样本中检测到这7个lncRNAs中的6个的表达水平：ANKRD20A5P、SMILO、PDLIM1P4、ENSG0000249988、ENSG0000254343和RHOT1P1（补充图5A）。ANKRD20A5P、SMILO、ensg0000254343和RHOT1P1的表达与amp（1q）的存在显著相关，而PDLIM1P4和ensg0000249988的表达与不同的MM基因群没有显著相关性（补充图3）。

为了评估MM特异性lncRNA的表达是否与MM患者的预后相关，根据每个lncRNA的表达水平分析了这些患者的PFS和OS，根据表达水平将病例分为两组（补充图5B）。进行了单变量统计生存分析，将MM患者分为两个危险因素组，观察到PDLIM1P4、ENSG0000249988和ENSG0000254343的表达与PFS相关（图4A-C；补充图6A-C）。在OS分析中，PDLIM1P4、SMILO和ENSG0000249988的表达显示出具有统计意义的结果（图4D-F；补充图6D-F）。

在单变量分析后，对单变量分析结果显著的lncRNAs进行了多变量统计分析，并对PFS和OS的不同临床和遗传改变进行了统计分析。检测到PDLIM1P4的高表达以及ISS的2期和3期、del（13q）、t（8，14）、TP53、性别男性，以及硼替佐米IMIDs和Carfilzomib IMIDs治疗导致PFS的统计学显著性（图4G）。在PFS分析中，硼替佐米联合IMIDs和卡非佐米联合IMIDs对MM患者有良好的预后。在OS分析中，显示PDLIM1P4和ENSG0000249988的高表达以及ISS的2期和3期、硼替佐米IMIDs治疗、年龄超过65岁、amp（1q）、del（13q）、del（17p）和性别男性可将MM患者分为不同的风险组（图4H）。ENSG0000249988的过度表达和硼替佐米与IMIDs联合使用与较长的OS相关。

最后，还进行了ANOVA检验，比较单独来自临床和遗传高危因素的模型，或与lncRNAs表达相结合的模型，发现PFS和OS的第二种情况都有显著改善。

四、结论

综上所述，本研究为MM的lncRNAs转录组提供了一个全面的图景，表明这些非编码元件在MM细胞中异质、动态、特异性表达，并且在某些情况下，在MM细胞中是重新激活的。此外，发现lncRNAs可能在MM的发病机制中起着重要的作用，它们可以作为预后生物标志物，甚至作为最终改善MM患者预后的治疗靶点。

lncRNA数据处理（get fasta>>>novel lncRNA))）

参考技术A 情况说明：我有两个实验，每个实验两个重复，我现在要鉴定到新的lncRNA之后进行差异分析。大体思路，我会先通过常规的鉴定novel lncRNA方法鉴定到新的转录本，将转录本与参考基因组的转录本合并，之后进行定量分析。因此，我之前处理好各个文件的接头序列之后比对基因组后，会将四个bam文件合并成一个，之后进行鉴定（通过length、fpkm、cov、编码能力预测）之后会将鉴定到的转录本合并（与基因组的转录本），使用RSEM定量

使用到的软件有--fastQC、fastp、STAR、stringtie、cuffcompare、cpc（cpc2）、CNCI、CPAT、

1.常规的数据处理步骤

    1）拿到数据后先进行质检（常用fastqc看测序质量、用fastp做处理去接头什么的）

        fastqc -t 4 ` ls * `

         fastp -i $file-r1.fastq.gz -o $file-f.r1.fastq.gz -I $file-r2.fastq.gz -O $file-f.r2.fastq.gz -Q --thread=5 -5 --detect_adapter_for_pe 2> $file.log

            ###$file是我用for循环做的变量（这里应该是你的测序文件，PE150双端150bp测序，所以会有两个文件，分别是r1和r2。fastp软件使用可以直接看manual ，比较杂，我主要用了其中的去接头处理，关闭质量过滤”Q” 2>log文件是使用该软件处理过程中生成的日志。

    2）开始进行比对（根据物种来选择一个合适的比对软件，植物多倍体推荐bwa，STAR，我用的是STAR，大体参数有这些

        STAR --runThreadN 20 \

#        --outSAMmapqUnique 60 \

#        --twopassMode Basic \

#        --outSAMtype BAM SortedByCoordinate \

#        --limitBAMsortRAM 200000000000 \

#        --outFilterMultimapNmax 1 \

#        --outSAMprimaryFlag AllBestScore \

#        --outReadsUnmapped Fastx \

#        --outFilterMatchNmin 50 \

#        --alignMatesGapMax 10000 \

#        --outSAMstrandField intronMotif \

#        --outFilterMismatchNoverLmax 0.03 \

#        --outFilterType BySJout \

#        --readFilesCommand zcat \

#        --alignSoftClipAtReferenceEnds No \

#        --alignEndsType Local \

#        --alignIntronMin 20 \

#        --alignIntronMax 15000 \

#        --outSAMattributes All \

#        --outFilterMismatchNoverReadLmax 0.02 \

这个参数比较多，但实际应用的时候可以相应减少，其实很多时候默认参数就可以，发文章的时候默认参数比较有说服力。

--runThreadN 20 --genomeDir $index \

        --readFilesIn $name_1P.fq.gz $name_2P.fq.gz \

        --outReadsUnmapped Fastx \

        --readFilesCommand zcat \

        --outSAMtype BAM SortedByCoordinate \

        --outFileNamePrefix ./$name_star/$name

###这是我用的参数，仅供参考，我是为了之后找lncRNA

    3)比对完之后要进行组装，有cufflinks stringtie等，此处用的是stringtie

"""stringtie -p 4 -G ref-genome.gtf -o output.gtf input_file """

    ###有的是在这里会将多个组装好的gtf文件merge到一起，可以使用stringtie 的merge

    4）使用cuffconpare进行分类

""" cuffcompare -r ref-genome.gtf  -o output input_file """

       ###会从中得到好多文件，一般我们会从*.tmap文件中进行筛选

"""awk 'if($7>=0.5 && $10 > 1 && $11 >200) print $0' 自己命名的.tmap |awk 'if ($3=="u" || $3=="x" || $3=="i" || $3=="j" || $3=="o")print $0' > newfile.gtf 

       ###其中$7是fpkm值、$10是cov值、$11是length ，使用管道直接筛选出$3是u、x、i、j、o的code，具体的分类值可以从cuffcompare的manual看到

    5）到了写脚本的时候啦，要做的有，从上一步筛选到的转录本的到它的转录本序列（我是得到它的每个转录本有多少个exon，位置是多少，从已知基因组中得到序列后将其合并到一起，形成转录本序列，之后进行coding potential calculator（CPC）CNCI等软件来做，CPC我是安装到liunx 本地的，