1-跟着science学习宏基因组-联合数据质控
数据质控
本小节数据已更新：https://github.com/taowenmicro/Megagenome_learing.
准备conda虚拟环境
我们创建一个虚拟环境，这里python环境为2.7版本，因为大部分软件都是2.7版本编写的，目前更新比较频繁的软件用了3以上的版本。
conda env create -n meta python=2.7
conda activate science
# 切换文件夹路径
cd Megagenome_learing/
step1-合并和重命名（merge_and_rename）
这部分对序列进行重整理，没有什么实际上的功能，只是要注意这种大文件用多线程压缩和解压缩工具pigz。
单个样本运行
# 软件安装
sudo apt-get install pigz
mkdir unpack
pigz -p 6 -dc ./rawdata_0.01/SUBERR793599_2.fq.gz | pigz -p 6 > ./unpack/SUBERR793599_2.fq.gz
pigz -p 6 -dc ./rawdata_0.01/SUBERR793599_1.fq.gz | pigz -p 6 > ./unpack/SUBERR793599_1.fq.gz
d 首先解压缩 压缩文件
c 将输出文件写入stdotu
p 线程数量
批量运行
mkdir -p unpack
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
pigz -p 6 -dc ./rawdata_0.01/${base}_1.fq.gz | pigz -p 6 > ./unpack/${base}_1.fq.gz
pigz -p 6 -dc ./rawdata_0.01/${base}_2.fq.gz | pigz -p 6 > ./unpack/${base}_2.fq.gz
done
step2 使用skewer去除barcode等非生物序列（rule skewer）
质控软件skewer进行处理。skewer质控我并没有见到很多的中文介绍。当然用法什么的就没有了。目前去重barcode也有好几种方式选择，下面就可以看到了。
skewer：是良好的barcode去除软件，用于双端测序序列barcode的去除。
这里去除的barcode是science给出的。
# 安装软件
conda install skewer# 构造结果存储文件夹
mkdir skewer
skewer -x AGATGTGTATAAGAGACAG -m head -1 -t 6 --quiet ./unpack/SUBERR793599_1.fq.gz 2>> skewer.head.log|skewer -x CTGTCTCTTATACACATCT -m tail -t 6 --quiet -1 - 2>> skewer.tail.log > ./skewer/SUBERR793599_1.fq
skewer -x AGATGTGTATAAGAGACAG -m head -1 -t 6 --quiet ./unpack/SUBERR793599_2.fq.gz 2>> skewer.head.log|skewer -x CTGTCTCTTATACACATCT -m tail -t 6 --quiet -1 - 2>> skewer.tail.log > ./skewer/SUBERR793599_2.fq
x barcode序列
m 模式 ，有两种，基于单条序列可选输入：
head（序列开头，默认5端开头）；
tail（序列结尾，3端结尾）；
any（一条序列的任何地方）。
paired-end — pe模式：
省略。
t 线程数，这里使用6个线程
—quiet屏幕上不用输出信息
批量运行barcode去除
mkdir -p skewer
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
skewer -x AGATGTGTATAAGAGACAG -m head -1 -t 6 --quiet ./unpack/${base}_1.fq.gz 2>> skewer.head.log|skewer -x CTGTCTCTTATACACATCT -m tail -t 6 --quiet -1 - 2>> skewer.tail.log > ./skewer/${base}_1.fq
skewer -x AGATGTGTATAAGAGACAG -m head -1 -t 6 --quiet ./unpack/${base}_2.fq.gz 2>> skewer.head.log|skewer -x CTGTCTCTTATACACATCT -m tail -t 6 --quiet -1 - 2>> skewer.tail.log > ./skewer/${base}_2.fq
done
使用结果做sickle pe质控
step3 sickle 质控pe模式
生信菜鸟团在17年的时候有一篇关于sickle的推送大家可以看看；质控原理：Sickle使用滑窗，连同碱基质量和长度阈值来确定何时质量足够低以去除3’末端reads，及何时足够高以去除5’末端reads。它利用质量值滑窗，窗口长度为read长度的0.1倍。
例如read为150bp，则窗口为15bp。若这个长度<1，则窗口长度等于read长度。若长度>1，则沿着质量值滑窗，直到窗口中的平均质量超过阈值，这时，算法确定质量在窗口的哪个位置开始升高，并在此处cut read，作为5’末端。接着滑窗，直到窗口中平均质量低于阈值，这时，算法确定是窗口中哪个位置开始下降的，并在此处cut read，作为3’末端。若此时长度小于指定的长度阈值，则丢弃这条reads。
软件安装
# 安装-建议conda安装检查依赖
conda install sickle
# 或者
sudo apt install sickle
代码运行：Sickle 双末端的read (sickle pe)；它可以使用两个双末端的fastq文件作为输入，然后输出两个已经修饰后的双末端fastq文件和一个“单个”文件
mkdir trimming
sickle pe -f ./skewer/SUBERR793599_1.fq -r ./skewer/SUBERR793599_2.fq -o ./trimming/SUBERR793599_1.fq.gz -p ./trimming/SUBERR793599_2.fq.gz -s ./trimming/SUBERR793599_unpaired.fq.gz -l 150 -q 30 -t sanger
f 双端序列前端
r 双端序列后端
o 输出前端质控过文件
p 输出后端质控过文件
s输出单个文件；
“单个”文件包含已经过滤后正向或反向的read中的其中一个。
-t sanger:sanger指的就是phred33
结果存于triming文件夹中
批量运行 质控
mkdir trimming
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
sickle pe -f ./skewer/${base}_1.fq -r ./skewer/${base}_2.fq -o ./trimming/${base}_1.fq -p ./trimming/${base}_2.fq -s ./trimming/${base}_unpaired.fq -l 150 -q 30 -t sanger
done
step4 使用trimmomatic质控
trimmomatic 质控是目前较为主流和常见的质控方式；这一步还去除了barcode。这一步不需要安装软件，直接调用和下载的可执行文件即可，注意环境中必需有java。
想必大家开始有疑问了，前面先去除barcode然后在使用sickle有什么不同呢？为什么还要做这一步呢？其实这篇science的流程部分都是单独发表的，不仅仅是软件的胡乱堆砌，在我梳理了整个流程之后，发现作者有多套逻辑在里面，等整个流程更新完成后，我将解释这个问题。
conda install trimmomaticmkdir trimmomatic
trimmomatic PE -threads 6 -phred33 ./unpack/SUBERR793599_1.fq.gz ./unpack/SUBERR793599_2.fq.gz ./trimmomatic/SUBERR793599_forward_paired.fq.gz ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz ./trimmomatic/SUBERR793599_reverse_paired.fq.gz ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz ILLUMINACLIP:../database/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:30 MINLEN:100 2> SUBERR793599.log
#--如果不好安装，可以直接调用jar文件运行，也是一样的
java -jar ~/sra/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 6 -phred33 ./unpack/SUBERR793599_1.fastq.gz ./unpack/SUBERR793599_2.fastq.gz ./trimmomatic/SUBERR793599_forward_paired.fq.gz ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz ./trimmomatic/SUBERR793599_reverse_paired.fq.gz ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz ILLUMINACLIP:../database/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:30 MINLEN:100 2> SUBERR793599.log
使用我们百分之一过滤的数据，几乎秒出
批量质控
mkdir trimmomatic
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
trimmomatic PE -threads 6 -phred33 ./unpack/${base}_1.fq.gz ./unpack/${base}_2.fq.gz ./trimmomatic/${base}_forward_paired.fq.gz ./trimmomatic/${base}_forward_unpaired.fq.gz ./trimmomatic/${base}_reverse_paired.fq.gz ./trimmomatic/${base}_reverse_unpaired.fq.gz ILLUMINACLIP:../database/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:30 MINLEN:100 2> ${base}.log
done
合并双端未匹配的文件
小工具介绍：使用zcat 做到不解压也可以合并文件
cat是用于查看普通文件的。zcat 是用于查看压缩的文件，
将前后端匹配上的序列直接合并到SUBERR793599_unpaired_combined.fq.gz中。
# -单个样本操作
zcat ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz | gzip -c > ./trimmomatic/SUBERR793599_unpaired_combined.fq.gz
# 批量操作
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
zcat ./trimmomatic/${base}_forward_unpaired.fq.gz ./trimmomatic/${base}_reverse_unpaired.fq.gz | gzip -c > ./trimmomatic/${base}_unpaired_combined.fq.gz
done
大家也不要有疑问，为什么合起来，要干什么呢？我可以提前透露给大家，其实就是为了在后面进一步利用序列，这篇文章的作者很好的将所有序列进行了深度的利用，我个人认为是一个典范。
step5 使用readstats.py工具统计原始序列
对原始序列进行序列统计。安装工具khmer,这也是一个质控工具。这里我们并没有使用khmer质控，仅仅使用其中一个小工具做了序列统计。统计序列数量，长度，bp数量。令我没想到的是作者使用的是单条序列竟然有300bp的长度，这与现在150bp的不同。
软件安装：
conda install khmer
统计原始序列 ：这里readstats.py脚本的使用方法也是十分简单，无论要统计多少个fq文件的信息，直接提供路径就好了。— csv 输出为csv格式文件
readstats.py ./unpack/SUBERR793599_1.fq.gz ./unpack/SUBERR793599_2.fq.gz --csv -o raw.readstat.csv 2> ./unpack/raw.readstat.log
统计Sickle质控后的序列
输出trimmed.readstat.csv文件，同时将log文件也输出trimmed.readstat.log
mkdir ./stats
readstats.py ./trimming/SUBERR793599_1.fastq.gz ./trimmomatic/SUBERR793599_forward_paired.fq.gz ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz ./trimmomatic/SUBERR793599_reverse_paired.fq.gz ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz --csv -o ./stats/trimmed.readstat.csv 2> ./stats/trimmed.readstat.log
查看序列统计信息输出
# 批量操作
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
readstats.py ./unpack/${base}_1.fq.gz ./unpack/${base}_2.fq.gz --csv -o raw.readstat.csv 2> ./unpack/raw.readstat.log
cat raw.readstat.csv
readstats.py ./trimming/${base}_1.fastq.gz ./trimmomatic/${base}_forward_paired.fq.gz ./trimmomatic/${base}_forward_unpaired.fq.gz ./trimmomatic/${base}_reverse_paired.fq.gz ./trimmomatic/${base}_reverse_unpaired.fq.gz --csv -o ./stats/trimmed.readstat.csv 2> ./stats/trimmed.readstat.log
cat ./stats/trimmed.readstat.csv
done
这里作者将输出的信息都保存在log文件中。
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。
团队工作及其成果 （点击查看）
了解 交流 合作
团队成员邮箱 袁军：
junyuan@njau.edu.cn；
文涛：
2018203048@njau.edu.cn
团队公众号：