初探mRNA、lncRNA联合分析之上游

本周就将尝试复现一篇mRNA、lncRNA联合分析的文章，内容比较简单，和常规分析流程查相差不多，比较适合我过度学习

RNA-Seq Comprehensive Analysis Reveals the Long Noncoding RNA Expression Profile and Coexpressed mRNA in Adult Degenerative Scoliosis(https://pubmed.ncbi.nlm.nih.gov/36035195/)

目的：随着世界范围内老龄化进程的加剧，成人退行性脊柱侧弯（ADS）的发病率正以惊人的速度增长。然而，与ADS病因相关的基因组研究在世界范围内很少报道。由于长链非编码RNA（lncRNA）在人类疾病进展中发挥着关键作用，本研究旨在通过RNA测序（RNA-seq）研究ADS相关信使RNA（mRNA）和lncRNA，并基于lncRNA-mRNA共表达网络和蛋白质-蛋白质相互作用（PPI）网络进行全面的生物信息学分析。方法：最初，从三名ADS和三名接受手术的非退行性腰椎创伤患者中获得六份全血（WB）样本，进行RNA-seq检测，以构建差异的mRNA和lncRNA表达谱。随后，进行定量RT-PCR（qRT-PCR），以验证从其他14名受试者（分别为7名ADS患者和7名非退行性腰椎创伤患者）的髓核（NP）组织中随机选择的3种差异表达的mRNA和lncRNA。结果：从RNA-Seq数据中共筛选出1651个上调和1524个下调的mRNA，147个上调和83个下调的lncRNA，构建了共表达网络，进一步研究了它们的调控相互作用。GO基因功能预测显示，lncRNA靶向基因可能通过参与多种生物学过程，如AMPK信号通路、溶酶体、泛素介导的蛋白水解以及细胞代谢过程，在ADS中发挥重要作用。此外，通过qRT-PCR分别验证了三种选定lncRNA和mRNA的表达水平，表明相对表达水平与RNA-seq数据一致。值得注意的是，失调的RNA AKT1、UBA52、PTPN12和CLEC16A在ADS WB样品中显著差异表达，可能成为未来研究的潜在调控基因。结论：本研究首次深入了解了与ADS相关的长链非编码RNA的转录组变化，为进一步探索这种鲜为人知的退行性疾病的临床生物标志物和分子调控机制铺平了道路。然而，ADS中这些候选lncRNA的详细生物学机制需要在未来的研究中进一步阐明。

基本分析流程：

• 上游mRNA、lncRNA定量
• mRNA、lncRNA差异表达分析
• 差异表达lncRNAs共表达mRNAs的富集分析
• 共表达网络、PPI

当然这篇文章是曾老师给我找的，曾老师认为文章中的分析可能存在问题，如热图比较奇怪：

同实验分组不同表达谱特征（红色区块）

我们接下来就从上游开始走一遍，学习这个简单联合分析流程，并尝试发现文章可能存在的问题

需要注意的是，因为篇幅较长，我把上下游拆开，本周将只涉及上游mRNA、lncRNA的定量，获取表达矩阵，下周再进行下游分享

我最近在学lncRNA包括上游的组装、预测和鉴定，当然这里没有鉴定新lncRNA，只对已知mRNA和lncRNA有个比对、定量

所以还是走这套流程：hisat2+stringtie

作者走的：hisat2+HTSeq

本文参考：

LncRNA从入门到精通 「生信技能树」LncRNA数据分析实战(https://www.bilibili.com/video/BV1Zg4y187ff/?spm_id_from=333.1007.top_right_bar_window_custom_collection.content.click) 明明PCA区分非常好，但是差异基因数量很少？ 使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？

下载fq文件；下载人的参考基因组fasta序列，并且构建hisat2的索引文件

发现和前面做单细胞从cellranger下载的参考基因组fa文件还是有区别

同时cellranger下载的还有gtf注释文件和STAR构建的索引

nohup hisat2-build -p 12 Homo_sapiens.GRCh38.dna.toplevel.fa human_hisat2 &

人类fa建索引还挺长时间的

注释文件：

https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gtf.gz

下载fq文件，获取srr_acc_list.txt

fastp

hisat2+samtools

stingtie

组装：

WARNING: no reference transcripts were found for the genomic sequences where reads were mapped! Please make sure the -G annotation file uses the same naming convention for the genome sequences.

(比较之前其他项目使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？stingtie组装结果)

可以发现少了reference_id、reference_gene_id、ref_gene

如果一条道走到黑:

合并

（相比）

定量

（相比）

解决问题·：WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!

Please make sure the -G annotation file uses the same naming convention for the genome sequences.

参考：

no reference transcripts were found for the genomic sequences where reads were mapped!(https://github.com/gpertea/stringtie/issues/287)

“I solved this issue by simply deleting every "chr" which shows chromosome number in the gtf file: mm10.ncbiRefSeq.gtf (saved as mm10.ncbiRefSeq_MODIFIED.gtf) . i.e., replaced "chr1 to chr22" with "1 to 22" stringtie -p 8 -G GTF_FILES/mm10.ncbiRefSeq_MODIFIED.gtf -o STRINGTIE/Trim_MT1_chr19.gtf -l MT1 SAMTOOLS/Trim_MT1_chr19.bam”

查看：

发现和ENSEMBL给的注释文件相比GENECODE每行开头染色体编号确实都带上了”chr“

修改：

测试：

成功

修改注释文件重新组装：

合并：

定量：

使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？ RNA-seq : Hisat2+Stringtie+DESeq2

后面我们会介绍使用gffcompare发现新转录本的流程，这里我们仍走这个流程，但只对已知转录本定量

可以发现如果直接基于bam文件定量，给到的 -G 参数就是一开始的参考基因组gtf文件，而不是需要多走一步组装后merge的gtf文件

提取表达矩阵：

nohup python prepDE.py3 -i ./5.quantity/sample_list &

明明PCA区分非常好，但是差异基因数量很少？

整理基因表达矩阵：

直接提取有GeneSymbol就行

• 提取出第一列通过”|“划分两部分，并且第二部分不是以”ENSG“开头的行
• 在这个基础上，去除空白行

整理转录本表达矩阵：

可见基本上就是ENSEMBL的人转录本ENST编号和”新转录本“MSTRG

看作者的表述

并没有都是ENSEMBL还是存在对应的gene symbol

我个人觉得应该是先拿这些转录本去分析，最后有意义的再map到对应的基因上，因为一个基因可能对应多个转录本

所以后面我们先不对转录本id进行转换，最后再map