Влияние меди на антибиотикорезистентность у синегнойной палочки

Задание

Скачайте данные SRR12905311 из SRA NCBI. Опишите, что за исследование проводилось, на каком организме. Можно воспользоваться GEO (GSM4862160).
Препроцессируйте риды (оценка качества, триммирование по необходимости, выравнивание на референс), и проведите поиск SNP.
Выберите 3-4 выявленные мутации. Проверьте по NCBI, каким генам они принадлежат, и предположите, как эти мутации могут влиять на жизнедеятельность организма.
Не забывайте, что ваши мутации не относятся к человеку.

Используемое ПО

sudo apt install fastp fastqc multiqc
sudo apt install spades; pip3 install quast
sudo apt install samtools bwa
sudo apt install kraken2 centrifuge
sudo apt install samtools bamtools freebayes bedtools tabix bcftools  # libvcflib-tools
sudo apt install busco augustus prokka igv
sudo apt install mafft raxml iqtree bedtools
sudo apt install snpeff genometools bedtools
sudo apt install ea-utils sickle seqtk bowtie2 muscle

bash project.sh

Bash-скрипт для обработки данных

#!/bin/bash -x

# [SRX9370194: GSM4862160: PAO1_MIC_Exp5_t=72 (DNA-seq); Pseudomonas aeruginosa; OTHER](https://www.ncbi.nlm.nih.gov/sra/SRR12905311)
# Study: Do Copper enhance antibiotic resistance (AR) in resistant bacteria and/or induce AR in sensitive bacteria? [DNA-seq]
# Organism: [Pseudomonas aeruginosa (Синегнойная палочка)](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=287)
# Reference genome: [ASM676v1](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000006765.1)
# GEO Accession: [GSM4862160](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4862160)
# SRA: [SRP288731](https://trace.ncbi.nlm.nih.gov/Traces/?view=study&acc=SRP288731)
# BioProject: [PRJNA672287](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA672287)
# GEO: [GSE160188](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160188)

ORGANISM=PAeruginosa
GENOME=GCF_000006765.1
RUN=SRR12905311

ABRICATE=~/Apps/abricate/bin/abricate
QUALIMAP=~/Apps/qualimap_v2.3/qualimap


RESULT_DIR=/mnt/c/Users/Asus/BioGen
[ -d ${RESULT_DIR} ] && rm -rf ${RESULT_DIR} && mkdir ${RESULT_DIR}

#mkdir -p project && cd project

mkdir -p reads ref

# Скачивание ридов в формате FASTQ
[ -f reads/${RUN}_1.fastq.gz ] || wget -P reads ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR129/011/SRR12905311/SRR12905311_1.fastq.gz
[ -f reads/${RUN}_2.fastq.gz ] || wget -P reads ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR129/011/SRR12905311/SRR12905311_2.fastq.gz

# Скачивание референсного генома для Pseudomonas aeruginosa в формате FASTA и GFF3
[ -f ${GENOME}.zip ] || curl -OJX GET "https://api.ncbi.nlm.nih.gov/datasets/v2alpha/genome/accession/GCF_000006765.1/download?include_annotation_type=GENOME_FASTA,GENOME_GFF,RNA_FASTA,CDS_FASTA,PROT_FASTA,SEQUENCE_REPORT&filename=GCF_000006765.1.zip" -H "Accept: application/zip"
unzip ${GENOME}.zip
mv ncbi_dataset/data/${GENOME}/${GENOME}_ASM*_genomic.fna ref/${ORGANISM}.fasta
mv ncbi_dataset/data/${GENOME}/genomic.gff ref/${ORGANISM}.gff
rm -rf ncbi_dataset README.md

rm -rf quality_check && mkdir -p quality_check
# Тримминг
mkdir -p trimmed
fastp --overrepresentation_analysis --trim_poly_g \
    -f 12 -t 2 \
    -h quality_check/${RUN}.fastp.html -j quality_check/${RUN}.fastp.json --thread `nproc` \
    -i reads/${RUN}_1.fastq.gz -I reads/${RUN}_2.fastq.gz -o trimmed/${RUN}_1.fastq.gz -O trimmed/${RUN}_2.fastq.gz
firefox quality_check/${RUN}.fastp.html

FASTQ_DIR=trimmed  # reads || trimmed
# Оценка качества ридов
fastqc -o quality_check $FASTQ_DIR/*.fastq.gz ../${RUN}.fastq.gz
firefox quality_check/${RUN}_*_fastqc.html
rm -rf multiqc_data multiqc_report.html
multiqc quality_check $FASTQ_DIR -o quality_check
firefox quality_check/multiqc_report.html || mv quality_check ${RESULT_DIR}

#mkdir -p assembly
# Creating a genome assembly
#spades.py -o assembly/spades-original --careful -1 $FASTQ_DIR/${RUN}_1.fastq.gz -2 $FASTQ_DIR/${RUN}_2.fastq.gz --threads `nproc`
# Assembly quality assessment
#quast.py -o assembly/quast assembly/spades-original/scaffolds.fasta --threads `nproc`
#firefox assembly/quast/report.pdf || mv assembly/quast/report.pdf ${RESULT_DIR}/assembly_report.pdf
#${ABRICATE} assembly/spades-original/scaffolds.fasta --threads `nproc` > abricate_report.tsv

# Выравнивание ридов на геном
mkdir -p index mappings
# generate index for mapping
bwa index ref/${ORGANISM}.fasta
#bowtie2-build ref/${ORGANISM}.fasta index/${ORGANISM}
# mapping reads
bwa mem ref/${ORGANISM}.fasta $FASTQ_DIR/${RUN}_1.fastq.gz $FASTQ_DIR/${RUN}_2.fastq.gz > mappings/${ORGANISM}.sam -t `nproc`
#bowtie2 -x index/${ORGANISM} -1 $FASTQ_DIR/${RUN}_1.fastq.gz -2 $FASTQ_DIR/${RUN}_2.fastq.gz -S mappings/${ORGANISM}.sam --threads `nproc`

# converting output from SAM to BAM
samtools view -b -o mappings/${ORGANISM}.bam mappings/${ORGANISM}.sam
rm mappings/${ORGANISM}.sam
# show statistics
samtools flagstat mappings/${ORGANISM}.bam
# sorting by BAM coordinate (good practice)
samtools sort -@ 8 -o mappings/${ORGANISM}.sorted.bam mappings/${ORGANISM}.bam

# building index for BAM and FASTA (to access the files faster)
samtools index mappings/${ORGANISM}.sorted.bam
samtools faidx ref/${ORGANISM}.fasta

# collect BAM statistics and generate plots
rm -rf bam_stats
samtools stats mappings/${ORGANISM}.sorted.bam > mappings/${ORGANISM}.sorted.bam.bc
plot-bamstats -p bam_stats/ mappings/${ORGANISM}.sorted.bam.bc
firefox bam_stats/index.html
${QUALIMAP} bamqc -bam mappings/${ORGANISM}.sorted.bam
firefox mappings/${ORGANISM}.sorted_stats/qualimapReport.html

# SNP calling
mkdir -p variants
#bcftools mpileup -f ref/${ORGANISM}.fasta mappings/${ORGANISM}.sorted.bam | bcftools call -mv -o calls.vcf && cat calls.vcf | grep -v ^## | head -5 && bcftools view -i '%QUAL>=20' calls.vcf
freebayes -p 1 -f ref/${ORGANISM}.fasta mappings/${ORGANISM}.sorted.bam > variants/${ORGANISM}.freebayes.vcf  # '-p 1' - specifies the ploidy level: E.Coli are haploid

# collect VCF statistics and generate plots
rm -rf vcf_stats
bgzip -f variants/${ORGANISM}.freebayes.vcf
tabix -p vcf variants/${ORGANISM}.freebayes.vcf.gz
bcftools stats -F ref/${ORGANISM}.fasta -s - variants/${ORGANISM}.freebayes.vcf.gz > variants/${ORGANISM}.freebayes.vcf.gz.stats
plot-vcfstats -p vcf_stats/ variants/${ORGANISM}.freebayes.vcf.gz.stats
firefox vcf_stats/summary.pdf
#bgzip -d variants/${ORGANISM}.freebayes.vcf.gz

# Отображение результатов
igv -g ref/${ORGANISM}.fasta mappings/${ORGANISM}.sorted.bam variants/${ORGANISM}.freebayes.vcf.gz ref/${ORGANISM}.gff

# Удаление ненужных данных
mv mappings/${ORGANISM}.sorted_stats ./qualimap_report
[ -d ${RESULT_DIR} ] && mv bam_stats vcf_stats qualimap_report ${RESULT_DIR}
mkdir -p data && mv ref/${ORGANISM}.fasta mappings/${ORGANISM}.sorted.bam* variants/${ORGANISM}.freebayes.vcf.gz GSE160188_core.vcf.gz ref/${ORGANISM}.gff data
rm -rf variants mappings index assembly trimmed ref

---

Результаты:

Аннотация:

• SRX9370194: GSM4862160: PAO1_MIC_Exp5_t=72 (DNA-seq); Pseudomonas aeruginosa; OTHER
• Исследование: Повышает ли медь устойчивость к антибиотикам у резистентных бактерий и/или индуцирует антибиотикорезистентность у чувствительных бактерий?
• Организм: Pseudomonas aeruginosa (Синегнойная палочка)
• Штамм: PAO1 (дикий тип)
• Референсный геном: ASM676v1
• GEO Accession: GSM4862160
• SRA: SRP288731
• BioProject: PRJNA672287
• GEO: GSE160188

Отчёты:

Результаты тримминга
Оценка качества (рид 1)
Оценка качества (рид 2)
Проверка качества (сводный отчёт)
Сборка ридов
Выравнивание ридов на референс
Проверка качества выравнивания
Поиск SNP
IGV: igv -g data/PAeruginosa.fasta data/PAeruginosa.sorted.bam data/PAeruginosa.freebayes.vcf.gz data/PAeruginosa.gff
ABRicate: Наличие генов устойчивости к противомикробным препаратам

Комментарий: Для улучшения качества ридов была проведена операция тримминга. В результате распределение ATGC в последовательностях (Per base sequence content) стало более равномерным. Распределение длины последовательности (Sequence Length Distribution) и содержание GC (Per sequence GC content) вполне сносное для дальнейшего анализа.

Гены и мутации:

Ген	Symbol	Locus Tag	IGV window	Пример мутации	Комментарий
aminoglycoside 3-N-acetyltransferase	aph	PA4119		Type: MNP Position: 4607781-4607788 Reference: GCAGATTT* Alternate: AGAGGTTG	Предположительно ген влияет на резистентность к KANAMYCIN. Множественные мутации в гене вероятно повлияли на резистентность.
chloramphenicol acetyltransferase	cat	PA0706		Type: MNP Position: 779812-779814 Reference: ACA* Alternate: GCC	Предположительно ген влияет на резистентность к CHLORAMPHENICOL. Несколько мутаций в гене могли повлиять на резистентность.
glutathione transferase FosA	-	PA1129		Type: SNP Position: 1222092 Reference: С* Alternate: T	Предположительно ген влияет на резистентность к FOSFOMYCIN. Маловероятно, что одна мутация серьёзно повлияла на резистентность.