| ARG | | read count | | gene length | | coverage | |
|---|
| argannot~~~(Bla)OXA-347~~~JN086160:1583-2407 | 41 | 825 | 2D822M1D |
| argannot~~~(Bla)cfxA5~~~AY769934:28-993 | 235 | 966 | 962M4D |
| argannot~~~(MLS)ErmB~~~M11180:714-1451 | 245 | 738 | 2D728M8D |
| argannot~~~(Tet)TetQ~~~Z21523:362-2287 | 1117 | 1974 | 3D1971M |
The final column of the GROOT report output is a CIGAR-like representation of the covered bases: M=covered D=absent. The read count and coverage (plus the coverage plots) can help us assess the confidence in the reported genes when we use `--lowCov`. We can also check if reported ARGs share classified reads (i.e. multimappers). To compare reads aligning to 2 ARGs, you could try some commands like this:
```bash
# split the bam by ARG
bamtools split -in out.bam -reference
# extract reads aligning to specific ARGs
bedtools bamtofastq -i ARG1.bam -fq reads-aligned-to-ARG1.fq
bedtools bamtofastq -i ARG27.bam -fq reads-aligned-to-ARG27.fq
# find matching reads
bbduk.sh in=reads-aligned-to-ARG1.fq ref=reads-aligned-to-ARG27.fq outm=matched.fq k=100 mm=f
```
Now let's try plotting the resistome profiles:
More to content to come...
---
## 5. ARG context
Now we have a set of ARGs we know are present in one or more samples, we might want to take a look at what is carrying these genes. We will use [metacherchant](https://www.ncbi.nlm.nih.gov/pubmed/29092015) to determine the genetic context of ARGs using subgraphs. We will also use [Bandage](https://rrwick.github.io/Bandage/) and [Kraken](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-3-r46) (not installed as part of our conda environment).
Download the full ARG-ANNOT set of genes and then index with samtools:
```bash
wget https://github.com/will-rowe/groot/raw/master/db/full-ARG-databases/arg-annot-db/argannot-args.fna
samtools faidx argannot-args.fna
```
Now we can extract all the genes present in our resistome profiles:
```bash
samtools faidx argannot-args.fna `cut -f1 combined-profiles.urban.tsv` > urban-args.fna
samtools faidx argannot-args.fna `cut -f1 combined-profiles.rural.tsv` > rural-args.fna
```
Remove any duplicate genes:
```bash
cat urban-args.fna | seqkit rmdup -s -o urban-args.fna
cat rural-args.fna | seqkit rmdup -s -o rural-args.fna
```
As a first example, let's look at the genomic context of one ARG (`(Bla)OXA-347`) in urban samples. To begin, get the samples that contain this ARG:
```bash
grep -l "argannot~~~(Bla)OXA-347~~~JN086160:1583-2407" *.lowCov.report | sed 's/.lowCov.report//' > samples-containing-blaOXA-347.list
```
We also need to store the reference sequence:
```bash
samtools faidx argannot-args.fna "argannot~~~(Bla)OXA-347~~~JN086160:1583-2407" > blaOXA-347.fna
```
Now we can run metacherchant:
```bash
while read -r sample
do
metacherchant.sh --tool environment-finder \
-k 31 \
--coverage=2 \
--maxradius 1000 \
--reads ${sample}.fastq \
--seq blaOXA-347.fna \
--output "./metacherchant/${sample}/output" \
--work-dir "./metacherchant/${sample}/workDir" \
-p 42 \
--trim
done < samples-containing-blaOXA-347.list
```
To classify the contigs assembled around the detected ARGs, try Kraken:
```bash
while read -r sample
do
kraken --threads 8 --preload --fasta-input --db /path/to/minikraken_20141208 metacherchant/${sample}/output/seqs.fasta | kraken-report --db /path/to/minikraken_20141208 > metacherchant/${sample}/${sample}-kraken.report
done < samples-containing-blaOXA-347.list
```