Linux做RNA-seq上游分析基本流程

本机处理器两核心两线程（多核心多线程在比对时占优），操作系统是Linux，发行版是deepin20.1社区版
参考

• http://www.bio-info-trainee.com/2218.html
• https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50177

一、安装miniconda

miniconda相对于anaconda更为便捷和容易上手，anaconda类似于一个百宝箱，除了conda和一个固定版本的python，还有众多包和科学计算工具等，功能齐全但也操作复杂，文件较大，且有大部分库和包用不上，本文择使用miniconda

 wget -c  https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh

或者在官网下载

 wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh

windows系统也可下载，一般电脑均为64位，选择Python3的64位版本即可，32位电脑可下载32位版本，https://docs.conda.io/en/latest/miniconda.html

chmod 777 Miniconda3-latest-Linux-x86_64.sh   sh 
#修改权限
Miniconda3-latest-Linux-x86_64.sh   
#运行

一直Enter下去然后yes，最后询问是否初始化Miniconda3时，输入no

conda
#验证是否能正常使用

二、更换国内源

在国内下载后续软件速度大家定会抓狂，可以更换清华镜像源

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 --set show_channel_urls yes
pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple

三、创建rna虚拟环境

为了避免污染Linux原生终端，创建虚拟环境。关闭终端再次进入会发现前面多了base字样，默认自动进入base环境，还会影响打开终端速度，可以关闭自动进入虚拟环境。

conda create -n rna python=3 
#创建名为rna的软件安装环境

conda info --envs 
#查看conda环境,带星号的是当前环境

source activate rna 
conda activate rna 
#激活conda的rna环境

 source deactivate 
#注销当前的rna环境

conda config --set auto_activate_base false
#关闭自动进入虚拟环境

四、下载常用软件

conda search sra 
#可以先搜索软件

conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc 
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon

五、下载数据

新建一个名为 SRR_Acc_List.txt的文档，将SRR号码保存在文档内，一个号码占据一行。

wkd=/home/xzw/project/airway/
cd /home/xzw/project/airway
#先设置工作目录

cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} &);done
#用prefetch下载数据，也可直接手动下载

文章数据所在网址
https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP029245&o=acc_s%3Aa

六、检测数据质量

FastQC是一款基于Java的软件，可以快速多线程地对测序原始数据进行质量评估，其官网地址：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

cd $wkd/

fastq-dump *.sra
#转换用到fastq-dump

fastqc *.fastq
#这里选择默认，直接生成对应质控报告

multiqc *.zip
#multiqc将各个样本的质控报告整合为一份报告

电脑配置略低且数据少时建议分开转换，后续操作同理参考（血的教训）

基本格式：

 fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

#-o：用来指定输出文件的所在目录，默认是不新建目录的
#–noextract：输出的结果是.zip文件，默认自动解压缩，加上命令则不解压缩分析报告
#-f：用来强制指定输入文件格式，默认会自动检测。
#-c：用来指定一个contaminant文件，fastqc会把overrepresented sequences往这个contaminant文件里搜索。
#-q：进入安静模式，则不会出现以下提示：

结果解读

打开单独的html文件，根据分析结果评判是否需要修剪

上图即网页中的Summary部分就是整个报告的目录，整个报告分成若干个部分。合格会有个绿色的对勾，警告是个黄色的叹号，不合格是个红色的叉子。
具体解读可以参考下这篇文章：https://www.jianshu.com/p/96d0dc0c08eb?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation

七、原始数据修剪

这里用的是trim_galore，该软件是Trim Galore是对FastQC和cutadapt的包装，适用于所有高通量测序，且包括Nextera和smallRNA测序平台的双端和单端数据。
参考https://www.jianshu.com/p/7a3de6b8e503

dir=‘/home/xzw/project/clean/’   
/home/miniconda3/bin/trim_galore -q 22 --phred33 --length 36 -e 0.1 --stringency 3 -o $dir /home/xzw/project/airway/SRR957677
/home/miniconda3/bin/trim_galore -q 22 --phred33 --length 36 -e 0.1 --stringency 3 -o $dir /home/xzw/project/airway/SRR957678
/home/miniconda3/bin/trim_galore -q 22 --phred33 --length 36 -e 0.1 --stringency 3 -o $dir /home/xzw/project/airway/SRR957679
/home/miniconda3/bin/trim_galore -q 22 --phred33 --length 36 -e 0.1 --stringency 3 -o $dir /home/xzw/project/airway/SRR957680

#我电脑配置不高为了稳定所以就单个单个执行了，后面比对也是

我这里中途出了问题，可能是一开始同时下载软件的时候报过错没留意，这里上官网又重新下载了一遍。https://anaconda.org/

conda install -c grst trim_galore 

这里其实还有另外两个常用软件，感兴趣的可以了解下

#Trimmomatic：针对Illumina高通量测序平台设计的接头去除和低质量reads清洗软件
#cutadapt：删除adapter,primer. polyA尾等序列，去除低质量序列

八、比对

关于一个问题就是在判断是否需要比对时，这里有一篇文章大家可以借鉴一下https://www.jianshu.com/p/681e02e7f9af
这里用的是hisat2，速度上相对比较快且二类错误较少，还有bowtie和tophat2大家感兴趣的也可以学习一下，还有一个subread，在本文章后面计数会用到。

hisat2 -x /home/xzw/project/1hg19/genome -p 8 -U /home/xzw/project/clean/SRR957677_trimmed.fq -S /home/xzw/project/aligned/SRR957677_trimmed.sam 
hisat2 -x /home/xzw/project/1hg19/genome -p 8 -U /home/xzw/project/clean/SRR957678_trimmed.fq -S /home/xzw/project/aligned/SRR957678_trimmed.sam 
hisat2 -x /home/xzw/project/1hg19/genome -p 8 -U /home/xzw/project/clean/SRR957679_trimmed.fq -S /home/xzw/project/aligned/SRR957679_trimmed.sam 
hisat2 -x /home/xzw/project/1hg19/genome -p 8 -U /home/xzw/project/clean/SRR957680_trimmed.fq -S /home/xzw/project/aligned/SRR957680_trimmed.sam

基本格式：

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]

-x 参考基因组索引文件的前缀。
-1 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-U 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。 –sra-acc
输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
-S 指定输出的SAM文件。

九、计算RNA表达量

这里我用的是featureCounts，须安装subread（种草安利此工具，本文用的是2.0.1版本，链接为https://sourceforge.net/projects/subread/files/subread-2.0.1/）

先下载注释文件

wget ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz

featureCounts -T 8 -a /home/xzw/project/gencode.v29.annotation.gtf -o 77.count -t exon  /home/xzw/project/aligned/SRR957677.bam
featureCounts -T 8 -a /home/xzw/project/gencode.v29.annotation.gtf -o 78.count -t exon  /home/xzw/project/aligned/SRR957678.bam
featureCounts -T 8 -a /home/xzw/project/gencode.v29.annotation.gtf -o 79.count -t exon  /home/xzw/project/aligned/SRR957679.bam
featureCounts -T 8 -a /home/xzw/project/gencode.v29.annotation.gtf -o 80.count -t exon  /home/xzw/project/aligned/SRR957680.bam

基本格式：

featureCounts [options] <input.file>

-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对双端数据
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数

此外也可以用htseq-count实现操作，随后得到的counts数据可用于R做差异分析了。