RNA-seq（一） rawdata分析

酵母（Saccharomyces cerevisiae ）细胞，单条线性染色体的RNA-seq数据分析

1 prefetch下载数据

使用：把sratoolkit.2.10.8/bin加入环境变量使用相对路径；调用prefetch的绝对路径。
如果数据过大，可加参数--max-size 100G

vim download_data.sh
#!/bin/bash
for i in `seq 4 1 9`
do
  prefetch-orig.2.10.8 `srapath-orig.2.10.8 SRR705970${i}`
done
#for循环，下载6个数据
qsub -N download -cwd download_data.sh -q g5.q 
#-N：提交的任务名，-cwd：在当前路径下提交，日志会输出在当前目录下；-q ：指定提交的队列；-o：重定向标准输出；-e：重定向标准错误信息输出（不设置-e，-o会生成提交任务id的.o和.e文件）

vi文本编译器，vim 是 vi 的增强版（vi Improved），与 vi 编辑器完全兼容，而且实现了很多增强功能。详细介绍：Linux vi/vim | 菜鸟教程 (runoob.com)
不要直接在服务器终端跑大型命令，把命令提交给服务器。
qsub命令：用来从登录节点上向计算节点进行任务投递。SGE(SGE, Sun Grid Engine）集群管理工具可以用来提交批处理作业，SGE支持单或多节点作业，它将用户投递的任务进行排队，然后将任务交给能够运行的结算节点执行，工作流程可以分为四步：接受用户提交的任务；任务运行前将任务放到一个存储区域；发送任务到一个执行设备，并监控任务的运行；运行结束写出结果并记录运行日志。
qstat：可查看提交任务的运行状态，"r"代表任务正在执行；"qw"代表任务在等待执行，一旦有计算资源了会马上运行；Eqw，投递任务出错；
qdel：删除提交的任务；qdel+ job_id

2 fastq-dump 格式转换

fastq-dump是sratoolkit.2.10.8的子工具，我们需要把下载的数据sra格式转换为fastq格式。
FASTQ格式是将核酸序列与其测序质量得分信息合并在一起的文本格式。质量得分是指该碱基的错误概率的对数值。
在FASTQ文件中，一个序列通常由四行组成：
第一行：@开头，后面是序列描述信息（与FASTA格式的描述行类似）
第二行：序列信息
第三行：+开头，后面是序列描述信息
第四行：碱基质量得分，与第二行的序列相对应，长度也与第二行相同
Sequence Read Archive (SRA)数据库，用于储存来自不同高通量测序平台的原始序列数据和比对信息。在提交测序数据到SRA数据库时，每一个样本会获得一个以SRR开头的数字编号：SRRxxx，数据格式为 .sra 。
生物信息学常见数据格式汇总(https://zhuanlan.zhihu.com/p/67402565)

vim sra2fastq.sh
#!/bin/bash
for i in `seq 4 1 9`
do
  fastq-dump-orig.2.10.8 --split-3 -O fastq --gzip SRR705970${i}/SRR705970${i}.sra
done
qsub -N fastq -cwd sra2fastq.sh

--split-3：将双端测序分为两份，放在不同的文件，但是对于一方有而一方没有的reads会单独放在一个文件夹里；（懒人做法，无论单端还是双端的SRA文件，都可加--split-3）。
--gzip：输出gz格式,，节省空间。

3 fastqc质控

vim fastqc.sh
#!bin/bash
for i in `seq 4 1 9`
do
  fastqc SRR705970${i}_1.fastq.gz SRR705970${i}_2.fastq.gz
done
qsub -N  fastqc  -cwd fastqc.sh

评价所获得的数据能不能用，这一步就叫做质量控制（Quality Control, QC）。详细讲解：RNASeq原始数据质量控制-FastQC - 知乎 (zhihu.com)

multiqc *fastqc.zip

multiqc 可整合QC质控结果，快速查看结果

4 trim_galore 过滤

可参阅Babraham Bioinformatics - Trim Galore!

vim trim_galore.sh
#!bin/bash
for i in `seq 4 1 9`
do
  trim_galore --output_dir clean_data -q 25 --phred33 --length 36 --paired SRR705970${i}_1.fastq.gz SRR705970${i}_2.fastq.gz
done
qsub -N trim_galore   -cwd trim_galore.sh

--length：设定输出reads长度阈值，小于设定值会被抛弃。
-q 25：低质量区域，若碱基质量值<25，后面的碱基被cut
Overview of removed sequences 记录、统计去掉的reads
--paired 如果序列对中的一个(或两个)小于设置的长度截断值，则移除整个序列对。 removes entire sequence pairs if one (or both) of the two reads became shorter than the set length cutoff.