DAY5学习内容：过滤和比对
过滤
通过DAY4_fastqc结果，我们开始对数据进行过滤，主要是去除接头和低质量碱基，一般使用两种软件可以完成这个过程：trim_galore和trimmoatic，而这两种软件遵循的原则是直接从不合格的位置向后全部切掉，而不是只挖掉这块，为了避免人为改变测序序列。
两款软件的学习资料：
trim_galore：https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
trimmomatic：https://github.com/mossmatters/KewHybSeqWorkshop/blob/master/FastQC_Trimmomatic.md
trim_galore
trim_galore在去除接头序列时依赖的是Cutadapt，一般接头出现在3'末端；
选项参数：
-q：设置碱基质量阈值
--phredxx：质量体系；phred33和phred64，不过现在绝大部分 Illumina 平台的产出数据也都转为使用 phred33 格式了。
如何区分phred33和phred64：fa第四行的碱基质量值如果大部分都是大写字母那就是phred64，如果存在特殊符号类似#、%、@就是phred33；或者看fastqc报告中basic statistics结果：Illumina 1.8版本以上的就是phred33
--length：表示测序片段长度阈值，例如--length 50，表示如果测序长度为150bp，如果左右切掉接头和低质量碱基后，长度低于了50bp，就直接把这整条去掉
--stringency：设置的数字表示：认为有几个碱基和接头序列是有重叠的，默认值为1，意思是从3'的对应位置检测，出现一个碱基与常用接头有重叠就把这个碱基以后的序列都去掉。
常用的接头前12-13bp序列
Illumina: AGATCGGAAGAGC
Small RNA: TGGAATTCTCGG
Nextera: CTGTCTCTTATA
--paired：双末端测序
-o：输出路径
以上参数及其他参数执行trim_galore -help均可查看
trim_galore -help
通过帮助文件了解软件使用后，我们开始用起它，两种情况：
-第一种情况：只对一组fq文件进行过滤
$dir=/YOUR_PATH/ #定义输出文件路径
$fq1=/YOUR_PATH/xxx.fq1 #将xxx.fq1文件路径及本身赋值给fq1
$fq2=/YOUR_PATH/xxx.fq2 #将xxx.fq2文件路径及本身赋值给fq2
trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 #-e 允许的最大错误率（错误数除以匹配区域的长度）（默认值：0.1）
-第二种情况：批量对多组fq文件进行过滤
首先要将所有的fq文件放在同一个文件中，基于前期学习我们这共有8组fq文件，如下：
如何快速的完成这个过程呢？
#进入fq的文件目录
cd raw
#找到1端测序文件，并保存到另一个文件中，例如这里叫做fastq_1
ls *_1* >fastq_1
#找到2端测序文件，并保存到另一个文件中，例如这里叫做fastq_2
ls *_2* >fastq_2
#接下来我们将这两个新生成的文件合并起来，使用paste命令
paste fastq_1 fastq_2
#屏幕显示所有fq文件
为了保证每个文件在任何时候和目录下都可以重复使用，我们需要添加其绝对路径，因此我们要对上面的工作再加上其路径，并将其合并结果放到一个新的文件中，命令如下：
ls `pwd`/*_1* >new_fastq_1
ls `pwd`/*_2* >new_fastq_2
paste new_fastq_1 new_fastq_2 >conf_fastq
#conf意思是配置文件的意思，然后#cat conf_fastq，显示加上了路径
对上面的conf_fastq进行trim_galore操作，同样使用循环命令，如下：
#将脚本写入script目录中，
cd script
#新建一个名字为filter.sh脚本，并添加内容
cat >filter.sh #或
vi filter.sh
#脚本内容
filter_data=~/raw
cat $filter_data/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
echo $i
# nohup trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o /hwfssz1/ST_AGRIC/F13ZHQYJSY1409_LP/USER/liujia/TCM/TD/clean $fq1 $fq2 &
done
还是先将主要执行的命令注释掉，看看循环可行不，如下：
没问题！释放执行脚本！
filter_data=~/raw
cat $filter_data/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
# echo $i
nohup trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o /hwfssz1/ST_AGRIC/F13ZHQYJSY1409_LP/USER/liujia/TCM/TD/clean $fq1 $fq2 &
done
在上述脚本中，一个新学的只是就是创建数组，进行数据提取
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
#使用（）创建数组，将每一个$i（一对PE测序文件）放入其中
#然后使用${[]}进行访问，[]表示访问第几列，从0开始计数
#每当读入一对fq时，${fqs[0]}表示将第一个测序文件也就是_1文件赋值给fq1，同样${fqs[1]}表示将第二个测序文件也就是_2文件赋值给fq2，此过程必须要带绝对路径。
脚本编辑完，我们就可以运行脚本了，bash filter.sh 或 sh filter.sh
sh和bash区别：
sh、bash都是linux中的解释器，但有的sh是没有数组array解释器的，因此当你使用sh去运行整个命令时，有可能会产生报错：
Syntax error: "(" unexpected (expecting "done")
但这并不是说代码写错了，而是选择错了解释器。
trimmomatic：
该软件有两种过滤模式，分别对应SE和PE测序数据，同时支持gzip和bzip2压缩文件。Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关
过滤步骤：
- ILLUMINACLIP；
- SLIDINGWINDOW；
- MAXINFO；
- LEADING；
- TRAILING；
- CROP；
- HEADCROP；
- MINLEN；
- AVGQUAL；
- TOPHRED33；
- TOPHRED64
参数选择：
对于双末端测序模式：
PE模式
输入输出文件：
PE 模式的两个输入文件：sample_R1.fastq sample_R2.fastq以及四个输出文件：sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq
SLIDINGWINDOW：滑窗剪切，统计滑窗口中所有碱基的平均质量值，如果低于设定的阈值，则切掉窗口。
SLIDINGWINDOW:<windowSize>:<requiredQuality>
widowSize:设置窗口大小
requiredQuality：设置窗口碱基平均质量阈值
SLIDINGWINDOW:5:20 设置 5bp 窗口，碱基平均质量值阈值 20
LEADING：从 reads 的起始端开始切除质量值低于设定的阈值的碱基，直到有一个碱基其质量值达到阈值。
LEADING:<quality>
quality:设定碱基质量值阈值，低于这个阈值将被切除。
LEADING:5
TRAILING：从 reads 的末端开始切除质量值低于设定阈值的碱基，直到有一个碱基质量值达到阈值。Illumina 平台有些低质量的碱基质量值被标记为 2 ，所以设置为 3 可以过滤掉这部分低质量碱基。TRAILING:<quality>
quality:设定碱基质量值阈值，低于这个阈值将被切除。
TRAILING:5
MINLEN：设定一个最短 read 长度，当 reads 经过前面的过滤之后，如果保留下来的长度低于这个阈值，整条 reads 被丢弃。被丢弃的 reads 数会被统计在 Trimmomatic 日志的 dropped reads 中。
MINLEN:<length>
length：可被保留的最短 read 长度
MINLEN:20
-trimlog：输出Trimmomatic 日志
TOPHRED33/64：
此选项可以将过滤之后的 Fastq 文件中质量值这一行转为 phred-33 格式 或 phred-64 格式或。
更多关于Trimmomatic的学习资料请看考NGS 数据过滤之 Trimmomatic 详细说明或其--help
了解Trimmomatic软件使用后，我们开始用起它，脚本如下，参数如上：
mkdir $rna/clean/trimmomatic && cd "$_"
cat >filter-trim.sh
rna=/YOUR_PATH/rnaseq
cat $rna/raw/conf_fastq | while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
tri1=`basename $fq1`
tri2=`basename $fq2`
nohup trimmomatic PE -phred33 \
-trimlog trim.logfile\
$fq1 $fq2 \
clean.$tri1 unpaired.$tri1 \
clean.$tri2 unpaired.$tri2 \
SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:20 &
done
bash filter-trim.sh
任务已经提交，那我们怎么看我们这个任务跑的是否成功呢
利用htop工具（进程管理器）功能和top一样，只是它是彩色的，信息详细
安装学习链接：
conda安装看这里：https://anaconda.org/conda-forge/htop
官网安装看这里：https://hisham.hm/htop/
安装后，直接输入htop回车就可以查看当前进程状态：CPU、内存、正在跑的进程
如果你只想看你自己的进程状态：htop -u your_ID
过滤好的数据我们还要进行质控，看看接头和低质量的碱基是否被处理掉，主要看接头！
比对
基于之前工作得到的cleandata，开始进行比对过程
比对所需的文件：
-过滤并质控合格后的数据
-参考基因组和注释文件
如何生成小数据用于测试流程
生成小数据，说白了就是从你的数据里提取出一部分，用于测试，这里我们使用seqtk软件，https://github.com/lh3/seqtk，用conda安装即可，
首先我们先整理一下之前生成的cleandata，方法如上：
# 先配置cleandata数据集文件
cd clean
ls `pwd`/*_1*gz >fastqc_1.clean
ls `pwd`/*_2*gz >fastqc_2.clean
paste fastqc_1.clean fastqc_2.clean >fastq_clean.conf
cat fastq_clean.conf
在项目总的目录下新建一个test目录
#创建test目录并进入
mkdir test && cd $_
#vi编译名字为sample_fq.sh脚本
vi sample_fq.sh
#脚本内容
clean_data=~/clean
cat $clean_data/fastq_clean.conf | while read i;
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
file=`basename $fq1`
#echo $file
surname=${file%%_*}
#echo $fq1 $fq2 $surname
# 随机选了10000条reads
seqtk sample $fq1 10000 >test.${surname}_1.fastq
seqtk sample $fq2 10000 >test.${surname}_2.fastq
done
bash sample_fq.sh
以上脚本可以通过echo $file，可得知file做了什么，basename就是数据这个路径下最底层的名字，输出如下：
echo $surname又做了什么？输出如下
$fq1：~/clean/SRR1039508_1_val_1. fq.gz
$fq2：~/clean/SRR1039508_2_val_2. fq.gz
$surname：SRR1039508
file%%_*：%%就是将file中_符号后面的全部删除，也就是把_2_val_2. fq.gz删除，只留下SRR1039508
生成小数据结果如下：
比对第一课：hisat2比对
比对又叫mapping，mapping的目的就是找到reads在基因组上的位置。
hisat2官网：http://ccb.jhu.edu/software/hisat2/index.shtml
选项参数：http://ccb.jhu.edu/software/hisat2/manual.shtml
文章推荐：NC的文章：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
更多背景知识请自行查找，开始对上一步生成的小数据进行mapping
mapping第一步：构建基因组索引
构建索引的目的：快速在基因组上进行定位，因为reads一般都是几千万条，基因组也有30亿碱基，肯定不能一条reads在基因组遍历一遍，再进行下一条reads。
mkdir ~/align/hisat2 && cd "$_"
vi index.sh #创建index.sh脚本
#脚本内容
rna=/YOUR_PATH/rnaseq
genome=$rna/ref/hg19.fa #DAY4下载的参考基因组
hisat2-build -p 10 $genome hg19 #建立索引
mapping第二步：开始比对
建立好索引后，我们开始进行比对（其实第二步在mapping之后还包含sam转bam和排序），脚本如下，改脚本在~/align/hisat2路径下执行：
# 小数据集合的路径
ls ~/test/*_1* >hisat2_1.test
ls ~/test/*_2* >hisat2_2.test
paste 1.test 2.test >test_hisat2.conf
vi hisat2_aln.sh
cat ~/test/test_hisat2.conf
| while read i
do
fqs=($i)
fq1=${fqs[0]}
fq2=${fqs[1]}
sample=`basename $fq1 | cut -d '_' -f1`
hisat2 -p 10 -x hg19 -1 $fq1 -2 $fq2 | samtools sort -T PREFIX.nnnn.bam -O bam \
-@ 10 -o ${sample}_hisat.sorted.bam
samtools index -b ${sample}_hisat.sorted.bam
done
参数信息：
-p：进程数，越打越快
-x：参考基因组索引文件的前缀
-1：双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2：双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
samtools sort：使用samtools工具对生成的sam文件进行排序，并转换成二进制bam文件
-T：添加一个临时文件是必须的，格式：PREFIX.nnnn.bam
-O：设置输出格式包括：bam，sam或cram
-@：除主线程外还要使用的BAM压缩线程数 10
samtools index：对bam进行排好队后，需要给大家一个编号，就是index
-b：建立索引后将产生后缀为.bai的文件。
以上脚本工作流程如下：hisat2将reads比对到了参考基因组，理论上会生成sam文件，但是这个文件太大，因此通过samtools sort -O bam直接转换成二进制的bam文件并排序，与sam内容一样而且节省空间。最后为了下一步导入进行可视化并且为了定量的需要，我们又对bam构建了索引（samtools index）。
对于samtools安装，conda install samtools=1.8，一定要设置1.8安装版本，不然会报错
关于更多samtools请学习：处理SAM、BAM你需要Samtools；英文资料：http://www.htslib.org/doc/
比对后生成文件，mapping率在nohup.out下看（less nohup.out）如下：
