宏基因组数据组装及注释
测序文件-组装（megahit，fq1 fq2 - A1.contigs.fa）- 基因预测（prokka，A1.contigs.fa - faa ffn）-
去冗余（cd-hit； faa ffn - all.cds.cdhit/unigene序列文件 ）
组装
## 单组数据组装
singularity exec ../../software/MetaGenome.sif megahit \
-1 ../../dataFQ/A1_1.fq.gz \ # 输入，fq1
-2 ../../dataFQ/A1_2.fq.gz \ # 输入，fq2
--min-contig-len 1000 \ # contig最小长度
--tmp-dir ./ \ # 设置tmp目录
--memory 0.4 \ # 内存占用
--num-cpu-threads 4 \ # 线程数
--out-dir A1_megahit \ # 输出目录
--out-prefix A1 # 输出前缀
## 多组数据组装, 输入数据逗号分隔
singularity exec ../../software/MetaGenome.sif megahit \
-1 ../../dataFQ/A1_1.fq.gz,../../dataFQ/A2_1.fq.gz \
-2 ../../dataFQ/A1_2.fq.gz,../../dataFQ/A2_2.fq.gz \
--min-contig-len 1000 \
--tmp-dir ./ \
--memory 0.4 \
--num-cpu-threads 4 \
--out-dir multi_megahit \
--out-prefix multi
输出结果文件: A1_megahit/A1.contigs.fa
使用 metaspades 进行组装：
## 单组数据组装
singularity exec ../../software/MetaGenome.sif spades.py \
--meta \ # 宏基因组模式
-t 4 \ # 线程
-k 21,33 \ # kmer
-1 ../../dataFQ/A1_1.fq.gz \ # 输入，fq1
-2 ../../dataFQ/A1_2.fq.gz \ # 输入，fq2
-o A1_metaspades # 输出目录
## 多组数据组装
singularity exec ../../software/MetaGenome.sif spades.py --meta \
-t 4 \
-k 21,33,55,77 \
--pe-1 1 ../../dataFQ/A1_1.fq.gz \ #输入，第1组fq1
--pe-2 1 ../../dataFQ/A1_2.fq.gz \ #输入，第1组fq2
--pe-1 2 ../../dataFQ/A2_1.fq.gz \ #输入，第2组fq1
--pe-2 2 ../../dataFQ/A2_2.fq.gz \ #输入，第2组fq2
-o multi_metaspades #输出目录
输出结果文件：A1_metaspades/scaffolds.fasta
组装结果可以使用 quast 进行汇总统计：
# 准备组装结果
ln -s ../01.megahit/A1_megahit/A1.contigs.fa
ln -s ../01.metaspades/A1_metaspades/scaffolds.fasta ./A1.spades.fa
# 使用quast进行汇总统计
singularity exec ../../software/MetaGenome.sif quast.py ./*.fa
输出结果为网页版报告：report.html
基因预测及去冗余
# 准备输入文件
ln -s ../01.megahit/A1_megahit/A1.contigs.fa
# 基因预测
singularity exec ../../software/MetaGenome.sif prokka \
--outdir A1_prokka \ # 输出目录
--prefix A1 \ # 输出文件前缀
--addgenes \ # 添加gene信息
--addmrna \ # 添加mrna信息
--locustag A1 \ # 基因ID添加tag
--kingdom Bacteria \ #物种分类Archaea|Bacteria|Mitochondria|Viruses
--metagenome \ # 宏基因组模式
--cpus 8 \ # 线程数
./A1.contigs.fa # 输入文件
主要结果：
A1.faa : 预测的氨基酸序列文件
A1.ffn ：预测出的转录本序列
A1.fna ：输入的基因组文件
A1.fsa ：添加Sequin tags的基因组文件
A1.gbf ：Genbank格式注释结果
A1.gff ：gff3格式注释结果
A1.sqn ：ASN1 格式 Sequin文件，修改后可用于提交Genbank
A1.tbl : Feature Table, 创建sequin文件使用
A1.tsv : features: locus_tag,ftype,len_bp,gene,EC_number,COG,product
A1.txt : 注释结果统计信息
对多个样品基因预测结果去冗余，得到 unigene:
## 将多样品转录本和氨基酸序列分别合并
cat ../02.prokka/*_prokka/*.ffn > all.trans.fa
cat ../02.prokka/*_prokka/*.faa > all.pep.fa
## 提取编码序列ID
sed -n "s/ˆ>\(\S\+\).*$/\1/p" all.pep.fa > all.cds.id
## 基于ID提取CDS序列
singularity exec ../../software/MetaGenome.sif seqtk subseq \
all.trans.fa all.cds.id > all.cds.fa
## 对CDS进行聚类
singularity exec ../../software/MetaGenome.sif cd-hit-est \
-i all.cds.fa \
-o all.cds.cdhit \
-c 0.95 \
-aS 0.9 \
-M 10000 \
-T 4
输出结果文件：
all.cds.cdhit ： unigene 序列文件 all.cds.cdhit.clstr：聚类cluster信息
由于 unigene 序列中使用原始的序列名称，这里我们可以将其统一修改 成 clusterID：
cp all.cds.cdhit unigene_cds.fasta
# 提取原始unigene ID
awk '$1 ~/ˆ>/{print $1}' all.cds.cdhit | \
sed 's/ˆ>//' > unigene.id
# 基于unigene ID 提取其氨基酸序列文件
singularity exec ../../software/MetaGenome.sif seqtk subseq \
all.pep.fa unigene.id > unigene_pep.fasta
# 提取cluster ID和原始unigeneID 对应关系
awk '{if($1~/ˆ>/){printf $1 $2} else if($NF ~ /*$/){print "\t"$3}}' \
all.cds.cdhit.clstr | \
sed 's/>//g; s/...$//'| \
awk '{print $2"\t"$1}' > map_id.txt
# 对原始unigeneID进行替换(原文件修改)
../script/map_fasta_ids map_id.txt unigene_cds.fasta
../script/map_fasta_ids map_id.txt unigene_pep.fasta
输出结果：
unigene_cds.fasta ：unigene CDS序列文件
unigene_pep.fasta ：unigene 氨基酸序列文件