个人电脑也做做宏基因组玩玩

宏基因组分析，首先需要的结果是获得物种分类信息，前面我们提到，宏基因组有两种分析方式分别是基于序列比对和组装的，组装对电脑硬件的要求是超级高的，不过比对，就轻松多了，得益于算法的优化，有的软件可以实现在个人电脑上进行分析。最近的“AMD YES”也很给力，把普通个人电脑推向了8核16G RAM这种配置，终于追上了手机（虽然二者性能不可同日而语）。算力有了，充分利用起来呀，来个宏基因组玩玩呀！

软件安装和环境准备

人生苦短，用conda吧！软件安装这事，能不折腾就不折腾。顺便提一句，我是使用win10自带的内置linux子系统WSL进行的，虽然还不完美（有人遭遇了权限问题），也算挺方便了。

# 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 提取码:xn3a

# 下载人类基因组
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/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的运行时间过后，你就得到了你的样本物种分类结果。

参考了几篇博客的内容。