Processing for Data_Tam_DNAseq_2025_ATCC19606

Leave a reply

  1. Download the raw data.

    see README_Downloads

  2. Call variant calling using snippy

     ln -s ~/Tools/bacto/db/ .;
 ln -s ~/Tools/bacto/envs/ .;
 ln -s ~/Tools/bacto/local/ .;
 cp ~/Tools/bacto/Snakefile .;
 cp ~/Tools/bacto/bacto-0.1.json .;
 cp ~/Tools/bacto/cluster.json .;

 mkdir raw_data; cd raw_data;

 # Note that the names must be ending with fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W2/W2_1.fq.gz W2_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W2/W2_2.fq.gz W2_R2.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W3/W3_1.fq.gz W3_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W3/W3_2.fq.gz W3_R2.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W4/W4_1.fq.gz W4_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/W4/W4_2.fq.gz W4_R2.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y2/Y2_1.fq.gz Y2_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y2/Y2_2.fq.gz Y2_R2.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y3/Y3_1.fq.gz Y3_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y3/Y3_2.fq.gz Y3_R2.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y4/Y4_1.fq.gz Y4_R1.fastq.gz
 ln -s ../RSMD00304/X101SC24065637-Z01/X101SC24065637-Z01-J002/01.RawData/Y4/Y4_2.fq.gz Y4_R2.fastq.gz

 #download CP059040.gb from GenBank
 mv ~/Downloads/sequence\(1\).gb db/CP059040.gb
 #setting the following in bacto-0.1.json

     "fastqc": false,
     "taxonomic_classifier": false,
     "assembly": true,
     "typing_ariba": false,
     "typing_mlst": true,
     "pangenome": true,
     "variants_calling": true,
     "phylogeny_fasttree": true,
     "phylogeny_raxml": true,
     "recombination": false, (due to gubbins-error set false)

     "genus": "Acinetobacter",
     "kingdom": "Bacteria",
     "species": "baumannii",  (in both prokka and mykrobe)
     "reference": "db/CP059040.gb"

 mamba activate /home/jhuang/miniconda3/envs/bengal3_ac3
 (bengal3_ac3) /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds

  3. Download complete_genome_470_ena_taxon_descendants_count.fasta (Done in ~/REFs)

     #Acinetobacter baumannii Taxonomy ID: 470
 #esearch -db nucleotide -query "txid470[Organism:exp]" | efetch -format fasta -email j.huang@uke.de > genome_470_ncbi.fasta
 #python ~/Scripts/filter_fasta.py genome_470_ncbi.fasta complete_genome_470_ncbi.fasta  #

 # ---- Download related genomes from ENA ----
 https://www.ebi.ac.uk/ena/browser/view/470
 #Click "Sequence" and download "Counts" (13059) and "Taxon descendants count" (16091) if there is enough time! Downloading time points is 28.02.2025.
 python ~/Scripts/filter_fasta.py  ena_470_sequence.fasta complete_genome_470_ena_taxon_descendants_count.fasta  #16091-->920
 #python ~/Scripts/filter_fasta.py ena_470_sequence_Counts.fasta complete_genome_470_ena_Counts.fasta  #xxx, 5.8G

  4. Run vrap

     #replace --virus to the specific taxonomy (e.g. Acinetobacter baumannii) --> change virus_user_db --> specific_bacteria_user_db
 cp trimmed/Y2_trimmed_P*.fastq .
 cp trimmed/W2_trimmed_P*.fastq .
 gzip *.fastq
 ln -s ~/Tools/vrap/ .
 mamba activate /home/jhuang/miniconda3/envs/vrap
 (vrap) vrap/vrap.py  -1 Y2_trimmed_P_1.fastq.gz -2 Y2_trimmed_P_2.fastq.gz -o vrap_Y2 --bt2idx=/home/jhuang/REFs/genome --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/REFs/complete_genome_470_ena_taxon_descendants_count.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
 (vrap) vrap/vrap.py  -1 W2_trimmed_P_1.fastq.gz -2 W2_trimmed_P_2.fastq.gz -o vrap_W2 --bt2idx=/home/jhuang/REFs/genome --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/REFs/complete_genome_470_ena_taxon_descendants_count.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g

  5. Using spandx calling variants (almost the same results to the one from viral-ngs!)

     #DONE
 mkdir ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/CP059040
 cp db/CU459141.gb  ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/CP059040/genes.gbk
 vim ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config
 #CU459141.chromosome : CU459141
 /home/jhuang/miniconda3/envs/spandx/bin/snpEff build -c ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config -genbank CU459141 -v   #-d
 #00:00:02 Saving sequences for chromosome 'CU459141.1 ...
 mamba activate /home/jhuang/miniconda3/envs/bengal3_ac3
 /home/jhuang/miniconda3/envs/bengal3_ac3/bin/snpEff build -c ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config -genbank CP059040 -v

 ~/Scripts/genbank2fasta.py CP059040.gb
 mv CP059040.gb_converted.fna CP059040.fasta    #rename "CP059040.1 xxxxx" to "CP059040" in the fasta-file

 ln -s /home/jhuang/Tools/spandx/ spandx
 mamba activate /home/jhuang/miniconda3/envs/spandx
 (spandx) nextflow run spandx/main.nf --fastq "trimmed/*_P_{1,2}.fastq" --ref CP059040.fasta --annotation --database CP059040 -resume  # 386 records in Outputs/Phylogeny_and_annotation/All_SNPs_indels_annotated.txt

 #DEBUG the results of spandx
 cd Outputs/Master_vcf
 (spandx) cp -r ../../snippy/Y2/reference .
 (spandx) cp ../../spandx/bin/SNP_matrix.sh ./
 #Adapt "snpEff eff -no-downstream -no-intergenic -ud 100 -formatEff -v ${variant_genome_path} out.vcf > out.annotated.vcf" to
 "/home/jhuang/miniconda3/envs/bengal3_ac3/bin/snpEff eff -no-downstream -no-intergenic -ud 100 -formatEff -c reference/snpeff.config -dataDir . ref out.vcf > out.annotated.vcf" in SNP_matrix.sh
 (spandx) bash SNP_matrix.sh CP059040 .

 python3 ~/Scripts/summarize_snippy_res.py snippy
 grep -v "None,,,,,,None,None" summary_snps_indels.csv > summary_snps_indels_.csv  # 29 records

 diff snippy/id2.txt Outputs/Master_vcf/id2.txt

  6. Draw the BRIC plots of Y2, Y3, Y4, W2, W3, W4 on reference.

     mkdir ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606
 cp CP059040.gb ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/W2/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/W2.fasta
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/W3/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/W3.fasta
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/W4/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/W4.fasta
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/Y2/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/Y2.fasta
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/Y3/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/Y3.fasta
 cp ~/DATA/Data_Tam_DNAseq_2025_ATCC19606/shovill/Y4/contigs.fa ~/brigWorkspace_Tam_DNAseq_2025_ATCC19606/gen/Y4.fasta

 #!!!!IMPORTANT_CRITICAL_DEBUG!!!!: we should use jdk1.6, otherwise results in ERROR!
 ~/Tools/Java1.6/jdk1.6.0_45/bin/java -Xmx25000M -jar BRIG.jar
 #DEBUG2
 Output-file should be "/home/jhuang/CP059040" to avoid to cover CP059040.gb.

 mv W2.fastaVsCP059040.gb.fna.tab W2.fastaVsCP059040.gb.fna.tab_
 mv W3.fastaVsCP059040.gb.fna.tab W3.fastaVsCP059040.gb.fna.tab_
 mv W4.fastaVsCP059040.gb.fna.tab W4.fastaVsCP059040.gb.fna.tab_
 mv Y2.fastaVsCP059040.gb.fna.tab Y2.fastaVsCP059040.gb.fna.tab_
 mv Y3.fastaVsCP059040.gb.fna.tab Y3.fastaVsCP059040.gb.fna.tab_
 mv Y4.fastaVsCP059040.gb.fna.tab Y4.fastaVsCP059040.gb.fna.tab_
 python3 ~/Scripts/filter_blast.py W2.fastaVsCP059040.gb.fna.tab_ W2.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95
 python3 ~/Scripts/filter_blast.py W3.fastaVsCP059040.gb.fna.tab_ W3.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95
 python3 ~/Scripts/filter_blast.py W4.fastaVsCP059040.gb.fna.tab_ W4.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95
 python3 ~/Scripts/filter_blast.py Y2.fastaVsCP059040.gb.fna.tab_ Y2.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95
 python3 ~/Scripts/filter_blast.py Y3.fastaVsCP059040.gb.fna.tab_ Y3.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95
 python3 ~/Scripts/filter_blast.py Y4.fastaVsCP059040.gb.fna.tab_ Y4.fastaVsCP059040.gb.fna.tab --evalue 1e-10 --min_length 200 --min_identity 95

 #Rerun without 'Re-do BLAST'

  7. Detect the gap positions

     blastn -query gen/Y4.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out Y4_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py Y4_vs_ref.tsv Y4_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' Y4_vs_ref.tsv_ | sort -k13,13n > Y4_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' Y4_sorted.tsv

 Gap in reference: 17584 to 17854, size: 271 nt
 Gap in reference: 737225 to 741666, size: 4442 nt * (adeI and adeJ)
 Gap in reference: 864814 to 867580, size: 2767 nt * (H0N29_04000)
 Gap in reference: 1326848 to 1327449, size: 602 nt
 Gap in reference: 1328557 to 1328922, size: 366 nt
 Gap in reference: 1481945 to 1482157, size: 213 nt
 Gap in reference: 1615852 to 1616118, size: 267 nt
 Gap in reference: 1641267 to 1641756, size: 490 nt
 Gap in reference: 1679492 to 1692554, size: 13063 nt
 Gap in reference: 1725690 to 1806693, size: 81004 nt
 Gap in reference: 1817484 to 1817867, size: 384 nt
 Gap in reference: 1851794 to 1852056, size: 263 nt
 Gap in reference: 1864443 to 1864741, size: 299 nt
 Gap in reference: 2123729 to 2123929, size: 201 nt
 Gap in reference: 2146439 to 2147201, size: 763 nt
 Gap in reference: 2269547 to 2270147, size: 601 nt
 Gap in reference: 2382428 to 2383507, size: 1080 nt
 Gap in reference: 2451609 to 2452242, size: 634 nt
 #between H0N29_13300 and H0N29_13660
 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *
 Gap in reference: 2899740 to 2900190, size: 451 nt

 Gap in reference: 3010504 to 3010847, size: 344 nt
 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3860147, size: 25556 nt
 Gap in reference: 3862342 to 3862920, size: 579 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

 blastn -query gen/Y3.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out Y3_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py Y3_vs_ref.tsv Y3_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' Y3_vs_ref.tsv_ | sort -k13,13n > Y3_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' Y3_sorted.tsv

 Gap in reference: 737225 to 741666, size: 4442 nt *
 Gap in reference: 792096 to 792573, size: 478 nt

 Gap in reference: 1679492 to 1724344, size: 44853 nt
 Gap in reference: 1725690 to 2190107, size: 464418 nt

 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *

 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3881864, size: 47273 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

 blastn -query gen/Y2.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out Y2_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py Y2_vs_ref.tsv Y2_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' Y2_vs_ref.tsv_ | sort -k13,13n > Y2_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' Y2_sorted.tsv
 Gap in reference: 737225 to 741666, size: 4442 nt *
 Gap in reference: 794508 to 1327449, size: 532942 nt
 Gap in reference: 1679492 to 1724344, size: 44853 nt
 Gap in reference: 1725690 to 2190107, size: 464418 nt
 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *
 Gap in reference: 3288370 to 3642512, size: 354143 nt
 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3881863, size: 47272 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

 blastn -query gen/W4.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out W4_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py W4_vs_ref.tsv W4_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' W4_vs_ref.tsv_ | sort -k13,13n > W4_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' W4_sorted.tsv

 Gap in reference: 737225 to 741666, size: 4442 nt *
 Gap in reference: 792438 to 793713, size: 1276 nt
 Gap in reference: 1679492 to 1724344, size: 44853 nt
 Gap in reference: 1725690 to 1817596, size: 91907 nt
 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *
 Gap in reference: 3288370 to 3642512, size: 354143 nt
 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3881863, size: 47272 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

 blastn -query gen/W3.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out W3_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py W3_vs_ref.tsv W3_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' W3_vs_ref.tsv_ | sort -k13,13n > W3_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' W3_sorted.tsv

 Gap in reference: 737225 to 741666, size: 4442 nt *
 Gap in reference: 792096 to 792573, size: 478 nt
 Gap in reference: 1679492 to 1724344, size: 44853 nt
 Gap in reference: 1725690 to 2190107, size: 464418 nt
 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *
 Gap in reference: 3288370 to 3642512, size: 354143 nt
 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3853883, size: 19292 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

 blastn -query gen/W2.fasta -subject scratch/CP059040.gb.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out W2_vs_ref.tsv
 python3 ~/Scripts/filter_blast.py W2_vs_ref.tsv W2_vs_ref.tsv_ --evalue 1e-10 --min_length 200 --min_identity 95
 awk '{start=($9<$10)?$9:$10; end=($9>$10)?$9:$10; print $0, start, end}' W2_vs_ref.tsv_ | sort -k13,13n > W2_sorted.tsv
 awk '
 {
 if (NR > 1) {
         gap_size = $13 - prev_end - 1
         if (gap_size > 200) {
         print "Gap in reference: " prev_end+1 " to " $13-1 ", size: " gap_size " nt"
         }
 }
 prev_end = $14
 }' W2_sorted.tsv

 Gap in reference: 737225 to 741666, size: 4442 nt *
 Gap in reference: 1679492 to 1724344, size: 44853 nt
 Gap in reference: 1725690 to 2190107, size: 464418 nt
 Gap in reference: 2810875 to 2813304, size: 2430 nt *
 Gap in reference: 2814695 to 2858885, size: 44191 nt *
 Gap in reference: 2860666 to 2862228, size: 1563 nt *
 Gap in reference: 2862967 to 2863469, size: 503 nt *
 Gap in reference: 3288370 to 3642512, size: 354143 nt
 Gap in reference: 3833923 to 3834384, size: 462 nt
 Gap in reference: 3834592 to 3881863, size: 47272 nt
 Gap in reference: 3892085 to 3913973, size: 21889 nt
 Gap in reference: 3914242 to 3942704, size: 28463 nt

Leave a Reply

Your email address will not be published. Required fields are marked *