---
# System prepended metadata

title: Yellow Cardinal (Gubernatrix cristata) reference genome assembly (PacBio HiFi reads)

---

# Yellow Cardinal (Gubernatrix cristata) reference genome assembly (PacBio HiFi reads)

**Genome assembly naming** (according to Rhie et al. 2021, and corroborated in CGP species list - page 23): *bGubCri1.1*

![](https://i.imgur.com/FMXuVyb.jpg)
![](https://i.imgur.com/2k8MGb3.jpg)
*Photos of the individual being sequenced for species reference genome.*

## Background knowledge

Genome Assemblies of close species:

![](https://i.imgur.com/Nb4Qugb.png)

## Sample Collection
Sample ID: 5697
Sample type: **blood in 96% ethanol**
Sex: **female**
Collectors: Dr. Bettina Mahler and Lic. Melina Atencio in the rescue center Fundacion Temaiken (CRET)
Origin: confiscated, in Rescue Center since 16.10.2019. Assigned to the management unit of the east, where there is no hybridization (markers: mtDNA + 10 microsatellites),
Date of collection: **19.04.2021**
Date of sample delivery from Argentina to  Germany: **19.04.2021**
Date of sample arrival to Uni Potsdam: **10.05.2021**

## DNA isolation
Date of DNA extraction: **17.05.2021**
Extraction kit for MHW DNA: Circulomics Nanobind CBB Big DNA
Protocol: [here](https://drive.google.com/file/d/1h-WXobb3CII9CpjwkXuXini435InkVqo/view?usp=sharing) but before 2X washes with PBS to remove ethanol:

80 ul of Ethanol/blood mix pipetted in Protein LoBind tubes.
1. Spin 5 min 2200 g 4°C to remove ethanol
2. Add 500 ul of 1x PBS to wash. Spin 5 min 2200 g 4°C and remove supernatant. Repeat for a total of 2 washes.
3. Pellet is directly dissolved in 200 ul of 1x PBS.

I sent 3 extractions from same individual to be considered for sequencing.

## Sequencing
DRESDEN-concept Genome Center (MPI-Dresden)
Contact: Dr. Sylke Winkler (winkler@mpi-cbg.de)
Date of sample delivery to sequencing center: **20.07.2021**
Date of sample reception in the sequencing center: **21.07.2021**
Template Specimen Form: [here](https://drive.google.com/file/d/1UtUPy8uKe2T8J_iR-3ePsoxTK7461igF/view?usp=sharing)
Date of sequencing completion: 20.08.2021
Link to Data (mail Pippel): https://bds.mpi-cbg.de/pippelshare/bGubCri_qpw6S5CzT2zAUpy11I5kF2HfPAxgHk2s8GDy5YKMbn/
QC in sequencing center (extraction N3 chosen): [here](https://drive.google.com/file/d/1FUG46828MQJY2LbqZ3P4qJbRvu_qaztq/view?usp=sharing)

## PacBio HiFi reads pre-processing
1. **bam to fq**: bam2fastx using SMRTtool (also tested with bedtools 2.9)
[submit_bam2fastx_newHPC.slurm](https://drive.google.com/file/d/1lSYmctsZnRZmFoMAxNEUQFu0KblVdWWT/view?usp=sharing)
```
#1. index
module purge
module load lang/Anaconda3/2020.11
source activate pbbam-env

pbindex /path/m64046_210816_140027.ccs.bam

#2. convert bam to fq
module load lang/Anaconda3/2020.11
source activate bam2fastx-env

bam2fasta -u -o /path/bam2fastx/ /path/m64046_210816_140027.ccs.bam

bam2fastq -u -o /path/bam2fastx/ /path/m64046_210816_140027.ccs.bam

conda deactivate
```
*I used bam2fasta 1.3.1 and bam2fastq 1.3.1*


**2. Trimming of PacBio SMRT adapters**  from reads (cutadapt) --> 4.5 % reads removed
[submitcutadaptnewHPC.slurm](https://drive.google.com/file/d/1_Ni22w8yN3suxcHC0J6GH7A7IbtARzht/view?usp=sharing)
```
cutadapt -j 10 -o /path/reads_noAdapt.fastq /path/reads.fastq -b ATCTCTCTCAACAACAACAACGGAGGAGGAGGAAAAGAGAGAGAT -b ATCTCTCTCTTTTCCTCCTCCTCCGTTGTTGTTGTTGAGAGAGAT --discard-trimmed
```
# Reads Statistics
**1. Fastqc on trimmed reads**
[submit_fastqc_trimmed_reads.slurm](https://drive.google.com/file/d/154yKNeLz8_4fHKZwD1vbQ0DNk19yD4uW/view?usp=sharing)
```
fastqc -t 6 /path/reads_noAdapt.fastq -o /path/fastqc_TrimmedReads
```
* link to full fastqc report **before** trimming sequencing adapters from reads [here](https://drive.google.com/file/d/1ppck5pY7ZYKCswjB4ArVzNu7WIARzEJw/view?usp=sharing)
* link to full fastqc report on **trimmed** reads [here](https://drive.google.com/file/d/1kO6I8tKme3-2Zy6u2LoTLcQKX13QKpDj/view?usp=sharing), main results:

![](https://i.imgur.com/gwG3BS9.png)

**2. Plot reads length**
```
#not working in HPC, I try in local computer Alan (I limit the transfer to 1Mb/sec so it does not get stalled)

#Transfer from HPC to local (and then from local to Alan):
cd /home/Marisol

scp dominguez@n-hpc-login1.rz.uni-potsdam.de:/work/dominguez/YelCar_GenomeAssembly_MPI/reads/bam2fastx/reads_noAdapt.fasta .

#it was interrupted by lost connection
#resume transfer:

rsync -P --rsh=ssh dominguez@n-hpc-login1.rz.uni-potsk/dominguez/YelCar_GenomeAssembly_MPI/reads/bam2fastx/reads_noAdapt.fasta ./reads_noAdapt.fasta

#From local to Alan:
scp ./reads_noAdapt.fasta alan@141.89.108.18:/home/alL/reads_trimmed/

python plot_fasta_length_v3.py /home/alan/Desktop/MARISOL/reads_trimmed/reads_noAdapt.fasta /home/alan/Desktop/MARISOL/reads_trimmed/reads_noAdapt.fasta.length.png
```
![](https://i.imgur.com/rMCh5LP.png)

**3. Reads statistics**
Count reads on .fasta file (N = 2 489 222)
```
grep ">" /path_to_reads/reads_noAdapt.fasta | wc -l
```
Count reads on .fastq file (N = 2 489 222)
```
expr $(cat /path_to_reads/reads_noAdapt.fastq | wc -l) / 4
```
Using asmstats script:
[submit_trimmed_asmstats.slurm]([https:/](https://drive.google.com/file/d/1Rjb8n7jA5L60mSb1NNQAwoh7TOv_cSKK/view?usp=sharing)/)
```
module purge
module load lang/Anaconda3/2020.11
source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats /path_to_reads/reads_noAdapt.fastq  > /path/reads_statistics/reads_noAdapt.fastq.stats

```
Outfile reads_noAdapt.fastq.stats [here](https://drive.google.com/file/d/1_jhrUgOdlufE5bDcxJqbYsdVCj0q6FYI/view?usp=sharing).
![](https://i.imgur.com/F74vQEg.png)

**4.  Kmer analysis on reads**
Assuming k21:
[submit_Jellyfish_trimmed_k21.slurm](https://drive.google.com/file/d/1id-VkCkCUNTB2nJIh49_3OmBZHUdoHnS/view?usp=sharing)
*I used jellyfish version 2.2.10*
```
module purge
module load lang/Anaconda3/2020.11
source activate jellyfish-env

#jellyfish count
jellyfish count -C -m 21 -s 1000000 -t 20 -o /path/reads_statistics/YelCar_reads_trim.k21.jf /path/reads/reads_noAdapt.fastq

#jellyfish histo
jellyfish histo /path/reads_statistics/YelCar_reads_trim.k21.jf > /path/reads_statistics/YelCar_reads_trim.k21.histo
```

Assuming k31:
[submit_Jellyfish_trimmed_k31.slurm](https://drive.google.com/file/d/1_GhZVkBATwnBDOngLZN4yhSVSpv1drkh/view?usp=sharing)
```
module purge
module load lang/Anaconda3/2020.11
source activate jellyfish-env

jellyfish count -C -m 31 -s 1000000 -t 20 -o /path/reads_statistics/YelCar_reads_trim.k31.jf /path/reads/reads_noAdapt.fastq

#jellyfish histo
jellyfish histo /path/reads_statistics/YelCar_reads_trim.k31.jf > /path/reads_statistics/YelCar_reads_trim.k31.histo
```

Open the histo files in [genomescope2 website](http://qb.cshl.edu/genomescope/genomescope2.0/) 

[k21.histo](https://drive.google.com/file/d/1yvNxjOaxRKsiGfDcAfLmUArRstqzcmfZ/view?usp=sharing) 

![](https://i.imgur.com/814sXYx.png)
Link to genomescope2 results [here](http://qb.cshl.edu/genomescope/genomescope2.0/analysis.php?code=5d0pXxIZINtIM6igWHcg)

Kmer = 21
Size: 1 057 467 553 bp ~ 1.06 Gb
Herterozygosity: 1.1%
Repts (100-uniq): 11.8%

Parameters for purge_dups:
- Estimated genome haploid size: 1057467553 bp
- Transition Parameter (kmercov*1.5): 21.045 -> 21 reads (integer)
- Maximum read depth (Transition Par*3): 63.135 -> 63 reads (integer)


[k31.histo](https://drive.google.com/file/d/19zt3RybTbybSlnzphghbwDwRIS2cB8yV/view?usp=sharing) 

![](https://i.imgur.com/zxIeYwI.png)

Link to genomescope2 results [here](http://qb.cshl.edu/genomescope/genomescope2.0/analysis.php?code=VtMACah2eEURIERJjRmW)

Kmer = 31
Size:  1 058 655 883 bp ~ **1.06 Gb**
Herterozygosity: 0.99 %
Repts (100-uniq): 9.9 %

Parameters for purge_dups:
- Estimated genome haploid size: 1 058 655 883 bp
- Transition Parameter (kmercov*1.5): 20.7 -> 21 reads (integer)
- Maximum read depth (Transition Par*3): 62.1 -> 62 reads (interger)


**Reads Coverage:  30.3X**
(number of bases sequenced -sum value in .stats file-) / (genome size in genomescope2 according to kmer analysis)
(32036044402/1058655883)

The **estimated genome size** of both kmers is around **1.06 Gb**, this is also close to the **zebra finch** (Taeniopygia guttata, 1.056 Gb), and other birds like **small tree finc**h (Camarynchus parvulus, 1.064 Gb), and **medium ground finch** (Geospiza fortis, 1.065 Gb).

![](https://i.imgur.com/kc8d6Ig.png)

## Genome Assemblies
1. **Hifiasm** (https://github.com/chhylp123/hifiasm) with different parameters (l0,l1,l2 and nok, k21, k31)
*I used Hifiasm version 0.16.1-r375*


- scripts for assemblies *(change parameters -l and -k accordingly)*
```
module purge
module load lang/Anaconda3/2020.11
source activate hifiasm-env

hifiasm -t 28 -l2 -- primary -o /path/asm/hifiasm.trim.l2.k31.asm /path/reads/reads_noAdapt.fastq -k 31

#took Real time: 27160.206 sec; CPU: 922008.523 sec; Peak RSS: 54.377 GB requesting 38 cpus-per-task with 16G mem-per-cpu.
```

- script for statistics for the asm:
```
/path/gfa2fa /path/hifiasm.trim.l2.k31.asm.bp.p_ctg.gfa > /path/YelCar_hifiasm_trimmed_l2_k31.asm.p_ctg.fa

module purge
module load lang/Anaconda3/2020.11
source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats /path/YelCar_hifiasm_trimmed_l2_k31.asm.p_ctg.fa > /path/hifiasm.trim.l2.k31.asm.p_ctg.fa.stats
```
### Primary contigs viewed by Bandage:
Primmary assembly:
YelCar_hifiasm_trim_l2.k31.PRI.gfa 
![](https://i.imgur.com/FecY4fY.png)
*Looks good, most of the contigs must do not have 2 branches.*

Zoom-in of all contigs I found with  bubbles:
![](https://i.imgur.com/clbmCWd.png)


*(hopefully these will be solved after purging - I can not corroborate since purgedups does not provide a .gfa file)*

Alternative assembly:
YelCar_hifiasm_trimmed_l2_k31.asm.a_ctg.gfa
![](https://i.imgur.com/wR9bb2A.png)
*It is more fragmented, and much more bubbles.*


- scripts for purging - 1st ROUND (purge PRI p1 to obtain **c1**)


```
module purge
module load lang/Anaconda3/2020.11
source activate purge_dups_env


## Change this: ##########################################################
OUT_DIR="/path/hifiasm/purge_dups"
NAME="c1.l2.k31"
ASM="/path/YelCar_hifiasm_trimmed_l2_k31.asm.p_ctg.fa"
ASM_NAME=$(basename $ASM .fa)
HIFI_DIR="/path_to_reads"
TRANS=21
MAX_DEPTH=62
CUTOFFS=$TRANS\_$MAX_DEPTH
###############################################################################

if [ ! -d $OUT_DIR/$NAME ]
then
  mkdir -p $OUT_DIR/$NAME
fi

##
# Align HiFi reads and generate BAM files, then calculate read depth histogram and base-level read depth
cd $OUT_DIR/$NAME

for i in $(ls $HIFI_DIR/*_noAdapt.fastq.gz)
do
	READ_NAME=$(basename $i .fastq.gz)
	minimap2 -xasm20 -I 4G -t 24 $ASM $i | gzip -c - > $OUT_DIR/$NAME/$READ_NAME\.paf.gz
done
## -xmap-hifi is not in latest release. See https://github.com/lh3/minimap2/issues/739
##  we will use the Reads mapped to contigs pairwise mapping format (PAF) file for calculating some statistics required in a later stage
## purge_dups initially produces a read-depth histogram from base-level coverages. This information is used for estimating the coverage cutoffs, taking into account that collapsed haplotype contigs will lead to reads from both alleles mapping to those contigs, whereas if the alleles have assembled as separate contigs, then the reads will be split over the two contigs, resulting in half the read-depth (Roach et al. 2018).

pbcstat *.paf.gz 
# this produces PB.base.cov and PB.stat files
calcuts PB.stat > cutoffs_default 2>calcults_default.log
python3 /path/hist_plot.py -c cutoffs_default PB.stat PB_default.png

### VGP suggested cutoffs:
calcuts -m$TRANS -u$MAX_DEPTH PB.stat > cutoffs_$CUTOFFS
python3 /path/hist_plot.py -c cutoffs_$CUTOFFS PB.stat PB_$CUTOFFS\.png

##
# Split an assembly and do a self-self alignment
cd $OUT_DIR/$NAME
split_fa $ASM > $ASM_NAME\.split
minimap2 -t 18 -xasm5 -DP $ASM_NAME\.split $ASM_NAME\.split | pigz -p 6 -c - > $ASM_NAME\.split.self.paf.gz

##### 

# Purge haplotigs and overlaps:
cd $OUT_DIR/$NAME
purge_dups -2 -T cutoffs_default -c PB.base.cov $ASM_NAME\.split.self.paf.gz > dups_default.bed 2> purge_dups_default.log
purge_dups -2 -T cutoffs_$CUTOFFS -c PB.base.cov $ASM_NAME\.split.self.paf.gz > dups_$CUTOFFS\.bed 2> purge_dups_$CUTOFFS\.log

##
# Get purged primary and haplotig sequences from draft assembly:
cd $OUT_DIR/$NAME
get_seqs -e dups_default.bed $ASM
mv purged.fa purged_default.fa
mv hap.fa hap_default.fa
get_seqs -e dups_$CUTOFFS\.bed $ASM
mv purged.fa purged_$CUTOFFS\.fa
mv hap.fa hap_$CUTOFFS\.fa

############################################################################
module purge
module load lang/Anaconda3/2020.11
source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats /path/hifiasm/purge_dups/c1.l2.k31/hap_21_62.fa > /path/hifiasm/purge_dups/c1.l2.k31/k31_hap_21_62.fa.stats
asmstats /path/hifiasm/purge_dups/c1.l2.k31/purged_21_62.fa > /path/hifiasm/purge_dups/c1.l2.k31/k31_purged_21_62.fa.stats
asmstats /path/hifiasm/purge_dups/c1.l2.k31/hap_default.fa > /path/hifiasm/purge_dups/c1.l2.k31/k31_hap_default.fa.stats
asmstats /path/hifiasm/purge_dups/c1.l2.k31/purged_default.fa > /path/hifiasm/purge_dups/c1.l2.k31/k31_purged_default.fa.stats
########################################################################

# "cutoffs" file has six numbers: the first is low coverage, any contig with an average coverge less than the number is a junk contig, the forth number is diploid (middle) coverage, the last number is high coverage for repeats, contigs with an average coverage over this number are usually repeats 

# Next step: merge hap.fa and alt_asm and redo the above steps to get the haplotig set.

#job was ran with --cpus-per-task=24, #SBATCH --mem=64G. It took 15 min.
```


- scripts for purging - 2nd ROUND (purge ALT c2p2 to obtain **q2**)
- - using VGP suggested cutoffs:
```
module purge
module load lang/Anaconda3/2020.11
source activate purge_dups_env

## Change this: ##########################################################
OUT_DIR="/path/hifiasm/purge_dups"
NAME="c1.l2.k31"
ASM="/path/YelCar_hifiasm_trimmed_l2_k31.asm.p_ctg.fa"
ASM_NAME=$(basename $ASM .fa)
ALT="/path/YelCar_hifiasm_trimmed_l2_k31.asm.a_ctg.fa"
HIFI_DIR="/path_to_reads"
TRANS=21
MAX_DEPTH=62
CUTOFFS=$TRANS\_$MAX_DEPTH
###############################################################################

if [ ! -d $OUT_DIR/$NAME/q2_VGPcutoffs ]
then
  mkdir -p $OUT_DIR/$NAME/q2_VGPcutoffs
fi

#Merge hap_.fa and ALT 
##c2 is before purging so it is in another folder
cd $OUT_DIR/$NAME/q2_VGPcutoffs
cat $ALT >> c2p2.fa; cat /path/hifiasm/purge_dups/c1.l2.k31/hap_21_62.fa >> c2p2.fa

#align HiFi reads to c2p2 and generate BAM files, then calculate read depth histogram and base-level read depth
cd $OUT_DIR/$NAME/q2_VGPcutoffs
for i in $(ls $HIFI_DIR/*_noAdapt.fastq.gz)
do
        READ_NAME=$(basename $i .fastq.gz)
        minimap2 -xasm20 -I 4G -t 24 c2p2.fa $i | gzip -c - > $OUT_DIR/$NAME/q2_VGPcutoffs/$READ_NAME\.paf.gz
done
pbcstat *.paf.gz

### VGP suggested cutoffs:
calcuts -m$TRANS -u$MAX_DEPTH PB.stat > cutoffs_$CUTOFFS
python3 /path/hist_plot.py -c cutoffs_$CUTOFFS PB.stat PB_$CUTOFFS\.png
##


# Split an assembly and do a self-self alignment
cd $OUT_DIR/$NAME/q2_VGPcutoffs
split_fa c2p2.fa > c2p2.split
minimap2 -t 18 -xasm5 -DP c2p2.split c2p2.split | pigz -p 6 -c - > c2p2.split.self.paf.gz

#############
# Purge haplotigs and overlaps
cd $OUT_DIR/$NAME/q2_VGPcutoffs
purge_dups -2 -T cutoffs_$CUTOFFS -c PB.base.cov c2p2.split.self.paf.gz > dups_$CUTOFFS\.bed 2> purge_dups_$CUTOFFS\.log


# Get purged and haplotig sequences from draft assembly
cd $OUT_DIR/$NAME/q2_VGPcutoffs
get_seqs -e dups_$CUTOFFS\.bed c2p2.fa
mv purged.fa purged_$CUTOFFS\.fa
mv hap.fa hap_$CUTOFFS\.fa

#########################################################################
module purge
module load lang/Anaconda3/2020.11
source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats $OUT_DIR/$NAME/q2_VGPcutoffs/purged_$CUTOFFS\.fa > $OUT_DIR/$NAME/q2_VGPcutoffs/purged_$CUTOFFS\.fa.stats
asmstats $OUT_DIR/$NAME/q2_VGPcutoffs/hap_$CUTOFFS\.fa > $OUT_DIR/$NAME/q2_VGPcutoffs/hap_$CUTOFFS\.fa.stats

## requested --cpus-per-task=24, --mem=64G, and it took 15 min to run
```
Values based on **kmercov** from genomescope2: 
For k31:
- Transition Parameter (kmercov*1.5): 20.7 -> 21 reads (integer)
- Maximum read depth (Transition Par*3): 62.1 -> **62** reads (interger)
For k21:
- Transition Parameter (kmercov*1.5): ): 21.045  -> 21 reads (integer)
- Maximum read depth (Transition Par*3): 63.135 -> **63** reads (interger)

- - using default values
```
module purge
module load lang/Anaconda3/2020.11
source activate purge_dups_env

## change this: ##########################################################
OUT_DIR="/path/hifiasm/purge_dups"
NAME="c1.l2.k31"
ASM="/path/YelCar_hifiasm_trimmed_l2_k31.asm.p_ctg.fa"
ASM_NAME=$(basename $ASM .fa)
ALT="/path/YelCar_hifiasm_trimmed_l2_k31.asm.a_ctg.fa"
HIFI_DIR="/path_to_reads"
###############################################################################
if [ ! -d $OUT_DIR/$NAME/q2_default ]
then
  mkdir -p $OUT_DIR/$NAME/q2_default
fi

#Merge hap_.fa and ALT 
##c2 is before purging so it is in another folder
cd $OUT_DIR/$NAME/q2_default
cat $ALT >> c2p2.fa; cat /hifiasm/purge_dups/c1.l2.k31/hap_default.fa >> c2p2.fa

# align HiFi reads to c2p2 and generate BAM files, then calculate read depth histogram and base-level read depth
cd $OUT_DIR/$NAME/q2_default
for i in $(ls $HIFI_DIR/*_noAdapt.fastq.gz)
do
        READ_NAME=$(basename $i .fastq.gz)
        minimap2 -xasm20 -I 4G -t 24 c2p2.fa $i | gzip -c - > $OUT_DIR/$NAME/q2_default/$READ_NAME\.paf.gz
done
pbcstat *.paf.gz
calcuts PB.stat > cutoffs_default 2>calcults_default.log
python3 /path/hist_plot.py -c cutoffs_default PB.stat PB_default.png

# Split an assembly and do a self-self alignment
cd $OUT_DIR/$NAME/q2_default
split_fa c2p2.fa > c2p2.split
minimap2 -t 18 -xasm5 -DP c2p2.split c2p2.split | pigz -p 6 -c - > c2p2.split.self.paf.gz

# Purge haplotigs and overlaps
cd $OUT_DIR/$NAME/q2_default
purge_dups -2 -T cutoffs_default -c PB.base.cov c2p2.split.self.paf.gz > dups_default.bed 2> purge_dups_default.log


# Get purged and haplotig sequences from draft assembly
cd $OUT_DIR/$NAME/q2_default
get_seqs -e dups_default.bed c2p2.fa
mv purged.fa purged_default.fa
mv hap.fa hap_default.fa
#########################################################################
module purge
module load lang/Anaconda3/2020.11
source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats $OUT_DIR/$NAME/q2_default/purged_default.fa > $OUT_DIR/$NAME/q2_default/purged_default.fa.stats
asmstats $OUT_DIR/$NAME/q2_default/hap_default.fa > $OUT_DIR/$NAME/q2_default/hap_default.fa.stats

## requested --cpus-per-task=24, --mem=64G, and it took 15 min to run
```
2. IPA
https://github.com/PacificBiosciences/pbipa
*I used ipa (wrapper) version=1.3.1*

* script for the assembly of the primmary (PRI) and alternative (ALT) assemblies

```
module purge
module load lang/Anaconda3/2020.11
source activate snakemake-env2
source activate ipa-env

ipa local --nthreads 60 --njobs 1 -i /path/reads/filename_ccs_reads.fa

# It was ran requesting 50G of mem and lasted 2 days to finish.
```

- - 1st ROUND of PURGING (purge PRI p1 to obtain **c1**)
*I used purgedups v1.2.5*
https://github.com/dfguan/purge_dups

```
module purge
module load lang/Anaconda3/2020.11
source activate purge_dups_env

## Change this: ##########################################################
OUT_DIR="/path/IPA/purge_dups"
NAME="IPA_c1.purged"
ASM="/path/IPA/final.p_ctg.fasta"
ASM_NAME=$(basename $ASM .fa)
HIFI_DIR="/path_to_reads"
###############################################################################

if [ ! -d $OUT_DIR/$NAME ]
then
  mkdir -p $OUT_DIR/$NAME
fi

# Align HiFi reads and generate BAM files, then calculate read depth histogram and base-level read depth
cd $OUT_DIR/$NAME

for i in $(ls $HIFI_DIR/*_noAdapt.fastq.gz)
do
	READ_NAME=$(basename $i .fastq.gz)
	minimap2 -xasm20 -I 4G -t 24 $ASM $i | gzip -c - > $OUT_DIR/$NAME/$READ_NAME\.paf.gz
done
## -xmap-hifi is not in latest release. See https://github.com/lh3/minimap2/issues/739
## we will use the Reads mapped to contigs pairwise mapping format (PAF) file for calculating some statistics required in a later stage
## purge_dups initially produces a read-depth histogram from base-level coverages. This information is used for estimating the coverage cutoffs, taking into account that collapsed haplotype contigs will lead to reads from both alleles mapping to those contigs, whereas if the alleles have assembled as separate contigs, then the reads will be split over the two contigs, resulting in half the read-depth (Roach et al. 2018).

pbcstat *.paf.gz #(produces PB.base.cov and PB.stat files)
calcuts PB.stat > cutoffs_default 2>calcults_default.log
python3 /path/hist_plot.py -c cutoffs_default PB.stat PB_default.png

# Split an assembly and do a self-self alignment:
cd $OUT_DIR/$NAME
split_fa $ASM > $ASM_NAME\.split
minimap2 -t 18 -xasm5 -DP $ASM_NAME\.split $ASM_NAME\.split | pigz -p 6 -c - > $ASM_NAME\.split.self.paf.gz

# Purge haplotigs and overlaps:
cd $OUT_DIR/$NAME
purge_dups -2 -T cutoffs_default -c PB.base.cov $ASM_NAME\.split.self.paf.gz > dups_default.bed 2> purge_dups_default.log

# Get purged primary and haplotig sequences from draft assembly:
cd $OUT_DIR/$NAME
get_seqs -e dups_default.bed $ASM
mv purged.fa purged_default.fa
mv hap.fa hap_default.fa

############################################################################
module purge
module load lang/Anaconda3/2020.11

source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats $OUT_DIR/$NAME/hap_default.fa > $OUT_DIR/$NAME/IPA_hap_default.fa.stats

asmstats $OUT_DIR/$NAME/purged_default.fa > $OUT_DIR/$NAME/IPA_purged_default.fa.stats

#Next step: Merge hap.fa and alt_asm and redo the above steps to get the haplotig set.

# It was ran requesting --cpus-per-task=24, --mem=64G. It took 10 min.
```

- - 2nd ROUND of PURGING (purge ALT c2p2 to obtain **q2**)

```
module purge
module load lang/Anaconda3/2020.11
source activate purge_dups_env

## Change this: ##########################################################
OUT_DIR="/path/reads_trimmed/IPA/purge_dups/"
NAME="IPA_c1.purged"
ASM="/path_to_PRIasm/final.p_ctg.fasta"
ASM_NAME=$(basename $ASM .fa)
ALT="/path_to_ALTasm/final.a_ctg.fasta"
HIFI_DIR="/path_to_reads"
###############################################################################

if [ ! -d $OUT_DIR/$NAME/q2 ]
then
  mkdir -p $OUT_DIR/$NAME/q2
fi

# Merge hap_.fa and ALT 
##c2 is before purging so it is in another folder
cd $OUT_DIR/$NAME/q2
cat $ALT >> c2p2.fa; cat /path/IPA/purge_dups/IPA_c1.purged/hap_default.fa >> c2p2.fa

# Align HiFi reads to c2p2 and generate BAM files, then calculate read depth histogram and base-level read depth:
cd $OUT_DIR/$NAME/q2
for i in $(ls $HIFI_DIR/*_noAdapt.fastq.gz)
do
        READ_NAME=$(basename $i .fastq.gz)
        minimap2 -xasm20 -I 4G -t 24 c2p2.fa $i | gzip -c - > $OUT_DIR/$NAME/q2/$READ_NAME\.paf.gz
done

pbcstat *.paf.gz

calcuts PB.stat > cutoffs_default 2>calcults_default.log
python3 /path/hist_plot.py -c cutoffs_default PB.stat PB_default.png

# Split an assembly and do a self-self alignment:
cd $OUT_DIR/$NAME/q2
split_fa c2p2.fa > c2p2.split
minimap2 -t 18 -xasm5 -DP c2p2.split c2p2.split | pigz -p 6 -c - > c2p2.split.self.paf.gz

# Purge haplotigs and overlaps:
cd $OUT_DIR/$NAME/q2
purge_dups -2 -T cutoffs_default -c PB.base.cov c2p2.split.self.paf.gz > dups_default.bed 2> purge_dups_default.log

# Get purged and haplotig sequences from draft assembly:
cd $OUT_DIR/$NAME/q2
get_seqs -e dups_default.bed c2p2.fa
mv purged.fa purged_default.fa
mv hap.fa hap_default.fa
#########################################################################

module purge
module load lang/Anaconda3/2020.11

source activate seqtk-env

export PATH=$PATH:/work/dominguez/YelCar_GenomeAssembly_MPI/

asmstats $OUT_DIR/$NAME/q2/purged_default.fa > $OUT_DIR/$NAME/q2/purged_default.fa.stats
asmstats $OUT_DIR/$NAME/q2/hap_default.fa > $OUT_DIR/$NAME/q2/hap_default.fa.stats

# It was ran with --cpus-per-task=24, #SBATCH --mem=24G and took 10 min.
```

## Genome Evaluation Pipeline (GEP)
https://git.imp.fu-berlin.de/cmazzoni/GEP

Developers: James Sullivan and Dr. Camila Mazzoni.

STEP 1: building the meryl kmer database
```
module load lang/Anaconda3
source activate GEP

nohup snakemake --profile SUBMIT_CONFIG/slurm --cores 20 &

```
 STEP 2: run the genome asm evaluation
 ```
module load lang/Anaconda3
source activate GEP
nohup snakemake --profile SUBMIT_CONFIG/slurm --cores 80 --jobs 5 &

```
 
 
| Assembly               | hifiasm_l0_nok   | purgedups_l0_nok_default   | purgedups_l0_nok_VGPcutoffs   | hifiasm_l0_k21   | purgedups_l0_k21_default   | purgedups_l0_k21_VGPcutoffs   | hifiasm_l0_k31   | purgedups_l0_k31_default   | purgedups_l0_k31_VGPcutoffs   | hifiasm_l1_nok   | purgedups_l1_nok_default   | purgedups_l1_nok_VGPcutoffs   | hifiasm_l1_k21   | purgedups_l1_k21_default   | purgedups_l1_k21_VGPcutoffs   | hifiasm_l1_k31   | purgedups_l1_k31_default   | purgedups_l1_k31_VGPcutoffs   | hifiasm_l2_nok   | purgedups_l2_nok_default   | purgedups_l2_nok_VGPcutoffs   | hifiasm_l2_k21   | purgedups_l2_k21_default   | purgedups_l2_k21_VGPcutoffs   | hifiasm_l2_k31   | purgedups_l2_k31_default   | purgedups_l2_k31_VGPcutoffs   | IPA           | IPA_purged    |
|:-----------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:-----------------|:---------------------------|:------------------------------|:--------------|:--------------|
| total_num_bases        | 2,224,168,269    | 1,351,576,889              | 1,188,367,148                 | 2,222,471,140    | 1,346,836,325              | 1,182,863,508                 | 2,223,097,156    | 1,355,784,751              | 1,182,862,056                 | 1,332,099,590    | 1,176,726,747              | 1,176,647,996                 | 1,331,886,907    | 1,180,543,672              | 1,180,465,568                 | 1,330,877,578    | 1,179,270,444              | 1,177,972,303                 | 1,245,686,462    | 1,165,793,779              | 1,165,793,779                 | 1,244,971,205    | 1,172,156,874              | 1,172,156,874                 | 1,243,890,104    | 1,164,967,086              | 1,164,967,086                 | 1,080,623,443 | 1,070,029,210 |
| size_estimate          | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261    | 1,058,734,261              | 1,058,734,261                 | 1,058,734,261 | 1,058,734,261 |
| max_heterozygosity     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%        | 0.997831%                  | 0.997831%                     | 0.997831%     | 0.997831%     |
| GC_Content             | 43.4626          | 43.1192                    | 43.3178                       | 43.4481          | 43.1054                    | 43.2881                       | 43.4551          | 43.0657                    | 43.2764                       | 43.9681          | 43.2168                    | 43.2159                       | 43.9665          | 43.2379                    | 43.2376                       | 43.9685          | 43.2309                    | 43.2214                       | 43.7381          | 43.1243                    | 43.1243                       | 43.7162          | 43.1699                    | 43.1699                       | 43.7021          | 43.1111                    | 43.1111                       | 42.4761       | 42.3778       |
| total_num_scaffolds    | 2,757            | 1,215                      | 1,083                         | 2,781            | 1,245                      | 1,115                         | 2,754            | 1,205                      | 1,073                         | 1,261            | 798                        | 794                           | 1,311            | 830                        | 829                           | 1,288            | 793                        | 786                           | 972              | 616                        | 616                           | 1,005            | 650                        | 650                           | 978              | 597                        | 597                           | 575           | 490           |
| total_num_contigs      | 2,757            | 1,215                      | 1,083                         | 2,781            | 1,245                      | 1,115                         | 2,754            | 1,205                      | 1,073                         | 1,261            | 798                        | 794                           | 1,311            | 830                        | 829                           | 1,288            | 793                        | 786                           | 972              | 616                        | 616                           | 1,005            | 650                        | 650                           | 978              | 597                        | 597                           | 590           | 501           |
| number_of_gaps         | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 0                | 0                          | 0                             | 15            | 11            |
| Longest_scaffold       | 46,214,770       | 46,214,770                 | 46,214,770                    | 46,214,784       | 46,214,784                 | 46,214,784                    | 46,214,776       | 46,214,776                 | 46,214,776                    | 50,262,241       | 50,262,241                 | 50,262,241                    | 50,372,431       | 50,372,431                 | 50,372,431                    | 50,428,791       | 50,428,791                 | 50,428,791                    | 64,130,911       | 64,130,911                 | 64,130,911                    | 65,868,658       | 65,868,658                 | 65,868,658                    | 69,566,365       | 69,566,365                 | 69,566,365                    | 33,552,962    | 33,552,962    |
| Scaffold_N50_length    | 2,955,405        | 7,258,968                  | 8,379,364                     | 2,998,866        | 7,383,908                  | 8,477,060                     | 2,998,866        | 7,258,958                  | 8,379,360                     | 10,709,157       | 14,151,982                 | 14,151,982                    | 11,908,915       | 14,446,381                 | 14,446,381                    | 12,368,032       | 14,424,108                 | 15,226,941                    | 17,289,722       | 17,435,227                 | 17,435,227                    | 17,291,129       | 17,435,228                 | 17,435,228                    | 17,435,235       | 21,431,074                 | 21,431,074                    | 8,893,208     | 8,893,208     |
| Scaffold_NG50_length   | 10,045,395       | 10,045,395                 | 10,045,395                    | 10,097,027       | 10,072,371                 | 9,722,664                     | 10,111,749       | 10,111,749                 | 10,045,389                    | 15,882,063       | 15,882,063                 | 15,882,063                    | 16,399,956       | 16,399,956                 | 16,399,956                    | 16,399,925       | 16,399,925                 | 16,399,925                    | 21,425,933       | 21,425,933                 | 21,425,933                    | 21,431,078       | 21,431,078                 | 21,431,078                    | 22,844,882       | 22,844,882                 | 22,844,882                    | 8,901,910     | 8,901,910     |
| Scaffold_N95_length    | 232,223          | 362,570                    | 374,275                       | 228,491          | 340,168                    | 317,809                       | 231,040          | 364,671                    | 364,465                       | 312,400          | 596,648                    | 596,648                       | 276,936          | 542,153                    | 542,153                       | 299,635          | 557,635                    | 559,069                       | 445,326          | 891,889                    | 891,889                       | 407,217          | 733,833                    | 733,833                       | 436,642          | 879,043                    | 879,043                       | 533,164       | 660,386       |
| Scaffold_NG95_length   | 3,485,820        | 2,473,102                  | 1,574,909                     | 3,496,061        | 2,473,082                  | 1,604,414                     | 3,554,867        | 2,531,255                  | 1,574,891                     | 2,730,996        | 2,153,208                  | 2,153,208                     | 2,698,482        | 2,219,116                  | 2,219,116                     | 2,698,460        | 2,219,116                  | 2,219,116                     | 3,399,398        | 2,986,878                  | 2,986,878                     | 3,320,541        | 3,236,925                  | 3,236,925                     | 3,386,191        | 3,231,517                  | 3,231,517                     | 868,503       | 868,503       |
| Longest_contig         | 46,214,770       | 46,214,770                 | 46,214,770                    | 46,214,784       | 46,214,784                 | 46,214,784                    | 46,214,776       | 46,214,776                 | 46,214,776                    | 50,262,241       | 50,262,241                 | 50,262,241                    | 50,372,431       | 50,372,431                 | 50,372,431                    | 50,428,791       | 50,428,791                 | 50,428,791                    | 64,130,911       | 64,130,911                 | 64,130,911                    | 65,868,658       | 65,868,658                 | 65,868,658                    | 69,566,365       | 69,566,365                 | 69,566,365                    | 33,552,962    | 33,552,962    |
| Contig_N50_length      | 2,955,405        | 7,258,968                  | 8,379,364                     | 2,998,866        | 7,383,908                  | 8,477,060                     | 2,998,866        | 7,258,958                  | 8,379,360                     | 10,709,157       | 14,151,982                 | 14,151,982                    | 11,908,915       | 14,446,381                 | 14,446,381                    | 12,368,032       | 14,424,108                 | 15,226,941                    | 17,289,722       | 17,435,227                 | 17,435,227                    | 17,291,129       | 17,435,228                 | 17,435,228                    | 17,435,235       | 21,431,074                 | 21,431,074                    | 8,893,208     | 8,893,208     |
| Contig_NG50_length     | 10,045,395       | 10,045,395                 | 10,045,395                    | 10,097,027       | 10,072,371                 | 9,722,664                     | 10,111,749       | 10,111,749                 | 10,045,389                    | 15,882,063       | 15,882,063                 | 15,882,063                    | 16,399,956       | 16,399,956                 | 16,399,956                    | 16,399,925       | 16,399,925                 | 16,399,925                    | 21,425,933       | 21,425,933                 | 21,425,933                    | 21,431,078       | 21,431,078                 | 21,431,078                    | 22,844,882       | 22,844,882                 | 22,844,882                    | 8,901,910     | 8,901,910     |
| Contig_N95_length      | 232,223          | 362,570                    | 374,275                       | 228,491          | 340,168                    | 317,809                       | 231,040          | 364,671                    | 364,465                       | 312,400          | 596,648                    | 596,648                       | 276,936          | 542,153                    | 542,153                       | 299,635          | 557,635                    | 559,069                       | 445,326          | 891,889                    | 891,889                       | 407,217          | 733,833                    | 733,833                       | 436,642          | 879,043                    | 879,043                       | 527,794       | 625,785       |
| Contig_NG95_length     | 3,485,820        | 2,473,102                  | 1,574,909                     | 3,496,061        | 2,473,082                  | 1,604,414                     | 3,554,867        | 2,531,255                  | 1,574,891                     | 2,730,996        | 2,153,208                  | 2,153,208                     | 2,698,482        | 2,219,116                  | 2,219,116                     | 2,698,460        | 2,219,116                  | 2,219,116                     | 3,399,398        | 2,986,878                  | 2,986,878                     | 3,320,541        | 3,236,925                  | 3,236,925                     | 3,386,191        | 3,231,517                  | 3,231,517                     | 856,071       | 856,071       |
| qv_score               | 61.9987          | 62.0147                    | 61.6227                       | 61.945           | 62.0304                    | 61.624                        | 62.0604          | 62.1082                    | 61.6803                       | 60.8608          | 61.6912                    | 61.6909                       | 60.8165          | 61.7131                    | 61.7129                       | 60.8343          | 61.617                     | 61.635                        | 60.5357          | 61.342                     | 61.342                        | 60.389           | 61.3572                    | 61.3572                       | 60.4856          | 61.3088                    | 61.3088                       | 56.3843       | 57.2844       |
| kmer_completeness      | 99.3468          | 85.5233                    | 82.6667                       | 99.3354          | 85.4438                    | 82.561                        | 99.3428          | 85.6085                    | 82.5495                       | 84.2954          | 82.507                     | 82.5067                       | 84.3931          | 82.6174                    | 82.6171                       | 84.4225          | 82.6217                    | 82.6059                       | 83.2772          | 82.4928                    | 82.4928                       | 83.3569          | 82.6469                    | 82.6469                       | 83.3613          | 82.5778                    | 82.5778                       | 81.3829       | 81.2586       |
| Complete_Buscos        | 96.2%            | 95.8%                      | 95.8%                         | 96.3%            | 95.8%                      | 95.9%                         | 96.3%            | 95.8%                      | 95.7%                         | 96.0%            | 95.8%                      | 95.8%                         | 96.0%            | 95.9%                      | 95.9%                         | 96.0%            | 95.8%                      | 95.8%                         | 96.0%            | 95.9%                      | 95.9%                         | 96.0%            | 95.9%                      | 95.9%                         | 96.1%            | 96.0%                      | 96.0%                         | 95.6%         | 95.7%         |
| Complete_single_buscos | 16.8%            | 78.6%                      | 94.8%                         | 16.9%            | 78.7%                      | 94.7%                         | 16.8%            | 78.4%                      | 94.7%                         | 87.1%            | 94.7%                      | 94.7%                         | 87.5%            | 94.8%                      | 94.8%                         | 87.2%            | 94.8%                      | 94.8%                         | 93.6%            | 94.9%                      | 94.9%                         | 93.7%            | 94.9%                      | 94.9%                         | 93.8%            | 95.0%                      | 95.0%                         | 94.6%         | 94.7%         |

PDF report for each asm evaluated [here](https://drive.google.com/file/d/1emPu1VEXbFs84HnUCY7okc1EiCcOWlEU/view?usp=sharing) 
*Achtung: the K-mer Multiplicity PRI (stacked) plot is wrong - the plot is actually showing the flattened pri+alt. The correct plot is in the merquey folder: asmnamemerqOutput.spectra-cn.st.png*

* Notes on [merqury](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02134-9) results:

-**Merqury copy number spectrum plots** (similar to KAT): 
tracks the multiplicity of each k-mer found in the reads and colors it by the number of times it is found in a given assembly 
->when using hifiasm.l0, the plots for the primmary asm have also a blue peak (denoting kmers appearing twice), which should have not been there since the asm is already phased. For hifiasm.l1 too, but after purging the blue is gone.
->no asms showed a bar at the origin of the plot representing kmers found only in the asm
-> the higher fraction of k-mers missing in the assembly (black peak, read-only) seen in all asm plots probably reflects kmers specific to the other haplotig! This is corroborated when looking at the plot flattened pri+alt, in which the grey peak of read-only kmers disappears. 

-**kmer completeness**: first filter out low-copy k-mers that are likely to represent sequencing errors, and then define k-mer completeness as the fraction of reliable k-mers in the read set that are also found in the assembly.

-**qv-score**:(Phred quality score or kmer survival rate) assumes that kmers found only in the asm are errors. It depends on the coverage,  quality of the read set, and kmer size chosen. Higher QVs indicate a more accurate consensus.

[FullTable_asms_comparisons.xlsx](https://docs.google.com/spreadsheets/d/1RD4zSMCuin1e0ovnP5k5UIenwJLs5cFM/edit?usp=sharing&ouid=108013701707137496295&rtpof=true&sd=true)

* Criteria for "good" assembly:


| Criterion | Best assembly | Best value |
| -------- | -------- | -------- |
| Contig N50>10Mb     | **hifiasm.l2.k31.purged** (both with default or VGP suggested values)     |   21 Mb  |
| Busco completness>95%     | all, but best: hifiasm.l0 (before purging), and **hifiasm.l2.k31.purged**  (both with default or VGP suggested values)  |   96%  |
| Busco duplicated rate<3%     | **all**, except hifiasm.l0 and l1 before purging, and l0 after purging with default values     |   1%  |
| < Number of contigs    | IPA, IPApurged, **hifiasm.l2.k31.purged** (default and 
| > qv-score    |  hifiasm.l0.k31.purged.default     |   62.11  |
| > kmer completness    |  hifiasm.l0.k31     |   99.34  |
| copy number spectrum plots    |  all hifiasm.l1.purged, hifiasm.l2 and all **hifiasm.l2.purged**     |   no blue peak, no bar of kmers found only in the asm  |

copy number spectrum plots

## Best assembly
All the QC metrics are better for asm **hifiasm.l2.k31.purged** except the qv-score and kmer completness.

![](https://i.imgur.com/ZhLbjR1.png)
![](https://i.imgur.com/dNhtrng.jpg)



*****
### Other assemblers I tried: Hicanu
 https://github.com/marbl/canu/releases

```
module purge
module load lang/Anaconda3/2020.11
source activate java11-env
export PATH=/work/dominguez/YelCar_GenomeAssembly_MPI/canu-2.1.1/bin:$PATH

canu -assemble -d /path/hicanu_FatNode/ -p YelCar_hicanu_assemble genomeSize=1.06g -pacbio-hifi /path_ro_reads/reads_noAdapt.fastq useGrid=false

#change genomeSize for what was predicted by genomescope2
#Requested --mem=150G, --cpus-per-task=8, and took xx to finish.
```

*****
### Other assemblers I tried: LJA
https://github.com/AntonBankevich/LJA
```

module purge
module load lib/zlib/1.2.11-GCCcore-10.2.0
module load lang/Anaconda3
source activate cmake-env

cd /path/LJA/
/work/dominguez/LJA/bin/lja -o /path/LJA_fatnode_Cintel/ --reads /path_to_reads/reads_noAdapt.fastq --threads 8 --diploid
```
The job to assemble the hifi reads seems not to use much memory over most of its runtime. It seems that it is only at the end that it performs a huge allocation. And it is killed because of that.  So, I then ran it in our "fat" node (where I can request --mem=1500G). In this case, the job failed after 14 days of running due to a time limitation. I, then, tested it with a shorter version of the infile, requesting 24 CPUs (--mem-per-cpu=16G). It took 1h and then also failed in the step of uncompressing homopolymers in contigs. It produced the expected output files except from the final fasta assembly and final gfa. 
 
*when fatnode is free, try: submitLJAtrimmedFatNodeintelv2.slurm*
 
*****
### Other assemblers: Rust-mdbg
https://github.com/ekimb/rust-mdbg
I did not try it given that it does not resolve both haplotypes.
*****

### MUMMER
https://github.com/mummer4/mummer
https://mummer4.github.io/manual/manual.html

**Compare my best asm  with a reference from a related species genome.**

There are no VGP genomes in the same family (Thraupidae).
I could use another family, like:
Zebra finch genome VGP([GCF_003957565.2_bTaeGut1.4.pri](https://www.ncbi.nlm.nih.gov/assembly/GCF_003957565.2/))...
or a species of the same family (but not so good asm), like  **small tree finch **([GCA_902806625.1_Camarhynchus_parvulus_V1.1](https://https://www.ncbi.nlm.nih.gov/assembly/GCA_902806625.1/)), and medium ground finch ([GCF_000277835.1_GeoFor_1.0](https://www.ncbi.nlm.nih.gov/assembly/GCF_000277835.1/)).

The genome asm for the small tree finch (Camarhynchus parvulus) is slightly better (genome coverage 31X, less gaps).


 all-vs-all comparison of nucleotide sequences contained in multi-Fasta data files. It is best used for similar sequences that may have large rearrangements.

download Camarhynchus parvulus (Campar) asm from ncbi (in terminal):
```
wget -O Campar_genome1.1.fna.gz ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/902/806/625/GCA_902806625.1_Camarhynchus_parvulus_V1.1/GCA_902806625.1_Camarhynchus_parvulus_V1.1_genomic.fna.gz
```

job script: [submit_mummer_Hifiasm.l2.k31.p.PRI_vs_Campar.slurm](https://drive.google.com/file/d/1ABDOC2OYKV0spl6pCpQ-FVUueXqVAAML/view?usp=sharing) (Ran with --cpus-per-task=24, --mem=64G, and took 3hs).
```
module purge
module load lang/Anaconda3/2020.11
source activate mummer-env

#unzip it: 
gunzip /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Campar_genome1.1.fna.gz

#run mummer:
nucmer -p mummer.l2.k31.p.Pri_Campar /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Campar_genome1.1.fna /work/dominguez/YelCar_GenomeAssembly_MPI/asm/reads_trimmed/hifiasm/purge_dups/c1.l2.k31/purged_default.fa
```
*I used NUCmer (NUCleotide MUMmer) version 3.1*

Dotplot between two assemblies


PartB
[submit_mummer_Hifiasm.l2.k31.p.PRI_vs_Campar_2.slurm
](https://drive.google.com/file/d/1oXEq6G959TbFwNN3hhswXHYqfi5d4JVZ/view?usp=sharing)
```
#run delta filter:
delta-filter -i 89 -l 10000 /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/mummer.l2.k31.p.Pri_Campar.delta > /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/mummer2.l2.k31.p.Pri_Campar.filtered.delta

#Run show-coords:
show-coords /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/mummer2.l2.k31.p.Pri_Campar.filtered.delta> /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/mummer2.l2.k31.p.Pri_Campar.filtered.delta.coords

mummerplot --large --filter --fat /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/mummer2.l2.k31.p.Pri_Campar.filtered.delta -p mummer2.l2.k31.p.Pri_CamPar.filtered -t png
```
This produced a .rplot, gp, and a fplot files in /work/dominguez/. I mv them to /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/

It fails producing the png ("Could not find gnuplot, plot will not be rendered"). 

I installed gnuplot with conda.
```
gnuplot /work/dominguez/YelCar_GenomeAssembly_MPI/Camarhynchus_parvulus_Asm/Mummer_filtered10000/mummer2.l2.k31.p.Pri_CamPar.filtered.gp > mummer2.l2.k31.p.Pri_CamPar.filtered.png
```
![](https://i.imgur.com/Af1UcEx.png)
Change in the .gp file towards the bottom this lines:
```
set linestyle 1  lt 1 lw 0.3 pt 6 ps 0.1
set linestyle 2  lt 3 lw 0.3 pt 6 ps 0.1
set linestyle 3  lt 2 lw 0.3 pt 6 ps 0.1
```
![](https://i.imgur.com/yyB33BV.png)

I can do a circos plot with the mummer .coords file.

### NETX STEPS:
1. assembly mitochondrial genome
2. Repeats annotation
3. Genes Annotation
4. Gene family expansion/constraction analysis considering other species (CAFE)


## Under construction:



## Contamination check ##
#https://kat.readthedocs.io/en/latest/
module purge
module load lang/Anaconda3/2020.11
source activate kat-env2

kat gcp -o /path/reads_statistics/filename_ccs_fasta.kat -t 40 -p png /path/reads/filename_ccs_reads.fa 


#Blobtools?
#Bandage (before purging, but not possible after purging because there is no gfa)


## Tutorials / resources
 
kmer evaluation: https://github.com/marbl/merqury/wiki/2.-Overall-k-mer-evaluation


https://eukaryotic-genome-assembly.github.io/handsOn_1/

https://training.galaxyproject.org/training-material/topics/assembly/tutorials/vgp_genome_assembly/tutorial.html#post-assembly-processing

https://ucdavis-bioinformatics-training.github.io/2020-Genome_Assembly_Workshop/ONT_Assembly/nanopore

https://canu.readthedocs.io/en/latest/tutorial.html#outputs


https://github.com/dfguan/purge_dups

