B站:RNA-seq转录组数据分析入门实战
1 linux常用命令
touch text.txt #新建文件
rm -rf /var/log/httpd/access #将会删除/var/log/httpd/access目录以及其下所有文件、文件夹
rm -f *html #删除所有html格式文件
rm -f *zip #删除所有zip格式文件
tar zxvf #解压tar.gz文件
tar jxvf samtools-1.11.tar.bz2 #解压tar.bz2文件
#gzip压缩后的格式为:*.gz,这种压缩方式不能保存原文件;
gzip buodo #压缩
gunzip buodo.gz #解压
.bashrc中可以添加一些常用快捷方式,也可以将软件添加到环境变量
如何用vi编辑和保存文件
vi ~/.bashrc #打开.bashrc,然后按a进行编辑
alias用来设置快捷方式
比如:
alias ..='cd ..' 返回上一级
alias ..='cd ../..' 返回上上级
设置完成后按ESC,然后按:wq保存退出。
. .bashrc 启动bashrc文件,刚才设置的就可以用了。
2 conda的下载和软件的安装
conda的下载
复制miniconda3的lunux版本链接地址,-c是支持断点下载,就是网断了再次连接可以继续下载。
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh #下载好后,进行解压安装
vi ~/.bashrc #然后进行查看,是否已经将miniconda3的bin文件添加到了环境变量中去。
source .bashrc #执行,将miniconda3添加到环境变量
导入重要的channel,因为后期会用到
conda config --add channels r
conda config --add channels conda--forge
conda config --add channels bioconda
image.png
软件安装,升级和卸载
conda isntall fastqc #安装fastqc
fastqc -help #判断是否安装成功
conda install samtools=1.4.1 #先安装一个旧版本
which samtools #查看安装的位置
conda update samtools #查看最新版本并进行安装
conda remove samtools #卸载软件
编译安装软件
image.png
先下载一个samtools安装包,链接从官网复制
wget -c https://github.com/samtools/samtools/releases/download/1.11/samtools-1.11.tar.bz2
tar jxvf samtools-1.11.tar.bz2 #解压
cd samtools-1.11/ #cd到当前软件目录下,可以看到一个configure文件,然后进行配置。
1. 配置(在root下新建一个samtools,最后将软件安装在这里)
./configure --prefix='/root/samtools'
报错处理:
configure: error: no acceptable C compiler found in $PATH 显示错误主要是没有C编译器。
安装C编译器 :yum -y install gcc
2. 编译
make
报错处理: (问题都是一样,包没有安装,后续百度安装即可)
2.1 CentOS编译安装软件过程中遇到zlib.h: No such file or directory,使用命令:yum install zlib-devel 解决问题。
2.2 curses.h: No such file or directory 错误解决办法yum install ncurses-devel ncurses
2.3 centos编译时报错:error: bzlib.h: No such file or directory
解释:bzlib就是bzip2包 .
首先查看你yum源中的bzip2包,命令为:yum search bzip2。
然后安装你yum源中的bzip2包,执行命令:yum install bzip2-devel.x86_64
3. 安装
make install
添加环境变量
安装完成后,如果不在当前路径运行samtools报错,那就需要添加环境变量
vi ~/.bashrc
image.png
. ~/.bashrc #加载环境变量
3 数据预处理
3.1准备工作
image.png
1.构建项目目录 用tree将rnaseq下的文件进行树状展示
image.png
2. 参考序列的下载 参考基因组fasta,注释信息gtf/gff
Genecode数据库中下载人类的参考基因组和注释信息文件,然后上传到服务器的00ref文件中。
image.png
image.png
gunzip *gz #然后在00ref中对两个gz文件进行解压
3. 原始数据的上传
原始数据通过WINscp软件进行上传。
image.png
然后,通过md5值检测数据的完整性
#整个操作就在01raw_data文件下进行
md5sum *gz >md5.txt #给所有的gz文件生成md5值,储存到md5.txt中。
cat md5.txt #进行查看
md5sum -c md5.txt #比对gz文件中已有的md5值和生成的md5是否匹配,匹配返回OK.
image.png
3.2 质量控制
image.png
fastqc sample1_R1.fastq.gz #单一样本质控
fastqc sample*gz #批量质控开头为sample,gz结尾的文件
for i in 'ls *gz';do fastqc $i;done #批量质控,遍历该文件夹下所有gz结尾的文件,do fastqc,然后done结束
并行批量质控
#xargs命令是将上一个命令的输出作为下一个命令的输入
#nonhup &是将运行程序挂在后台,这样就可以把每个都挂在后台,达到批量的目的。
#>fastqc.sh是将这个应用保存在shell脚本
ls *gz |xargs -I[] echo 'nohup fastqc []&' >fastqc.sh
cat fastqc.sh #可以看到就是将四个程序分别挂到后台运行
#nohup fastqc sample1_R1.fastq.gz&
#nohup fastqc sample1_R2.fastq.gz&
#nohup fastqc sample2_R1.fastq.gz&
#nohup fastqc sample2_R2.fastq.gz&
bash fastqc.sh #运行shell脚本
top -c #查看运行进程,ctrl+推出查看
在运行的结果中,可以看到每个样本有一个html(保存质控信息)和zip(保存的文本内容)文件
image.png
多个质控结果的查看
首先运用conda安装multiqc,multiqc支持多种结果的整合,同时展示多个结果文件,可以自动检测当前文件下可以合并的文件。
multiqc ./
image.png
然后将multiqc_report.html 下载到本地进行查看。
3.3 质量过滤
image.png
image.png
接头序列的选择:图中可以看到是Truseq adapter ,所以就根据单端或者双端选择Truseq3。
另外去接头的参数选择TURE,保留下来的reads较多。
image.png
image.png
用trimmomatic进行质量过滤,在存在fastq.gz的文件下进行运行
#在01raw_data下运行
# PE是双端测序,threads 是线程,../02clean_data/是生成文件的保存位置
#sample1_paired_clean是生成的成对read的文件,sample1_unpair_clean是生成的去除了一条read的不成对的文件
#ILLUMINACLIP接头的去除,先说明位置,在说明怎么去(2:30:10:1:true)
#SLIDINGWINDOW:4:20 滑窗设置的是4bp,如果低于20则被去除掉
#MINLEN:50 将read小于50的去除
#TOPHRED33 将reads的碱基质量体系转为phred-33
trimmomatic PE -threads 4 \
sample1_R1.fastq.gz sample1_R2.fastq.gz \
../02clean_data/sample1_paired_clean_R1.fastq.gz \
../02clean_data/sample1_unpair_clean_R1.fastq.gz \
../02clean_data/sample1_paired_clean_R2.fastq.gz \
../02clean_data/sample1_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/root/miniconda3/share/trimmomatic-0.39-1/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
80%的reads被留下
接下来就用paired文件进行后的序列比对和定量分析
image.png
4 序列比对
基于基因组和基于转录本的序列比对
转录组和基因组比对区别
image.png
image.png
新建文件夹 mkdir arab_STAR_genome保存索引文件
先解压00ref里面的基因组和注释文件用gunzip,然后建立索引
整个操作是在rnaseq文件夹下完成
建立索引
# 工作模式是genomeGenerate,arab_STAR_genome生成的索引文件的位置
#genomeFastaFiles和sjdbGTFfile分别表示参考基因组和注释信息位置
#sjdbOverhang是可变剪切的预处理,默认值是100。也可以用read长度减1的方式来做,我们的reads是150.
STAR --runThreadN 6 --runMode genomeGenerate \
--genomeDir arab_STAR_genome \
--genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \
--sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
--sjdbOverhang 149
STAR比对
#genomeDir索引的位置,zcat对所有gz文件解压,readFilesIn数据清洗后要读入的进行比对数据
#outFileNamePrefix输出位置,并给输出文件加上sample2_ 前缀
#outSAMtype BAM SortedByCoordinate将SAM格式转换为二进制的BAM格式,并进行排序
#outBAMsortingThreadN 5 排序采用5线程
#quantMode TranscriptomeSAM GeneCounts将基因组比对结果转换为转录组结果,并计数每个基因的reads数。
#TranscriptomeSAM 为了使用RSEM进行定量分析做准备。
STAR --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \
--readFilesIn 02clean_data/sample2_paired_clean_R1.fastq.gz \
02clean_data/sample2_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/sample2_ \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN 5 \
--quantMode TranscriptomeSAM GeneCounts
输出比对bam文件
5 表达定量
image.png
两种方式都需要先比对再定量
5.1 STAR+RSEM进行定量
mkdir arab_RSEM #在rnaseq下新建文件夹arab_RSEM保存resm参考文件
conda install RSEM #先安装RSEM,需要时间比较久
#rsem prepare reference
#需要基因组和注释文件,提取参考基因组中每一个转录组的fasta信息
#arab_RSEM/arab_rsem 保存在arab_RSEM,加arab_rsem前缀。
rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
00ref/TAIR10_Chr.all.fasta \
arab_RSEM/arab_rsem
mkdir 04rsem_out#新建RSEM的输出文件
#表达定量
#因为数据是双端测序的所以用paired-end,不需要输出bam文件,所以no-bam-output
#-p 5线程是5,-q 是输入bam文件的保存位置
#arab_RSEM/arab_rsem \之前准备工作的索引位置
#输出到04rsem_out,加sample2_rsem前缀
rsem-calculate-expression --paired-end --no-bam-output \
--alignments -p 5 \
-q 03align_out/sample2_Aligned.toTranscriptome.out.bam \
arab_RSEM/arab_rsem \
04rsem_out/sample2_rsem
less sample2_rsem.genes.results |head #查看基于基因的定量结果
image.png
5.2 STAR+featurecounts进行定量
可以STAR后,直接用featurecounts进行定量,不做RSEM,将得到的out_counts.txt文件下载到本地,然后进行后续处理
image.png
就在STAR后的结果文件03align_out中运行
conda install subread #安装featurecounts
#-p就是双端测序,-a后面跟基因组文件,-o后面跟输出的文件
#-T是线程,-t根据外显子进行定量,-g输出gene
#sample*_Aligned.sortedByCoord.out.bam 对当前文件中所有sample开头和bam结尾的文件进行定量
featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 6 -t exon -g gene_id sample*_Aligned.sortedByCoord.out.bam
out_counts.txt文件中包含了counts数
out_counts.txt
6 将RSEM定量后的结果转换为表达矩阵
image.png
对RSEM定量分析后得到的sample1_rsem.genes.results 文件进行合并,构建表达矩阵
image.png
#对所有当前目录下_rsem.genes.results结尾的文件进行合并,输出到output.matrix
rsem-generate-data-matrix *_rsem.genes.results >output.matrix
head output.matrix #查看
wc -l output.matrix #查看文件行数