RosaBioinfo

Computational pipeline for treatment of long reads, including masking host reads and assembly of a Wolbachia genome

Blandine Trouche, Alice Brunner, Julie Reveillaud 2025-12-08

DWK=/Users/btrouche/Desktop/Bioinformatique/RosaBioinfo/Celg3_Long_reads

0- Basecalling of raw fast5 data into fastq

This script was run on i-TROP by doing sbatch script.sh

#!/bin/bash

##### Configuration slurm #####
# define job name
#SBATCH --job-name=guppy_basecall_SUP
# define target partition
#SBATCH -p gpu
# gpu specific options
#SBATCH -A gpu_group
#SBATCH --gres=gpu:1
#SBATCH -c 8
# email alert: where and when 
#SBATCH --mail-user=blandine.trouche@ird.fr
#SBATCH --mail-type=ALL
###############################

# create working directory in /scratch
cd /scratch 
mkdir trouche-$SLURM_JOB_ID
cd trouche-$SLURM_JOB_ID

# get data from san
scp -r san:/projects/large/RosaLind/Long_Read_Sequencing/Cel_g3_20240626/20230302_1616_MC-113479_FAU44798_43dcecde/fast5_pass /scratch/t
rouche-$SLURM_JOB_ID

# load software to use
module load guppy-gpu/6.5.7

# analysis: running basecalling 
guppy_basecaller -c dna_r9.4.1_450bps_sup.cfg -i fast5_pass --recursive \
      -s SUP_basecall_no_adapters --num_callers 8 --gpu_runners_per_device 8 \
      --device auto --min_qscore 7 --compress_fastq \
      --detect_adapter --trim_adapters

# transfer results back to nas
scp -r SUP_basecall_no_adapters san:/projects/large/RosaLind/Long_Read_Sequencing/Cel_g3_20240626/20230302_1616_MC-113479_FAU44798_43dce
cde

# remove data on scratch (after transfer)
cd /scratch
rm -rf trouche-$SLURM_JOB_ID

1- Concatenating fastqs and quality plots

Prepare data by concatenating fastqs

# prepare data for analysis
cd $DWK

mkdir -p 00_RAW_DATA

# unzip
gunzip 00_RAW_DATA/SUP_basecall_no_adapters/pass/*.gz

# concatenate all files in 1
cat 00_RAW_DATA/SUP_basecall_no_adapters/pass/*.fastq > 00_RAW_DATA/Cel_g3_SUP.fastq

# rezip data for space
gzip 00_RAW_DATA/SUP_basecall_no_adapters/pass/*.fastq
gzip 00_RAW_DATA/Cel_g3_SUP.fastq

Produce summary of quality with Nanoplot

conda activate nanoplot

NanoPlot -o 00_RAW_DATA/Nanoplot_out --tsv_stats \
          --fastq 00_RAW_DATA/Cel_g3_SUP.fastq.gz \
          --loglength

2- Preparing the data: masking Culex genome, removing lambda phage DNA, Culex DNA

Masking Culex pipiens molestus genome with wPipMol

cd $DWK

# Symlink to reference data
ln -s /Users/btrouche/Desktop/Bioinformatique/00_Reference_wPip/01_FASTA 02_FASTA_REF
ln -s /Users/btrouche/Desktop/Bioinformatique/00_Reference_wPip/02_CONTIGS 03_CONTIGS_REF

Masking done on Migale, with script /work_projet/culex/Blandine/01_Metagenomics/03_Cxq/scripts/S1b_mask_Cxm.sh

#!/bin/bash
source /usr/local/genome/Anaconda3/etc/profile.d/conda.sh

# Shell to use for job execution
#$ -S /bin/bash

# Job name
#$ -N bbmask_Cxm

# Nom de la queue
#$ -q short.q

# Export all environment variables
#$ -V

# Run command from current repertory
#$ -cwd

# CPUs number
#$ -pe thread 6

# Logs
#$ -o /work_projet/culex/Blandine/01_Metagenomics/03_Cxq/scripts/logs 
#$ -e /work_projet/culex/Blandine/01_Metagenomics/03_Cxq/scripts/logs

# activate conda environment
conda activate bbmap-38.93

ref_folder=/work_projet/culex/Blandine/00_Ref_data
wPipMol_folder=$ref_folder/01_FASTA/wPipMol

# make output folder
mkdir -p $ref_folder/Cxm_masked

# unzip ref genomes 
gzip $ref_folder/CulexMol_GCA_024516115.1.fna
gunzip $wPipMol_folder/wPipMol-contigs-prefix-formatted-only.fa.gz

# use bbmask to mask the Culex molestus genome and avoid future problems when removing it by mapping
# 1) prepare wPipMol for mapping
shred.sh in=$wPipMol_folder/wPipMol-contigs-prefix-formatted-only.fa \
        out=$ref_folder/Cxm_masked/wPipMol-contigs_shredded.fa \
        length=80 minlength=70 overlap=40

# 2) map Culex molestus to wPipMol to identify similar sequences
bbmap.sh ref=$ref_folder/CulexMol_GCA_024516115.1.fna.gz \
        in=$ref_folder/Cxm_masked/wPipMol-contigs_shredded.fa \
        outm=$ref_folder/Cxm_masked/Cxm_to_wPipMol_mapped.sam \
        minid=0.85 maxindel=2 \
        threads=6

# 3) mask similar sequences, zones with low entropy and tandem-repeated kmers
bbmask.sh in=$ref_folder/CulexMol_GCA_024516115.1.fna.gz \
        out=$ref_folder/Cxm_masked/Cxm_masked.fa \
        entropy=0.7 \
        sam=$ref_folder/Cxm_masked/Cxm_to_wPipMol_mapped.sam \
        threads=6

gzip $ref_folder/Cxm_masked/*.fa
gzip $wPipMol_folder/*.fa

# deactivate conda environment
conda deactivate

Sort phage reads by mapping them on the reference lambda genome

# mapping of reads to find lambda phage reads
# this NEEDS the option to not output unaligned reads, otherwise the following filter removes everything
conda activate pomoxis

mkdir -p 04_MAPPING

minimap2 -ax map-ont -a /Users/btrouche/Desktop/Bioinformatique/Long_Read_Sequencing/MinION_seq_runs/00_DATA/Lambda_phage_ref.fasta.gz \
        --sam-hit-only \
        00_RAW_DATA/Cel_g3_SUP.fastq.gz \
        > 04_MAPPING/Cel_g3_lambda_phage_alignment.sam

# sorting lambda phage reads from mosquito + Wolbachia reads
conda activate anvio-8
samtools view 04_MAPPING/Cel_g3_lambda_phage_alignment.sam | cut -f 1 > ids_to_remove_lambda.txt

iu-remove-ids-from-fastq -i 00_RAW_DATA/Cel_g3_SUP.fastq.gz -l ids_to_remove_lambda.txt -d " "
# 32,230 reads removed out of 309,559
# ~90% of reads would be from the organism!!

mv 00_RAW_DATA/Cel_g3_SUP.fastq.gz.survived 00_RAW_DATA/Cel_g3_SUP.target.fastq
mv 00_RAW_DATA/Cel_g3_SUP.fastq.gz.removed 00_RAW_DATA/Cel_g3_SUP.phage.fastq

gzip 00_RAW_DATA/*.fastq

Sort Culex reads by mapping them on the reference Culex pipiens molestus genome

Here I will use the genome of Culex Molestus masked against wPipMol

cd $DWK

# mapping
conda activate pomoxis
minimap2 -ax map-ont -a /Users/btrouche/Desktop/Bioinformatique/00_Ref_data/Cxm_masked/Cxm_masked.fa.gz \
          --sam-hit-only \
          00_RAW_DATA/Cel_g3_SUP.target.fastq.gz > 04_MAPPING/Cel_g3_CulexMol_alignment.sam

# separating Culex and bacterial reads
conda activate anvio-8 
samtools view 04_MAPPING/Cel_g3_CulexMol_alignment.sam | cut -f 1 > ids_to_remove_Culex.txt

iu-remove-ids-from-fastq -i 00_RAW_DATA/Cel_g3_SUP.target.fastq.gz -l ids_to_remove_Culex.txt -d " "
# processed: 277,329 reads, removed: 229,608 reads

mv 00_RAW_DATA/Cel_g3_SUP.target.fastq.gz.survived 00_RAW_DATA/Cel_g3_SUP_noCulex.fastq
rm 00_RAW_DATA/Cel_g3_SUP.target.fastq.gz.removed
gzip 00_RAW_DATA/*.fastq

3- Assembly without lambda phage reads and wihtout Culex reads

cd $DWK

conda activate flye

mkdir -p 02_ASSEMBLY_onlybac

# assembly of long reads : Flye (for metagenomes?)
flye --nano-raw 00_RAW_DATA/Cel_g3_SUP_noCulex.fastq.gz \
     --out-dir 02_ASSEMBLY_onlybac --meta 
# Raw reads N50/N90: 3614 / 909

4- Create associated Anvi’o databases with only bac assembly

cd $DWK

conda activate anvio-8

mkdir -p 03_CONTIGS
mkdir -p 00_LOGS

anvi-gen-contigs-database -f 02_ASSEMBLY_onlybac/assembly.fasta \
              -o 03_CONTIGS/onlybac_assembly-contigs.db \
              --project-name Celg3_onlybac >> 00_LOGS/assembly_only_bac-anvi_gen_contigs_database.log 2>&1

anvi-run-hmms -c 03_CONTIGS/onlybac_assembly-contigs.db >> 00_LOGS/assembly_only_bac-anvi-run-hmms.log 2>&1

anvi-run-scg-taxonomy -c 03_CONTIGS/onlybac_assembly-contigs.db \
      >> 00_LOGS/assembly_only_bac-anvi_run_scg_taxonomy.log 2>&1

####
anvi-estimate-scg-taxonomy -c 03_CONTIGS/onlybac_assembly-contigs.db \
                           --metagenome-mode

anvi-get-sequences-for-hmm-hits -c 03_CONTIGS/onlybac_assembly-contigs.db --hmm-source Bacteria_71
# mostly in contig 2

anvi-run-kegg-kofams -c 03_CONTIGS/onlybac_assembly-contigs.db   

anvi-run-ncbi-cogs -c 03_CONTIGS/onlybac_assembly-contigs.db

# check taxonomy of contigs with centrifuge : 
anvi-get-sequences-for-gene-calls -c 03_CONTIGS/onlybac_assembly-contigs.db \
                                  -o 03_CONTIGS/only_bac-gene_calls.fa \
                                  >> 00_LOGS/export_gene_calls_only_bac.log 2>&1

# second step has to be done on a cluster with centrifuge installed (and more particularly its database)
# first send your data on the server (scp or filezilla)
# then make a script with the following command
centrifuge -f -x /home/shared/db/centrifuge_db/ncbi_nr/nt \
              only_bac-gene_calls.fa \
              -S out_centrifuge-only_bac.tsv \
              --report-file centrifuge-only_bac.tsv \
              --threads 3 >> 00_LOGS/centrifuge-only_bac.log 2>&1

# get the results back locally (scp or filezilla)
anvi-import-taxonomy-for-genes -c 03_CONTIGS/onlybac_assembly-contigs.db \
                                -i 03_CONTIGS/centrifuge-only_bac.tsv \
                                03_CONTIGS/out_centrifuge-only_bac.tsv \
                                -p centrifuge >> 00_LOGS/anvi-import-taxonomy-for-genes_centrifuge_only_bac.log 2>&1

mkdir -p 04_MAPPING/only_bac

# mapping raw reads on assembled contigs
conda activate pomoxis
minimap2 -ax map-ont -a 02_ASSEMBLY_onlybac/assembly.fasta \
                      00_RAW_DATA/Cel_g3_SUP.fastq.gz \
                      --sam-hit-only > 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned.sam
                      
conda activate anvio-8                   
samtools view -bS 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned.sam \
              -o 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned-RAW.bam  

anvi-init-bam 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned-RAW.bam \
              -o 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned.bam  

anvi-profile -i 04_MAPPING/only_bac/Celg3_SUP_contigs_aligned.bam \
              -c 03_CONTIGS/onlybac_assembly-contigs.db --cluster-contigs  \
              -o 05_ANVIO_PROFILE/only_bac --sample-name Celg3  \
              --overwrite-output-destinations --min-contig-length 0 >> 00_LOGS/anvi-profile_only_bac.log 2>&1

anvi-interactive -p 05_ANVIO_PROFILE/only_bac/PROFILE.db -c 03_CONTIGS/onlybac_assembly-contigs.db

anvi-summarize -p 05_ANVIO_PROFILE/only_bac/PROFILE.db \
              -c 03_CONTIGS/onlybac_assembly-contigs.db \
               -C default \
               -o 06_SUMMARY/only_bac

5- Proovframe

Apply proovframe and check how it affects completion and redundancy of MAG. FYI, with super-accuracy basecalled reads, we do not seem to run as much into the problem of broken ORFs and overestimated redundancy.

mkdir -p 07_PROOVFRAME

DIAMOND_db=/Users/btrouche/opt/miniconda3/envs/anvio-7.1/lib/python3.6/site-packages/anvio/data/misc/COG/COG20/DB_DIAMOND/COG.dmnd

conda activate anvio-8

### proovframe on target assembly (phage reads filtered out + binning done)
/Users/btrouche/Desktop/proovframe/bin/proovframe map -d $DIAMOND_db \
                                                      -o 07_PROOVFRAME/Wolbachia_onlybac_SUP_corrected.tsv \
                                                      06_SUMMARY/only_bac/bin_by_bin/Bin_1/Bin_1-contigs.fa

/Users/btrouche/Desktop/proovframe/bin/proovframe fix -o 07_PROOVFRAME/Wolbachia_onlybac_SUP_corrected.fa \
                                                      06_SUMMARY/only_bac/bin_by_bin/Bin_1/Bin_1-contigs.fa \
                                                      07_PROOVFRAME/Wolbachia_onlybac_SUP_corrected.tsv 

# get completeness
anvi-gen-contigs-database -f 07_PROOVFRAME/Wolbachia_onlybac_SUP_corrected.fa \
              -o 07_PROOVFRAME/Wolb_onlybac_corrected-contigs.db \
              --project-name Celg3_onlybac

anvi-run-hmms -c 07_PROOVFRAME/Wolb_onlybac_corrected-contigs.db

anvi-estimate-genome-completeness -c 07_PROOVFRAME/Wolb_onlybac_corrected-contigs.db

6- Visualizing genome and mapping in IGV

mkdir -p 07_PROOVFRAME/04_MAPPING

# mapping raw reads on assembled contigs
conda activate pomoxis
minimap2 -ax map-ont -a 07_PROOVFRAME/Wolbachia_onlybac_SUP_corrected.fa \
                      00_RAW_DATA/Cel_g3_SUP.fastq.gz \
                      --sam-hit-only > 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected.sam
                      
conda activate anvio-8    

samtools view -bS 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected.sam \
              -o 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected-RAW.bam  

anvi-init-bam 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected-RAW.bam  \
              -o 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected.bam   

anvi-get-sequences-for-gene-calls -c 07_PROOVFRAME/Wolb_onlybac_corrected-contigs.db \
            -o 07_PROOVFRAME/Wolb_onlybac_corrected_genecalls.gff --export-gff3 \
            --annotation-source COG20_FUNCTION
            
anvi-run-ncbi-cogs -c 07_PROOVFRAME/Wolb_onlybac_corrected-contigs.db

Extracting reads of interest for blasting

conda activate samtools
samtools view 07_PROOVFRAME/04_MAPPING/Celg3_SUP-Wolb_onlybac_corrected.bam | grep -w "dee6845e-22e6-40a2-b5cb-317536e63ab6"