宏基因組分析，首先須要的結果是得到物種分類信息，前面咱們提到，宏基因組有兩種分析方式分別是基於序列比對和組裝的，組裝對電腦硬件的要求是超級高的，不過比對，就輕鬆多了，得益於算法的優化，有的軟件能夠實如今我的電腦上進行分析。最近的「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源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。算法