ONT_pilot_w_GDB_samples

This is the story of how harmless curiosity for Oxford Nanopore Technology lead to an awesome learning experience!

Documented here for those who wish to follow in our footsteps..

It started with a conversation titled, "what is the current state of the ONT data from 2018 and 2019?"

And here, the journey begins with me (SBB) trying to understand how different assembly methods involving short-read and/or long-read sequences influence the number and kind of proteins recovered.
And if one is an improvement over the other

Chapter I - Running clustering analyses using Prodigal output

CCL built an amazing and comprehensive Snakefile
Everything from assemblies, to gene calls, to genome coverages were estimated using short-reads and long-reads from a Generous donor B (GDB) sample collected in 2018 and in 2019
The ensuing story is being told from my perspective including the minor role I will have played in the downstream analyes

The story starts with first working in my own folders

cd /scratch/users/sbusi/ONT
# symlinking cedric's data to mine and then modifying files such as the snakefile
ln -sd /scratch/users/claczny/ont/fecal_pilot_testing/results/annotation/proteins/ cedric_ont_proteins

# included a rule to run cd-hit
rule cdhit:
   input:
     db1=expand("{datadir}/flye/lr/merged/barcode07/assembly.faa", datadir=DATA_DIR),
     db2=expand("{datadir}/megahit/NEB2_MG_S17/final.contigs.faa", datadir=DATA_DIR),
     db3=expand("{datadir}/metaspades_hybrid/lr_barcode07-sr_NEB2_MG_S17/contigs.faa", datadir=DATA_DIR)
   output: g12=expand("{resultsdir}/flye_megahit_novel.fasta", resultsdir=RESULTS_DIR),g23=expand("{resultsdir}/megahit_metaspades_novel.fasta", resultsdir=RESULTS_DIR),g13=expand("{resultsdir}/flye_metaspades_novel.fasta", resultsdir=RESULTS_DIR),g21=expand("{resultsdir}/megahit_flye_novel.fasta", resultsdir=RESULTS_DIR),g32=expand("{resultsdir}/metaspades_megahit_novel.fasta", resultsdir=RESULTS_DIR),g31=expand("{resultsdir}/metaspades_flye_novel.fasta", resultsdir=RESULTS_DIR),gout=expand("{resultsdir}/cd-hit.done", resultsdir=RESULTS_DIR)
   conda:  "/scratch/users/sbusi/cd-hit.yml"
   shell:  """
           date
           cd-hit-2d -i {input.db1} -i2 {input.db2} -o {output.g12} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           cd-hit-2d -i {input.db2} -i2 {input.db3} -o {output.g23} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           cd-hit-2d -i {input.db1} -i2 {input.db3} -o {output.g13} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           cd-hit-2d -i {input.db2} -i2 {input.db1} -o {output.g21} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           cd-hit-2d -i {input.db3} -i2 {input.db2} -o {output.g32} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           cd-hit-2d -i {input.db3} -i2 {input.db1} -o {output.g31} -c 0.9 -n 5 -d 0 -M 16000 -T 4
           touch {output.gout}
           date
           """

# ran the file as follows:
./run_snakemake.sh 

# subsequently used auxilliary scripts to merge the cluster (output) files
cd /scratch/users/sbusi/ONT/cd_hit_output

# making plots for all .clstr files (http://weizhongli-lab.org/lab-wiki/doku.php?id=cd-hit-user-guide#cd-hit-2d)
for file in *.clstr 
do      
    echo "$file"
    plot_len1.pl "$file"  \
           1,2-4,5-9,10-19,20-49,50-99,100-299,500-99999   \
                 10-59,60-149,150-499,500-1999,2000-999999
done >> cluster_plots

# merged multiple clstr files
clstr_merge.pl flye_megahit_novel.fasta.clstr flye_metaspades_novel.fasta.clstr megahit_metaspades_novel.fasta.clstr > merged.clstr
# plotting
plot_len1.pl test.clstr  \
           1,2-4,5-9,10-19,20-49,50-99,100-299,500-99999   \
                 10-59,60-149,150-499,500-1999,2000-999999

however the need to re-run the same samples by swapping was not only cumbersome, but depending on the orientation of the samples, obtained non-overlapping sets of fasta list

So, I decided to pursue mmseqs2 instead

The output of this is now incorporated into the updated_SNAKEFILE that will be used for the analyses going forward.
see also the following:

And this is where the fun begins..

Like any good research story - one question always leads to another, and one answer opens up the rabbit hole of unending possibilities
And this is where the ensuing chapter arrives. Let's call it "guppy"!

Chapter II - Running `methylation-aware` basecalling to improve the otherwise dismal sequencing error rates from ONT

Installed guppy as follows:

# Installation
cd /home/users/sbusi/apps 
wget https://mirror.oxfordnanoportal.com/software/analysis/ont-guppy_research_rle_0.0.0_linux64.tar.gz
tar -xzvf ont-guppy_research_rle_0.0.0_linux64.tar.gz
export PATH=$PATH:"/home/users/sbusi/apps/ont-guppy/bin/"

Basecalling with methylation-aware latest 'guppy' update

# first one has to make the appropriate directories
cd /scratch/users/sbusi/ONT/ 
mkdir cedric_ont_basecalling/

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling
cp /scratch/users/claczny/ont/fecal_pilot_testing/cluster.json .
cp /scratch/users/claczny/ont/fecal_pilot_testing/config.yaml .
cp /scratch/users/claczny/ont/fecal_pilot_testing/Snakefile .
cp /scratch/users/claczny/ont/fecal_pilot_testing/src/* .
mkdir src && mv *.sh src/.

# edited the above files as needed with correct paths (OR) do the following
ln -sd /scratch/users/claczny/ont/fecal_pilot_testing/data data
ln -sd /scratch/users/claczny/ont/fecal_pilot_testing/dbs dbs
ln -sd /scratch/users/claczny/ont/fecal_pilot_testing/envs envs
ln -sd /scratch/users/claczny/ont/fecal_pilot_testing/scripts scripts

# basecalling run as follows after modifying the "rule all" inputs for basecalling.done
./src/snakemake_run_use_conda.sh

# individual snakemake rules can also be run separately as follows:
snakemake -p --use-conda BASECALL_MERGE_QC
snakemake -p --use-conda COVERAGE_OF_REFERENCES

# LOG can be checked as follows:
less default.log

Chapter III - What about 'non-methlyation-aware' basecalling, a.k.a. `non-mod` basecalling?

The question arose as to what effect this might have
What better way to answer this, than test it?!

#########################
### NO_MOD_basecalled ###
#########################
- Running the basecalled data that was generated `without` any "modifications" or "methylation-awareness"
- To do this, we first have to move the "results" folder to a new name
- then create a new "results" and keep the "basecalled_NO_MOD", but rename this folder to "basecalled"
- Can be done via script as follows:

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/
[./move_results_folder.sh](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/scripts/move_results_folder.sh)

# removing snakemake 'touched' files to re-run freshly the "non-mod-basecalled" data
rm annotate.done
rm assemble_and_coverage.done
rm basecall_merge_qc.done
rm coverage_of_references.done

# due to QosGrpCPUlimit, replaced the 'batch' in "cluster.json" to 'bigmem'
# Running the 'Snakefile' using the wrappers
./src/snakemake_run_use_conda.sh BASECALL_MERGE_QC
./src/snakemake_run_use_conda.sh COVERAGE_OF_REFERENCES
./src/snakemake_run_use_conda.sh ASSEMBLE_AND_COVERAGE
./src/snakemake_run_use_conda.sh TEST
./src/snakemake_run_use_conda.sh TEST_PRODIGAL
./src/snakemake_run_use_conda.sh TEST_DIAMOND
./src/snakemake_run_use_conda.sh POLISH_AND_COVERAGE 
./src/snakemake_run_use_conda.sh DB_MMSEQ2
./src/snakemake_run_use_conda.sh COMPARE_MMSEQ2
./src/snakemake_run_use_conda.sh CONVERT_MMSEQ2
./src/snakemake_run_use_conda.sh PREPARE_FILES
./src/snakemake_run_use_conda.sh PLOT_MMSEQ2

Chapter IV - Metagenome-assembled genomes (MAGs)

When you're having so much fun, why give up so soon?!
We had assemblies, gene calls, coverages, and comparisons between gene calls
So, we decided to go the full monty, a.k.a., "can we bin MAGs and will there be differences (if any)?"
Note: the differences were expected, but we are scientists - we need data to validate!!

################### 
##### Binning #####
###################
- We decided to bin using [MaxBin2](https://sourceforge.net/projects/maxbin2/) and [MetaBAT2](https://github.com/songweizhi/Katana_cmds/wiki/MetaBAT2)
- Next, [DASTool](https://github.com/cmks/DAS_Tool) was used to identify 'high-quality, non-redundant' bins
- However, before we get to the fun parts, we first had to map the short- and long-reads to the "hybrid" assembly using [MetaSPADES](https://github.com/ablab/spades)
- The rules for [MAPPING](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/rules/MAPPING_RULES) and [BINNING](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/rules/BINNING_RULES) provde our strategy to tackle this
- Note: we decided to test both [bwa-mem2](https://github.com/bwa-mem2/bwa-mem2url) and [minimap2](https://github.com/lh3/minimap2url) as mappers to see how they perform/affect our analyses

#################### 
##### Taxonomy #####
####################
- Now that we had MAGs, the logical thing to do was test their completion and contamination among other metrics
- [CheckM](https://github.com/Ecogenomics/CheckM) was used for testing completion
- Then we assigned taxonomy using the [GTDBtk](https://github.com/Ecogenomics/GTDBTk) toolkit
- See the rule here: [TAXONOMY_RULES](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/rules/TAXONOMY_RULES)

Chapter V - Let's throw in some METATRANSCRIPTOMICS

Since we had corresponding metaT data for the 2018-GDB sample, we decided to map the metaT reads to the different assemblies
See METAT_RULES

#########
# metaT #
#########
# Running the metaT processing and mapping as follows:
./src/snakemake_run_use_conda.sh METAT

# calculating average coverage of the assemblies via metaT was done as indicated below
cat sr/GDB_2018_metaT_reads-x-NEB2_MG_S17-megahit_contigs.avg_cov.txt | awk '{x+=$2; next} END{print x/NR}'
#5.52716
cat sr/GDB_2018_metaT_reads-x-lr_barcode07_sr_NEB2_MG_S17-metaspades_hybrid_contigs.avg_cov.txt | awk '{x+=$2; next} END{print x/NR}'
#3.02177
cat lr/GDB_2018_metaT_reads-x-barcode07-flye_contigs.avg_cov.txt | awk '{x+=$2; next} END{print x/NR}'
#30.7486

Chapter VI - the `big kahuna`

This part of the story is very much serendipitous
Having a 1200 line Snakefile can be amazing; unless, for some unknown reason the entire run is triggered even when files are present
Due to this slight inconvenience of time, we decided to break down the Snakefile into "modular" sections
All the rules were broken down to rules and workflows

####################
# MODULAR WORKFLOW #
####################
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling
snakemake -np -s [updated_SNAKEFILE](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/updated_SNAKEFILEurl)

Chapter VII - Validating the analyses and workflows on a different machine

In the wetlab, one repeats experiments, especially a protocol which is tested across several labs
So too, in the computational world, it is ideal to validate the working nature of the workflow on a different machine
In light of the this, we call this part of the story, "The Terminal Transfer"

Transfer to ESB

Transferred the files to ESB to work there due to 'batch' and 'bigmem limitations

# Transferring the files from IRIS to esb-compute-01
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/
si10     # interactive session on IRIS for 10 hours, set in the bash_profile
rsync --rsh='ssh -p 8022' -avp MODULAR_SNAKEFILE susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/
rsync --rsh='ssh -p 8022' -avp results susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/
rsync --rsh='ssh -p 8022' -avp /scratch/users/claczny/fecal_pilot_testing/data \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/.
rsync --rsh='ssh -p 8022' -avp /scratch/users/claczny/ont/fecal_pilot/data/raw \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/data/.
rsync --rsh='ssh -p 8022' -avp /scratch/users/claczny/ont/fecal_pilot_testing/dbs \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/.
rsync --rsh='ssh -p 8022' -avp /home/users/sbusi/apps/DAS_Tool-master/db \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/home/susheel.busi/databases/ 
rsync --rsh='ssh -p 8022' -avp /scratch/users/sbusi/ONT/cedric_ont_basecalling/metaT \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/
rsync --rsh='ssh -p 8022' -avp /home/users/sbusi/apps/db/gtdbtk/release89 \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/home/susheel.busi/databases/
rsync --rsh='ssh -p 8022' -avp /mnt/isilon/projects/ecosystem_biology/NOMIS/DIAMOND/new_nr.dmnd \
  susheel.busi@lcsb-cdc-esb-compute-01.uni.lux:/mnt/data/sbusi/MODULAR_SNAKEFILE/dbs/

# on ESB - # This was a pain - have to download diamond, gtdbtk database, dastool database, paths etc
cd /mnt/data/sbusi/MODULAR_SNAKEFILE

mv /home/susheel.busi/databases/db /home/susheel.busi/databases/DASTool_db
mv /home/susheel.busi/databases/release89 /home/susheel.busi/databases/gtdbtk_release89

# edited the launcher files to remove the "cluster" part
cp src/snakemake_run_use_conda_FINAL.sh src/snakemake_run_on_ESB.sh # edited the 'snakemake_run_on_ESB.sh' file 

cd /mnt/data/sbusi/MODULAR_SNAKEFILE
conda activate snakemake
./src/snakemake_run_on_ESB.sh     (OR)    snakemake -j 4 -p -s updated_SNAKEFILE --rerun-incomplete --use-cond

# Running into diamond database issues
- So extracted sequences and re-made the diamond database
- updated the config/CONFIG.yaml file
.snakemake/conda/4fa897e7/bin/diamond getseq -d /mnt/data/sbusi/MODULAR_SNAKEFILE/dbs/new_nr.dmnd | \
  .snakemake/conda/4fa897e7/bin/diamond makedb -d /mnt/data/sbusi/MODULAR_SNAKEFILE/dbs/diamond_db.dmnd

# 'activate' was missing so exported the PATH as follows:
export PATH="/home/susheel.busi/miniconda3/bin:$PATH"
./src/snakemake_run_on_ESB.sh       # to run the snakefile
grep 'steps' default.log          # to check the progress quickly without having to 'less' the file each time

# Error in checkM databases. Missing!!
# setting up the checkM databaases
cd /home/susheel.busi/.checkm
curl -L -O curl -L -O https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar xzf checkm_data_2015_01_16.tar.gz 
/mnt/data/sbusi/MODULAR_SNAKEFILE/.snakemake/conda/0dfd8768/bin/checkm data setRoot /home/susheel.busi/.checkm/

# re-running only the Taxonomy part of the pipeline since everything up to DASTool was done
# did this by editing the config file
# Finally, everything ran beautifully - MACHINE VALIDATION COMPLETE!!

Chapter VIII - Are we there yet? Nope, we've still got 2019

While all around 2020 has turned out to be a crappy year thanks to SARS-CoV2
2019 was a great year for GDB, who donated sufficient sample to run metaG, metaT, metaP (proteins) and also ONT (long-read) sequencing
This was also the first time in our lab that we ran a MinION flowcell to its exhaustion, acquiring ~31GB of data, thanks to RH's efforts
This section of the narrative will be titled, "the shit storm that wasn't!"

##### 2019 GDB DATA #####
# since the new guppy-3.6 was released with a better accuracy, I downloaded and installed the same locally
# downloaded from: https://community.nanoporetech.com/downloads
cd /home/users/sbusi/apps
wget https://mirror.oxfordnanoportal.com/software/analysis/ont-guppy_3.6.0_linux64.tar.gz
tar -xzvf ont-guppy_3.6.0_linux64.tar.gz     # exported to PATH to be used later

# Preparing the data folders etc.
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling
mkdir 2019_GDB
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB

mkdir data
mkdir data/metaT
cp -vrf ../MODULAR_SNAKEFILE/config .
cp -vrf ../MODULAR_SNAKEFILE/cluster.json .
cp -vrf ../MODULAR_SNAKEFILE/envs .
cp -vrf ../MODULAR_SNAKEFILE/dbs . 
cp -vrf ../MODULAR_SNAKEFILE/rules . 
cp -vrf ../MODULAR_SNAKEFILE/scripts . 
cp -vrf ../MODULAR_SNAKEFILE/src . 
cp -vrf ../MODULAR_SNAKEFILE/workflows .
cp -vrf ../MODULAR_SNAKEFILE/workflows_SNAKEFILE .
cp -vrf ../MODULAR_SNAKEFILE/updated_SNAKEFILE .

mkdir data/raw
mkdir data/raw/short_reads

# symlinking to the raw ONT reads
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/raw

ln -sd /scratch/users/claczny/ont/fecal_pilot-Aug2019_run/data/raw/DonerB_fecal \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/raw/.       # just in case

cp -vrf /scratch/users/claczny/ont/fecal_pilot-Aug2019_run/data/raw/DonerB_fecal/S1_SizeSelected/20190828_0614_MN22103_FAK04421_24db0054/fast5/ \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/multifast5/S1_SizeSelected
cp -vrf /scratch/users/claczny/ont/fecal_pilot-Aug2019_run/data/raw/DonerB_fecal/S3_Gtube/20190828_1842_MN22103_FAK04421_7a557a8b/fast5/ \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/multifast5/S3_Gtube

# copying the metaG reads
cp -vrf /mnt/isilon/projects/lcsb_sequencing/transfer/bioecosystem/Rashi/2019/Nov/MG/fastq/*.gz \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/raw/short_reads/.

# copying the metaT reads
cp -vrf /mnt/isilon/projects/lcsb_sequencing/transfer/bioecosystem/Rashi/2019/Nov/MT/fastq/*Full* \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/metaT/.
cp -vrf /mnt/isilon/projects/lcsb_sequencing/transfer/bioecosystem/Rashi/2019/Nov/MT/fastq/*Half*   \
  /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/metaT/.

# Since there were several metaT samples to choose from (processed with different library preparations)
# Counted the number of reads in the metaT fastq.gz files to decide which sample to use
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/data/metaT
for file in *.fastq.gz
do 
  echo "$file"; zcat "$file" | echo $((`wc -l`/4))
done
# FastSelectFull1_MT_Rashi_S14_R1_001.fastq.gz
# 40105481
# FastSelectFull1_MT_Rashi_S14_R2_001.fastq.gz
# 40105481
# FastSelectFull2_MT_Rashi_S15_R1_001.fastq.gz
# 38354199
# FastSelectFull2_MT_Rashi_S15_R2_001.fastq.gz
# 38354199
# FastSelectHalf1_MT_Rashi_S12_R1_001.fastq.gz
# 40251171
# FastSelectHalf1_MT_Rashi_S12_R2_001.fastq.gz
# 40251171
# FastSelectHalf2_MT_Rashi_S13_R1_001.fastq.gz
# 39937608
# FastSelectHalf2_MT_Rashi_S13_R2_001.fastq.gz
# 39937608

# Chose the "FastSelectHalf1" sample based on reads

## Edited the "config/CONFIG.yaml" files with the appropriate file paths
vi config/CONFIG.yaml
# since there are only 2 runs: S1_SizeSelected, S3_Gtube 
# Also, had to adjust the rules/ASSEMBLY_ANNOTATION_RULES, and the workflows/assembly.smk files
# Adjusted the short-read "name" in the rules files from "NEB2_MG_S17" to "ONT3_MG_xx_Rashi_S11"
# Adjusted the config file to run the first line before going on to the second line and further
# Included the workflows as steps in the [CONFIG.YAML](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/master/config/CONFIG.yamlurl) file

steps: "assembly_annotation metaT mapping mmseq binning taxonomy"

# Basecalling taking too long, so running separately
# created a "basecalling_snakefile" and an associated launcher script "run_basecalling_snakemake.sh"
cp src/snakemake_run_use_conda_FINAL.sh run_basecalling_snakemake.sh
./run_basecalling_snakemake.sh

##### Running non-methylation-aware-basecalled results #####
mv -vf results methylation_aware_results
mkdir results 
cp -vrf methylation_aware_results/basecalled_NO_MOD results/basecalled
rm *.done

# re-running the whole workflow for the non-methylation aware basecalled results
./src/snakemake_run_use_conda_FINAL.sh    

# issues with the "FLYE ASSEMBLY"
# Error: https://github.com/fenderglass/Flye/issues/61
# creating a new conda environment with the latest flye
conda create -n flye flye
conda activate flye
conda env export > envs/flye_v2_7.yaml
sed -i 's/=[^=]*//2g' envs/flye_v2_7.yaml

# adjusted "rules/ASSEMBLY_ANNOTATION_RULES" to look for the new flye conda env
# re-running with new 'flye' environment
./src/snakemake_run_use_conda_FINAL.sh

Chapter IX - If it ain't broke, why fix it? Well, to beautify it!!

2019_GDB data does not have any barcodes, unlike 2018 data which was multiplexed
So, it is essential to create a generic workflow that can be used with others that have barcodes
To improve this, we created the "WORKING_checkpoint_snakefile"
this includes the guppy basecalling, merging fastq(s) and creating "dummy" folders rules
Added new 'rule' and 'workflow' to the respective folders as follows, by adding the updated 'checkpoint' rules
files created are as below

################################
##### Checkpoint Snakefile #####
################################
cp updated_SNAKEFILE checkpoint_SNAKEFILE
# edited the paths to the 'rules' and 'workflows' in the "checkpoint_SNAKEFILE", to include the above files 
# test run the "checkpoint_SNAKEILF"
snakemake -np -s checkpoint_SNAKEFILE --use-conda     # test-run seemed to work fine

Chapter X - Is the short-read metaspades part of the hybrid assembly contributing to the increase in protein counts compared to long-read alone?

A random Tuesday morning conversation, led to the conclusion that we should also use the MetaSPADES assembler and assemble the short-reads
This part of the now, super-long story is called, "sr assembly with metaspades"

###############################
# SR assembly with metaspades #
###############################
# based on cedric's suggestion - to contrast the mmseq results and eliminate the bias of the assembler
# using metaspades to assemble short-reads
# then will do mmseq, and compare megahit, metaspades and metaspades_hybrid
# saved the below analyses method in the "scripts" folders as "cdhit_comparisons.sh"

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
# edited the following files
1. rules/ASSEMBLY_ANNOTATION_rule         # line 783 == added rule "assemble_short_reads_w_metaspades"
2. workflows/assembly_annotation.smk      # line 126 == added "metaspades" as a new target
                                          # line 90 == added 'metaspades' in TEST prodigal
3. cluster.json     # addded parameters for the new rule
4. rules/MMSEQ_RULES          # need to edit the "prepare_plot_files.sh" to include the new additions
5. workflows/mmseq.smk

MMSEQ

mmseqs2 comparisons were performed as part of the snakefile
Required adjustment of the plotting prep and R scripts to account for the new "metaspades" group

cp mmseq_plots.R [mmseq_plots_w_metaspades.R](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/scripts/mmseq_plots_w_metaspades.R)
cp prepare_plot_files.sh [prepare_plot_files_w_metaspades.sh](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/checkpoint_snakefile/2019_GDB/scripts/prepare_plot_files_w_metaspades.sh)

CD-HIT

Using cd-hit compared metaspades_sr and metaspades_hybrid proteins to see which contributes more to the mmseq comparisons

cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
mkdir cd-hit && cd cd-hit
ln -s ../results/annotation/proteins/metaspades/ONT3_MG_xx_Rashi_S11/final.contigs.faa metaspades.faa
ln -s ../results/annotation/proteins/metaspades_hybrid/lr_no_barcode-sr_ONT3_MG_xx_Rashi_S11/contigs.faa metaspades_hybrid.faa

# getting the count of partially-called genes based on Prodigal output
grep -E -o "partial.{0,3}" metaspades.faa | grep '10\|01\|11' | wc -l > metaspades_partial_counts.txt
grep -E -o "partial.{0,3}" metaspades_hybrid.faa | grep '10\|01\|11' | wc -l > metaspades_hybrid_partial_counts.txt

# getting the counts into one file
for file in *.txt; do echo $file; cat $file; done >> counts
cat counts | paste - - > partial_gene_counts.txt

# editing fasta headers to make it easier after cd-hit clustering
si    # interactive
conda activate bbmap
rename.sh in=metaspades.faa out=spades.faa prefix=spades ignorejunk=t
rename.sh in=metaspades_hybrid.faa out=hybrid.faa prefix=hybrid ignorejunk=t

conda activate cd-hit
cd-hit-2d -i spades.faa -i2 hybrid.faa -o spades_hybrid -c 0.9 -n 5 -d 0 -M 16000 -T 8
cd-hit-2d -i hybrid.faa -i2 spades.faa -o hybrid_spades -c 0.9 -n 5 -d 0 -M 16000 -T 8

# determining number of unique sequences
# according to http://weizhongli-lab.org/lab-wiki/doku.php?id=cd-hit-user-guide#cd-hit-2d
# CD-HIT-2D outputs two files: a fasta file of proteins in "db2" that are not similar to db1 and a text file that lists similar sequences between db1 & db2
# non-mod-basecalled
grep -c '>' spades_hybrid   # db2 == hybrid
# 63911 == no.unique in hybrid
grep -c '>' hybrid_spades   # db2 == spades
# 27526 == no.unique in spades

# "methylation-aware"_basecalled
grep -c '>' spades_hybrid   # db2 == hybrid
# 61828 == no.unique in hybrid
grep -c '>' hybrid_spades   # db2 == spades
# 30807 == no.unique in spades

# using auxilliary scripts to merge the cluster (output) files
# making plots for all .clstr files (http://weizhongli-lab.org/lab-wiki/doku.php?id=cd-hit-user-guide#cd-hit-2d)
for file in *.clstr 
do      
    echo "$file"
    plot_len1.pl "$file"  \
           1,2-4,5-9,10-19,20-49,50-99,100-299,500-99999   \
                 10-59,60-149,150-499,500-1999,2000-999999
done >> cluster_plots

cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/
mv -vf cd-hit cd-hit_non_mod_basecalled   # renamed the folder since it was run with the "non_mod_basecalled" data
# Repeated the cd-hit analyses for the "methylation-aware" basecalled data
# The data shown above is for the "methylation-aware" run juxtaposed against the 'non-methylation' basecalled data

Chapter XI - The miscellaneous or nearly-forgotten side projects

Due to the multifaceted nature of the beast, i.e. the modular workflow, we tested several aspects separately
For example: since we used two mappers bwa-mem and minimap for the reads, we binned each sample separtely based on the mapper
Additionally, we needed to merge bam files for the "hybrid"-binning, so we compared bins using sourmash and compared assemblies using quast

QUAST - 2019_GDB

Like any good experiment, sometimes the basics are key evaluators of the data
To address this, used (meta)QUAST to assess how the assemblies look

# Testing assembly stats using quast on the 2019_GDB "methylation-aware" data
srun -p interactive -t 4:00:00 --pty bash -i    # interactive session
conda activate quast
cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/
mkdir quast_methylation_aware
cd quast_methylation_aware
ln -s  ../results/assembly/*.fa .
ln -s ../results/assembly/metaspades/ONT3_MG_xx_Rashi_S11/final.contigs.fna metaspades.fa
rm bwa* mmi*      # removing non-essential files

metaquast.py --max-ref-number 0 --threads 24 *fa -o methylation_2019_quast_results

# Testing assembly stats using quast on the 2019_GDB "non_mod" basecalled data
srun -p interactive -t 4:00:00 --pty bash -i    # interactive session
conda activate quast
cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/
mkdir quast_non_mod
cd quast_non_mod
ln -s ../non_methylation_aware_results/assembly/*.fa .
ln -s ../non_methylation_aware_results/assembly/metaspades/ONT3_MG_xx_Rashi_S11/final.contigs.fna metaspades.fa
rm bwa* mmi*      # removing non-essential files

metaquast.py --max-ref-number 0 --threads 24 *fa -o non_mod_2019_quast_results

## comparing methylation vs non-methylation assemblies ##
# running quast to compare methylation vs non-methylation assemblies
conda activate quast
cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/
mkdir quast_methylation_v_non_comparison
cd quast_methylation_v_non_comparison

# getting the methylation assemblies
ln -s  ../results/assembly/*.fa .
ln -s ../results/assembly/metaspades/ONT3_MG_xx_Rashi_S11/final.contigs.fna metaspades.fa
# rename files
for file in *.fa; do mv "$file" meth_"$file"; done

# getting the non-methylation assemblies
ln -s ../non_methylation_aware_results/assembly/*.fa .
ln -s ../non_methylation_aware_results/assembly/metaspades/ONT3_MG_xx_Rashi_S11/final.contigs.fna metaspades.fa
rm bwa* mmi*      # removing non-essential files

metaquast.py --max-ref-number 0 --threads 24 *fa -o quast_comparision_meth_v_non_results

# downloaded all the folders to desktop at ~/Documents/Nanopore_ONT/.
rsync --rsh='ssh -p 8022' -auvp \
  sbusi@access-iris.uni.lu:/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/quast_methylation_v_non_comparison/quast_comparision_meth_v_non_results \
    ~/Downloads/.

rsync --rsh='ssh -p 8022' -auvp \
  sbusi@access-iris.uni.lu:/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/quast_non_mod/non_mod_2019_quast_results \
    ~/Downloads/.

rsync --rsh='ssh -p 8022' -auvp \
  sbusi@access-iris.uni.lu:/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/quast_methylation_aware/methylation_2019_quast_results \
    ~/Downloads/.

SOURMASH

To check if the bins across different metaspades_hybrid 'bam' files are similar
to decide whether to go with sr_bam or merged_bam files
therefore, comparing ONLY the {mapper}_{reads}metaspades_hybrid vs {mapper}{reads}_metaspades_hybrid bins

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/non_mod_basecalled_results
cd Binning
mkdir sourmash
mkdir bins
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/non_mod_basecalled_results/Binning/sourmash/bins

for file in `cat list`
do
  for i in `ls /scratch/users/sbusi/ONT/cedric_ont_basecalling/non_mod_basecalled_results/Binning/"$file"/dastool_output/"$file"_DASTool_bins/*.fa`
    do 
      j=`echo $i | cut -d'.' -f2,3`
      ln -s "$i" "$file"_"$j"
    done
done

# renaming files that do not have the '.fa' extension but are in reality fasta files
for f in *_sub; do 
    mv -- "$f" "${f%_sub}_sub.fa"
done

# Running sourmash
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/non_mod_basecalled_results/Binning/sourmash
si    # interactive
conda activate sourmash
sourmash compute -k 31 bins/*.fa
sourmash compare *.sig -o ont_hybrid_bins_sourmash_comp
saving labels to: ont_hybrid_bins_sourmash_comp.labels.txt
# saving distance matrix to: ont_hybrid_bins_sourmash_comp
# sourmash compare *.sig --csv ont_hybrid_bins_sourmash_comp.csv
sourmash plot ont_hybrid_bins_sourmash_comp --labels
sourmash plot ont_hybrid_bins_sourmash_comp --labels --pdf
# downloaded the .csv file to the "Desktop" and used an additional Rscript for plotting

dREP

Binning the metaspades_hybrid assembly was a conundrum
Binners require a bam file (typically) to estimate coverage
So, the question arose which bam file should/could we use?
We had the following options to choose from: bwa_sr_metaspades_hybrid, bwa_lr_metaspades_hybrid, bwa_merged_metaspades_hybrid, mmi_sr_metaspades_hybrid, mmi_lr_metaspades_hybrid, and mmi_merged_metaspades_hybrid
Here above, sr==short-reads mapped to assembly, lr==long-reads, and merged==merging of sr- and lr- bam files
If you haven't caught on to how we do things around here - we do them all
This is the part of the story, I like to call "BinCHECK" using drep

# checking to see if the bins reconstructed from different 'bam' files are similar
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
mkdir bins
cp /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/Binning/checkm_output/sample_list bins/.
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/bins

# symlink files by renaming incrementally
export i=1
for f in `cat sample_list`
do
  for file in /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/Binning/"$f"/dastool_output/"$f"_DASTool_bins/*.fa
  do
     # Copy file to new location, safer than move in case you make a mistake
     ln -s "$file" /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/bins/"$i"_"$f".fa

   # Increment
   let i=i+1
   # Reset if exceeds 12
   if [ $i -gt 500 ]
   then
      export i=1
   fi
  done
done

# removing some broken links
rm 20_metaspades.fa 21_metaspades_hybrid.fa

# Creating a new conda environment for drep
srun -p interactive -t 4:00:00 --pty bash -i
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
conda create -n drep -c bioconda drep
conda activate drep
conda install -c bioconda fastani
cd ~/apps/miniconda3/envs/drep/
wget https://ani.jgi.doe.gov/download_files/ANIcalculator_v1.tgz
tar -xzvf ANIcalculator_v1.tgz
cp ANIcalculator_v1/ANIcalculator bin/.
# testing drep installation
dRep bonus testDir --check_dependencies

# Running drep on the bins to see similarities
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/bins
dRep compare drep_output -g ./*.fa

Assessing MAGs

Although expected, as stated somewhere above, "we are scientists and so, we need evidence!"
I decided to confirm that the bins recovered via each of the assembly methods, and additionally the mapping methods would be different
To get there, first we needed to link the CheckM and GTDBtk outputs to the binned MAGs

# Working directory
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/Binning/checkm_output

# making a list of all the samples
ls -1 *.txt | sed 's/_output.txt//g' > sample_list

# collating the checkm_output to a readable format
for file in `cat sample_list`
do
  sed '1,/Bin Id/d;/INFO/,$d' "$file"_output.txt | \
  sed '/^-----/d' | \
    awk '{print $1"\t"$13"\t"$14"\t"$15}' | \
      sed $'1 i\\\nbin\tcompletion\tcontamination\tstrain_heterogeneity' > "$file"_checkm.txt 
done

# collecting the gtdbtk information
for file in `cat sample_list`
do
  awk '{print $1"\t"$2}' ../gtdbtk_output/"$file"/classify/gtdbtk.bac120.summary.tsv | sed 's/user_genome/bin/g' > "$file"_gtdbtk.txt
done

# merging files based on a common first column
for file in `cat sample_list`
do
  awk -F"\t" 'BEGIN{OFS="\t"} {if (NR==FNR) {a[$1]=$2"\t"$3"\t"$4; next} if ($1 in a) {print $1, $2, a[$1]}}' \
    "$file"_checkm.txt "$file"_gtdbtk.txt > "$file"_taxonomy_contamination.txt 
done

# keeping only those bins with completion over 60% and contamination below 10% (this number can be adjusted as needed)
for file in `cat sample_list`
do
  awk '$3>=60 && $4 <10 {print}' "$file"_taxonomy_contamination.txt > "$file"_HQ_bins.txt
done

# The "$file"_HQ_bins.txt files were subsequently analysed on the Desktop in R using the "MAGS_2019_GDB_analyses.R" script

VizBin

The binners did great for the sr- and hybrid assemblies, however, fared poorly with the lr- (flye) assembly
A sign of a good researcher is that they never trust their own data, and validate it twelve-ways-to-Sunday
Not saying that I'm good, but I like to validate, for lack-of-a-better-word in my vocabulary and at my discretion, the "heck" out of something
And since "seeing is believing", I took the assemblies to our very own CCL's VizBin implemented via BusyBee

# Running flye assemblies on the BusyBee web server
# Downloaded the "methylation-aware" flye assembly to the Desktop
# And uploaded it to the BusyBee server
# Changed the minimum contig length to 1500
# Run ID: b9d043cf-9a55-4ba8-9d1b-2bef4861cd61

QUAST

Testing assembly stats using quast

si     # interactive session on hpc-IRIS
conda activate quast
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/Binning
metaquast.py --max-ref-number 0 --threads 24 *fna -o non_mod_basecalled_quast_results
# downloaded the folder to desktop at ~/Documents/Nanopore_ONT/.
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/mod_basecalled_results/assembly 
ln -sd flye/lr/merged/barcode07/assembly.fasta flye.fa
ln -sd megahit/NEB2_MG_S17/final.contigs.fa megahit.fa
ln -sd metaspades_hybrid/lr_barcode07-sr_NEB2_MG_S17/contigs.fasta metaspades.fa

metaquast.py --max-ref-number 0 --threads 24 *fa -o methylation_basecalled_quast_results
# downloaded the folder to desktop at ~/Documents/Nanopore_ONT/.

GITLAB

copying all changed files to gitlab: https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling
git pull https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab.git

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/ont_pilot_gitlab
cp ../Snakefile .

git checkout work
git add Snakefile 
git commit -m "any message of note"
git push 
# go to the provided URL, and then submit the merge request, for it to be approved later.

SLURM files

moving all slurm files into one folder for 'cleaner' appearance

mkdir slurm_files
mv slurm* slurm_files/.

Chapter XII - Work(s) in Progress

CIRCULARY POTENTIAL

This is where the brain doesn't know when to shut off, and we went on a tanget
The question being - does flye or metaspades_hybrid allow for circulisation of the contigs, especially for binning

############
# CIRCULAR #           "TO-DO"
############
# checking how many contigs in each bin
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/results/Binning
# made a list with the following. excluded 'flye' 'cos bins are not real
bwa_lr_metaspades
bwa_merged_metaspades
bwa_sr_metaspades
mmi_lr_metaspades
mmi_merged_metaspades
mmi_sr_metaspades
megahit

for i in `cat list`
do 
    for file in "$i"/dastool_output/"$i"_DASTool_bins/*.fa
    do
        echo "$file"
        grep -c '>' "$file"
    done >> "$i"_bins_contigs
done 
# None of the bins came from a single contig

# checking the circularity of the bins
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/
mkdir circular
cd circular

conda activate [circlator](https://github.com/sanger-pathogens/circlator)
conda install [canu](https://github.com/marbl/canu)
circlator all ../results/assembly/flye/lr/merged/barcode07/assembly.fna ../results/basecalled/merged/barcode07.fna flye_circlator_output

# failing in 10-hour interactive so submitted as bash script
sbatch -p bigmem -t 2-00:00:00 -N 1 -n 24 run_circlator.sh (OR) sbatch -p bigmem -t 0-16:00:00 -N 1 -n 24 run_circlator.sh      # failed after a while, so going to manual methods below

# Checking graph with [Bandage](https://github.com/rrwick/Bandage)
conda activate bandage
Bandage info ../results/assembly/metaspades_hybrid/lr_barcode07-sr_NEB2_MG_S17/assembly_graph_with_scaffolds.gfa

# removing singletons == DOES NOT WORK yet..
cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/circular
ml swenv/bioinfo-env/v0.1-20170602-production
ml swenv/default-env/v0.1-20170602-production
ml bio/BioPerl/1.7.1-foss-2017a-Perl-5.24.1

perl -0076 -ne ‘{chomp; m/(.*):(.*);/; print “>$_” if $1 ne $2};’ < /scratch/users/sbusi/ONT/cedric_ont_basecalling/results/assembly/metaspades_hybrid/lr_barcode07-sr_NEB2_MG_S17/assembly_graph.fastg > megahit_singleton_removed.fastg

GraphBin

https://github.com/Vini2/GraphBin - for binning from graphs
downloaded the "methylation-aware" flye assembly and .gfa files to desktop
also downloaded the flye_maxbin_scaffolds2bin.tsv and saved as .csv - working on the desktop (MAC)

cd ~/Documents/Nanopore_ONT
git clone https://github.com/Vini2/GraphBin.git
cd GraphBin/
conda env create -f environment.yml
conda activate graphbin

python graphbin.py --assembler flye --graph assembly_graph.gfa --paths assembly.fasta \
  --binned flye_maxbin.scaffolds2bin.csv --output flye_methylation_graphbin --max_iteration 10

# downloaded the spades contigs.fasta, contigs.paths, assembly_graph_with_scaffolds.gfa, and the scaffolds2bin.tsv file from bwa_sr_metaspades_hybrid to desktop
cd ~/Documents/Nanopore_ONT/GraphBin/spades_input
conda activate graphbin

# removed the "maxbin_output." from the .tsv file and saved as .csv
cd /Users/susheel.busi/Documents/Nanopore_ONT/GraphBin
python graphbin.py --assembler spades --graph spades_input/assembly_graph_with_scaffolds.gfa \
  --paths spades_input/contigs.paths --binned spades_input/bwa_sr_metaspades_hybrid_maxbin.scaffolds2bin.csv \
    --output spades_graphbin_output

AGB

Trying AGB for circularity analyses

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
conda create -c almiheenko -c bioconda -n AGB agb
chmod -R 775 /home/users/sbusi/apps/miniconda3/envs/AGB/

si
conda activate AGB
agb.py -i methylation_aware_results/assembly/flye/lr/merged/no_barcode/ -a Flye -o flye_methylation_agb_output

MAPPABILITY

Assessing how "mappable" the assemblies are with the short, long and merged bam files
i.e. the number of reads mapping back to the assemblies

cd /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/mapping
conda create -n analysis -c bioconda samtools
conda activate analysis
conda env export > /scratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/envs/analysis.yaml
srun -p interactive -t 4:00:00 --pty bash -i    # interactive session

# getting the "mappability" using a launcher script for the following:
sbatch -t 1:00:00 -N 1 -n 4 [get_mappability.sh](https://git-r3lab.uni.lu/susheel.busi/ont_pilot_gitlab/-/blob/master/2019_GDB/scripts/get_mappability.sh)
>>>
# metaG
samtools flagstat sr/bwa_mem/ONT3_MG_xx_Rashi_S11_reads-x-ONT3_MG_xx_Rashi_S11-megahit_contigs.bam > sr_mappability.txt
samtools flagstat lr/merged/no_barcode/no_barcode_reads-x-no_barcode-flye_contigs.bam > lr_mappability.txt

samtools flagstat bwa_merged_metaspades_hybrid/bwa_merged_metaspades_hybrid.bam > hybrid_bwa_merged.txt 
samtools flagstat bwa_sr_metaspades_hybrid/bwa_sr_metaspades_hybrid.bam > hybrid_bwa_sr.txt
samtools flagstat bwa_lr_metaspades_hybrid/bwa_lr_metaspades_hybrid.bam > hybrid_bwa_lr.txt

samtools flagstat mmi_merged_metaspades_hybrid/mmi_merged_metaspades_hybrid.bam > hybrid_mmi_merged.txt
samtools flagstat mmi_sr_metaspades_hybrid/mmi_sr_metaspades_hybrid.bam > hybrid_mmi_sr.txt
samtools flagstat mmi_lr_metaspades_hybrid/mmi_lr_metaspades_hybrid.bam > hybrid_mmi_lr.txt

# metaT
samtools flagstat ./metaT/sr/FastSelectHalf1_MT_Rashi_S12_reads-x-lr_no_barcode_sr_ONT3_MG_xx_Rashi_S11-metaspades_hybrid_contigs.bam > hybrid_metaT_sr.txt
samtools flagstat ./metaT/sr/FastSelectHalf1_MT_Rashi_S12_reads-x-ONT3_MG_xx_Rashi_S11-megahit_contigs.bam > megahit_metaT_sr.txt
>>> 

# collating all the results in one place
for file in *.txt; do echo $file; egrep 'total | mapped' $file | head -n 2; done | paste - - - | \
   sed $'1 i\\\nsample\ttotal\tmapped\tpercent_mapped' | sed 's/ + 0 in total (QC-passed reads + QC-failed reads)//g' | \
     sed s'/ + 0 mapped//g' | sed 's/ : N\/A)//g' | sed s'/(//g' | sed s'/.txt//g' | \
       awk '{print$0=$1"\t"$2"\t"$3"\t"$4}' > mappability_index.tsv


#### Question: will the hybrid make a difference if the smaller contigs are removed
# Answer: lets disregard the smaller contigs, starting at 1000bp and then get the "mappability_index"
# getting read counts for only contigs > (greater than) 1Kbp
# samtools idxstats output:: contig_name, contig_length, mapped_reads, unmapped_reads; note: * at the end indicates number of unmapped reads

# metaG
samtools idxstats sr/bwa_mem/ONT3_MG_xx_Rashi_S11_reads-x-ONT3_MG_xx_Rashi_S11-megahit_contigs.bam > megahit_counts.txt
samtools idxstats lr/merged/no_barcode/no_barcode_reads-x-no_barcode-flye_contigs.bam > flye_counts.txt

samtools idxstats bwa_merged_metaspades_hybrid/bwa_merged_metaspades_hybrid.bam > hybrid_bwa_merged_counts.txt 
samtools idxstats bwa_sr_metaspades_hybrid/bwa_sr_metaspades_hybrid.bam > hybrid_bwa_sr_counts.txt
samtools idxstats bwa_lr_metaspades_hybrid/bwa_lr_metaspades_hybrid.bam > hybrid_bwa_lr_counts.txt

samtools idxstats mmi_merged_metaspades_hybrid/mmi_merged_metaspades_hybrid.bam > hybrid_mmi_merged_counts.txt
samtools idxstats mmi_sr_metaspades_hybrid/mmi_sr_metaspades_hybrid.bam > hybrid_mmi_sr_counts.txt
samtools idxstats mmi_lr_metaspades_hybrid/mmi_lr_metaspades_hybrid.bam > hybrid_mmi_lr_counts.txt

# metaT
samtools idxstats ./metaT/sr/FastSelectHalf1_MT_Rashi_S12_reads-x-lr_no_barcode_sr_ONT3_MG_xx_Rashi_S11-metaspades_hybrid_contigs.bam > hybrid_metaT_sr_counts.txt
samtools idxstats ./metaT/sr/FastSelectHalf1_MT_Rashi_S12_reads-x-ONT3_MG_xx_Rashi_S11-megahit_contigs.bam > megahit_metaT_sr_counts.txt

# collating the sums for contigs larger than 1000 bp
for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 1000 {sum += $3} END {print sum}' "$file"  # sum of number of reads (3rd column) if 2nd column is greater/equal to 1000
  awk '{sum+=$3+$4} END {print sum}' "$file"            # sum of 3rd and 4th columns to get total reads
done | paste - - - | awk '$4=100*$2/$3' | \             # getting the percentage from columns 2 and 3 and writing to column-4
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_1000bp.txt

# collating the sums for contigs larger than 1500 bp
for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 1500 {sum += $3} END {print sum}' "$file"
  awk '{sum+=$3+$4} END {print sum}' "$file"
done | paste - - - | awk '$4=100*$2/$3' | \
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_1500bp.txt

# collating the sums for contigs larger than 2000 bp
for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 2000 {sum += $3} END {print sum}' "$file"
  awk '{sum+=$3+$4} END {print sum}' "$file"
done | paste - - - | awk '$4=100*$2/$3' | \
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_2000bp.txt

# collating the sums for contigs larger than 5000 bp
for file in *counts.txt
do 
  echo "$file"
  awk '$2 >= 5000 {sum += $3} END {print sum}' "$file"      
  awk '{sum+=$3+$4} END {print sum}' "$file"
done | paste - - - | awk '$4=100*$2/$3' | \
  sed $'1 i\\\nsample\tmapped_reads\ttotal\tpercent_mapped' > mappability_over_5000bp.txt

## Assessing the number of "unique_mapping_reads" from BAM files ##
-based on information here: https://bioinformatics.stackexchange.com/questions/508/obtaining-uniquely-mapped-reads-from-bwa-mem-alignment/4542#4542
samtools view -h megahit/megahit.bam | grep -v "XT:A:U"| wc -l
# 51826913
# total = 
# percent_mapped =
samtools view flye/flye.bam | grep -v "XT:A:U" | uniq | wc -l
# 8351240
# total = 
# percent_mapped =
samtools view bwa_merged_metaspades_hybrid/bwa_merged_metaspades_hybrid.bam | grep -v "XT:A:U" | wc -l
# 60114498
# total = 
# percent_mapped =

NOTE: all of the above can be run using the **rules/ANALYSIS_RULES**

CRISPR

A key feature of the long-read technology is the ability to sequence repeat regions
So, CRISPRS are exactly that, and why don't we look for this tiny little pesky repeats that give short-reads such a bad name
We decided to go with CASC and minced
Rob Edwards from UCSD suggested trying the CRISPRDetect, but @SBB got stuck in perl-hell, so we decided to put that away for another time
NOTE: rather than following all of this long-winded way, use the "workflows/analysis.smk" for the CRISPR detection and quantification

cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB
mkdir crispr
cd crispr
srun -p interactive -t 4:00:00 --pty bash -i
conda activate analysis
conda install -c anaconda perl
# downloading and installing CASC: https://github.com/dnasko/CASC
git clone https://github.com/dnasko/CASC.git

# installing locally
cd CASC
PERL Makefile.PL PREFIX=/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr
make
make test
make install

# updating the PATHs
export PATH=$PATH:/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/bin
export PERL5LIB=/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/lib/site_perl

# testing with 'flye'
casc -i 
/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/assembly/flye.fa \
  -o /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/flye_crispr_output \
    -n 12 --conservative

# running in a loop for all samples since it only takes 20 minutes per sample
for file in /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/assembly/*.fa
do 
  casc -i "$file" \
  -o /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/$(basename "$file")_crispr_output -n 12 --conservative
done

# collating the number of CRISPR spacers
awk '$8 == "true" {sum += $2} END {print sum}' flye_crispr_output/flye.results.txt
# 28
awk '$8 == "true" {sum += $2} END {print sum}' bwa_lr_metaspades_hybrid.fa_crispr_output/bwa_lr_metaspades_hybrid.results.txt
# 61
awk '$8 == "true" {sum += $2} END {print sum}' bwa_merged_metaspades_hybrid.fa_crispr_output/bwa_merged_metaspades_hybrid.results.txt
# 61

for i in `cat files`
do 
  echo "$i"
  awk '$8 == "true" {sum += $2} END {print sum}' "$i".fa_crispr_output/"$i".results.txt
done | paste - - | sed $'1 i\\\nassembly\tnumber_CRISPR_spacers'  > casc_CRISPR_output.txt

##### MinCED #####
# https://github.com/ctSkennerton/minced
git clone https://github.com/ctSkennerton/minced.git
cd minced 
conda install -c bioconda java-jdk
make 

for file in /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/results/assembly/*.fa
do
  ln -s $file .
done 

for file in *.fa
do 
  ./minced "$file" "$file".txt "$file".gff
  echo "$file" >> minced_CRISPR_output.txt
  grep -c 'CRISPR' "$file".txt >> minced_CRISPR_output.txt
done

## Based on a discussion with Rob Edwards from UCSD, decided to try out CRISPRdetect-parser
# from here: https://github.com/hwalinga/crisprdetect-parser
cd /mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr
conda activate emboss
conda install -c bioconda viennarna
conda install -c anaconda perl

git clone https://github.com/hwalinga/crisprdetect-parser.git
cd crisprdetect-parser/CRISPRDetect_2.4_grep_patched
export PATH=$PATH:/home/users/sbusi/apps/miniconda3/envs/cd-hit/bin     # for cd-hit
export PATH=$PATH:/home/users/sbusi/apps/miniconda3/envs/anvio6/bin     # for blastn
export PATH=$PATH:/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/bin    # for perl
export PERL5LIB=/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/lib/site_perl      # for perl

wget http://search.cpan.org/CPAN/authors/id/Y/YA/YANICK/Parallel-ForkManager-1.19.tar.gz
tar xvzf Parallel-ForkManager-1.19.tar.gz
cd Parallel-ForkManager-1.19
perl Makefile.PL PREFIX=/mnt/lscratch/users/sbusi/ONT/cedric_ont_basecalling/2019_GDB/crispr/crisprdetect-parser/CRISPRDetect_2.4_grep_patched/Parallel-ForkManager-1.19
make test && make install

####### INSTALLATION NOT WORKING ########

PLASMIDS

LDN is running PathoFact on the ONT samples, so obtained the plasmid data through the MGE/plasmid/$i_plasmids_selected.tsv file generated in PathoFact

# checking the number of plasmids in the assemblies
# grep -hnr "plasmid" $i_plasflow_prediction_final.tsv > $i_plasmids_selected.tsv

for file in /scratch/users/ldenies/ONT/PathoFact/PathoFact_results/MGE/plasmid/*final.tsv
do 
  echo $file
  grep -c 'plasmid.' "$file"
done
# /scratch/users/ldenies/ONT/PathoFact/PathoFact_results/MGE/plasmid/flye_plasflow_prediction_final.tsv
# 280
# /scratch/users/ldenies/ONT/PathoFact/PathoFact_results/MGE/plasmid/megahit_plasflow_prediction_final.tsv
# 5186
# /scratch/users/ldenies/ONT/PathoFact/PathoFact_results/MGE/plasmid/metaspades_hybrid_plasflow_prediction_final.tsv
# 3211
# /scratch/users/ldenies/ONT/PathoFact/PathoFact_results/MGE/plasmid/metaspades_plasflow_prediction_final.tsv
# 5718

@SBB: is PlasFlow affected by the "fragmentation or lack thereof" of the assemblies. Why are very few plasmids detected in the FLYE results?
# so checked the number of contigs used by plasflow for the analyses
# number of contigs used as input:
#   2533982 flye_filtered.part_001.fna
#    975951 megahit_filtered.part_001.fna
#    829650 megahit_filtered.part_002.fna
#    677918 megahit_filtered.part_003.fna
#   2483519 total
#   1992017 metaspades_filtered.part_001.fna
#    316220 metaspades_filtered.part_002.fna
#    185613 metaspades_filtered.part_003.fna
#   2493850 total
#   2847904 metaspades_hybrid_filtered.part_001.fna
#    306145 metaspades_hybrid_filtered.part_002.fna
#      7023 metaspades_hybrid_filtered.part_003.fna
#   3161072 total

###############

Troubleshooting - SNAKEMAKE wanting to run everything

DIAMOND - RE-RUN

###############

When trying to run a single step., i.e. "TEST_DIAMOND" from the "ASSEMBLY_ANNOTATION_RULES"
SNAKEMAKE wanted to re-run everything including assemblies
using the -r flag in a dry-run we saw that some files had a newer timestamp

TODO: figure out why certain folders have a newer timestamp

The "WORKAROUND" for this is presented below

# Looking at timestamp with "stat" and editing "timestamp" with touch
1. stat results/annotation/proteins/flye/lr/merged/no_barcode/assembly.faa
2. # saved the original fasta file
cp results/assembly/flye/lr/merged/no_barcode/assembly.fasta results/assembly/flye/lr/merged/no_barcode/assembly.fasta.orig
3. # replacing the timestamp
touch -d '14 May 2020 05:17' results/assembly/flye/lr/merged/no_barcode/assembly.fasta
4. # dry-run
snakemake -npr -s workflows/assembly_annotation.smk results/annotation/diamond/flye/lr/merged/no_barcode/assembly.tsv

# created a new bash script to RE-RUN "DIAMOND with the "nr_uniprot_trembl" database
./src/run_diamond_only.sh "results/annotation/diamond/metaspades_hybrid/lr_no_barcode-sr_ONT3_MG_xx_Rashi_S11/contigs.tsv \
  results/annotation/diamond/megahit/ONT3_MG_xx_Rashi_S11/final.contigs.tsv \
    results/annotation/diamond/flye/lr/merged/no_barcode/assembly.tsv"

Notes - for 2019_GDB data

Since diamond was inadvertently run on the "new_nr.dmnd (BLAST nucleotide database), the folder was renamed mv diamond diamond_new_nr

Notes - for 2018_GDB data

Flye does not produce metabat bins. Reason == only 2 contigs (less than 1500 bp) have coverages over 1.
Should we adjust the "minCV" parameter for metabat?
Added the snakemake "--keep-going" or "-k" flag to the launcher to finish the independent jobs, and restarted
Need to fix checkm and gtdbktk runs
CheckM was throwing file missing erros, so did the following: "cp /home/users/sbusi/apps/checkm/* /home/users/sbusi/.checkm/."
Use "checkm data setRoot" to specify the location of CheckM database files.

ONT_pilot_w_GDB_samples

This is the story of how harmless curiosity for Oxford Nanopore Technology lead to an awesome learning experience!

It started with a conversation titled, "what is the current state of the ONT data from 2018 and 2019?"

Chapter I - Running clustering analyses using Prodigal output

The story starts with first working in my own folders

So, I decided to pursue mmseqs2 instead

And this is where the fun begins..

Chapter II - Running methylation-aware basecalling to improve the otherwise dismal sequencing error rates from ONT

Chapter III - What about 'non-methlyation-aware' basecalling, a.k.a. non-mod basecalling?

Chapter IV - Metagenome-assembled genomes (MAGs)

Chapter V - Let's throw in some METATRANSCRIPTOMICS

Chapter VI - the big kahuna

Chapter VII - Validating the analyses and workflows on a different machine

Transfer to ESB

Chapter VIII - Are we there yet? Nope, we've still got 2019

Chapter IX - If it ain't broke, why fix it? Well, to beautify it!!

Chapter X - Is the short-read metaspades part of the hybrid assembly contributing to the increase in protein counts compared to long-read alone?

MMSEQ

CD-HIT

Chapter XI - The miscellaneous or nearly-forgotten side projects

QUAST - 2019_GDB

SOURMASH

dREP

Assessing MAGs

VizBin

QUAST

GITLAB

SLURM files

Chapter XII - Work(s) in Progress

CIRCULARY POTENTIAL

GraphBin

AGB

MAPPABILITY

CRISPR

PLASMIDS

Troubleshooting - SNAKEMAKE wanting to run everything

DIAMOND - RE-RUN

TODO: figure out why certain folders have a newer timestamp

Notes - for 2019_GDB data

Notes - for 2018_GDB data

Chapter II - Running `methylation-aware` basecalling to improve the otherwise dismal sequencing error rates from ONT

Chapter III - What about 'non-methlyation-aware' basecalling, a.k.a. `non-mod` basecalling?

Chapter VI - the `big kahuna`