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

宏基因组分析,首先需要的结果是获得物种分类信息,前面我们提到,宏基因组有两种分析方式分别是基于序列比对和组装的,组装对电脑硬件的要求是超级高的,不过比对,就轻松多了,得益于算法的优化,有的软件可以实现在个人电脑上进行分析。最近的“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的运行时间过后,你就得到了你的样本物种分类结果。

参考了几篇博客的内容。

发表评论