宏基因组SOAPdenovo组装

mkdir 0.raw_data
cd 0.raw_data
gzip -dc ../0.initial_data/A_1.fq.gz > A.1.fastq
gzip -dc ../0.initial_data/A_2.fq.gz > A.2.fastq
gzip -dc ../0.initial_data/B_1.fq.gz > B.1.fastq
gzip -dc ../0.initial_data/B_2.fq.gz > B.2.fastq
 
for i in `ls *.1.fastq`
do
    i=${i/.1.fastq/}
    echo $i
done > ../samples.txt
cd ..

mkdir 1.sequencing_data_quality_control
cd 1.sequencing_data_quality_control
 
# prepare adaptor cluster
echo "AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
	CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT
	AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
	AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
	CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT
	AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTG
	TTTTTTTTTTAATGATACGGCGACCACCGAGATCTACAC
	TTTTTTTTTTCAAGCAGAAGACGGCATACGA
	TACACTCTTTCCCTACACGACGCTCTTCCGATCT
	GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
	TACACTCTTTCCCTACACGACGCTCTTCCGATCT
	AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
	GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
	AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC" > primer_adaptor.txt
 
# remove low quality reads
for i in `cat ../samples.txt`
do
    echo "/opt/biosoft/NGSQCToolkit_v2.3.3/QC/IlluQC_PRLL.pl -pe ../0.raw_data/$i.1.fastq ../0.raw_data/$i.2.fastq primer_adaptor.txt 1 -l 60 -s 13 -c 8 -o $i"
done  > command.remove_low_quality_reads.list
ParaFly -c command.remove_low_quality_reads.list -CPU 4
 
# remove non-ATCG reads
for i in `cat ../samples.txt`
do
    echo "/opt/biosoft/NGSQCToolkit_v2.3.3/Trimming/AmbiguityFiltering.pl -i ./$i/$i.1.fastq_filtered -irev ./$i/$i.2.fastq_filtered -p 5"
done  > command.remove_non-ATCG_reads.list
ParaFly -c command.remove_non-ATCG_reads.list -CPU 4
cd ..

mkdir 2.genome_assembling/SOAPdenovo -p
cd 2.genome_assembling/SOAPdenovo
for i in `cat ../../samples.txt`
do
    echo "ln -s ../../1.sequencing_data_quality_control/$i/$i.1.fastq_filtered_trimmed $i.1.fastq"
    echo "ln -s ../../1.sequencing_data_quality_control/$i/$i.2.fastq_filtered_trimmed $i.2.fastq"
done  > command.copy_data.list
sh command.copy_data.list 
 
echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
q1=$HOME/2.genome_assembling/SOAPdenovo/A.1.fastq
q2=$HOME/2.genome_assembling/SOAPdenovo/A.2.fastq" > A_config.txt
 
mkdir A
SOAPdenovo-127mer all -s A_config.txt -o ./A/A -K 81 -d 1 -M 3 -p 8 -R -u -F &> A_soapdenovo.log
 
cd A
fasta_no_blank.pl A.scafSeq > A.scafSeq.noblank && cp A.scafSeq.noblank A.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" A.Scaftigs
rename_fasta_identifiers.pl A.Scaftigs > A.Scaftigs.rename
genome_seq_clear.pl --seq_prefix A_Scaftigs_500 --min_length 500 A.Scaftigs.rename > A_Scaftigs_500.fasta
 
cd ..
 
#mapping to get unpaired reads
/opt/biosoft/SOAPaligner/soap2.21release/2bwt-builder A.Scaftigs.rename
/opt/biosoft/SOAPaligner/soap2.21release/soap -a ../../../2.genome_assembling/SOAPdenovo/A.1.fastq -b ../../../2.genome_assembling/SOAPdenovo/A.2.fastq -D A.Scaftigs.rename.index -o A_PE_output -u A_unmap_output -2 A_SE_output -m 200 -p 20
 
echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
q1=$HOME/2.genome_assembling/SOAPdenovo/B.1.fastq
q2=$HOME/2.genome_assembling/SOAPdenovo/B.2.fastq" > B_config.txt
 
mkdir B
SOAPdenovo-127mer all -s B_config.txt -o ./B/B -K 81 -d 1 -M 3 -p 8 -R -u -F &> B_soapdenovo.log
 
cd N
fasta_no_blank.pl B.scafSeq > B.scafSeq.noblank && cp B.scafSeq.noblank B.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" B.Scaftigs
rename_fasta_identifiers.pl B.Scaftigs > B.Scaftigs.rename
genome_seq_clear.pl --seq_prefix B_Scaftigs_500 --min_length 500 B.Scaftigs.rename > B_Scaftigs_500.fasta
 
#mapping to get unpaired reads
/opt/biosoft/SOAPaligner/soap2.21release/2bwt-builder B.Scaftigs.rename
/opt/biosoft/SOAPaligner/soap2.21release/soap -a ../../../2.genome_assembling/SOAPdenovo/B.1.fastq -b ../../../2.genome_assembling/SOAPdenovo/B.2.fastq -D B.Scaftigs.rename.index -o B_PE_output -u B_unmap_output -2 B_SE_output -m 200 -p 20
 
cd ..

mkdir Mix
cd Mix
cat $HOME/2.genome_assembling/SOAPdenovo/A/A_unmap_output $HOME/2.genome_assembling/SOAPdenovo/B/B_unmap_output > Mix_seq
echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
p=$HOME/2.genome_assembling/SOAPdenovo/Mix_seq" > Mix_config.txt
 
SOAPdenovo-127mer all -s Mix_config.txt -o ./Mix -K 81 -d 1 -M 3 -p 8 -R -u -F &> Mix_soapdenovo.log
 
fasta_no_blank.pl Mix.scafSeq > Mix.scafSeq.noblank && cp Mix.scafSeq.noblank Mix.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" Mix.Scaftigs
rename_fasta_identifiers.pl Mix.Scaftigs > Mix.Scaftigs.rename
genome_seq_clear.pl --seq_prefix Mix_Scaftigs_500 --min_length 500 Mix.Scaftigs.rename > Mix_Scaftigs_500.fasta

-s  string      配置文件
-p  int         程序运行时设定的cpu线程数，默认值[8]
-o  string      输出的文件名前缀 
-K  int         输入的K-mer值大小，K-mer值必须是奇数，默认值[23]
-a  int         假设初始化内存，为了避免进一步的再分配，默认[0G]
-d  int         为了减少错误测序的影响，去除kmers频数不大于该值（kmerFreqCutoff）的k-mer，默认值[0] 
-R  (optional)  利用read鉴别短的重复序列，默认值[NO]
-D  int         为了减少错误测序的影响，去除频数不大于该值的由k-mer连接的边，默认值[1]
-M  int         连contiging时，合并相似序列的强度，最小值0，最大值3，默认值为[1]
-e  int         两边缘之间的弧的权重大于该值，将被线性化，默认值[0]
-m  int         最大的kmer size（max 127）用于multi-kmer，默认[NO]
-E  (optional)  在iterate迭代之前，合并clean bubble功能，仅在当使用multi-kmer时且设置-M参数，默认[NO]
-F  (optional)  利用read对scaffold中的gap进行填补，默认[NO]
-u  (optional)  构建scaffolding前不屏蔽高/低覆盖度的contigs，高频率覆盖度指平均contig覆盖深度的2倍。默认[屏蔽]
-w  (optional)  在scaffold中，保持contigs弱连接于其他contigs，默认[NO]
-G  int         估计gap的大小和实际补gap的大小的差异值，默认值[50bp]
-L  int         用于构建scaffold的contig的最短长度(minContigLen)，默认为：[Kmer参数值+2]
-c  float       在scaffolding之前，contigs小于100bp，且覆盖率小于该最小contig覆盖率（c*avgCvg），将被屏蔽，除非参数-u已设定，默认值[0.1]
-C  float       在scaffolding之前，contigs覆盖率大于该最大contigs覆盖率（C*avgCvg），或者contigs小于100bp且覆盖率大于0.8*（C*avgCvg），将被屏蔽，除非参数-u已设定，默认值[2]
-b  float       当处理contigs间的pair-end连接时，如果参数>1，该插入片段的上限值（b*avg_ins）将被用来作为大的插入片段（>1000）的上限，默认值[1.5]
-B  float       如果两个contigs的覆盖率都小于bubble覆盖率（bubbleCoverage）乘以contigs平均覆盖率（bubbleCoverage*avgCvg），则去除在bubble结构中的低覆盖率contig，默认值[0.6]
-N  int         统计基因组大小，默认值[0]
-V  (optional)  输出相关信息(Hawkeye)为了可视化组装，默认[NO]

max_rd_len=150        #read的最大长度
[LIB]                 #文库信息以此开头
avg_ins=300           #文库平均插入长度
reverse_seq=0         #序列是否需要被反转，插入片段大于等于2k的采用了环化，所以对于插入长度大于等于2k文库，序列需要反转，reverse_seq＝1，小片段设为0
asm_flags=3           #设为1：reads只用于只构建contig；设为2：reads只用于构建scaffold；设为3：reads同时用于构建contig和scaffold；设为4：reads只用于补洞gap
rank=1                #rank该值取整数，值越低，reads越优先用于构建scaffold
pair_num_cutoff=3     #该参数规定了连接两个contig或者是pre-scaffold 的可信连接的阈值，即，当连接数大于该值，连接才算有效。小片段文库(<2k)默认值[3]，大片段文库默认值[5]。
map_len=32            #该参数规定了在map过程中reads和contig的比对长度必须达到该值才能作为一个可信的比对。小片段文库一般设32，大片段文库一般设35，默认值[K＋2]。
q1=$DATA_PATH/read_1.fq
q2=$DATA_PATH/read_2.fq #双端fastq序列，写绝对路径
f1=./pathfasta_read_1.fa
f2=./pathfasta_read_2.fa #双端fasta序列，写绝对路径
q=./pathfastq_read_single.fq #单向fastq
f=./pathfasta_read_single.fa #单向fasta
p=./pathpairs_in_one_file.fa #双向测序得到的一个fasta格式的序列文件