| RNA | 您所在的位置：网站首页 › featurecounts结果解读 › RNA |
RNA
|
RNA-seq数据分析
简介1 生物基础1.1 中心法则1.2 RNA-seq Protocol1.3 RNA-seq总的路线图
2 数据分析2.1 前期准备2.1.1 创建目录&安装conda2.1.2 常用文件格式简介
2.2 软件安装2.2.1 conda安装软件2.2.2 预编译版本软件安装2.2.3 源码安装
2.3 数据下载
下载完成之后一定要检查数据完整性，不然分析白做使用awk命令提取md5值，生成md5文件，运行以下命令检查文件完整性2.5 质控（QC）2.6 去接头2.7 hisat2比对2.8 定量和转录本组装
参考
简介
基因表达是功能基因组学研究的一个重要领域。基因表达与基因信息从基因组DNA到功能蛋白产物的流动有关。RNA-seq已成为一种标准的基因表达检测方法，尤其用于检测转录本相对丰度和多样性。已有研究表明，其可靠性可以与其他成熟的方法如定量聚合酶链式反应(qPCR)相媲美。 声明： 文中如有不足请留言讨论！ 1 生物基础吐槽： 系统性总结真的好困难，刚开始阅读文献也很吃力。但是行动是知识特有的果实，慢慢积累吧。 1.1 中心法则中心法则对于学习生物的人来说再熟悉不过了，基因信息的流动，蛋白的产生等一系列生物学事件都基于中心法则，是基因表达分析的一条主线。遗传信息从双链基因组DNA模板到翻译后修饰蛋白的流动是每个阶段至关重要的分子特征。RNA-seq通常针对成熟的mRNA分子。
典型的RNA-seq流程包括从感兴趣的样本中分离RNA、建库、高通量测序产生数以百万计的reads(长度一般为30-300bp)、比对到参考基因组或转录组，差异表达分析、发现转录本亚型和其它的下游分析。下图很好的概括了RNA-seq数据产生的过程。
这张图引自A survey of best practices for RNA-seq data analysis。包括了三个板块，分别是预分析，核心分析和高级分析。很好的概括了实验设计和拿到测序reads后的数据分析，并介绍了不同的分析路线和每一步的数据分析方法，可以说是一个很好的大纲。 注意： 这里省略了进入相关文件目录的步骤，大家分析时要注意。本次数据是小鼠早期胚胎测序得到的双末端数据，单末端相关参数可以使用–help查看帮助文档。但分析总体流程相同。 我使用的是HISAT2+featureCounts+StringTie流程。 2.1 前期准备 2.1.1 创建目录&安装conda创建目录 # 首先使用cd命令需要进入自己的目录 cd ~/ # 创建文件夹，用于存储参考基因组 mkdir -p /mnt/data/hwb/database/mouse_genome_esml_84/grcm38_tran/ # 建立软件安装目录 mkdir biosoft # 此外还需要建立项目分析的目录以及备份文件，方便查找conda安装 Conda 是一种通用包管理系统，旨在构建和管理任何语言的任何类型的软件。通常与 Anaconda (集成了更多软件包，https://www.anaconda.com/download/#download) 和 Miniconda(只包含基本功能软件包, https://conda.io/miniconda. html) 一起分发。 # 从官网下载最新版Miniconda3安装包，但速度较慢 wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh # 安装Miniconda3：安装过程只需要输入 yes 或者按 Enter bash Miniconda3-latest-Linux-x86_64.sh # miniconda3安装成功，并成功配置好环境变量 source ~/.bashrc # 镜像设置 # 下面这四行配置清华大学的bioconda的channel地址，国内用户推荐 conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/ # 配置镜像成功 conda config --set show_channel_urls yes conda config --set channel_priority flexible # 查看配置文件 cat ~/.condarc
conda常用命令 conda create -y -n trans python=3 # 创建小环境并成功安装python3 conda info --e # 查看当前conda环境 conda-env list # 列出所有小环境 conda activate trans # 激活小环境 conda deactivate # 退出小环境 2.1.2 常用文件格式简介1. SRA：NCBI SRA数据库存放格式 SRA（Sequence Read Archive）：SRA是一个数据库，NCBI为了解决高通量数据庞大的存储压力而设计的一种数据压缩方案。 一般使用fastq-dump或fasterq-dump来将其转换为fastq格式的数据，才能做后续分析。 2. FASTQ：高通量数据存储格式 FASTQ格式将序列和Phred质量存储在一个文件中。序列和质量得分皆由单个ASCII字符编码。最早由Sanger Institute开发使用，目前已经是高通量测序结果的标准格式。 Phred质量值简介： fasta序列格式如下： 必须包含的3列： (1). chrom - 染色体名字(e.g. chr3，chrY, chr2_random)或scafflold 的名字(e.g. scaffold0671 )。 (2). chromStart - 染色体或scaffold的起始位置，染色体第一个碱基的位置是0。 (3). chromEnd - 染色体或scaffold的结束位置，染色体的末端位置没有包含到显示信息里面。例如，首先得100个碱基的染色体定义为chromStart =0 .chromEnd=100, 碱基的数目是0-99。 9 个额外的可选列: (4). name - 指定BED行的名字，这个名字标签会展示在基因组浏览器中的bed行的左侧。 (5). score - 0到1000的分值，如果在注释数据的设定中将原始基线设置为１，那么这个分值会决定现示灰度水平（数字越大，灰度越高），下面的这个表格显示GenomeBrowser
这里介绍三种软件安装方法。 2.2.1 conda安装软件 # 可以一次安装多个软件 conda install -y -c hcc aspera-cli conda install -y sra-tools conda install -y fastqc trim-galore hisat2 multiqc samtools featureCounts # 运行以下命令，检查软件是否安装成功 hisat2 --help which ascp # 查看软件安装路径 2.2.2 预编译版本软件安装 # 以aspera为例 wget -c https://d3gcli72yxqn2z.cloudfront.net/connect_3_9_1_171801_ga/bin/ibm-aspera-connect-3.9.1.171801-linux-g2.12-64.tar.gz # 解压并运行脚本 tar -xzf ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh # 添加环境变量 export PATH=~/.aspera/connect/bin:\$PATH >>~/.bashrc # source使变量生效 source ~/.bashrc 2.2.3 源码安装 # 首先使用wget下载 # 安装分三步： # 配置 ./configure --prefix=想要指定的路径 # 编译 make # 安装 make install # 具体的安装在下载之后可以查看readme文件，一般有详细的介绍 2.3 数据下载我这里使用aspera下载sra数据。 首先我们通过阅读文章获取数据的GEO accession，在GEO数据库中找到对应的BioProject，这里建议在ebi数据库下载对应的SRR号的相关信号，并获取aspera下载链接以及md5值等相关信息。 # ------------------------------------------------------ # GEO accession：GSE70605 # ------------------------------------------------------ # 下载得到tsv文件首先需要使用awk命令提取文件内aspera链接，保存为srr.aspra文件。之后运行下面脚本 创建一个脚本,命名为aspera.sh： ```bash # 找到openssh文件位置，替换以下代码文件路径 find /home/hwb/ -name asperaweb_id_dsa.openssh #!/bin/bash cat srr.aspera | while read id do ascp -QT -l 300m -P33001 -i ~/miniconda3/envs/aspera/etc/asperaweb_id_dsa.openssh era-fasp@${id} ./ done # 运行脚本： ./aspera.sh 下载完成之后一定要检查数据完整性，不然分析白做 使用awk命令提取md5值，生成md5文件，运行以下命令检查文件完整性md5sum -c md5 md5文件如下图：  ## 2.4 sra2fastq 将sra文件转换为fastq文件 ```bash 创建脚本文件，个人觉得这样很方便。方便下次使用 #!bin/bash cat srr.list | while read id do nohup fastq-dump --split-3 $id -O ./ & done 2.5 质控（QC）这一步生成的是HTML文件，可以用浏览器打开查看。其中，各种颜色是各项标准分析结果：绿色代表"PASS"；黄色代表"WARN"；红色代表"FAIL"。 ls ../raw/*.fastq|xargs fastqc -t 10 -o ./ multiqc ./# 整合质控结果
去接头完成之后，一定要再次质控，确保数据tidy 同样，写个简单的脚本 #!/bin/bish # paired sequence # trim adapter paste ../clean/trim.log & done 2.7 hisat2比对我这里使用的是HISAT2，index是在HISAT2官网下载的。所以省去了建index的步骤。hisat2的index有八个文件。这里建议新手将参考基因组的fa文件，注释gtf文件和index文件放在同一文件夹下，不然报错很难找。 |
图中对于中心法则做了系统的总结，感兴趣可以click here。
RNA-seq的应用： 检测所有基因在特定条件下的表达水平(发育阶段、不同组织、正常与疾病、药物治疗)。 发现新基因，可变剪切，基因突变和基因融合等。
更多的细节可以通过阅读文章中的引用进一步了解。
在FASTQ文件中，一个序列通常由四行组成： 第一行由@开头，之后为序列的标识符和描述信息； 第二行为序列信息； 第三行以+开头，之后可以再次加上序列的标识和描述信息； 第四行为质量的分信息，与第二行的序列相对应，其长度必须与第二行相同。
3. FASTA：记录信息的格式 对于每条序列，首行以“>”开头，其之后是注释； 在首行之后，是以单字母标准编码表达的实际序列数据 FASTA的序列表达： 1） 核算编码：A,C,G,T,N,U,R,Y,K,M,S,W,B,D,H,V 2） 氨基酸编码：[A-Z],*（代表终止密码子）
4. SAM/BAM：高通量数据比对存放格式 SAM文件是由比对产生的以tab建分割的文件格式，BAM是SAM文件的二进制压缩版本。使用samtools view -S -b -o my.bam my.sam可以将SAM文件转换为BAM文件。 
5. BED：基因组浏览器常用格式 BED 文件格式提供了一种灵活的方式来定义的数据行，以用来描述注释信息。BED行有3个必须的列和9个额外可选的列。每行的数据格式要求一致。
(6). strand - 定义链的方向，’’+” 或者”-”。 (7). thickStart - 起始位置（The starting position atwhich the feature is drawn thickly）(例如，基因起始编码位置）。 (8). thickEnd - 终止位置（The ending position at whichthe feature is drawn thickly）（例如：基因终止编码位置）。 (9). itemRGB - 是一个RGB值的形式, R, G, B (eg. 255, 0,0), 如果itemRgb设置为’On”, 这个RBG值将决定数据的显示的颜色。 (10). blockCount - BED行中的block数目，也就是外显子数目。 (11). blockSize - 用逗号分割的外显子的大小, 这个item的数目对应于BlockCount的数目。 (12). blockStarts - 用逗号分割的列表, 所有外显子的起始位置，数目也与blockCount数目对应。 6. GFF/GTF GFF(General Feature Format)是基于Sanger GFF2的一种格式。GFF有9个必需字段，这些字段必须用制表符分隔。如果字段用空格而不是制表符分隔，则将不能正确显示。GTF (Gene Transfer Format, GTF2.2)是GFF的一种扩展格式。 
【本文地址】
| 今日新闻 |
| 推荐新闻 |
| 专题文章 |