Skip to content
This repository was archived by the owner on Dec 18, 2020. It is now read-only.

Project 1

Jump to bottom
klimke edited this page Nov 25, 2020 · 19 revisions

Analysis of data from Cartography of opportunistic pathogens and antibiotic resistance genes in a tertiary hospital environment

Step 0: Log into your VM and authorize

You will need to log into your newly created VM by clicking the 'SSH' link.

To use google cloud utilities you will need to give permissions to the google cloud tools access to your account. The following command gives the VM full access. We're doing this for simplicity in this workshop; you may want to grant more limited permissions you grant for your own work.

Paste in the following command and follow the instructions.

gcloud auth login

Create a bucket for the workshop

gsutil mb gs://asm-ngs-$USER

Step 1: Get the data

We've done this task for you because it is complex, time consuming, and illustrates a problem with the way data is included from many papers today. Little thought is given to how the data can be accessed even when it is "publicly available". To get the dataset to do further analysis, we've created a separate page listing what we did Project 1: Get IDs. For this workshop you will just copy the processed list of accessions into a work directory

# create work directories
mkdir -p ~/project1/step1
cd ~/project1/step1
    wget https://raw.githubusercontent.com/evolarjun/workshop/main/project1/singapore_samples_srafull.txt

Step 2: Assemble isolates using SKESA

We've used SKESA to assemble the isolates because it is a fast and high-accuracy assembler. You could use another, more thorough assembler such as SPAdes to gain some additional sequence, but at the cost of some expense and accuracy.

We're going to set you up to do an assembly in parallel using dsub. In the interests of saving some compute we'll set up for computation of just two assemblies, but the same setup was used to do all 106 assemblies from the paper.

Create a working directory

cd ~/project1
mkdir step2
cd step2

Link in the sample accession file from the previous step

    # Link in the sample accession file (from the previous step)
    ln -s ../step1/singapore_samples_srafull.txt

Create the submit_list.tsv file which is used to tell dsub what parameters to parallelize

    # Small script to create the submit_list.tsv file for dsub
    cat <<END > create_submit_list.sh
cut -f 3 singapore_samples_srafull.txt \
| perl -ne 'BEGIN { print "--env ERR\t--output ASM\n" }
        chomp;
        print join("\t",
                "\$_",
                "gs://asm-ngs-$USER/project1/step2/asm/\$_.asm.fa"), "\n"
        '
END
    bash create_submit_list.sh > submit_list.tsv
 
    # For the workshop we want to just run two assemblies
    head -3 submit_list.tsv > submit_list_two.tsv

** ARJUN AS I SAID I WANT THE DSUB COMMAND TO USE FASTQ TO SHOW USERS HOW TO ASSEMBLE DATA NOT IN SRA

Create the dsub command to perform the assembly

    # The assembly command. I like to create a script for easier editing and record-keeping
        cat > assemble <<END
dsub \\
  --provider google-cls-v2 \\
  --project ncbi-pd-amr \\
  --regions us-central1 \\
  --logging gs://asm-ngs-$USER/project1/logs \\
  --image staphb/skesa \\
  --tasks submit_list_two.tsv \\
  --command   'skesa --sra_run \${ERR} --contigs_out \${ASM} --cores 6 ' \\
  --min-ram 24 \\
  --min-cores 6 \\
  --name skesa \\
  --wait
END

Run the assembly command

    # Run the assemblies in parallel on new VMs using dsub
    # You should see new VMs created and running in the VM instances list
    time bash -x assemble 2>&1 | tee assemble.out

You should have to wait 5-10 minutes for this to finish. You might want to open a new SSH window and run the status check command (see the output of dsub)

Copy the assemblies you just made from the bucket to your VM disk

    # copy the assemblies from the bucket onto your local disk.
    mkdir asm
    gsutil -m cp gs://asm-ngs-$USER/acinetobacter/asm/* asm/

You should now be able to see the FASTA assemblies you just created.

    ls -l asm/*

In the interests of time we have pre-computed the rest of the assemblies.

    # copy pre-computed assemblies
    gsutil -m cp gs://asm-ngs/project1/step2/asm/* asm

Step 3: Run AMRFinderPlus on the assemblies

mkdir ~/project1/step3
cd ~/project1/step3

Create the submit_list.tsv file

        echo $'--input ASM\t--output AMR' > submit_list.tsv
        for asm in `gsutil ls gs://asm-ngs/project1/step2/asm`
        do
            acc=`basename $asm .asm.fa`
            echo $asm$'\t'"gs://asm-ngs-$USER/project1/step3/amrfinder/$acc.amrfinder"
        done >> submit_list.tsv

        # just run two in parallel
        head -3 submit_list.tsv > submit_list_two.tsv

Create the dsub command

        cat > run_amrfinder <<END
dsub \\
  --provider google-cls-v2 \\
  --project ncbi-pd-amr \\
  --regions us-central1 \\
  --logging gs://asm-ngs-$USER/project1/logs \\
  --image staphb/ncbi-amrfinderplus \\
  --tasks submit_list_two.tsv \\
  --command   'amrfinder --plus -O Acinetobacter_baumannii -n \${ASM} -o \${AMR}' \\
  --min-ram 3 \\
  --min-cores 2 \\
  --name amrfinder \\
  --wait
END

Run dsub

        time bash -x run_amrfinder 2>&1 | tee run_amrfinder.out

Copy amrfinder files from bucket to local disk

    gsutil -m cp -r gs://asm-ngs/project1/step3/amrfinder .

Switch to dsub

mkdir ~/project1/step3
cd ~/project1/step3
    # link in assembly directory
    ln -s ../asm .

    # for output
    mkdir amrfinder

    for asm in asm/*
    do
        base=`basename $asm .asm.fa`
        amrfinder -O Acinetobacter_baumannii -n $asm --plus --threads 8 > amrfinder/$base.amrfinder
    done

Step 4: Get resistance phenotypes from paper

We have written a short script to go from the "A. baumannii" tab of Supplementary Table 6 from the paper to an easy to work with "tidy" data format. We have already converted the Excel spreadsheet to tab-delimited format by copying and pasting it into a text file.

mkdir ~/project1/step4
cd ~/project1/step4
    
    # get data
    wget https://raw.githubusercontent.com/evolarjun/workshop/main/project1/spreadsheet-phenotypes.tab

    # get conversion script
    wget https://raw.githubusercontent.com/evolarjun/workshop/main/project1/spreadsheet2tidy.pl

    # run conversion
    perl spreadsheet2tidy.pl > spreadsheet-phenotypes.tidy.tab

Step 5: AMR genotype to Phenotype prediction

This step is not trivial for many taxa. For some drugs especially, Acinetobacter baumannii is particularly difficult. For Acinetobacter we have codified what we do know in a small script, but even then, the presence of chromosomal carbapenemases and efflux pumps complicates interpretation.

We wrote a short script to go from AMRFinderPlus reports to predictions of resistance. Note that chromosomal efflux pumps and beta-lactamases play a big role in resistance in Acinetobacter, but it is not well understood which variants affect resistance. For that reason we are not predicting resistance based on efflux pumps, and the presence of beta-lactamases with known resistance profiles may not cause phenotypic resistance.

Note: This small prediction script was created just for this exercise and should not be relied upon to predict phenotypes in any other context, especially not a clinical one.

Set up directory and get script

mkdir ~/project1/step5
cd ~/project1/step5

    # link in data from other steps
    ln -s ../step4/spreadsheet-phenotypes.tidy.tab .
    ln -s ../step3/amrfinder .
    ln -s ../step1/singapore_samples_srafull.txt .

    # get script to convert AMRFinderPlus results to phenotype predictions for comparison
    wget https://raw.githubusercontent.com/evolarjun/workshop/main/project1/acinetobacter_amrfinder2phenotype.pl
    wget https://raw.githubusercontent.com/evolarjun/workshop/main/project1/subclass2phenotype.tab
    chmod 755 acinetobacter_amrfinder2phenotype.pl

Convert AMRFinderPlus output files to phenotypic predictions

    for file in amrfinder/*.amrfinder
    do
        base=`basename $file .amrfinder`
        ./acinetobacter_amrfinder2phenotype.pl $base $file
    done > predicted-phenotypes.tab

Step 6: Identify discrepancies between predicted and measured phenotype

Merge tables

mkdir ~/project1/step6
cd ~/project1/step6
    
    # link in files from previous steps
    ln -s ../step5/predicted-phenotypes.tab . 
    ln -s ../step4/spreadsheet-phenotypes.tidy.tab .
    ln -s ../step1/singapore_samples_srafull.txt .
    
    # Merge tables using 'q' (https://github.com/harelba/q/)
        cat > merge_query.sql <<END
select trans.c1 samplename, trans.c3 accession, pheno.c2 drug, pheno.c3 phenotype,
       pred.c3 predicted, pred.c4 reason
  from spreadsheet-phenotypes.tidy.tab pheno
  join singapore_samples_srafull.txt trans
    on pheno.c1 = trans.c1
  left join predicted-phenotypes.tab pred
    on pred.c1 = trans.c3
   and pred.c2 = pheno.c2
END
    q -t -T -O -q merge_query.sql | awk 'BEGIN {IFS=OFS="\t"} { sub(/^$/, "S", $5)} 1' > all_phenotypes_predictions.tab

Find disagreements

    # find disagreements
    awk '$4 != $5 && !($4 == "I" && $5 == "R") { print }' all_phenotypes_predictions.tab > disagreements.tab

Clone this wiki locally