【RNA-Seq 实战】三、fastq下载及过滤数据
这里是佳奥！继续转录组分析的学习。
一般我们都是从.fq文件开始的，但是如果是自己寻找数据的话，要多一步骤。
1 获取fastq数据
1.1 sratoolkit
我们看文章的时候，找到对应的SRR号。
进入环境，下载sratoolkit软件，后续就可以直接在Linux内操作了。
最好的方式还是从官网下载（Github），把压缩包传到Linux，解压（conda下载的无法使用，版本过旧）
添加文件夹到环境变量即可
$ export PATH="$PATH:/home/kaoku/biosoft/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin"
配置软件：这个界面是可以鼠标点击的，设置两个路径到root/ncbi即可
$ vdb-config --interactive
安装成功后，使用prefetch便开始下载
prefetch SRR1039508
prefetch批量下载：SRR存在一个id.txt内
cat id | while read id; do prefetch $id; done
挂在后台下载：
cat id | while read id; do (prefetch $id &); done
使用top查看下载进程。
想关掉全部进程：
ps -ef | grep prefe | awk '{print $2}'| while read id; do kill $id; done
下载完成后，我们会看到一列的SRR1039508.sra文件，转化成fastq文件
fastq-dump --gzip --split-3 -o ./ /SRR1039508.sra
便开始生成SRR1039508_1.fastq.gz文件。
1.2 直接wget
这是一个研究对象为拟南芥的文章，所有的fastq数据存放于此， ID为E-MTAB-5130。
先获取.txt文件，再提取出URL，wget下载。
wget http://www.ebi.ac.uk/arrayexpress/files/E-MTAB-5130/E-MTAB-5130.sdrf.txt
head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI
tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}
随后下载参考基因组TAIR10 ver.24。
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gff3.gz &
nohup wget ftp://ftp.ensemblgenomes.org/pub/plants/release-28/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.28.gtf.gz &
后续就是按照流程一步一步了。
查看下载日志：
less raw/nohup.out
grep failed raw/nohup.out
2 fastqc查看测序质量
获取数据后，用fastqc查看数据质量。
fastqc SRR957678.fastq.gz
批量处理.fastqc.gz：一共10个文件
ls *gz | xargs fastqc -t 10
用multiqu批量查看.html结果
使用conda安装：
conda install -c bioconda multiqc
进入结果目录，直接运行：
multiqc ./
即可出现multiqc_report.html文件
下面我用数据跑一遍流程：
$ ls
SRR957677.fastq.gz SRR957678.fastq.gz SRR957679.fastq.gz SRR957680.fastq.gz
全部跑一遍fastqc
$ ls *gz | xargs fastqc -t 10
运行multqc
$ multiqc ./
/// MultiQC | v1.13.dev0
| multiqc | Search path : /home/kaoku/rnaseq/biotree_human
| searching | ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 14/14
| fastqc | Found 4 reports
| multiqc | Compressing plot data
| multiqc | Report : multiqc_report.html
| multiqc | Data : multiqc_data
| multiqc | MultiQC complete
$ ls 查看该multiqc_report.html
multiqc_data
multiqc_report.html
SRR957677_fastqc.html
SRR957677_fastqc.zip
SRR957677.fastq.gz
SRR957678_fastqc.html
SRR957678_fastqc.zip
SRR957678.fastq.gz
SRR957679_fastqc.html
SRR957679_fastqc.zip
SRR957679.fastq.gz
SRR957680_fastqc.html
SRR957680_fastqc.zip
SRR957680.fastq.gz
image.png
整体质量还不错。
查看了数据以后，我们就需要过滤测序数据。
3 过滤测序结果（双端测序结果）
如果Adapter Content是错误的，说明需要对reads进行修剪处理。
查看Adapter Content ：（fastqc软件测出的Adapter Content 序列后数个碱基是TCTTCTGCTTG）
zcat SRR957677_fastqc.zip | grep TCTTCTGCTTG
使用TrimGalore软件
在GitHub下载.tar.gz，解压
添加到环境变量（要进入rnaseq环境，依赖cutadapt）
$ export PATH="$PATH:/home/kaoku/biosoft/trim_galory/TrimGalore-0.6.7"
设置参数：（--length 30-40、--stringency 3-5）
dir='/home/kaoku/rnaseq/biotree_human/clean'（工作目录）
fq1='/home/kaoku/rnaseq/biotree_human/SRR957677_1_fastqc.zip'（举例）
fq2='/home/kaoku/rnaseq/biotree_human/SRR957677_2_fastqc.zip'（举例）
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2
可以继续调整参数（引自lncRNA组装流程的软件介绍之trim-galore）
--quality：设定Phred quality score阈值，默认为20。分析时可改成25，稍微严格一些。
--phred33：：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。具体怎么选择，看你用什么测序平台；例如：-phred33对应(Sanger/Illumina 1.9+ encoding)，-phred64对应(Illumina 1.5 encoding)
--adapter：输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
--stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length：设定输出reads长度阈值，小于设定值会被抛弃。
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
--output_dir：输入目录。需要提前建立目录，否则运行会报错。
--trim-n : 移除read一端的reads
批量处理的话：
ls *_1.fastqc.gz >1
ls *_2.fastqc.gz >2
paste 1 2 > config
dir='/home/kaoku/rnaseq/biotree_human/clean'
cat config | while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
echo $dir $fq1 $fq2
nohub trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
运行完成后也会有ERR1698194_2_val_2_fastqc.html报告，查看清洗前后差异。
拿到干净的测序数据后，就可以开始后面的比对分析流程了。
那么我们，下篇再见！
