2-跟着science学习宏基因组-去除宿主-评估测序质量是否足够
去除宿主徐序列bowtie2
本小节数据已更新：https://github.com/taowenmicro/Megagenome_learing.
bowtie2算是目前去除宿主的主流脚本之一了，使用的非常多，发展到目前，被许多组合型宏基因组软件所依赖。值得学习。许多组合型宏基因组软件所依赖。值得学习。我们一定要放开视野，如果测序样本中可能有多个宿主，那就下载多个宿主都建索引，进行去除操作。如果不确定是哪种物种，就找亲缘比较近的物种进行去除。多找几个都没关系的。这个道理也是这份science教程教给我们的。
构造宿主基因组索引
首先下载宿主基因组
网站
菜豆基因组下载位置：这里我选择的toplevel。做转录组的就比较清楚了，toplevel使用的最多了。这是全基因组。菜豆基因组 这个十分难下载；但是挂上vpn之后就简单了。构建索引解压基因组文件，然后构建索引；占用内存不大，但是消耗时间mkdir index
cd index
# 将下载的序列放到index文件夹中
#构建索引
bowtie2-build Phaseolus_vulgaris.PhaVulg1_0.dna_rm.toplevel.fa GCA_000499845.1_PhaVulg1_0_genomic
去除宿主序列
开始去除宿主
对trimmomatic的全部输出序列进行宿主去除。首先这个步骤作者使用的是trimmomatic后得到的匹配的两个序列文件去宿主。
首先是对trimmomatic后得到的匹配的两个序列文件进行比对。输出：SUBERR793599_fr_mapped_and_unmapped.sam和SUBERR793599_fr_unmapped_pairs.1和SUBERR793599_fr_unmapped_pairs.2bowtie2参数详解mkdir log
mkdir host_filtering
bowtie2 --very-sensitive -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -1 ./trimmomatic/SUBERR793599_forward_paired.fq.gz -2 ./trimmomatic/SUBERR793599_reverse_paired.fq.gz -S ./host_filtering/SUBERR793599_fr_mapped_and_unmapped.sam --un-conc ./host_filtering/SUBERR793599_fr_unmapped_pairs 2>> ./log/host_removal_SUBERR793599.log
参数解释：
1 前端质控过序列
2 后端质控后序列
p:线程数目
x:参考基因组索引
S:得到匹配上的序列文件，输出sam格式文件
—un-conc：
将不能比对的paired-end reads写入两个文件
下面继续对trimmomatic得到未匹配的双端序列 去宿主，trimmomatic得到未匹配的双端序列这里也需要做一个比对；对trimmomatic后未能匹配上的序列做比对。对未匹配上的序列做去除宿主操作：bowtie2 --very-sensitive -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz -U ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz -S ./host_filtering/SUBERR793599_s_mapped_and_unmapped.sam --un ./host_filtering/SUBERR793599_s_unmapped 2>> ./log/host_removal_SUBERR793599.log
输出文件：C1_s_mapped_and_unmapped.sam为比对上的文件和C1_s_unmapped是未必对上的文件，注意，合并了，只有一个。
补充参数解释：
U:逗号分隔的包含未匹配的双端序列文件读取参数；
这里读取正反两个文件。
—un：
未匹配的输出到同一个文件夹中。
过滤宿主序列# 创建一个文件夹用于过滤序列
mkdir ./host_filtering/raw/
#对前端序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/raw/SUBERR793599_R1_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/SUBERR793599_forward_paired.fq.gz -S ./host_filtering/SUBERR793599_fw_mapped_and_unmapped.sam 2>> ./log/host_filtering_SUBERR793599.log
#对后端序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/raw/SUBERR793599_R2_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/SUBERR793599_reverse_paired.fq.gz -S ./host_filtering/SUBERR793599_rev_mapped_and_unmapped.sam 2>> ./log/host_filtering_SUBERR793599.log
# 对未匹配上的序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/SUBERR793599_unpaired_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/SUBERR793599_forward_unpaired.fq.gz -U ./trimmomatic/SUBERR793599_reverse_unpaired.fq.gz -S ./host_filtering/SUBERR793599_fr_unpaired_mapped_and_unmapped.sam 2>> ./log/host_filtering_${base}.log
使用和宿主序列比对好的SUBERR793599_fw_mapped_and_unmapped.sam文件对trimmomatic/SUBERR793599_forward_paired.fq.gz和./trimmomatic/SUBERR793599_reverse_paired.fq.gz还有未匹配的双端测序文件进行过滤；（文件进行过滤。分别对前端，后端还有未匹配的双端序列进行操作）
输出三个文件，前两个是pair的两个过滤文件，存到raw文件夹中，第三个 ./host_filtering/SUBERR793599_unpaired_filtered.fastq文件直接存到host_filtering中。
过滤后文件对齐：为防止本来pair的序列没有对齐，这里进行操作，将其对齐
本工具将去除宿主后的序列进行对齐并输出三个文件：前两个分别为匹配上的双端文件，第三个为合并后的未匹配序列。这个工具是基于py2的 首先对这些工具进行下载：https://github.com/enormandeau/Scripts存放到一个位置并指定好，这个命令似乎是不管输出命令的，直接做了输入文件夹输出 这里我手动修改了输出文件名称：mkdir reference_filtered
python2 ./Scripts-master/fastqCombinePairedEnd.py ./host_filtering/raw/SUBERR793599_R1_filtered.fastq ./host_filtering/raw/SUBERR793599_R2_filtered.fastq ./host_filtering/SUBERR793599_R1_paired_filtered.fastq ./host_filtering/SUBERR793599_R2_paired_filtered.fastq ./reference_filtered/SUBERR793599_R1R2_singular_filtered.fastq
以上流程就做好了去除宿主的全部操作。得到两个文件：SUBERR793599_R1_paired_filtered.fastq 和SUBERR793599_R2_paired_filtered.fastq；
批量运行# ${base}
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
#-前端序列宿主比对
bowtie2 --very-sensitive -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -1 ./trimmomatic/${base}_forward_paired.fq.gz -2 ./trimmomatic/${base}_reverse_paired.fq.gz -S ./host_filtering/${base}_fr_mapped_and_unmapped.sam --un-conc ./host_filtering/${base}_fr_unmapped_pairs 2>> ./log/host_removal_${BASE}.log
#后端序列宿主比对
bowtie2 --very-sensitive -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/${base}_forward_unpaired.fq.gz -U ./trimmomatic/${base}_reverse_unpaired.fq.gz -S ./host_filtering/${base}_s_mapped_and_unmapped.sam --un ./host_filtering/${base}_s_unmapped 2>> ./log/host_removal_C1.log
# 创建一个文件夹用于过滤序列
mkdir ./host_filtering/raw/
#对前端序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/raw/${base}_R1_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/${base}_forward_paired.fq.gz -S ./host_filtering/${base}_fw_mapped_and_unmapped.sam 2>> ./log/host_filtering_${base}.log
#对后端序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/raw/${base}_R2_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/${base}_reverse_paired.fq.gz -S ./host_filtering/${base}_rev_mapped_and_unmapped.sam 2>> ./log/host_filtering_${base}.log
# 对未匹配上的序列去宿主
bowtie2 --very-sensitive --un ./host_filtering/${BASE}_unpaired_filtered.fastq -p 6 -x ./index/GCA_000499845.1_PhaVulg1_0_genomic -U ./trimmomatic/${base}_forward_unpaired.fq.gz -U ./trimmomatic/${base}_reverse_unpaired.fq.gz -S ./host_filtering/${base}_fr_unpaired_mapped_and_unmapped.sam 2>> ./log/host_filtering_${base}.log
# 序列对齐
python2 ./Scripts-master/fastqCombinePairedEnd.py ./host_filtering/raw/${base}_R1_filtered.fastq ./host_filtering/raw/${base}_R2_filtered.fastq ./host_filtering/${base}_R1_paired_filtered.fastq ./host_filtering/${base}_R2_paired_filtered.fastq ./reference_filtered/${base}_R1R2_singular_filtered.fastq
done
step7 质控后数据后处理
sickle质控后为匹配的双端序列拆分
将双端文件拆分开来
sickle 过滤后得到的文件位于trimming中，输出过滤后的对齐文件和双端在一起的未对齐文件。下面将未对齐的文件拆分出来两个。
切换到qiime1运行：这里的时候出现了问题，我从NCBI下载的数据除了文件名之外名没有在内容中区分前后端序列，所以这里报错了，我修改前端标识和后端标识后才重新成功；conda activate qiime1
extract_reads_from_interleaved_file.py -i ./trimming/SUBERR793599_unpaired.fq -o ./trimming/SUBERR793599_unpaired/
# 批量运行
# ${base}
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
extract_reads_from_interleaved_file.py -i ./trimming/${base}_unpaired.fq -o ./trimming/${base}_unpaired/
done
conda activate science
序列统计trimming文件夹中文件unpaired序列数量
这里统计前端unpaired序列数量 printf SUBERR793599'\t' >> ./trimming/unpaired_forward.count.txt && cat ./trimming/SUBERR793599_unpaired/forward_reads.fastq | printf $((`wc -l`/4)) >> ./trimming/unpaired_forward.count.txt && printf '\tforward\n' >> ./trimming/unpaired_forward.count.txt
# ${base}
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
printf ${base}'\t' >> ./trimming/unpaired_forward.count.txt && cat ./trimming/${base}_unpaired/forward_reads.fastq | printf $((`wc -l`/4)) >> ./trimming/unpaired_forward.count.txt && printf '\tforward\n' >> ./trimming/unpaired_forward.count.txt
done
统计后端unpaired序列数量printf SUBERR793599'\t' >> ./trimming/unpaired_reverse.count.txt && cat ./trimming/SUBERR793599_unpaired/reverse_reads.fastq | printf $((`wc -l`/4)) >> ./trimming/unpaired_reverse.count.tx && printf '\treverse\n' >> ./trimming/unpaired_reverse.count.txt
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
printf ${base}'\t' >> ./trimming/unpaired_reverse.count.txt && cat ./trimming/${base}_unpaired/reverse_reads.fastq | printf $((`wc -l`/4)) >> ./trimming/unpaired_reverse.count.tx && printf '\treverse\n' >> ./trimming/unpaired_reverse.count.tx
done
按照处理合并质控好的序列文件
按照处理合并序列文件
作者在这里可选两种方式，一种使用去除宿主的序列合并，另一个使用trimm后的序列合并。这里我选择的是去除宿主后的序列进行合并，一共合并三次，分别对对齐的两个文件和未对齐的文件进行合并。
思考：作者提供了这两种可选方式：如果不做宿主去除，使用sickle质控后的结果进行后续分析，如果进行宿主去除操作，则使用trimmati···软件的质控结果，这里可以理解为，如果宏基因组未包含宿主基因，使用后呢sickle质控会更好一些，如果包含，使用trim···质控然后去除宿主更好一些？#按照处理合并序列，建立一个文件夹
mkdir ./treatment
#合并同一个处理所有前端序列
cat ./host_filtering/raw/SUBERR793599_R1_filtered.fastq_pairs_R1.fastq ./host_filtering/raw/SUBERR793600_R1_filtered.fastq_pairs_R1.fastq ./host_filtering/raw/SUBERR793601_R1_filtered.fastq_pairs_R1.fastq > ./treatment/C_forward.fastq
# 合并同一个处理所有后端序列
cat ./host_filtering/raw/SUBERR793599_R2_filtered.fastq_pairs_R2.fastq ./host_filtering/raw/SUBERR793600_R2_filtered.fastq_pairs_R2.fastq ./host_filtering/raw/SUBERR793601_R2_filtered.fastq_pairs_R2.fastq > ./treatment/C_reverse.fastq
#合并双端未匹配上的序列
cat ./host_filtering/raw/SUBERR793599_R1_filtered.fastq_singles.fastq ./host_filtering/raw/SUBERR793600_R1_filtered.fastq_singles.fastq ./host_filtering/raw/SUBERR793601_R1_filtered.fastq_singles.fastq > ./treatment/C_unpaired.fastq
# 下面这个注释的代码错误了，你们能看出来哪里错误了吗？是空格不对
# cat ./host_filtering/raw/SUBERR793599_R1_filtered.fastq_singles.fastq ./host_filtering/raw/SUBERR793600_R1_filtered.fastq_singles.fastq ./host_filtering/raw/SUBERR793601_R1_filtered.fastq_singles.fastq > ./treatment/C_unpaired.fastq
使用nonpareil工具评估测序深度
该工具用来评估宏基因组序列覆盖度和多样性；通过读段冗余来估计平均覆盖率，并预测实现“几乎完全覆盖率”（定义为≥95％或≥99％的平均覆盖率）所需的序列数量。详情nonpareil本来想看看汉字的，却没有发现中文教程或者介绍。
安装：地址:https://nonpareil.readthedocs.io/en/latest/installation.html conda install nonpareil
nonpareil工具使用
一下参数均为science推荐，尚未完全习得参数内含，故不做修改。nonpareil -b ./nonpareil/nonpareil_forward -s ./treatment/nonpareil_forward.fastq -f fastq -t 6 -R 400000 -L 50 -T kmer
在R中运行可视化
注意在R语言中，脚本需要自己下载，我去github上搜的。source("./script/Nonpareil.R")
png("./nonpareil/nonpareil_forward.png")
np <- Nonpareil.curve("./nonpareil/nonpareil_forward.npo")
legend('bottomright', legend = c(paste("Coverage: ", round(np$C*100,digits=2), "%"),paste("Actual effort =", round(np$LR/1000000,digits=2), "Mbp"), paste("Required effort for 95% coverage=", round(np$LRstar/1000000,digits=2), "Mbp"), paste("Diversity =", round(np$diversity,digits=2))),bty = "n")
dev.off()
可以看到需要测140~G才能达到覆盖≥95％。
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。
团队工作及其成果 （点击查看）
团队关注
团队文章成果
团队成果-EasyStat专题
ggClusterNet专题
袁老师小小组
了解 交流 合作
团队成员邮箱 袁军：
junyuan@njau.edu.cn；
文涛：
2018203048@njau.edu.cn
团队公众号：