初步大纲
我三代数据量不多，只有30x，师姐说用二代补。她说先组装scaffold，然后polish，然后scaffolding，然后gapcloser。
说实话 不是很懂为什么这个顺序。跟chatgpt对话的缺点是一旦你不懂他说的，质疑他一下，他就会根据你说的改他的答案。
目前阶段我就是师姐让干啥我就干啥...以后写文章再细查吧...
1.shasta拼接ont（或者使用soapdenovo，有时间可以两种方法一起都做，看哪个质量好）
运行命令：
# Download the executable for the latest release.
curl -O -L https://github.com/paoloshasta/shasta/releases/download/0.11.1/shasta-Linux-0.11.1
# Grant execute permissions.
chmod ugo+x shasta-Linux-0.11.1
mkdir 01.shasta
~/01.Software/shasta-Linux-0.10.0 --input all.fa --config Nanopore-May2022
生成文件在 北京集群03.result/01.shasta/ShastaRun中。主要就是Assembly.fasta
Nanopore-May2022是默认运行参数
2.pilon初步组装后的polish
工作路径：mkdir 02.pilon
教程参考：https://www.jianshu.com/p/cceeb7d1f413
bwa index -p STR hj.dmo.cns.fa #smartdenovo的拼装质量好一些所以换了fa文件
bwa mem -t 20 STR /PATH/TO/CL100047426_L02_12_1.clean.fq.gz /PATH/TO/CL100047426_L02_12_2.clean.fq.gz> bwa.run.o 2> bwa.run.e &
samtools sort -@ 10 -O bam -o align.bam < bwa.run.o
#分别对12，13，14三个文件夹的WGS数据fq进行比对，然后进行samtools merge，得到merge.bam
samtools merge merge.bam align12.bam align13.bam align14.bam
samtools index merge.bam
sambamba markdup -t 10 merge.bam merge_markdup.bam
samtools view -@ 10 -q 30 merge_markdup.bam > align_filter.bam
samtools index align_filter.bam(未成功)(其实输出的是个sam，要转一下格式，把merge_markdup.bam的header给了merge.filter)cat
java -Xmx500G -jar /~/01.Software/pilon-1.24.jar --genome hj.dmo.cns.fa --frags merge_filter_withheader.bam --fix snps,indels --threads 48 --output pilon_polished --vcf（qsub提交，运行了差不多12小时。中间有一次内存不足，所以多设置了内存）
最后生成了pilon_polished.fasta，还有vcf文件
pilon_polished.fasta用于round2输入，重复上述步骤，共进行3次得到pilon_polished_r3.fasta
3.arks软件 10x数据scaffolding
mkdir 03.arks
分别安装arcs需要的软件，可参考：http://yangl.net/2015/11/12/abyss_install/
Boost (tested on 1.61)
GCC (6+)
Autotools (if cloning directly from repository)#解析安装包的
LINKS (tested on 1.8)
Google SparseHash #https://github.com/sparsehash/sparsehash
ABySS (if using long reads) #wget https://www.bcgsc.ca/downloads/abyss/abyss-1.3.7.tar.gz
btllib (1.4.3+)
最后找到了青岛集群的软件
先用longranger处理10x的fq文件
longranger basic \
> --id=Hc\
> --fastqs=longranger/10xfqs #fq的存放路径\
> --sample=Hc
输出一直提示路径不对，但是检查了很多次都是没问题的
是因为我的10x数据是BGI seq500测的，需要转格式（/03.10x_scaffold/step02.conversion.sh）输出到10x_out文件夹，后续在青岛集群进行分析。
/bin/longranger-2.2.0/longranger basic --id=Hc --fastqs=/data/10x_out \--sample=CL100065457 --jobmode=local --localcores=8 --localmem=60 --noexit
下一步使用ark的perl脚本处理生成的barcode-multiplicities
参照https://protocols.faircloth-lab.org/en/latest/protocols-computer/assembly/assembly-scaffolding-with-arks-and-links.html#purpose
运行arks.sh脚本中间报错说缓存不够了，单独把一行脚本提出来arkerr.sh单跑，然后注释掉01.bin/anaconda/envs/scaffolding/bin/tigmint-make里面的那行，生成文件后运行剩下代码。
最后生成了文件Hcari.contigs.tigmint.fa，busco检验
4.使用二代数据（WGS）gapcloser
三代可以用TGSGapcloser，但是minimap2总出问题，待会再试
师兄用的是软件gapcloser
需要soapdenovo构建lconfig文件
在/04.gapclose/01.run运行。
configfile
#输入配置文件参数
max_rd_len=100 #最大的read长度，可以根据fastqc得到。
[LIB]
avg_ins=350 #文库平均插入长度(或取插入大小分布图中的峰值位置)，根据建库情况设置。
reverse_seq=0 #序列是否需要反转,PE(插入片段小于1kb)测序不需要，设为0，SE(插入片段在2kb-5kb)测序需要，设为1.
asm_flags=4 #1:表示只组装contig，2:表示只组装scaffold，3:表示同时组装contig和scaffold，4:表示只补gap。
rd_len_cutoff=100 #作用跟max_rd_len相同，大于该长度的序列会被切除至该长度，一般该参数不用。
rank=1 #如果建库时候建了多个不同大小的插入片段文库，需设不同rank跑。对于短片段设为1。
pair_num_cutoff=3 #至少有几对reads支持才可连接contig构建scaffold,小片段为3、大片段为5。
map_len=32 #对于pair-end，默认为32，对于mate-pair，默认为35。
#一对fastq文件
q1=/12/CL100047426_L02_12_1.clean.fq.gz
q2=/12/CL100047426_L02_12_2.clean.fq.gz
q1=/13/CL100047426_L02_13_1.fq
q2=/13/CL100047426_L02_13_2.fq
q1=/14/CL100047426_L02_14_2.fq
q2=/14/CL100047426_L02_14_1.fq
gapclose.sh
bin/GapCloser -a Hcari.contigs.tigmint.fa -b config_file -o gapcloser.fasta -t 4
运行后生成gapcloser.fasta。应该是最后结果了。可以busco测一下。
等以后不忙了再理一下原理。
