analyzing high-frequency mutations in the flagellar genes of V. cholerae
Creators: name: Guozhong, Chen
Publisher: name: Zenodo
Year: 2025
Subjects: vibrio cholerae

description: High-frequency mutations in the flagellar genes of V. cholerae analysis

 

Step1: prodigal predictions

input folder-"/data/cgz"

for fasta file in "$input folder"/*.fna; do base name-$(basename "$fasta file" .fna)

prodigal -a "./$ fbase name}-P.faa"-d "./$ fbase name}-N.fna"-f gff\-o"./$ {base name}.gff'\p meta s "./Sfbase name}.stat"-i "Sfasta file"done

 

Step2:diamond 

diamond makedb --in *.fasta  --db name of libary

input_folder="/data/chengz/Cohort1/TotalBin/Bin_ALLcopy/Flitervc/prodigal"

diamond_db="./flagellumlibpro.dmnd"

for fasta_file in "$input_folder"/*.faa; do        base_name=$(basename "$fasta_file" .faa)        diamond blastp --db "$diamond_db" \                   --query "$fasta_file" \                   -e 0.00001 \                   --outfmt 6 \                   --more-sensitive \                   --max-target-seqs 1 \                   --threads 32 \                   --id 80 \                   --subject-cover 50 \                   --query-cover 50 \                   --out "./${base_name}.txt"done

 

 

metaWRAP for Cholera

 

Prefetch SRA Filesbash# Download SRA files based on SRR numbersprefetch SRRXXXXX

Move .sra Files to a Single Directorybashfind . -type f -name '*.sra' -exec mv {} /data/chengz/Cohort/CFS-1/CFS-1-fastq \;Check Folder Sizebashdu -sh *Convert .sra Files to FASTQ Formatbashfasterq-dump -p -e 80 --split-3 *.sraQuality Control with MetaWRAPbashsource activate metawrapmetawrap read_qc --skip-bmtagger -1 FLJ2_1.fq -2 FLJ2_2.fq -t 50 -o READ_QC/FLJ2--skip-bmtagger: Skips host contamination removal.Input files (-1 and -2) should be named as sample_1 and sample_2.-t 50: Use 50 threads.Merge FASTQ Filesbashmkdir ALLcat *_R1.fastq > ALL/ALL_READS_R1.fastqcat *_R2.fastq > ALL/ALL_READS_R2.fastqRemove Host ContaminationPaired-End Sequencingbashfor i in SRR26886987 SRR26886988 ...; do     bowtie2 -p 30 -x /mnt/database/humanindex/humangenome \        -1 ${i}_1.fastq.gz -2 ${i}_2.fastq.gz \        -S ${i}.sam --un-conc-gz /data/chengz/Cohort1/BangladeshControl/Fastq/${i}.fastq.gz; doneSingle-End Sequencingbashfor i in SRR26267438 SRR26267439 ...; do     bowtie2 -p 30 -x /mnt/database/humanindex/humangenome \        -U ${i}.fastq.gz -S ${i}.sam \        --un-conc-gz /data/chengz/Cohort1/IndianControl/Fastq/Output/${i}.fastq.gz; doneAssembly with MetaWRAPSingle-End Assemblybashmegahit -r SRR15408397_1.fastq --min-contig-len 800 -t 96 -o ASSEMBLYPaired-End Assemblybashfor i in *.fastq.1.gz; do     f=$(basename "$i" .fastq.1.gz)    metawrap assembly -1 $f.fastq.1.gz -2 $f.fastq.2.gz -t 50 -l 800 -o ASSEMBLY/${f}; doneBinning with MetaWRAPPrepare Filenamesbashfor file in input_dir/*; do    if [[ $file == *.fastq.1.gz ]]; then        mv $file $(echo $file | sed 's/\.fastq\.1\.gz/_1.fastq.gz/')    fidonefor file in input_dir/*; do    if [[ $file == *.fastq.2.gz ]]; then        mv $file $(echo $file | sed 's/\.fastq\.2\.gz/_2.fastq.gz/')    fidoneBinningbashfor sample in "${samples[@]}"; do     metawrap binning -o "${sample}_BINNING" -t 40 \        -a ASSEMBLY/${sample}/final_assembly.fasta \        --metabat2 --maxbin2 --concoct \        ${sample}_1.fastq ${sample}_2.fastqdoneBin Refinementbashfor dir in SRR*_BINNING; do     metawrap bin_refinement -o "${dir}_REBIN" -t 30 \        -A "${dir}/metabat2_bins/" -B "${dir}/maxbin2_bins/" -C "${dir}/concoct_bins/" \        -c 50 -x 10; doneExtract and Rename Refined Binsbashfind SRR*/metawrap_50_10_bins/ -type f -name "*.fa" -exec cp {} /path/to/Output/ \;Remove Redundant Bins with dRepTaxonomic Classification with GTDB-TkWorkflowbashsource activate gtdbtk?nohup gtdbtk classify_wf --genome_dir /data/chengz/Cohort1/TotalBin/BinALL --skip_ani_screen --out_dir /data/chengz/Cohort1/TotalBin/gtdbtkALL --extension fa --cpus 16 &

 

6N Label relative abundance analysis

source activate qiime2-2021.4nano vsearch.sh

!/bin/bash

mkdir -p merged_sequences

for r1 in *R1.clean.fq.gzdo        r2=${r1/_R1.clean.fq.gz/_R2.clean.fq.gz}       output="merged_sequences/$(basename "$r1" _R1.clean.fq.gz).merged.fastq"            vsearch --fastq_mergepairs "$r1" \            --reverse "$r2" \            --fastqout "$output" \            --label_suffix ";sample=$(basename "$r1"_R1.clean.fq.gz)\n"done

echo "Merge completed"

chmod +x vsearch.sh ./vsearch.sh

nano align2.py

import sysimport os

def rev_comp(s):    sc = ''    for i in (s):        if i is 'A':            sc += 'T'        elif i is 'C':            sc += 'G'        elif i is 'G':            sc += 'C'        elif i is 'T':            sc += 'A'    return sc

def count_line(s, pattern, pattern1):    count = 0    with open(s,'r') as input:        for line in input:            line = line.strip()            if (pattern in line) or (pattern1 in line):                count = count + 1    return count

sixbp_dir = '/ME4012-Vol01/lisc/chengz-seq/6N-label'input_dir = '/ME4012-Vol01/lisc/chengz-seq/merged_sequences'

for sixbp_file in os.listdir(sixbp_dir):    if sixbp_file.endswith(".fasta"):        sixbp_path = os.path.join(sixbp_dir, sixbp_file)        sixbp_name = os.path.splitext(sixbp_file)[0]

        for input_file in os.listdir(input_dir):            if input_file.endswith(".fastq"):                input_path = os.path.join(input_dir, input_file)                input_name = os.path.splitext(input_file)[0]

                output_file = f"{sixbp_name}-to-{input_name}.csv"                f = open(output_file, 'w')

                with open(sixbp_path, 'r') as sixbp_input:                    for line in sixbp_input:                        sixbp = line.strip()                        pattern = 'TGATCTGTTTGCTTGATGAAGATGG' + sixbp + 'TGTAGGCTGGAGCTGCTTCTA'                        pattern1 = rev_comp('TGATCTGTTTGCTTGATGAAGATGG') + rev_comp(sixbp) + rev_comp('TGTAGGCTGGAGCTGCTTCTA')                        count = count_line(input_path, pattern, pattern1)                        f.write(f"{sixbp},{count}\n")

                f.close()

python align2.py 

grep -c "^+$" *.fastq 

for file in *.csv; do     sum=$(awk -F',' 'FNR == 1 {next} {sum+=$2} END {print sum}' "$file");     echo "$file,$sum" >> file_sum-2.txt; done