宏基因组分析，首先须要的结果是得到物种分类信息，前面咱们提到，宏基因组有两种分析方式分别是基于序列比对和组装的，组装对电脑硬件的要求是超级高的，不过比对，就轻松多了，得益于算法的优化，有的软件能够实如今我的电脑上进行分析。最近的“AMD YES"也很给力，把普通我的电脑推向了8核16G RAM这种配置，终于追上了手机（虽然两者性能不可同日而语）。算力有了，充分利用起来呀，来个宏基因组玩玩呀！html
软件安装和环境准备
人生苦短，用conda吧！软件安装这事，能不折腾就不折腾。顺便提一句，我是使用win10自带的内置linux子系统WSL进行的，虽然还不完美（有人遭遇了权限问题），也算挺方便了。python
# http://blog.sciencenet.cn/blog-3334560-1115288.html
source ~/miniconda3/bin/activate
#conda create -n kraken2
conda activate kraken2
#conda install -y bracken kraken2
#conda install python=2.7 -y bowtie2 trimmomatic KneadData
下载数据库是一个难题，前面已经提到怎样使用腾讯云下载了，这里再也不重复。造福你们，我把数据库上传了百度云，听说如今宿主基因组，这里下载的g38。这个使用全平台多线程下载工具如motrix的下载速度挺快的，几M/s。直接下载软件准备好的数据库也是能够的。kraken2数据库 连接:https://pan.baidu.com/s/1sqI3Srw3YxRENpQbP80RKA 提取码:xn3alinux
# 下载人类基因组
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz
gunzip Homo_sapiens_assembly38.fasta.gz
#建索引
bowtie2-build --threads 8 Homo_sapiens_assembly38.fasta broad_hg38
#质控
kneaddata -i KM180115501/KM180115501.R1.fastq.gz -i KM180115501/KM180115501.R2.fastq.gz -o kneaddata_out -v \
-db /home/zd200572/Reference/broad_hg38 --trimmomatic /home/zd200572/miniconda3/envs/kraken2/share/trimmomatic-0.39-1/ --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
-t 8 --bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output
#汇总
kneaddata_read_count_table --input kneaddata_out --output kneaddata_read_counts.txt
# 清理宿主污染至指定目录备用
mkdir kneaddata_out/contam_seq
mv kneaddata_out/*_contam*fastq kneaddata_out/contam_seq
## 本质上只合并了1/2
concat_paired_end.pl -p 8 --no_R_match -o cat_reads kneaddata_out/*_paired_*.fastq
# 可选方法2：我感受是合并4个文件，应该重命名序列ID，不则双端序列重名字，shell命令合并单样品
for i in `tail -n+2 map.txt | cut -f 1`;do \
cat kneaddata_out/${i}_R1_kneaddata_paired_1.fastq kneaddata_out/${i}_R1_kneaddata_paired_2.fastq knead
#分类
kraken2 \
--db ~/Biosofts/minikraken2_v1_8GB \
--threads 8 \
--report report \
--paired \
kneaddata_out/KM180115501.R1_kneaddata_paired_1.fastq kneaddata_out/KM180115501.R1_kneaddata_paired_2.fastq
bracken -d ~/Biosofts/minikraken2_v1_8GB -t 8 -i report -o result.out
159s的运行时间事后，你就获得了你的样本物种分类结果。ios
参考了几篇博客的内容，若有错误，欢迎指正。web
本文分享自微信公众号 - 科技记者（kejijizhe）。
若有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一块儿分享。算法