組裝三代番木瓜基因組——by Serenity Fang

# 估算測序深度、reads數目、N50等值（自寫perl程序）：

$ perl ~/TangerScript/fqStat -i sunset.raw.subreads.fastq -g 372m

統計結果如下：

# 基因組組裝三步走1. Correction 2. Assembly 3. Polish 

## Step1: canu組裝（1. Correction 2. Assembly）

$ (nohup) canu -s spec.txt -p sunset -d sunset-auto genomeSize=400m -pacbio-raw sunset.raw.subreads.fastq &

$ cat spec.txt  注：spec文件為配置文件，根據不同服務器設置不同的參數。

   ### 組裝初步結果如下（自寫perl程序）：

$ cd /public1/home/Serenity/Sunset_Assembly/Canu-sunset-auto-201704

$ perl ~/perl_scripts/faSize.pl sunset.contigs.fasta

   ### 抽取unassembled.fasta中reads>5的contigs（自寫python程序）

$ python ~/python_scripts/extract_faread_filter.py sunset.unassembled.fasta

   ### 將上一步結果與 sunset.contigs.fasta合並

$ cat sunset.contigs.fasta sunset.unassembled.fastareadfilter > sunset.all.contigs.fasta

 

## Step2: 第一輪矯正（3. Polish）： quiver——取至少50x的三代數據做校正

$ cd /public1/home/Serenity/Sunset_Assembly/Canu-sunset-auto-201704/canu-quiver

$ ln -s ../sunset.all.contigs.fasta .

$ perl ~/TangerScript/runQuiver.pl -i sunset.all.contigs.fasta -d /public4/zhangxt/DATA/Papaya/sunset/baxh5 -t 16    注：run Quiver矯正，-t 設置節點數16-24個

$ for i in {1..27};do qsub script/script.${i}.pbs; done     注：結束后檢查outcmp里面的文件數目，檢查無誤后提交quiver.sh腳本

$ qsub quiver.sh   注：結束后得到consensus.fasta文件便是quiver校正后的基因組文件

  

## Step3: 第二輪矯正（3. Polish）： pilon——取至少50x的二代數據做校正

$ cd /public1/home/Serenity/Sunset_Assembly/sunset-reseq-raw-data

   ### 首先統計read長度、read數目、總鹼基數

$ zcat papaya_S1FR_CAGATC_L000_R1.fastq.gz | awk 'NR==2{a=length($1)}END{print "read length:"a"\nread num:"NR/4"\ntotal base:"a*NR/4*2"\n"}' > papaya_S1FR_CAGATC_L000_R1.fastq.gz.qstat.txt

$ cat papaya_S1FR_CAGATC_L000_R1.fastq.gz.qstat.txt   注：測序深度=total base/372000000

   ### bwa mem進行align

$ bwa index -a bwtsw consensus.fasta

$ bwa mem -t 24 -R '@RG\tID:S1FR_CAGATC\tSM:S1FR_CAGATC\tPL:Illumina\tLB:lib1' consensus.fasta papaya_S1FR_CAGATC_L000_R1.fastq.gz papaya_S1FR_CAGATC_L000_R2.fastq.gz > papaya_S1FR_CAGATC_L000.sam

$ samtools view -bS papaya_S1FR_CAGATC_L000.sam > papaya_S1FR_CAGATC_L000.bam

$ samtools sort papaya_S1FR_CAGATC_L000.bam -o papaya_S1FR_CAGATC_L000.sorted.bam

$ samtools index papaya_S1FR_CAGATC_L000.sorted.bam

$ qsub run_pilon.sh

$ cat run_pilon.sh    注：在本實驗室服務器指定13節點或者14節點，因為這兩個節點內存比較大，java設置內存300G，線程設置12以上 

 

   ### 組裝最終結果如下：

 $ perl ~/perl_scripts/faSize.pl sunset_pilon.fasta

  注：N50大概達到了1.2M，總基因組大小大概組裝到了330M