RosaBioinfo

Extraction of reads >9,000 nt aligning to the plasmid and verification ofadapter removal

Alice Brunner, Maxime Mahout 2025-01-30

Extraction of high quality reads of each contig

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

Creation of the .bam file with just high quality reads

samtools view -bq 2 Celg3_SUP-Wolb_onlybac_corrected.sam -o Celg3_SUP-Wolb_onlybac_quality.bam
samtools sort Celg3_SUP-Wolb_onlybac_quality.bam > Celg3_SUP-Wolb_onlybac_quality.bam-sorted.bam
samtools index Celg3_SUP-Wolb_onlybac_quality.bam-sorted.bam

The .bam file is located in /sequençage_LR_cel/seq_cel_3/15_anvio_quality/Celg3_SUP-Wolb_onlybac_quality.bam.

Creating of new profile.db

Creating of new profile.db from the new quality-filtered bam file

anvi-profile -i Celg3_SUP-Wolb_onlybac_quality.bam-sorted.bam -c onlybac_assembly-contigs.db --skip-SNV-profiling

# or, if it doesn't work : 

anvi-profile -i Celg3_SUP-Wolb_onlybac_quality.bam-sorted.bam -c onlybac_assembly-contigs.db --skip-SNV-profiling --cluster-contigs

Anvio’s Visualisation

anvi-interactive -c onlybac_assembly-contigs.db -p Celg3_SUP_Wolb_onlybac_quality_bam_sorted/PROFILE.db

# Observation of the mean_coverage of each contig interactively

Analyzing dimeric plasmid reads

Extracting long reads mapping to the contig 9

The following script, which we applied to all contigs, generate in particular a file called mapping_to_contig9_filtered.bam, corresponding to the plasmid.

for i in $(echo "1 2 3 5 6 7 8 9 10 13 14 15 16 17 18 19 20 21 22 23 24 26 27 38" | xargs); do 
    # extract reads mapping on plasmid contig (contig_${i}?)
    samtools view ../Reads_Alice/Celg3_SUP-Wolb_onlybac_corrected.bam "contig_${i}" -o mapping_to_contig${i}.bam

    # filter low quality mapping if you want        
    samtools view -bq 2 mapping_to_contig${i}.bam -o mapping_to_contig${i}_filtered.bam

    # get the read ids and get them from the original fastq
    mkdir -p 15_extracting_reads_qualfilter

    samtools view mapping_to_contig${i}_filtered.bam | cut -f 1 > 15_extracting_reads_qualfilter/ids_of_mapped_reads_contig${i}.txt

    iu-remove-ids-from-fastq -i ../Reads_Alice/Cel_g3_SUP_noCulex.fastq.gz \
                -l 15_extracting_reads_qualfilter/ids_of_mapped_reads_contig${i}.txt \
                -d " " --keep-ids

    # make the name more explicit and remove the file with sequences we are not interested in            
    mv ../Reads_Alice/Cel_g3_SUP_noCulex.fastq.gz.survived  15_extracting_reads_qualfilter/mapped_reads_contig${i}.fastq
    rm ../Reads_Alice/Cel_g3_SUP_noCulex.fastq.gz.removed

    # keep only sequence info (go from fastq to fasta)
    iu-fastq-to-fasta 15_extracting_reads_qualfilter/mapped_reads_contig${i}.fastq 
    mv 15_extracting_reads_qualfilter/mapped_reads_contig${i}.fastq-FASTA 15_extracting_reads_qualfilter/mapped_reads_contig${i}.fasta
done

Getting only reads of size over 9000 nt

This command gets the reads that map to contig_9, the plasmid and filter them to generate a .bam containing only those whose length is over 9000.

samtools view -h mapping_to_contig9_filtered.bam | awk '{if( ( $0 ~ /^@/ ) || ( length($10) > 9000) ) {print $0}}'| samtools view -bh -o over9000.bam

This command generates a FASTA file from the .bam, with the length and read identifier in the header of each sequence.

samtools view over9000.bam | awk '{print ">"$1"_"length($10)"\n"$10}' > sequences_of_bam.fasta

grep '>' files/sequences_of_bam.fasta

## >6119c259-f59d-4eaf-9a69-689c237fac78_15116
## >dee6845e-22e6-40a2-b5cb-317536e63ab6_19413
## >bf711754-5c4a-4ace-bc7d-4729c7dd7b1d_22496
## >ce5485f9-2c39-4a6c-ade4-448fe59cdf9f_18092
## >a88955c7-22b0-470b-af77-8410680ca062_9016
## >3215ad22-fc70-487a-919a-f20174d57d9a_9181
## >5e7757b8-f603-4965-9b9b-1828ac77338a_9210
## >95840707-71b8-44d6-b180-740da0e45301_9209
## >d5c2d544-e5f2-4278-81de-c7b2d67517dd_9194
## >e3621e35-8557-4ce0-b8a2-af6093995df9_9561
## >f6f5290f-9f7c-4c82-80b9-f35afacb449e_10204
## >fd312871-80fb-43e7-b391-35d57e9ac216_9009
## >bf711754-5c4a-4ace-bc7d-4729c7dd7b1d_14434
## >d3d74bac-291c-47eb-8fd6-33e62afb0716_9182
## >e5a05178-9539-4b77-9932-ef7f39e11bf5_10328
## >4f48d6ac-a63d-4a0c-90f8-8669dc04368e_10322
## >eefe0d2e-9e42-4988-a2e7-f30b596c1123_9263
## >55d77c18-3a14-43f8-8fb3-a5aba3aa3e39_9218
## >89b07288-ddc7-4142-b12c-9e3008b53682_9182
## >a2d06b4d-e27a-48c7-90c4-38829135c02b_9130
## >b60ecfc3-a042-4088-b678-50af8579e27f_9751
## >ef25d816-bfee-4b3b-9628-ac041ff2febb_9202
## >b5381dd8-8b27-49de-bf29-bbb53cfdd753_9341
## >5316a16c-b905-4281-8b2e-7d841ebb462c_9267
## >e0b43ef5-be84-43ff-9209-3d4b2f76b18a_9119

BLAST adapter test

Figure: Blast Contig 19 Figure: Blast Contig 19 Figure: Blast Contig 19 Figure: Blast Contig 19

Figure: Blast Contig 19

Figure: Blast Contig 19

Figure: Blast Contig 19

There are no adapters in the long reads extracted from contig_9, including for those that most support the idea of a dimeric plasmid.