RNA-Seq流程（包含部分Linux系统基础知识）
2019-04-05 本文已影响116人
大医院小牙医
- 写在前面之一：出于兴趣和知识整理的初衷，写了这篇博文，给网上众多RNA-Seq教程添砖加瓦，也给像我一点一点踩坑的小伙伴提供些许帮助，欢迎有兴趣的小伙伴批评指正讨论，由于 我的主业是👉▶▶ 口腔种植牙修复 ◀◀👈和👉▶▶ 阻生牙微创拔除 ◀◀👈，下班后还要看专业书和文献学习，时间和精力有限，无法帮助所有小伙伴解决生信分析问题，比起生信分析，还是欢迎大家咨询我的专业👉▶▶ 口腔种植牙修复 和 微创拔牙 ◀◀👈，我会尽量帮助大家解决牙齿的问题
- 写在前面之二：初学者一定要▶▶▶不求甚解◀◀◀，千万不要纠结于某个参数为什么那样写，举个🌰：不要纠结下面代码中 ~符号 的作用，就是这么规定的，programme language 也是一门语言，学习语言少问为什么，用多了就习惯了（类似英语中morning为什么是早上而不是晚上）
# RNA-Seq在得到表达矩阵后，通过R语言设计对比矩阵时，‘~‘符号的用法
# 构建实验设计矩阵
design<- model.matrix(~0+group)
所以...具有以下类似症状的小伙伴
▶没有任何软件/网站/脚本开发经验
▶傻傻分不清Windows系统分区/内存与硬盘
▶日常主要用电脑word/PPT/追剧/游戏(做过游戏进程注入修改数据的牛人请绕道)
按照本文代码+百度 / google 应该可以走一遍流程，不过会有很多坑等你踩...!@#$%^&*
但是，各位小伙伴▶要有信心◀，因为人生不就是重在折腾 奋斗吗！6个月前我也未真正接触过Linux系统，很多概念是经过实践后才逐渐由模糊到理解（比如Linux系统中一切皆文件的概念）
我的折腾 学习经历仅供参考，如果你也有类似的经历，恭喜你，直接开始生信分析，然后再补充点分子生物学基础理论知识：
▶很久~很久~以前考过C语言；
▶很久~很久~以前有脚本 折腾 开发经历：用VB+AutoIT脚本写过自动化交易策略+脚本程序；
▶校内网 用Flash ActionScript制作互动动画覆盖个人主页，从而实现动态博客的效果。校内网最初允许调用外部链接，后来就被禁了，再后来...
▶为了博士毕业 基于兴趣爱好，学了R和python；
▶为了拿基金 为了激发科研灵感，用Scrapy框架爬了NSFC网站。👏👏👏ps: 发现好多曾经的同学都已经手握国青了👍👍👍💪💪💪
- 写在前面之三：为了给
折腾自己学习知识的小伙伴建立整体框架感，按照本文代码完成后，RNA-Seq目录的结构如下：(有点进入正题的感觉了)
▶ RNA-Seq目录结构
teabio@tea-Bio:~/mount_Bio_disk/RNA-Seq$ tree -L 1
.
├── 1.QC
├── 2.mapping
├── 3.featureCounts
├── clean_data
├── index_file
├── PRJNA323422
├── ref
├── sra_files
▶ 如果有一定计算机编程基础，但是没有接触过Linux系统的小火伴，墙裂建议搞明白以下这些概念是什么鬼（如果真的是零基础小白，暂时忽略这些概念，因为从入门到放弃 是很容易做到的）：
- Linux文件系统结构
- Linux系统启动 / 用户登录 + 不同方式读取配置文件的顺序
- vim编辑器的简单使用
- source命令
- $PATH
👉实践中如果出现安装应用后无法启动，用 应用名称 --help 无法查看应用帮助文件 / which 应用名称 查不到路径，记得回来反复理解这些概念，等你哦😂😂😂👈
👉如果已经看不明白我在说什么了，说明已经开始说正事👈
- Linux文件系统结构
| 路径 | 作用 | 优先理解 |
|---|---|---|
| / | 是一切路径的起点 | ☆☆☆☆☆ |
| /usr | 存放与系统用户有关的文件和目录 | ☆☆ |
| /bin | bin是Binary的缩写,存放系统中最常用的可执行文件(二进制) | ☆☆ |
| /sbin | s就是Super User的意思,这里存放的是系统管理员使用的系统管理程序，如系统管理、目录查询等关键命令文件 | ☆☆ |
| /etc | 存放所有的系统管理所需要的配置文件和子目录 | ☆☆☆ |
| /home | 用户的主目录,Linux中,每个用户都有一个自己的目录,一般该目录名是以用户的账号命名的 | ☆☆☆☆☆ |
| /lib | 存放共享的库文件,包含许多被/bin和/sbin中程序使用的库文件 | ☆☆☆ |
| /mnt | 作为被挂载的文件系统的挂载点 | ☆☆ |
| /dev | dev是Device(设备),挂载硬盘要先查看这里 | ☆☆ |
| /boot | 这里存放的是Linux内核和系统启动文件,包括Grub、lilo启动器程序 | |
| /opt | 作为可选文件和程序的存放目录,主要被第三方开发者用来简易安装和卸载他们的软件 | |
| /proc | 一个虚拟的目录,它是系统内存的映射,通过直接访问此目录获取系统信息,比如cpuinfo存放cpu当前工作状态的数据 | |
| /root | 该目录为系统管理员,也称作超级权限者的用户主目录 | |
| /srv | 存放系统所提供的服务数据 | |
| /sys | 系统设备和文件层次结构 | |
| /tmp | 用来存放一些临时文件的,所有用户对此目录都有读写权限 |
- Linux系统启动 / 用户登录方式+不同方式读取配置文件顺序
Linux系统中，每个用户都有自己专属的运行环境，这个环境是由一组变量所定义，这些变量称之为环境变量。用户可以修改自己的环境变量以满足个性化需求，但是无权限修改其他用户的环境变量
（举个🌰：一栋楼里（一台服务器）相同户型的住户（权限一样的用户，一般是普通用户），每家每户都有一个唯一的门，进了这个门，就是进了自己家，也就是进了自己专属的运行环境，想怎么搞就怎么搞，不会破坏邻居家里的环境）
设置环境变量：export
显示某个环境变量：echo
令查看当前用户的环境变量：env
# 设置环境变量：
export 变量名="路径:$变量名" #临时变量，重启系统将失效。
#============================================#
#显示某个环境变量：
echo $变量名 # '$'符号在shell中是调用变量的含义
# 举个例子，显示路径变量（$PATH），这个变量也是初学Linux经常被坑的一个变量
echo $PATH
# 显示如下
/trainee/vip26/bin:/trainee/vip26/.local/bin:/trainee/vip26/miniconda2/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/trainee/vip26/.aspera/connect/bin
# 再举个例子，显示shell变量（$SHELL），运行/写 .sh脚本 时告诉系统调用哪个SHELL运行
echo $SHELL
# 显示如下
/bin/bash
#=============================================#
# 命令查看当前用户的环境变量
env
# 显示如下
XDG_SESSION_ID=32511
TERM=xterm-256color
SHELL=/bin/bash
HISTSIZE=3000
SSH_CLIENT=123.113.76.74 63488 22
SSH_TTY=/dev/pts/0
USER=vip26
MAIL=/var/mail/vip26
PATH=/trainee/vip26/bin:/trainee/vip26/.local/bin:/trainee/vip26/miniconda2/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/trainee/vip26/.aspera/connect/bin
PWD=/trainee/vip26
_=/usr/bin/env
当你启动（摁下电源开关）一台Linux电脑，登录一个用户，打开一个shell进行命令行操作，或者远程ssh登录一台已经开机（摁下电源+管理员登录）的服务器时，将会发生以下事情：
-
再次提醒：零基础的小白兔请跳过这部分内容，从“开始说正事”继续阅读此文，尽量避免 从入门到放弃 🚫⚠️
- 一台Linux计算机启动，首先是计算机系统启动，这个动作读取的环境配置文件是 → /etc/environment，是整个 计算机系统的环境配置文件
- 计算机启动，登录系统，获取最顶层shell，读取/etc/profile，是所有 用户的环境配置文件，这一步属于login shell
- 这一步注意区分两种不同登录方式
（1）login shell：用户通过终端，凭借用户名和密码登录控制台，这个的动作叫做login shell，也就是最终会调用login命令的操作都可称之为login shell，远程登录和本地登录都可以是login shell。
（2）non-login shell：用户在图形界面启动一个terminal，或者执行/bin/bash，/usr/bin/bash都属于non-login shell。
login shell与non-login shell的主要区别在于它们启动时会读取不同的配置文件，从而导致环境变量不一样。
login shell启动时首先读取/etc/profile全局配置，然后按顺序依次查找下面3个配置文件
- ~/.bash_profile
- ~/.bash_login
- ~/.profile
并且读取第1个被找到的并且可读的配置文件，后面的配置文件忽略不读取。login shell退出时读取并执行~/.bash_logout中的命令。
non-login shell启动时，交互式non-login shell启动时读取~/.bashrc资源文件。非交互式的non-login shell不读取上述所有配置文件，而是查找环境变量BASH_ENV，读取并执行BASH_ENV指向的文件中的命令。
- 通常我们要定制一些配置时，将配置写在 ~/.bashrc中，然后在~/.bash_profile中读取~/.bashrc，这样可以保证login shell和交互式non-login shell得到相同的配置（这句话是整理此文时才理解的❗️💪😂）。至于/etc/profile就不要轻易去改啦，毕竟会影响系统全局的配置。
通过以下2个表格来说明login shell 和 non-login shell 启动时配置文件读取顺序
- login shell
| 配置文件路径 | 启动顺序 |
|---|---|
| /etc/environment | 永远第1位 |
| /etc/profile | 2 |
| ~/.bash_profile | 3 |
| ~/.bash_login | 4 |
| ~/.profile | 5 |
| ~/.bashrc | 6 |
- non-login shell：从其父进程上继承过来环境变量
| 配置文件路径 | 启动顺序 |
|---|---|
| ~/.bashrc | 1 |
搞明白配置文件启动顺序了，以后找不到软件就知道去哪里排查问题了💪
👏👏👏开始说正事👏👏👏
坚持看到这里 或者 直接跳到这里的小伙伴们，基本具备了借助百度和谷歌独立完成以下内容的能力！
1. 安装软件
如果从未接触过Linux系统，或者使用OSX系统(苹果电脑操作系统)但未使用过终端(Terminal)的朋友，建议在实践中理解下面这个概念清单（不限于以下概念）
如果是一个全新的linux系统，可能要安装很多软件，常用软件在百度/科学上网查询均可找到安装方法，但由于不同Linux版本差异，安装不会完全一致，会碰到很多报错，初学者要有心理准备，以下举几个🌰
# 安装ruby
sudo apt-get install ruby-full
# 安装brew
/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
# miniconda安装
wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-4.5.11-Linux-x86_64.sh #下载
bash Miniconda2-4.5.11-Linux-x86_64.sh #安装
source ~/.bashrc #路径生效
#配置miniconda
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
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
2. 新建conda python2.7环境
- 在名为RNA-Seq的新建环境中安装软件，无论是使用脚本自动安装还是每个软件手动安装，一定要安装至RNA-Seq环境中。 理解conda管理环境 / 管理软件包
# 创建名为RNA-Seq的conda环境，python版本为2.7
conda create -n RNA-Seq python=2.7
- 安装软件
方法1：激活要安装软件的环境后（RNA-Seq），一个一个软件安装
# 激活RNA-Seq环境
source activate RNA-Seq
conda install 软件名
# 如果没有基础的同学在此处最好调用 '安装的程序 --help' 测试是否正确安装
# 举个栗子
conda install -y fastqc
# 显示：
(wes) vip26@VM-0-11-ubuntu:~$ fastqc --help
FastQC - A high throughput sequence QC analysis tool
SYNOPSIS
fastqc seqfile1 seqfile2 .. seqfileN
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam]
[-c contaminant file] seqfile1 .. seqfileN
DESCRIPTION
FastQC reads a set of sequence files and produces from each one a quality
control report consisting of a number of different modules, each one of
which will help to identify a different potential type of problem in your
data.
If no files to process are specified on the command line then the program
will start as an interactive graphical application. If files are provided
on the command line then the program will run with no user interaction
required. In this mode it is suitable for inclusion into a standardised
analysis pipeline.
The options for the program as as follows:
-h --help Print this help file and exit
方法2：通过sh脚本批量安装，不需要激活环境（RNA-Seq）
新建一个名为conda_auto_install_softwares.sh脚本，用文本编辑器打开，将如下内容拷贝至脚本中，保存退出。
#!/bin/bash
# 目的 → 实现软件自动安装至已经建立的conda的RNA-Seq环境下，RNA-Seq为通过 'conda create -n 环境名 python=2.7' 创建的环境
# 第一步 下载Miniconda
# OSX系统安装Miniconda：wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-latest-MacOSX-x86_64.sh
# Linux系统安装Miniconda：wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-latest-Linux-x86_64.sh
# 第二步 安装Miniconda （以ubuntu 18.04为例）
# bash Miniconda2-latest-Linux-x86_64.sh
# 第三步 配置镜像：
# conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
# 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中的RNA-Seq环境：conda create -n RNA-Seq python=2.7 #python版本根据不同需求设定
# 第五步 将conda_auto_install_softwares.sh 拷贝至~目录下
# 第六步 先输入：chmod 777 ~/conda_auto_install_softwares.sh #赋予权限
# 第七步 输入：~/conda_auto_install_softwares.sh #自动安装开始
# 第八步 输入：conda list -n RNA-Seq #查看以下软件是否正确安装至RNA-Seq环境中
conda install -y -n RNA-Seq sra-tools
conda install -y -n RNA-Seq fastqc
conda install -y -n RNA-Seq trim-galore
conda install -y -n RNA-Seq star
conda install -y -n RNA-Seq hisat2
conda install -y -n RNA-Seq bowtie2
conda install -y -n RNA-Seq subread
conda install -y -n RNA-Seq htseq
conda install -y -n RNA-Seq cutadapt
conda install -y -n RNA-Seq multiqc
conda install -y -n RNA-Seq samtools
conda_auto_install_softwares.sh 脚本获取权限
chmod 777 ~/conda_auto_install_softwares.sh
自动安装开始
~/conda_auto_install_softwares.sh
查看以下软件是否正确安装至RNA-Seq环境中
conda list -n RNA-Seq
验证软件安装是否成功同方法1
fastqc --help
3. 下载数据（下载sra文件）
以PRJNA323422中的4个小鼠测序文件为分析对象，下载4个sra文件:
SRR3589959
SRR3589960
SRR3589961
SRR3589962
勾选这4个sra文件，下载SRR_Acc_List.txt，放在~/RNA-Seq/PRJNA323422路径下。
注意：SRR_Acc_List.txt文件中的光标要在空行结束，如果在有内容的一行结束，则最后一行可能被忽略而不能被下载。
- 下载SRA数据，使用循环脚本提交至服务器后台下载：
nohup cat SRR_Acc_List.txt | while read num;do (prefetch -O ~/RNA_Seq/ $num); done &
- 将sra文件转换为fastq.gz格式
必须写 --split -3 参数 参数含义百度
此处实际用的是 fasterq-dump命令，比fastq-dump更好的一个命令，可以显示进度，可以多线程
fastq-dump --gzip --split-3 -O ~/RNA_Seq/ ./SRR3589959.sra ./SRR3589960.sra ./SRR3589961.sra ./SRR3589962.sra
# 此处实际用的是fasterq-dump命令, 命令格式和参数相似
4. 使用trim_galore去接头
- trim_galore安装：
方法1 ：在conda的RNA-Seq环境中安装(推荐初学者使用)
# 激活名为RNA-Seq的conda环境
conda activate RNA-Seq
#文件名一定要一模一样，否则在资源中查不到软件
conda resarch trim-galore
conda install -y trim-galore
# 验证是否安装成功
trim-galore --help
方法2：从以下网址下载至一个安装所有软件的路径下，后续建立软连接至环境变量指向的路径。
# 下载
wget -c https://github.com/FelixKrueger/TrimGalore/archive/0.6.0.tar.gz -O trim_galore.tar.gz #第二个参数 -O是大写的o
# 解压缩：
tar xvzf trim_galore.tar.gz
# 建立连接至环境变量
ln -s 软件实际路径 目标路径
使用trim_galore前需要确认需先安装fastqc和cutadapt是否已经安装，如果未安装，可以通过conda install 安装
trim_falore 参数：
-threads：设置线程数目
--phred33 使用ASCII+33质量得分作为Phred得分标准。默认选项
--fastqc 当剪切结束后用默认选项对结果文件进行fastqc分析
--fastqc_args "<ARGS>" 对fastQC设置参数。
--paired 设置双端序列
--dont_gzip 输出文件不压缩
--gzip 压缩输出文件，如果输入文件是压缩文件则自动压缩。
--length <INT> 设置低于数值长度的reads就抛弃掉，默认值20.
-q/--quality <INT> 切除质量得分低于设置值的序列。默认值20.
-a/--adapter <STRING> -a参数后为切除的接头序列
-a2/--adapter2 <STRING> 双端序列中read2所切除的接头序列，需配合'--paired'参数。
--rrbs 这是Trim Galore最独特的功能，也是我一开始找到使用这个软件的原因：特异性处理MspI所处理的RRBS样本数据(识别位点：CCGG)
trim_galore代码如下
nohup trim_galore -phred33 -q 26 -e 0.1 -length 36 --stringency 3 --paired -o ~/RNA-Seq ~/RNA_Seq/SRR3589959_1.fastq.gz ~/RNA_Seq/SRR3589959_2.fastq.gz ~/RNA_Seq/SRR3589960_1.fastq.gz ~/RNA_Seq/SRR3589960_2.fastq.gz ~/RNA_Seq/SRR3589961_1.fastq.gz ~/RNA_Seq/SRR3589961_2.fastq.gz ~/RNA_Seq/SRR3589962_1.fastq.gz ~/RNA_Seq/SRR3589962_2.fastq.gz &
5. 质控（Quality Control）
- Fastqc用法：
FastQC参数：
-o --outdir:输出路径
--extract：结果文件解压缩
--noextract：结果文件压缩
-f --format:输入文件格式.支持bam,sam,fastq文件格式
-t --threads:线程数
-c --contaminants：制定污染序列。文件格式 name[tab]sequence
-a --adapters：指定接头序列。文件格式name[tab]sequence
-k --kmers：指定kmers长度（2-10bp,默认7bp）
-q --quiet： 安静模式
# 使用循环脚本提交至服务器后台下载： -t：2线程; -o：输出路径
nohup fastqc -t 2 -o ~/RNA_Seq/1.QC/ ~/RNA_Seq/SRR3589959_1.fastq.gz ~/RNA_Seq/SRR3589959_2.fastq.gz ~/RNA_Seq/SRR3589960_1.fastq.gz ~/RNA_Seq/SRR3589960_2.fastq.gz ~/RNA_Seq/SRR3589961_1.fastq.gz ~/RNA_Seq/SRR3589961_2.fastq.gz ~/RNA_Seq/SRR3589962_1.fastq.gz ~/RNA_Seq/SRR3589962_2.fastq.gz &
fastqc_sequence_counts_plot.png
fastqc_per_base_sequence_quality_plot.png
fastqc_per_sequence_quality_scores_plot.png
fastqc_sequence_length_distribution_plot.png
fastqc_sequence_duplication_levels_plot.png
6. 比对（mapping）
采用Hisat2进行比对，需要自行下载index 文件，或者自己构建index文件，推荐下载，自己构建时间非常~非常~非常久，推荐使用迅雷直接下载人index / 小鼠index文件（hg38的index约4.1G）,下载完成后，新建index_file文件夹
mkdir ~/mount_Bio_disk/RNA-Seq/index_file
将下载的mm10.tar.gz移至~/mount_Bio_disk/RNA-Seq/index_file路径下，
mv mm10.tar.gz下载路径 ~/mount_Bio_disk/RNA-Seq/index_file
解压缩解包：
cd filepath #存放index文件的路径
tar -zxvf *.tar.gz #解压解包
新建文件夹 2.mapping,进行比对，代码如下：
# hisa2 -p 整数：线程数
# samtools sort @整数：线程数，默认单线程
for i in {59..62};do (hisat2 -p 50 -x ./index_file/mm10/genome -1 ./clean_data/SRR35899${i}.sra_1_val_1.fq -2 ./clean_data/SRR35899${i}.sra_2_val_2.fq | samtools sort -@ 50 -o ./2.mapping/SRR35899${i}.sort.bam);done
代码运行后显示：
29900227 reads; of these:
29900227 (100.00%) were paired; of these:
2356332 (7.88%) aligned concordantly 0 times
24105272 (80.62%) aligned concordantly exactly 1 time
3438623 (11.50%) aligned concordantly >1 times
----
2356332 pairs aligned concordantly 0 times; of these:
155034 (6.58%) aligned discordantly 1 time
----
2201298 pairs aligned 0 times concordantly or discordantly; of these:
4402596 mates make up the pairs; of these:
2756804 (62.62%) aligned 0 times
1165864 (26.48%) aligned exactly 1 time
479928 (10.90%) aligned >1 times
95.39% overall alignment rate
[bam_sort_core] merging from 20 files and 2 in-memory blocks...
52061584 reads; of these:
52061584 (100.00%) were paired; of these:
3872015 (7.44%) aligned concordantly 0 times
42501088 (81.64%) aligned concordantly exactly 1 time
5688481 (10.93%) aligned concordantly >1 times
----
3872015 pairs aligned concordantly 0 times; of these:
265752 (6.86%) aligned discordantly 1 time
----
3606263 pairs aligned 0 times concordantly or discordantly; of these:
7212526 mates make up the pairs; of these:
4560210 (63.23%) aligned 0 times
1887381 (26.17%) aligned exactly 1 time
764935 (10.61%) aligned >1 times
95.62% overall alignment rate
[bam_sort_core] merging from 34 files and 2 in-memory blocks...
35956411 reads; of these:
35956411 (100.00%) were paired; of these:
2675945 (7.44%) aligned concordantly 0 times
29017688 (80.70%) aligned concordantly exactly 1 time
4262778 (11.86%) aligned concordantly >1 times
----
2675945 pairs aligned concordantly 0 times; of these:
169988 (6.35%) aligned discordantly 1 time
----
2505957 pairs aligned 0 times concordantly or discordantly; of these:
5011914 mates make up the pairs; of these:
3072626 (61.31%) aligned 0 times
1343911 (26.81%) aligned exactly 1 time
595377 (11.88%) aligned >1 times
95.73% overall alignment rate
[bam_sort_core] merging from 24 files and 2 in-memory blocks...
42622324 reads; of these:
42622324 (100.00%) were paired; of these:
3261521 (7.65%) aligned concordantly 0 times
34871791 (81.82%) aligned concordantly exactly 1 time
4489012 (10.53%) aligned concordantly >1 times
----
3261521 pairs aligned concordantly 0 times; of these:
203644 (6.24%) aligned discordantly 1 time
----
3057877 pairs aligned 0 times concordantly or discordantly; of these:
6115754 mates make up the pairs; of these:
4114565 (67.28%) aligned 0 times
1442337 (23.58%) aligned exactly 1 time
558852 (9.14%) aligned >1 times
95.17% overall alignment rate
[bam_sort_core] merging from 26 files and 2 in-memory blocks...
查看比对、排序后的bam文件
# 查看比对、排序后的bam文件
samtools view -h ./SRR3589959.sort.bam
# 第一行显示如下：
@HD VN:1.0 SO:coordinate #未排序显示为 unsorted
7. 计数
- 建立 3.featureCounts 目录
# 建立 3.featureCounts 目录
mkdir 3.featureCounts
- 下载annotation
# 下载annotation
wget ftp://ftp.ensembl.org/pub/release-95/gtf/mus_musculus/Mus_musculus.GRCm38.95.chr.gtf.gz -o ~/mount_Bio_disk/RNA-Seq/ref/Mus_musculus.GRCm38.95.chr.gtf.gz
- 解压解包
# 解压解包
unzip ~/mount_Bio_disk/RNA-Seq/ref/Mus_musculus.GRCm38.95.chr.gtf.gz
- 计数，输出2个文件：all.id.txt 和 all.id.txt.summary，all.id.txt 是结果，用于R语言分析。
# 计数，输出2个文件：all.id.txt 和 all.id.txt.summary，all.id.txt 是结果，用于R语言分析。
(RNA-Seq) teabio@tea-Bio:~/mount_Bio_disk/RNA-Seq/3.featureCounts$ featureCounts -T 50 -p -t exon -g gene_id -a ~/mount_Bio_disk/RNA-Seq/ref/Mus_musculus.GRCm38.95.chr.gtf -o ./all.id.txt ../2.mapping/*sort.bam
###########################命令运行后显示内容如下##############################
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
v1.6.4
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 4 BAM files ||
|| P SRR3589959.sort.bam ||
|| P SRR3589960.sort.bam ||
|| P SRR3589961.sort.bam ||
|| P SRR3589962.sort.bam ||
|| ||
|| Output file : all.id.txt ||
|| Summary : all.id.txt.summary ||
|| Annotation : Mus_musculus.GRCm38.95.chr.gtf (GTF) ||
|| Dir for temp files : . ||
|| ||
|| Threads : 50 ||
|| Level : meta-feature level ||
|| Paired-end : yes ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
|| Chimeric reads : counted ||
|| Both ends mapped : not required ||
|| ||
\\============================================================================//
//================================= Running ==================================\\
|| ||
|| Load annotation file Mus_musculus.GRCm38.95.chr.gtf ... ||
|| Features : 818414 ||
|| Meta-features : 54752 ||
|| Chromosomes/contigs : 22 ||
|| ||
|| Process BAM file SRR3589959.sort.bam... ||
|| Paired-end reads are included. ||
|| Assign alignments (paired-end) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total alignments : 37290702 ||
|| Successfully assigned alignments : 22081740 (59.2%) ||
|| Running time : 0.35 minutes ||
|| ||
|| Process BAM file SRR3589960.sort.bam... ||
|| Paired-end reads are included. ||
|| Assign alignments (paired-end) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total alignments : 64097834 ||
|| Successfully assigned alignments : 38840169 (60.6%) ||
|| Running time : 0.58 minutes ||
|| ||
|| Process BAM file SRR3589961.sort.bam... ||
|| Paired-end reads are included. ||
|| Assign alignments (paired-end) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total alignments : 45141073 ||
|| Successfully assigned alignments : 26707618 (59.2%) ||
|| Running time : 0.48 minutes ||
|| ||
|| Process BAM file SRR3589962.sort.bam... ||
|| Paired-end reads are included. ||
|| Assign alignments (paired-end) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total alignments : 51938946 ||
|| Successfully assigned alignments : 31372932 (60.4%) ||
|| Running time : 0.44 minutes ||
|| ||
|| ||
|| Summary of counting results can be found in file "./all.id.txt.summary" ||
|| ||
\\============================================================================//
- 完成后，查看all_id.txt文件：
# Program:featureCounts v1.6.4; Command:"featureCounts" "-T" "50" "-p" "-t" "exon" "-g" "gene_id" "-a" "/home/teabio/mount_Bio_disk/RNA-Seq/ref/Mus_musculus.GRCm38.95.chr.gtf" "-o" "./all.id.txt" "../2.mapping/SRR3589959.sort.bam" "../2.mapping/SRR3589960.sort.bam" "../2.mapping/SRR3589961.sort.bam" "../2.mapping/SRR3589962.sort.bam"
Geneid Chr Start End Strand Length ../2.mapping/SRR3589959.sort.bam ../2.mapping/SRR3589960.sort.bam ../2.mapping/SRR3589961.sort.bam ../2.mapping/SRR3589962.sort.bam
ENSMUSG00000102693 1 3073253 3074322 + 1070 0 0 0 0
ENSMUSG00000064842 1 3102016 3102125 + 110 0 0 0 0
ENSMUSG00000051951 1;1;1;1;1;1;1 3205901;3206523;3213439;3213609;3214482;3421702;3670552 3207317;3207317;3215632;3216344;3216968;3421901;3671498 -;-;-;-;-;-;- 6094 15 16 17 6
ENSMUSG00000102851 1 3252757 3253236 + 480 0 0 0 0
ENSMUSG00000103377 1 3365731 3368549 - 2819 0 0 0 0
ENSMUSG00000104017 1 3375556 3377788 - 2233 0 0 0 0
ENSMUSG00000103025 1 3464977 3467285 - 2309 0 0 0 0
ENSMUSG00000089699 1;1 3466587;3513405 3466687;3513553 +;+ 250 0 0 0 0
ENSMUSG00000103201 1 3512451 3514507 - 2057 0 0 0 0
ENSMUSG00000103147 1 3531795 3532720 + 926 1 0 0 3
ENSMUSG00000103161 1 3592892 3595903 - 3012 0 0 0 0
ENSMUSG00000102331 1;1 3647309;3658847 3650509;3658904 -;- 3259 1 0 0 0
ENSMUSG00000102348 1 3680155 3681788 + 1634 0 0 0 0
ENSMUSG00000102592 1 3752010 3754360 + 2351 0 0 0 0
ENSMUSG00000088333 1 3783876 3783933 - 58 0 0 0 0
ENSMUSG00000102343 1;1;1;1;1 3905739;3984225;3985160;3985160;3986147 3906134;3984298;3985984;3985351;3986215 -;-;-;-;-
1364 0 0 0 3
ENSMUSG00000025900 1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1 3999557;4007656;4019070;4024736;4041888;4092617;4119668;4120015;4142612;4147812;4148612;4163855;4170205;4197534;4206660;4226611;4228443;4231053;4243133;4243417;4243543;4245031;4261527;4267469;4284766;4290846;4292926;4292981;4311270;4344146;4351910;4351910;4351910;4351910;4352202;4352202;4352202;4352202;4360200;4409170;4409170;4409170 3999617;4007737;4019148;4024890;4042107;4092780;4119712;4120073;4142766;4147963;4148744;4163941;4170404;4197641;4206837;4226823;4228619;4231144;4243262;4243448;4243619;4245106;4261605;4267620;4284898;4293012;4293012;4293012;4311433;4350091;4352081;4352081;4352081;4352081;4352837;4352837;4352837;4352837;4360314;4409187;4409241;4409241 -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;- 12311 1 0 0 1
ENSMUSG00000102948 1;1 4256234;4258847 4256427;4260519 -;- 1867 0 0 0 0
ENSMUSG00000104123 1 4363346 4364829 - 1484 0 0 0 0
ENSMUSG00000025902 1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1 4490931;4491250;4491390;4491713;4492457;4492458;4492465;4492467;4493100;4493100;4493100;4493100;4493768;4493772;4493772;4493772;4493772;4495136;4495136;4495136;4495136;4496291;4496291;4496291;4496291;4496291;4496291 4492668;4492663;4492668;4492668;4493604;4492668;4492668;4492668;4493490;4493735;4493490;4493466;4493863;4493863;4493863;4493863;4493863;4495942;4495198;4495942;4495942;4496396;4496330;4497354;4496363;4496757;4496413 -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;- 4772 55 63 58 91
ENSMUSG00000104238 1;1;1;1 4496551;4497474;4498007;4499470 4499378;4497654;4498211;4499558 +;+;+;+ 2917 0 0 0
0
ENSMUSG00000102269 1;1 4522905;4524446 4523603;4526737 +;+ 2991 3 9 3 1
ENSMUSG00000096126 1 4529017 4529123 + 107 0 0 0 0
ENSMUSG00000103003 1 4534837 4535286 - 450 0 0 0 0
ENSMUSG00000104328 1;1 4583129;4585937 4585585;4586252 -;- 2773 0 0 0 0
ENSMUSG00000102735 1 4610471 4611406 + 936 0 0 0 0
ENSMUSG00000098104 1 4687934 4689403 - 1470 141 255 146 206
8.R语言中进行分析和绘图
由于之前为了按时毕业 兴趣爱好做过几个完整做GEO数据分析(主要在R
中完成)，所以R语言的数据分析和绘图并不是本文的重点，会在后续整理中总结，敬请期待...
