RNA-seq第一部分(原始数据下载,提取fastq文件,fastqc质控)

原始数据下载地址:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50177

一.SRA数据下载

打开网页按以下步骤进行

将连接复制到迅雷下载,可以淘宝临时买个会员下载,速度较快,建议使用该方式下载。

也可复制地址用wget下载,速度贼慢,不建议。

1
wget 	https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR957677/SRR957677.1

数据下载后通过xftp上传到服务器

二.使用sratoolkit从sra文件提取fastq文件

1
fastq-dump --gzip --split-files SRR957677 #从SRR957677分别提取左端和右端reads

–gzip 生成压缩格式为gz的fastq文件

SRR957677是一个PE测序结果,所以,需要–split-files参数可以将其分解为两个fastq文件。如果不加该参数,则只有1个fastq文件(包含了两端测序的结果)

有四个文件分别是SRR957677,SRR957678,SRR957679,SRR957680,;可以一个一个地提取fastq文件,也可以批量提取,有多种方法批量提取。

1)for循环

检查当前文件夹是否有SRR文件,*是匹配任意字符

1
ls SRR*

检查循环是否正确,echo打印出脚本

1
for i in `ls SRR*`;do echo fastq-dump --gzip --split-files $i ;done 

将上述命令中的echo去掉直接运行即可

1
for i in `ls SRR*`;do fastq-dump --gzip --split-files $i ;done

2)while循环

新建文件列表这里命名SRR_list

1
ls SRR* >SRR_list

查看文件列表是否正确

1
cat SRR_list

检查循环是否正确,echo打印出脚本

1
cat SRR_list | while read id;do echo fastq-dump --gzip --split-files $id ;done

或者

1
ls SRR95* | while read id;do echo fastq-dump --gzip --split-files $id ;done#这里改为SRR95*是因为当前有个文件名为SRR_list,如果还是SRR*则会把这个文件也匹配进去

3)xargs

1
ls SRR95* | xargs -i echo echo fastq-dump --gzip --split-files {}

4)awk

1
awk '{print "nohup fastq-dump --gzip --split-files "$0 " &"}' SRR_list >fastq_dump.sh

加nohup & 是为了放后台,因为每一个子命令都有nohup &所以这些命令是并行,>fastq_dump.sh是将命令写入fastq_dump.sh中。

1
bash fastq_dump.sh

查看命令是否挂后台,top查看进程,-u 参数接用户名

1
top -u liyue

可以看到有四个fastq-dump命令运行中

还可以使用其他命令查询

1
ps -aux | grep "fastq_dump"

注意以上命令有个小问题,-split-files参数可以不加,因为这些测序文件是单端测序,所以即使加了该参数也只有一个fastq文件生成,如果是双端测序则会生成左端,右端两个fastq文件,不过加了也没影响

生成以下fastq文件

三.质控

1.质控使用fastqc软件

该软件可直接用conda安装,这里使用循环进行质控

1
ls *fastq.gz | awk '{print "nohup fastqc " $0  " &"}' >fastqc.sh

运行fastqc.sh脚本

1
bash fastqc.sh

每个fastq文件质控后生一个html和zip压缩格式文件,html可下载到本地用浏览器打开

2.multiqc

可以将当前各个fastqz文件的质控报告综合到一起,生成一个总的html命令multiqc接fastqc结果文件所在地址

1
multiqc ./

生成结果文件为multiqc_report.html,html可下载到本地用浏览器打开

3.去除低质量碱基

1.fastp

fastp去除低质量碱基的同时可以自动识别接头

单端测序:fastp -i 输入文件 -o 输出文件

双端测序:fastp - i reads1 -I reads2 -o reads1输出文件 -O reads2输出文件

1
2
3
4
5
ls *fastq.gz >fastq_list#生成fastq列表
ls *fastq.gz | xargs -i basename {} .1_1.fastq.gz#获得fastq.gz文件的前缀
ls *fastq.gz | xargs -i basename {} .1_1.fastq.gz | paste fastq_list - #生成第一列是fastq文件名,第二列是对应的前缀名
ls *fastq.gz | xargs -i basename {} .1_1.fastq.gz | paste fastq_list - | awk '{print "nohup fastp -i "$1" -o "$2 "_clean.fastq.gz &" }' >fastp.sh #在前一步的基础上用awk生成fastp脚本
bash fastp.sh#运行脚本

2.trim-galore

1
2
3
4
5
6
7
8
9
trim_galore [options] <filename>
--quality<int> #设定phred quality阈值。默认20(99%的read质量),如果测序深度较深,可以设定25
--phred33 #设定记分方式,代表Q+33=ASCII码的方式来记分方式。这是默认值。
--paired # 对于双端结果,一对reads中若一个read因为质量或其他原因被抛弃,则对应的另一个read也抛弃。
--output_dir #输出目录,需确保路径存在并可以访问
--length #设定长度阈值,小于此长度会被抛弃。这里测序长度是100我设定来75,感觉有点浪费
--strency #设定可以忍受的前后adapter重叠的碱基数,默认是1。不是很明白这个参数的意义
-e<ERROR rate> #设定默认质量控制数,默认是0.1,即ERROR rate大于10%的read 会被舍弃,如果添加来--paired参数则会舍弃一对reads
<filename> #如果是采用illumina双端测序的测序文件,应该同时输入两个文件。
1
2
3
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 \
--paired $dir/cmp/01raw_data/$fq1 $dir/cmp/01raw_data/$fq2 \
--gzip -o $input_data #双端测序

3.fastqc

再次用fastqc检测碱基质量,并与去除低质量碱基前的fastqc结果进行对比

客官打个赏咯.