Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

Code Block
name: flye
channels:
 - defaults
 - anaconda
 - bioconda
 - conda-forge
dependencies:
 - python=3.7
 - flye=2.9.1
 - seqkit=2.3.1
 - blast=2.13.0
 - nanoq=0.9.0
 - minimap2
 - samtools

Run the following command to generate the ‘flye’ conda environment:

...

Code Block
#!/bin/bash -l
#PBS -N FlyeAssembly
#PBS -l walltime=24:00:00
#PBS -l mem=16gb
#PBS -l ncpus=8

## Usage: qsub launch_flye_genome_assembly.pbs

cd $PBS_O_WORKDIR

################################################################################################################################
# USER DEFINED VARIABLES
################################################################################################################################
SAMPLEID=MyGenome
REFNAME=Accession_RefGenome
REF='/path/to/REF/genome.fasta'
ONT='/path/to/ONT/reads.fastq.gz'
################################################################################################################################

#activate flye conda environment 
conda activate flye

#ASSEMBLY

#STEP1: Run de novo genome assembly for >=Q20 data, use a combination of --nano-hq and --read-error 0.03
flye  --out-dir outdir_nano --threads 8 --read-error 0.03 --nano-hq $ONT  

#STEP2: Compare sequence similarity of assembled genome with reference sequence
blastn -query outdir_nano/assembly.fasta -subject $REF -evalue 1e-5 -out blastn_${REFNAME}_vs_${SAMPLEID}_assembly.txt -outfmt '6 qseqid sacc length pident mismatch gapopen qstart qend qlen sstart send slen evalue bitscore qcovhsp qcovs'

#STEP3: format output BLASTN table
echo "qseqid sacc length pident mismatch gapopen qstart qend qlen sstart send slen evalue bitscore qcovhsp qcovs" > header

#add header to blast output
cat header blastn_${REFNAME}_vs_${SAMPLEID}_assembly.txt > BLASTN_${REFNAME}_vs_${SAMPLEID}_assembly.txt 

#remove intermediate files
rm header blastn_${REFNAME}_vs_${SAMPLEID}_assembly.txt

#MAPPING

#generic mapping reads
minimap2 -a $REF $ONT > ${SAMPLEID}_aln.sam

#mapping noisy reads
#minimap2 -ax $REF $ONT > ${SAMPLEID}_aln2.sam  

#Samtools
samtools view -bt ${LIST} -o ${SAMPLEID}_aln.bam ${SAMPLEID}_aln.sam

samtools sort -T /tmp/aln.sorted -o ${SAMPLEID}_aln.sorted.bam ${SAMPLEID}_aln.bam

samtools index ${SAMPLEID}_aln.sorted.bam

submit the job

Code Block
qsub launch_flye_genome_assembly.pbs

...