这里是佳奥！继续转录组分析的学习。

一般我们都是从.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报告，查看清洗前后差异。

拿到干净的测序数据后，就可以开始后面的比对分析流程了。

那么我们，下篇再见！