不用linux转录组数据分析,RNA-seq转录组数据分析

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 #查看文件行数