酵母(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.
网友评论