conda install -c bioconda megahit

wget https://github.com/voutcn/megahit/releases/download/v1.2.8/MEGAHIT-1.2.8-Linux-x86_64-static.tar.gz
tar zvxf MEGAHIT-1.2.8-Linux-x86_64-static.tar.gz
cd MEGAHIT-1.2.8-Linux-x86_64-static/bin/
export PATH =<路径>:$PATH

# 常用参数：
-1/ -2 :分别输入双端序列
-m/--memory：SdBG 图构建过程中占用的最大内存字节数（如果设置为0-1，表示硬件内存的百分比。默认0.9）
--min-contig-len 设置拼接的contig最小长度
-t/--num-cpu-threads：使用的CPU线程数, 如果GPU允许激活，最少设置为2。
-o/--out-dir <字符型> 输出的文件夹
--out-prefix <字符型> 输出文件名的前缀
--k-list <奇数形式> 逗号分隔的kmer大小的列表，每个数必须是15-255的奇数,每两个相邻的kmer之间的增幅是小于28;默认：[21,29,39,59,79,99,119,141]
# 高级参数：
--no-mercy  Mercy k-mer用于宏基因组组装低覆盖序列。对于一般数据集>=30x，megahit可以使用-no-mercy选项生成更好的结果。
--low-local-ratio 定义低的局部覆盖的比率阈值 [0.2]
--max-tip-len  删除小于此值的端点[2*k]
--no-local  关闭局部组装模式

# 多样品组装
megahit -1 a1.fq,b1.fq -2 a2.fq,b2.fq -o final_megahit.out 
# 组装时候选组MEGHIT默认参数（MEGAHIT默认参数组装得到的contigs N50较⾼，质量较好）进⾏组装。也可以自己设置kemr的范围及最小Contig长度`--k-list 21,29,39,59,79,99,119,141 --min-contig-len 500`。

# 逐个单独组装
megahit -t 2 -m 0.95 --min-contig-len 500 -1 a1_.fq -2  a2_.fq -o  a_Contig.out
megahit -t 2 -m 0.95 --min-contig-len 500 -1 b1_.fq -2  b2_.fq -o  b_Contig.out

# 这里以样品a为例，首先根据它的Congtig序列构建索引
2bwt-builder final.contigs.fa  # 在a_Contig.out目录下

# 将质控后的序列比对至contig上
soap -p 6 -r 2 -m 200 -x 400 -M 4 -a a1_.fq -b a2_.fq -D final.contigs.fa.index -o a_PE.soap -2 a_SE.soap -u a_UN.fasta
# 参数介绍
-p  处理器数目，默认为1；
-r  比对到多个位置的序列如何输出，1表示不输出，2表示任意输出一次，3表示全部输出，默认为1；
-m  最小插入片段长度，默认400，single-end不需要设置该参数；
-x  最大插入片段长度，默认600，single-end不需要设置该参数；
-a  输入文件，pair-end时，输入其中一端文件；
-b  single-end不需要设置，pair-end输入另一端文件；
-u  没有比对上的reads；
-M  <int>   vi为每个read或read的seed部分匹配模式，不应该超过俩个误配。

1. 可以分别先将标识1 和 2的序列ID取出来,然后再根据id取序列
less a_UN.fasta |grep "1$" |sed 's/>//'> a1_unmapping.id
less a_UN.fasta |grep "2$" | sed 's/>//'> a2_unmapping.id
2. 根据id取对应的序列
seqtk subseq a_UN.fasta a1_unmapping.id  > a1_unmapping.fa
seqtk subseq a_UN.fasta a2_unmapping.id   > a2_unmapping.fa

####split_unmaping.py 
import sys
seq_file = {}
with open(sys.argv[1],"r") as old_fas: 
    for line in old_fas:
        line = line.strip()
        if line[0] == '>':
            seq_id = line
            seq_file[seq_id] = ''
        else:
            seq_file[seq_id] += line
new_fas1 = open(sys.argv[2],'w')
new_fas2 = open(sys.argv[3],'w')
for key,value in seq_file.items():
    if  '/1' in key:
        print(key + '\n' +value,file  = new_fas1)
    else:
        print(key + '\n'+ value, file = new_fas2)

# 同样的参数进行组装
megahit --min-contig-len 500 -1 a1_unmapping.fa,b1_unmapping.fa -2 a1_unmapping.fa,b2_unmapping.fa -o unmapping_megahit.out

# 下载文件
wget https://github.com/aquaskyline/SOAPdenovo2/archive/r241.tar.gz
tar xzvf r241.tar.gz
cd SOAPdenovo2-r241/
make
# conda安装也可以

##### 配置文件config
max_rd_len=151
[LIB]
avg_ins=379 #输出插入长度
reverse_seq=0
asm_flags=3
rd_len_cutoff=151
rank=1
q1=./a1.fq
q2=./a2.fq
########**配置文件说明**  
max_rd_len    表示read的最大长度；
[LIB]    表示文库信息标签；
avg_ins    文库中插入片段平均长度
reverse_seq    序列是否需要被反转，0(不反转)，1(反转)，一般插入片段大于等于2k文库，在建库是会将插入片段进行环化，此时须设置该参数为1：
asm_flags    表示reads用于组装哪个部分，可设为1,2,3, 1表示reads仅用于contig组装，2表示reads仅用于scaffold组装，3表示reads同时用于contig和scaffold组装；
rank    构建scaffold时，不同文库中reads的使用顺序，文库中reads序列越短，级别越高；
q1,q2    用于组装的双端fq文件。

SOAPdenovo-127mer all -s config -K 69 -M 3 -R -F -d 1 -u -o a.contig
## 参数
-s     config配置文件；
-K    k-mer的长度；
-o    输出文件前缀；
-d    [INT], kmerFreqCutoff, 去除频数小于等于该值的kmers，默认为0；
-M   [INT], mergeLevel连接contigs时, 合并相似序列的等级，默认为1，最小值为0，最大值为3，
-u    构建scaffold前屏蔽过高或过低覆盖度contigs，默认屏蔽；
-F     利用reads对scaffolds的gap进行填补，默认不执行；
-p    需要使用的cpu数目，默认8；

max_rd_len=151
[LIB]
avg_ins=379 #输出插入长度
reverse_seq=0
asm_flags=3
rd_len_cutoff=151
rank=1
q1=./a1_unmapping.fa
q2=./a2_unmapping.fa
q1=./b1_unmapping.fa
q2=./b2_unmapping.fa