RNA-seq鉴定LncRNA概述

发表于
RNA-seq鉴定LncRNA概述：lncRNA简介，分析pipline，分析软件，结果评估指标

lncRNA简介

  人类细胞中仅有少数转录生成的RNAs可作为蛋白质合成的模板。其余的RNA被称为非编码RNAs (ncRNAs)，它们位于编码蛋白的基因之间，其中长度超过200nt的ncRNAs称为lncRNAs(Long non-coding RNAs)

  LincRNA序列的获取有以下两种途径：①可参考LncRNA (Long non-coding RNA)序列的获取，只是在获得LncRNA的序列后需确认该LncRNA为LincRNA；②通过检索LincRNA的研究论文获取其序列，很多LincRNA的序列会在附件中包含或者文章提供链接地址。

即lncRNA是一大类，长度超过200nt的称为LncRNA，位于基因间区的称为LincRNA.

lncRNA特点

  • lincRNA典型的表现出显著的组织和细胞特异性表达
  • lincRNA可能与其共表达的蛋白编码基因具有共同的生物学功能
  • lincRNA能折叠成热力学稳定的二级或者更高级的结构，这是其发挥功能的基础
  • 相当大比例的lncRNA是编码基因的反义转录本或者与编码基因的外显子区正义方向部分重叠。

RNA-seq分析pipline

  1. RNA-Seq guide (RNA-seqlopedia)

转录组内的RNA，按照polyA形态：

  • 带polyA的RNA(mRNA和大部分的lncRNA)
  • 不带polyA的RNA(小RNA和小部分的lncRNA)

取样实例：晚期肝癌病人的肝组织（共四个）

  • 癌旁组织（N）
  • 原发灶（P）
  • 转移灶（M）
  • 门脉血栓转移灶（V）

  1. twbattaglia/RNAseq-workflow

  1. griffithlab/rnaseq_tutorial,  文章 ,   Lectures

Expression

重复文章 “RNA-seq的转录水平表达分析”

该流程不适用小样本数据：For small sample sizes (n < 4 per group), it is often better to perform regularization. This can be done using the limma package in Bioconductor

  1. RNACocktail：A comprehensive framework for accurate and efficient RNA-Seq analysis， 文章

  2. RNA-seq workshops and online tutorials

LncRNA注释GTF文件及fasta文件

  1. 参考基因组fa及GTF文件下载：链接
  • Homo_sapiens.GRCh38.dna.primary_assembly.fa
  • Homo_sapiens.GRCh38.92.gtf

参考基因组建立索引：

which hisat2_extract_exons.py $known_coding_gtf > genome.exon
which hisat2_extract_splice_sites.py $known_coding_gtf > genome.ss
time hisat2-build -p 8 $ref_hg38 –ss genome.ss –exon genome.exon genome_tran

  1. GENCODE    2. NONCODE

ps: NONCODE has updated to NONCODEv5. NONCODE2016 website has been moved to http://www.bioinfo.org/NONCODE2016

lncRNA建立索引：

which hisat2_extract_exons.py $lncRNA_gtf > NONCODE2016_human.exon
which hisat2_extract_splice_sites.py $lncRNA_gtf > NONCODE2016_human.ss
time hisat2-build -p 8 $lncRNA_ref –ss NONCODE2016_human.ss –exon NONCODE2016_human.exon NONCODE2016_index

  ENCODE计划揭示了大约76%的人类基因组转录产生非编码蛋白的RNA分子，其中包括大约10000lincRNAs。NOCODE最新版数据库收录了>30000个人的lincRNAs和 >20000 小鼠的lincRNA。

RNA-seq分析软件

推荐文献：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis

基因组和转录组比对

基因水平定量

官网：HTSeq: Analysing high-throughput sequencing data with Python

使用：使用HTSeq进行有参转录组的表达量计算（推荐）

扩展阅读：

  1. featureCounts or htseq-count?

  2. htseq-counts跟bedtools的区别

stringTie

NCBI-Hackathons/RNA_mapping

  RNA-Seq基因组比对工具HISAT2：HISAT2是TopHat2/Bowti2的继任者，使用改进的BWT算法，实现了更快的速度和更少的资源占用，作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。官网

gene annotation 一般选择 RefSeq 或者 Ensembl，这里，我们的参考基因组fasta文件和基因注释文件都选择 Ensembl的hg38版本（release-91），其ftp地址为：http://asia.ensembl.org/info/data/ftp/index.html

wget ftp://ftp.ensembl.org/pub/release-91/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-91/gtf/homo_sapiens/Homo_sapiens.GRCh38.91.chr.gtf.gz

tophat(基于python2.*)

tophat 一些小错误

#genome index
samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa
bowtie2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa hg38_genome

#transcriptome-index
tophat -G $ref_mRNA_gtf –transcriptome-index=$tophat_transcriptome_index $bowtieGenomeIndex

PS: Please note that it is highly recommended that a FASTA file with the sequence(s) the genome being indexed be present in the same directory with the Bowtie index files and having the name .fa. If not present, TopHat will automatically rebuild this FASTA file from the Bowtie index files

比对质量评估

RSeQC: An RNA-seq Quality Control Package

pip2.7 install RSeQC

转录本重构

转录本构建效果评估—评价指标：

  • 多外显子比率
  • 转录本长度
  • 转录本的可变剪切数目
  • 对已知基因的覆盖程度
基因数目 转录本数目 多外显子比率 多外显子转录本数目
5-8万 >10万 30%-50% 5万

对已知编码基因的覆盖程度：>60%

使用StringleTie的-个G模式构建转录本的输出文件主要有：

  1. .gtf文件:记录组装的转录本信息
  2. gene_abundances.tsv文件：以tab键分割的记录基因丰度信息
  3. 在使用-B 参数下，生成*.ctab文件：用于下游Ballgown软件做差异表达分析的输入文件

输出结果具体解释可查看StringTie官网的output部分

合并转录本

cole-trapnell-lab/cufflinks

可以直接下载解压预编译好的版本，省去安装 Boost C++ libraries的麻烦。在一下链接中选择linux版即可：cufflinks

事实上，对于多个样本构建的多套转录本，如何得到统一的一套转录本，有下面三种方法：

  • 1 在转录本拼接之前，把各样本的比对bam文件合并，然后用合并的bam跑cufflinks
  • 2 每个样本的比对bam文件分别单独跑cufflinks，各样本的转录本构建后，再用cuffcompare合并为一套转录本。
  • 3 每个样本的比对bam文件分别单独跑cufflinks，各样本的转录本构建后，再用cuffmerge合并为一套转录本。

这三种做法区别在于：

  • 第一种方法流程相对简单，所有的工作都抛给cufflinks一人完成，你都不需要知道cuffmerge、cuffcompare的用法。貌似是种完美解决方案。但很大的问题是：cufflinks能处理得了最终合并的bam吗？对于小物种的样本还可以，但对于人，若是7，8个样本合成的bam，cufflinks吃不销！
  • 第二和三种方法是类似的，都是在保留可变剪切结构的前提下，将转录本合并。不同的是，cuffcompare只有A、B两条转录本结构相同的时候，才将A、B合并。而cuffmerge是A、B某些部分互相overlap，就将它俩合并。事实上，cuffmerge再做合并的时候，是把overlap的transfrag重新调用了cufflinks，合成一个transfrag。
  • 第二和三种方法还有一个很大的不同是，cuffmerge可以带上参考注释有参考的进行合并，而cuffcompare不能如此。这是cuffmerge对cuffcompare的一个优势。

因此，我认为第三种方法比第二种方法是最接近于第一种方法，而第一种方法的可实现性较差，在现有条件，最完美的解决方案就是第三种方法：每个bam单独跑cufflinks，跑完的结果再用cuffmerge合并。

参考：转录组测序分析中cufflinks的使用及问题

鉴定novel lncRNA

  • step1: 对所有样品拼接得到的转录本使用cuffcompare软件进行合并，筛除链方向不明的转录本；
  • step2: 选择转录本长度>=200bp，Exon个数>=2的转录本；
  • step3: 通过cufflinks计算每条转录本的reads覆盖度，选择至少在一个样品中覆盖度>=3的转录本；
  • step4: 若该物种存在已知lncRNA数据，首先通过cuffcompare软件，将前一步得到的转录本与已知lncRNA进行比较，得到与已知lncRNA相同的转录本。这一部分转录本直接纳入最终的lncRNA集，不再进行后续筛选。之后，通过与该物种已知非lncRNA及非mRNA类型（rRNA,tRNA,snRNA,snoRNA,pre-miRNA,pseudogenes等）的转录本进行比较，筛除那些与以上已知转录本相似或相同的转录本。若无已知lncRNA数据，则直接进行与该物种已知非lncRNA及非mRNA类型转录本的比较；
  • step5: 通过与已知mRNA进行比较，并利用cuffcompare分析结果中的class_code(http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/index.html#transfrag-class-codes)信息筛选候选lincRNA，intronic lncRNA, anti-sense lncRNA类型的转录本。

class_code

Priority Code Description
1 = Complete match of intron chain
2 c Contained
3 j Potentially novel isoform (fragment): at least one splice junction is shared with a reference transcript
4 e Single exon transfrag overlapping a reference exon and at least 10 bp of a reference intron, indicating a possible pre-mRNA fragment.
5 i A transfrag falling entirely within a reference intron
6 o Generic exonic overlap with a reference transcript
7 p Possible polymerase run-on fragment (within 2Kbases of a reference transcript)
8 r Repeat. Currently determined by looking at the soft-masked reference sequence and applied to transcripts where at least 50% of the bases are lower case
9 u Unknown, intergenic transcript
10 x Exonic overlap with reference on the opposite strand
11 s An intron of the transfrag overlaps a reference intron on the opposite strand (likely due to read mapping errors)
12 . (.tracking file only, indicates multiple classifications)

转录水平本定量

lincRNA鉴定及功能注释

CPC

差异表达分析

RNA-Seq workflow: gene-level exploratory analysis and differential expression

cuffdiff -L 选项有两种写法：
偶数个样本：

奇数个样本：

ballgown安装

需要先安装libxml2

wget http://xmlsoft.org/sources/old/libxml2-2.7.1.tar.gz
untargz libxml2-2.7.1.targz
cd libxml2-2.7.1
./configure --prefix=/home/wangdong/softwares/libxml2
make && make install


echo 'export PATH=/home/wangdong/softwares/libxml2/bin:$PATH'>>~/.bashr    
echo 'export LD_LIBRARY_PATH=/home/wangdong/softwares/libxml2/lib:$LD_LIBRARY_PATH'>>~/.bashrc 
echo 'export PKG_CONFIG_PATH=/home/wangdong/softwares/libxml2/lib/pkgconfig:$PKG_CONFIG_PATH'>>~/.bashrc

ps:安装 dplyr包之前需要安装assertthat包

可参考：链接

使用可参考：链接

ps:pheno_data里面第一列样本名需要和ballgown下面的文件夹的样本名一样，且先后顺序保持一致，不然会报错

富集分析

Gene ontology or pathway analysis for differentially expressed genes, this step is only applicable for protein coding genes

参考

（1） Identification and function annotation of long intervening noncoding RNAs

（2）转录组分析工具哪家强？

-------------本文结束感谢您的阅读-------------