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

Project 1

Jump to bottom
Arjun Prasad 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

Step 1a

# 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.

Step 2a: Create a working directory

cd ~/project1
mkdir step2
cd step2

Step 2b: 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

Step 2c: 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

Step 2c: Create the dsub command to perform the assembly

We are going to use dsub to run skesa. Dsub is a command line utility that allows users a simple way to tie into cloud APIs to manage tasks. It can be run locally in order to prototype and quickly iterate. It uses Docker images for compute tasks.

    # The assembly command. I like to create a script for easier editing and record-keeping
        cat > assemble <<END
dsub \\
  --provider google-cls-v2 \\
  --project asm-ngs-workshop\\
  --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

** why are we doing this complicated step of writing a script?

Step 2d: 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)

Step 2e: 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/

Step 2f: You should now be able to see the FASTA assemblies you just created.

    ls -l asm/*

Step 2g: 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

Step 3a: Create directories

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

Step 3b: 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

Step 3c: 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

Step 3d: Run dsub

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

Step 3e: Copy amrfinder files from bucket to local disk

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

**** Step 3f: 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.

Step 4a: Create dirs, and copy data and files and run conversion

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.

Step 5a: 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

Step 5b: 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

Step 6a: 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

Step 6b: Find disagreements

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

Discuss findings

Clone this wiki locally